1 subroutine genlmass (x, shp,shgl) 2c 3 use pointer_data 4c 5 include "common.h" 6 include "mpif.h" 7c 8 real*8 x(numnp,nsd) 9c 10 real*8 shp(MAXTOP,maxsh,MAXQPT), 11 & shgl(MAXTOP,nsd,maxsh,MAXQPT) 12c 13 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 14 15c 16c gmass came in via pointer_data and will 17c be available wherever it is included (allocate it now). 18c 19 20 allocate (gmass(nshg)) 21 gmass=zero 22c 23c.... loop over the element-blocks 24c 25 do iblk = 1, nelblk 26 iel = lcblk(1,iblk) 27 lelCat = lcblk(2,iblk) 28 lcsyst = lcblk(3,iblk) 29 iorder = lcblk(4,iblk) 30 nenl = lcblk(5,iblk) ! no. of vertices per element 31 nshl = lcblk(10,iblk) 32 mattyp = lcblk(7,iblk) 33 npro = lcblk(1,iblk+1) - iel 34 inum = iel + npro - 1 35 ngauss = nint(lcsyst) 36c 37c 38c.... compute and assemble the residual and tangent matrix 39c 40 allocate (tmpshp(nshl,MAXQPT)) 41 allocate (tmpshgl(nsd,nshl,MAXQPT)) 42 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 43 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 44 45 46 call AsImass (x, tmpshp, 47 & tmpshgl, mien(iblk)%p, 48 & gmass) 49 50 deallocate ( tmpshp ) 51 deallocate ( tmpshgl ) 52c 53c.... end of interior element loop 54c 55 enddo 56 57 return 58 end 59 60 subroutine AsImass (x, shp, 61 & shgl, ien, 62 & gmass) 63c 64c---------------------------------------------------------------------- 65c 66c This routine computes and assembles the mass corresponding to the 67c each node. 68c 69c Ken Jansen, Winter 2000. (Fortran 90) 70c---------------------------------------------------------------------- 71c 72 include "common.h" 73c 74 real*8 x(numnp,nsd), 75 & shp(nshl,maxsh), shgl(nsd,nshl,ngauss), 76 & gmass(nshg) 77 78 integer ien(npro,nshl) 79 80c 81 real*8 xl(npro,nenl,nsd), WdetJ(npro), 82 & sgn(npro,nshl), shape(npro,nshl), 83 & locmass(npro,nshl), shg(npro,nshl,nsd), 84 & fmstot(npro), temp(npro), 85 & dxidx(npro,nsd,nsd), shdrv(npro,nsd,nshl) 86 87 integer aa 88c 89c 90c 91c.... gather the variables 92c 93c 94c.... get the matrix of mode signs for the hierarchic basis functions. 95c 96 if (ipord .gt. 1) then 97 call getsgn(ien,sgn) 98 endif 99 100 call localx(x, xl, ien, nsd, 'gather ') 101c 102c.... zero the matrices if they are being recalculated 103c 104 105 locmass=zero 106 fmstot=zero 107 108 do intp = 1, ngauss 109 110 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 111c 112c.... get the hierarchic shape functions at this int point 113c 114 call getshp(shp, shgl, sgn, 115 & shape, shdrv) 116 117c 118c.... ---------------------> Element Metrics <----------------------- 119c 120 call e3metric( xl, shdrv, dxidx, 121 & shg, WdetJ) 122 123c 124c get this quad points contribution to the integral of the square of the 125c shape function 126c 127 do aa = 1,nshl 128 locmass(:,aa)= locmass(:,aa) 129 & + shape(:,aa)*shape(:,aa)*WdetJ 130 enddo 131c 132c also accumulate this quad points contribution to the integral of the element 133c volume (integral Na^2 d Omega) 134c 135 fmstot= fmstot + WdetJ ! intregral d Omega 136c 137c.... end of integration loop 138c 139 enddo 140c 141c.... lumped mass if needed Note that the locmass factors accumulated 142c over integration points and weighted with WdetJ already. 143c 144 145c.... scale the LHS matrix contribution with special lumping weighting 146c 147c The first term we collect is the trace of integral Na^2 d Omega 148c 149 temp = zero 150 do aa = 1, nshl 151 temp = temp + locmass(:,aa) !reusing temp to save memory 152 enddo 153 154c 155c scale the diagonal so that the trace will still yield Omega^e (the volume 156c of the element) 157c 158 do aa = 1, nshl 159 locmass(:,aa) = locmass(:,aa) * fmstot / temp 160 enddo 161c 162c.... assemble the residual 163c 164 call local (gmass, locmass, ien, 1, 'scatter ') 165 166c 167c.... end 168c 169 return 170 end 171 172 subroutine lmassadd ( ac, res, 173 & rowp, colm, 174 & lhsK, gmass) 175c 176 include "common.h" 177c 178 real*8 ac(nshg,ndof), res(nshg,4), tmp,tmp1 179 real*8 lhsK(9,nnz_tot), gmass(nshg), rho(nshg) 180 integer rowp(nnz*nshg), colm(nshg+1) 181 integer n, k 182c 183 integer sparseloc 184c 185c 186 rho=datmat(1,1,1) ! needs to be generalized for VOF or level set 187 tmp1=flmpl*almi 188 if((flmpl.ne.0).and.(lhs.eq.1)) then 189c 190c.... Add lmass to diag of lhsK 191c 192 do n = 1, nshg 193 k = sparseloc( rowp(colm(n)), colm(n+1)-colm(n), n ) 194 & + colm(n)-1 195 tmp=gmass(n)*tmp1*rho(n) 196 lhsK(1,k) = lhsK(1,k) + tmp 197 lhsK(5,k) = lhsK(5,k) + tmp 198 lhsK(9,k) = lhsK(9,k) + tmp 199 enddo 200 endif 201 202 tmp1=flmpr 203 204 if(flmpr.ne.0) then 205 rho=rho*gmass*tmp1 ! reuse rho 206 res(:,1)=res(:,1)-ac(:,1)*rho(:) 207 res(:,2)=res(:,2)-ac(:,2)*rho(:) 208 res(:,3)=res(:,3)-ac(:,3)*rho(:) 209 endif 210 211 return 212 end 213 214 subroutine lmassaddSclr ( ac, res, 215 & rowp, colm, 216 & lhsS, gmass) 217c 218 include "common.h" 219c 220 real*8 ac(nshg), res(nshg), tmp, tmp1 221 real*8 lhsS(nnz_tot), gmass(nshg), rho(nshg) 222 integer rowp(nnz*nshg), colm(nshg+1) 223 integer n, k 224c 225 integer sparseloc 226c 227c 228 rho=datmat(1,1,1) ! needs to be generalized for VOF or level set 229 tmp1=flmpl*almi 230 if((flmpl.ne.0).and.(lhs.eq.1)) then 231c 232c.... Add lmass to diag of lhsK 233c 234 do n = 1, nshg 235 k = sparseloc( rowp(colm(n)), colm(n+1)-colm(n), n ) 236 & + colm(n)-1 237 tmp=gmass(n)*tmp1*rho(n) 238 lhsS(k) = lhsS(k) + tmp 239 enddo 240 endif 241 242 tmp1=flmpr 243 if(flmpr.ne.0) then 244 rho=rho*gmass*tmp1 ! reuse rho 245 res(:)=res(:)-ac(:)*rho(:) 246 endif 247 248 return 249 end 250