1 subroutine e3q (yl, dwl, shp, shgl, 2 & xl, ql, rmassl, 3 & xmudmi, sgn ) 4c 5c---------------------------------------------------------------------- 6c 7c This routine computes the element contribution to the 8c diffusive flux vector and the lumped mass matrix. 9c 10c input: 11c yl (npro,nshl,ndof) : Y variables 12c shp (nen,ngauss) : element shape-functions 13c shgl (nsd,nen,ngauss) : element local-grad-shape-functions 14c xl (npro,nshl,nsd) : nodal coordinates at current step 15c 16c output: 17c ql (npro,nshl,idflx) : element RHS diffusion residual 18c rmassl (npro,nshl) : element lumped mass matrix 19c 20c---------------------------------------------------------------------- 21c 22 include "common.h" 23c 24 dimension yl(npro,nshl,ndof), dwl(npro,nenl), 25 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 26 & xl(npro,nenl,nsd), 27 & ql(npro,nshl,idflx), rmassl(npro,nshl), 28 & xmudmi(npro,ngauss) 29c 30c local arrays 31c 32 dimension g1yi(npro,nflow), g2yi(npro,nflow), 33 & g3yi(npro,nflow), shg(npro,nshl,nsd), 34 & dxidx(npro,nsd,nsd), WdetJ(npro), 35 & rmu(npro) 36c 37 dimension qdi(npro,idflx),alph1(npro),alph2(npro) 38c 39 dimension sgn(npro,nshl), shape(npro,nshl), 40 & shdrv(npro,nsd,nshl), shpsum(npro) 41 42 real*8 tmp(npro) 43c 44c.... for surface tension 45c 46 dimension g1yti(npro), g2yti(npro), 47 & g3yti(npro) 48 integer idflow 49c 50c.... loop through the integration points 51c 52 53 54 alph1 = 0.d0 55 alph2 = 0.d0 56 57 do intp = 1, ngauss 58 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 59c 60 call getshp(shp, shgl, sgn, 61 & shape, shdrv) 62 63c 64c.... initialize 65c 66 qdi = zero 67c 68c 69c.... calculate the integration variables necessary for the 70c formation of q 71c 72 73 call e3qvar (yl, shdrv, 74 & xl, g1yi, 75 & g2yi, g3yi, shg, 76 & dxidx, WdetJ ) 77c 78 idflow = 9 ! we ALWAYS save space for tau_{ij} in q_i 79 ! even if idiff is not greater than 1 80 81 if(idiff >= 1) then !so taking care of all the idiff=1,3 82c 83c.... compute diffusive fluxes 84c 85c.... compute the viscosity 86c 87 call getdiff(dwl, yl, shape, xmudmi, xl, rmu, tmp) 88c 89c.... diffusive flux in x1-direction 90c 91 qdi(:,1) = two * rmu * g1yi(:,2) 92 qdi(:,4) = rmu * (g1yi(:,3) + g2yi(:,2)) 93 qdi(:,7) = rmu * (g1yi(:,4) + g3yi(:,2)) 94c 95c.... diffusive flux in x2-direction 96c 97 qdi(:,2) = rmu * (g1yi(:,3) + g2yi(:,2)) 98 qdi(:,5) = two * rmu * g2yi(:,3) 99 qdi(:,8) = rmu * (g2yi(:,4) + g3yi(:,3)) 100c 101c.... diffusive flux in x3-direction 102c 103 qdi(:,3) = rmu * (g1yi(:,4) + g3yi(:,2)) 104 qdi(:,6)= rmu * (g2yi(:,4) + g3yi(:,3)) 105 qdi(:,9)= two * rmu * g3yi(:,4) 106c 107c 108c.... assemble contribution of qdi to ql,i.e., contribution to 109c each element node 110c 111 do i=1,nshl 112 ql(:,i,1 ) = ql(:,i,1 )+ shape(:,i)*WdetJ*qdi(:,1 ) 113 ql(:,i,2 ) = ql(:,i,2 )+ shape(:,i)*WdetJ*qdi(:,2 ) 114 ql(:,i,3 ) = ql(:,i,3 )+ shape(:,i)*WdetJ*qdi(:,3 ) 115 116 ql(:,i,4 ) = ql(:,i,4 )+ shape(:,i)*WdetJ*qdi(:,4 ) 117 ql(:,i,5 ) = ql(:,i,5 )+ shape(:,i)*WdetJ*qdi(:,5 ) 118 ql(:,i,6 ) = ql(:,i,6 )+ shape(:,i)*WdetJ*qdi(:,6 ) 119 120 ql(:,i,7 ) = ql(:,i,7 )+ shape(:,i)*WdetJ*qdi(:,7 ) 121 ql(:,i,8 ) = ql(:,i,8 )+ shape(:,i)*WdetJ*qdi(:,8 ) 122 ql(:,i,9 ) = ql(:,i,9 )+ shape(:,i)*WdetJ*qdi(:,9 ) 123 124 enddo 125c 126c.... compute and assemble the element contribution to the lumped 127c mass matrix 128c 129c 130c.... row sum technique 131c 132 if ( idiff == 1 ) then 133 do i=1,nshl 134 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 135 enddo 136 endif 137c 138c.... "special lumping technique" (Hughes p. 445) 139c 140 if ( idiff == 3 ) then 141 shpsum = zero 142 do i=1,nshl 143 shpsum = shpsum + shape(:,i)*shape(:,i) 144 rmassl(:,i)=rmassl(:,i)+shape(:,i)*shape(:,i)*WdetJ 145 enddo 146 alph1 = alph1+WdetJ 147 alph2 = alph2+shpsum*WdetJ 148 endif 149 endif ! end of idiff=1 .or. 3 150c 151 if(isurf .eq. 1) then 152c 153c.... initialize 154c 155 g1yti = zero 156 g2yti = zero 157 g3yti = zero 158c 159c.... calculate the integration variables necessary for the 160c formation of q 161c 162c.... compute the global gradient of Yt-variables, assuming 6th entry as 163c.... the phase indicator function 164c 165c Yt_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(int) Yta) 166c 167 do n = 1, nshl 168 g1yti(:) = g1yti(:) + shg(:,n,1) * yl(:,n,6) 169 g2yti(:) = g2yti(:) + shg(:,n,2) * yl(:,n,6) 170 g3yti(:) = g3yti(:) + shg(:,n,3) * yl(:,n,6) 171 enddo 172c 173c computing N_{b}*N_{a,x_i)*yta*WdetJ 174c 175 do i=1,nshl 176 ql(:,i,idflow+1) = ql(:,i,idflow+1) 177 & + shape(:,i)*WdetJ*g1yti 178 ql(:,i,idflow+2) = ql(:,i,idflow+2) 179 & + shape(:,i)*WdetJ*g2yti 180 ql(:,i,idflow+3) = ql(:,i,idflow+3) 181 & + shape(:,i)*WdetJ*g3yti 182 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 183 enddo 184 endif !end of the isurf 185c 186c.... end of the loop over integration points 187c 188 enddo 189c 190c.... normalize the mass matrix for idiff == 3 191c 192 if ( idiff == 3 ) then 193 do i=1,nshl 194 rmassl(:,i) = rmassl(:,i)*alph1/alph2 195 enddo 196 endif 197 198 199c 200c.... return 201c 202 return 203 end 204 205 206 subroutine e3qSclr (yl, dwl, shp, shgl, 207 & xl, ql, rmassl, 208 & sgn ) 209c 210c---------------------------------------------------------------------- 211c 212c This routine computes the element contribution to the 213c diffusive flux vector and the lumped mass matrix. 214c 215c---------------------------------------------------------------------- 216c 217 include "common.h" 218c 219 dimension yl(npro,nshl,ndof), dwl(npro,nshl), 220 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 221 & xl(npro,nenl,nsd), 222 & ql(npro,nshl,nsd), rmassl(npro,nshl) 223c 224c local arrays 225c 226 dimension gradT(npro,nsd), 227 & dxidx(npro,nsd,nsd), WdetJ(npro) 228c 229 dimension qdi(npro,nsd),alph1(npro),alph2(npro) 230c 231 dimension sgn(npro,nshl), shape(npro,nshl), 232 & shdrv(npro,nsd,nshl), shpsum(npro) 233 234 real*8 diffus(npro) 235c 236c.... loop through the integration points 237c 238 239 240 alph1 = 0.d0 241 alph2 = 0.d0 242 243 do intp = 1, ngauss 244 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 245c 246 call getshp(shp, shgl, sgn, 247 & shape, shdrv) 248 249c 250c.... initialize 251c 252 qdi = zero 253c 254c 255c.... calculate the integration variables necessary for the 256c formation of q 257c 258 call e3qvarSclr (yl, shdrv, 259 & xl, gradT, 260 & dxidx, WdetJ ) 261 262c 263c.... compute diffusive flux vector at this integration point 264c 265 call getdiffsclr(shape, dwl, yl, diffus) 266 267c 268c.... diffusive flux 269c 270 qdi(:,1) = diffus * gradT(:,1) 271 qdi(:,2) = diffus * gradT(:,2) 272 qdi(:,3) = diffus * gradT(:,3) 273c 274c 275c.... assemble contribution of qdi to ql,i.e., contribution to 276c each element node 277c 278 do i=1,nshl 279 ql(:,i,1 ) = ql(:,i,1 )+ shape(:,i)*WdetJ*qdi(:,1 ) 280 ql(:,i,2 ) = ql(:,i,2 )+ shape(:,i)*WdetJ*qdi(:,2 ) 281 ql(:,i,3 ) = ql(:,i,3 )+ shape(:,i)*WdetJ*qdi(:,3 ) 282 283 enddo 284c 285c.... compute and assemble the element contribution to the lumped 286c mass matrix 287c 288c 289c.... row sum technique 290c 291 if ( idiff == 1 ) then 292 do i=1,nshl 293 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 294 enddo 295 endif 296c 297c.... "special lumping technique" (Hughes p. 445) 298c 299 if ( idiff == 3 ) then 300 shpsum = zero 301 do i=1,nshl 302 shpsum = shpsum + shape(:,i)*shape(:,i) 303 rmassl(:,i)=rmassl(:,i)+shape(:,i)*shape(:,i)*WdetJ 304 enddo 305 alph1 = alph1+WdetJ 306 alph2 = alph2+shpsum*WdetJ 307 endif 308c 309c.... end of the loop over integration points 310c 311 enddo 312c 313c.... normalize the mass matrix for idiff == 3 314c 315 if ( idiff == 3 ) then 316 do i=1,nshl 317 rmassl(:,i) = rmassl(:,i)*alph1/alph2 318 enddo 319 endif 320c 321c.... return 322c 323 return 324 end 325 326