1 subroutine elmStatsLhs( x, iBC, iper, ilwork ) 2c----------------------------------------------------------------------- 3c compute the necessary terms for the statistics projection 4c matrices. 5c----------------------------------------------------------------------- 6 use stats 7 use pointer_data 8 9 include "common.h" 10 11 real*8 x(numnp,3) 12 integer iBC(nshg), iper(nshg), ilwork(nlwork) 13 14 real*8, allocatable :: xl(:,:,:) 15 real*8, allocatable :: lStsVec(:,:,:) 16 17c 18c.... loop over element blocks 19c 20 stsVec = zero 21 22 do iblk = 1, nelblk 23 iel = lcblk(1,iblk) 24 lcsyst = lcblk(3,iblk) 25 nenl = lcblk(5,iblk) ! no. of vertices per element 26 nshl = lcblk(10,iblk) 27 ndofl = lcblk(8,iblk) 28 npro = lcblk(1,iblk+1) - iel 29 30 allocate ( xl(npro,nenl,3) ) 31 allocate ( lStsVec(npro,nshl,nResDims) ) 32c 33c.... localize needed data 34c 35 call localx ( x, xl, mien(iblk)%p, nsd, 'gather ' ) 36c 37c.... form the Lhs 38c 39 call e3StsLhs( xl, lStsVec ) 40c 41c.... assemble 42c 43 call local (stsVec, lStsVec, mien(iblk)%p, 44 & nResDims, 'scatter ' ) 45 46 deallocate ( xl ) 47 deallocate ( lStsVec ) 48c 49c.... end loop over element blocks 50c 51 enddo 52 53 if (numpe > 1) then 54 call commu (stsVec, ilwork, nResDims , 'in ') 55 endif 56c 57c.... local periodic boundary conditions (no communications) 58c 59 do j = 1,nshg 60 if (btest(iBC(j),10)) then 61 i = iper(j) 62 stsVec(i,:) = stsVec(i,:) + stsVec(j,:) 63 endif 64 enddo 65c 66 do i = 1,nshg 67 stsVec(i,:) = stsVec(iper(i),:) 68 enddo 69 if (numpe > 1) then 70 call commu (stsVec, ilwork, nResDims , 'out') 71 endif 72 73 return 74 end 75 76c----------------------------------------------------------------------- 77c Assemble the residual for the statistics 78c----------------------------------------------------------------------- 79 subroutine elmStatsRes( y, ac, x, shp, shgl, 80 & shpb, shglb, iBC, BC, 81 & iper, ilwork, rowp, colm, 82 & lhsK, lhsP ) 83 84 use stats 85 include "common.h" 86 87 88 real*8 y(nshg,ndof), ac(nshg,ndof), x(numnp,nsd), 89 & shp(MAXTOP,maxsh,MAXQPT), shgl(MAXTOP,nsd,maxsh,MAXQPT), 90 & shpb(MAXTOP,maxsh,MAXQPT), 91 & shglb(MAXTOP,nsd,maxsh,MAXQPT), 92 & BC(nshg,ndofBC), lhsK(9,nnz_tot), 93 & lhsP(4,nnz_tot), res(nshg,ndof) 94 95 integer iBC(nshg), iper(nshg), 96 & ilwork(nlwork), rowp(nshg,nnz), 97 & colm(nshg+1) 98 99 100 lhs = 0 101 stsVec = zero 102 103 stsResFlg = 1 104 ierrcalctmp=ierrcalc ! save the current value of ierrcalc 105 ierrcalc=0 ! set it to zero so that we don't calculate error 106 ! here (which would overflow memory around rjunk) 107 call ElmGMR (y, ac, x, 108 & shp, shgl, iBC, 109 & BC, shpb, shglb, 110 & res, iper, ilwork, 111 & rowp, colm, lhsK, 112 & lhsP, rjunk ) 113 stsResFlg = 0 114 ierrcalc=ierrcalctmp ! reset it back to original value 115 return 116 end 117 118