1 subroutine e3q (ycl, 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 ycl (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 ycl(npro,nshl,ndof), 25 & shp(nshl,ngauss), 26 & shgl(nsd,nshl,ngauss), 27 & xl(npro,nenl,nsd), 28 & ql(npro,nshl,idflx), rmassl(npro,nshl), 29 & xmudmi(npro,ngauss) 30c 31c local arrays 32c 33 dimension g1yi(npro,nflow), g2yi(npro,nflow), 34 & g3yi(npro,nflow), shg(npro,nshl,nsd), 35 & dxidx(npro,nsd,nsd), WdetJ(npro), 36 & T(npro), cp(npro), 37 & u1(npro), u2(npro), 38 & u3(npro), rmu(npro), 39 & rlm(npro), rlm2mu(npro), 40 & con(npro), rho(npro) 41c 42 dimension qdi(npro,idflx),alph1(npro),alph2(npro) 43c 44 dimension sgn(npro,nshl), shape(npro,nshl), 45 & shdrv(npro,nsd,nshl), shpsum(npro) 46c 47c.... for surface tension 48c 49 dimension g1yti(npro), g2yti(npro), 50 & g3yti(npro) 51 integer idflow 52c 53c.... loop through the integration points 54c 55 56 ttim(33) = ttim(33) - secs(0.0) 57 58 alph1 = 0.d0 59 alph2 = 0.d0 60 61 do intp = 1, ngauss 62 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 63c 64c.... create a matrix of shape functions (and derivatives) for each 65c element at this quadrature point. These arrays will contain 66c the correct signs for the hierarchic basis 67c 68 call getshp(shp, shgl, sgn, 69 & shape, shdrv) 70c 71c.... initialize 72c 73 qdi = zero 74c 75c 76c.... calculate the integration variables necessary for the 77c formation of q 78c 79 80 call e3qvar (ycl, shape, shdrv, 81 & rho, xl, g1yi, 82 & g2yi, g3yi, shg, 83 & dxidx, WdetJ, T, 84 & cp, u1, u2, 85 & u3 ) 86c 87c.... compute diffusive flux vector at this integration point 88c 89c 90c.... get material properties 91c 92 call getDiff (T, cp, rho, ycl, 93 & rmu, rlm, rlm2mu, con, shape, 94 & xmudmi, xl) 95 96 idflow = 0 97 if(idiff >= 1) then !so taking care of all the idiff=1,2,3 98 idflow = idflow+12 99c 100c.... compute diffusive fluxes 101c 102c.... diffusive flux in x1-direction 103c 104 qdi(:,1) = rlm2mu * g1yi(:,2) 105 & + rlm * g2yi(:,3) 106 & + rlm * g3yi(:,4) 107 qdi(:,2) = rmu * g1yi(:,3) 108 & + rmu * g2yi(:,2) 109 qdi(:,3) = rmu * g1yi(:,4) 110 & + rmu * g3yi(:,2) 111 qdi(:,4) = rlm2mu * u1 * g1yi(:,2) + rmu * u2 * g1yi(:,3) 112 & + rmu * u3 * g1yi(:,4) 113 & + rmu * u2 * g2yi(:,2) + rlm * u1 * g2yi(:,3) 114 & + rmu * u3 * g3yi(:,2) + rlm * u1 * g3yi(:,4) 115 & + con * g1yi(:,5) 116 117c 118c.... diffusive flux in x2-direction 119c 120 qdi(:,5) = rmu * g1yi(:,3) 121 & + rmu * g2yi(:,2) 122 qdi(:,6) = rlm * g1yi(:,2) 123 & + rlm2mu * g2yi(:,3) 124 & + rlm * g3yi(:,4) 125 qdi(:,7) = rmu * g2yi(:,4) 126 & + rmu * g3yi(:,3) 127 qdi(:,8) = rlm * u2 * g1yi(:,2) + rmu * u1 * g1yi(:,3) 128 & + rmu * u1 * g2yi(:,2) + rlm2mu * u2 * g2yi(:,3) 129 & + rmu * u3 * g2yi(:,4) 130 & + rmu * u3 * g3yi(:,3) + rlm * u2 * g3yi(:,4) 131 & + con * g2yi(:,5) 132c 133c.... diffusive flux in x3-direction 134c 135 qdi(:,9 ) = rmu * g1yi(:,4) 136 & + rmu * g3yi(:,2) 137 qdi(:,10) = rmu * g2yi(:,4) 138 & + rmu * g3yi(:,3) 139 qdi(:,11) = rlm * g1yi(:,2) 140 & + rlm * g2yi(:,3) 141 & + rlm2mu * g3yi(:,4) 142 qdi(:,12) = rlm * u3 * g1yi(:,2) + rmu * u1 * g1yi(:,4) 143 & + rlm * u3 * g2yi(:,3) + rmu * u2 * g2yi(:,4) 144 & + rmu * u1 * g3yi(:,2) + rmu * u2 * g3yi(:,3) 145 & + rlm2mu * u3 * g3yi(:,4) 146 & + con * g3yi(:,5) 147 148c 149c 150c.... assemble contribution of qdi to ql,i.e., contribution to 151c each element node 152c 153 do i=1,nshl 154 ql(:,i,1 ) = ql(:,i,1 )+ shape(:,i)*WdetJ*qdi(:,1 ) 155 ql(:,i,2 ) = ql(:,i,2 )+ shape(:,i)*WdetJ*qdi(:,2 ) 156 ql(:,i,3 ) = ql(:,i,3 )+ shape(:,i)*WdetJ*qdi(:,3 ) 157 ql(:,i,4 ) = ql(:,i,4 )+ shape(:,i)*WdetJ*qdi(:,4 ) 158 ql(:,i,5 ) = ql(:,i,5 )+ shape(:,i)*WdetJ*qdi(:,5 ) 159 ql(:,i,6 ) = ql(:,i,6 )+ shape(:,i)*WdetJ*qdi(:,6 ) 160 ql(:,i,7 ) = ql(:,i,7 )+ shape(:,i)*WdetJ*qdi(:,7 ) 161 ql(:,i,8 ) = ql(:,i,8 )+ shape(:,i)*WdetJ*qdi(:,8 ) 162 ql(:,i,9 ) = ql(:,i,9 )+ shape(:,i)*WdetJ*qdi(:,9 ) 163 ql(:,i,10) = ql(:,i,10)+ shape(:,i)*WdetJ*qdi(:,10) 164 ql(:,i,11) = ql(:,i,11)+ shape(:,i)*WdetJ*qdi(:,11) 165 ql(:,i,12) = ql(:,i,12)+ shape(:,i)*WdetJ*qdi(:,12) 166 enddo 167c 168c.... compute and assemble the element contribution to the lumped 169c mass matrix 170c 171c 172c.... row sum technique 173c 174 if ( idiff == 1 ) then 175 do i=1,nshl 176 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 177 enddo 178 endif 179c 180c.... "special lumping technique" (Hughes p. 445) 181c 182 if ( idiff == 3 ) then 183 shpsum = zero 184 do i=1,nshl 185 shpsum = shpsum + shape(:,i)*shape(:,i) 186 rmassl(:,i)=rmassl(:,i)+shape(:,i)*shape(:,i)*WdetJ 187 enddo 188 alph1 = alph1+WdetJ 189 alph2 = alph2+shpsum*WdetJ 190 endif 191 endif ! end of idiff=1 .or. 3 192c 193 if(isurf .eq. 1) then 194c 195c.... initialize 196c 197 g1yti = zero 198 g2yti = zero 199 g3yti = zero 200c 201c.... calculate the integration variables necessary for the 202c formation of q 203c 204c.... compute the global gradient of Yt-variables, assuming 6th entry as 205c.... the phase indicator function 206c 207c Yt_{,x_i}=SUM_{a=1}^nshl (N_{a,x_i}(int) Yta) 208c 209 do n = 1, nshl 210 g1yti(:) = g1yti(:) + shg(:,n,1) * ycl(:,n,6) 211 g2yti(:) = g2yti(:) + shg(:,n,2) * ycl(:,n,6) 212 g3yti(:) = g3yti(:) + shg(:,n,3) * ycl(:,n,6) 213 enddo 214c 215c computing N_{b}*N_{a,x_i)*yta*WdetJ 216c 217 do i=1,nshl 218 ql(:,i,idflow+1) = ql(:,i,idflow+1) 219 & + shape(:,i)*WdetJ*g1yti 220 ql(:,i,idflow+2) = ql(:,i,idflow+2) 221 & + shape(:,i)*WdetJ*g2yti 222 ql(:,i,idflow+3) = ql(:,i,idflow+3) 223 & + shape(:,i)*WdetJ*g3yti 224 rmassl(:,i) = rmassl(:,i) + shape(:,i)*WdetJ 225 enddo 226 endif !end of the isurf 227c 228c.... end of the loop over integration points 229c 230 enddo 231c 232c.... normalize the mass matrix for idiff == 3 233c 234 if ( idiff == 3 ) then 235 do i=1,nshl 236 rmassl(:,i) = rmassl(:,i)*alph1/alph2 237 enddo 238 endif 239 240 ttim(33) = ttim(33) + secs(0.0) 241 242c 243c.... return 244c 245 return 246 end 247