1 subroutine e3qvar (yl, shgl, 2 & xl, g1yi, 3 & g2yi, g3yi, shg, 4 & dxidx, WdetJ ) 5c 6c---------------------------------------------------------------------- 7c 8c This routine computes the variables at integration point 9c necessary for the computation of the diffusive flux vector. 10c 11c input: 12c yl (npro,nshl,ndof) : primitive variables 13c shgl (npro,nsd,nshl) : element local-grad-shape-functions 14c xl (npro,nenl,nsd) : nodal coordinates at current step 15c 16c output: 17c g1yi (npro,ndof) : grad-y in direction 1 18c g2yi (npro,ndof) : grad-y in direction 2 19c g3yi (npro,ndof) : grad-y in direction 3 20c shg (npro,nshl,nsd) : element global grad-shape-functions 21c dxidx (npro,nsd,nsd) : inverse of deformation gradient 22c WdetJ (npro) : weighted Jacobian 23c u1 (npro) : x1-velocity component 24c u2 (npro) : x2-velocity component 25c u3 (npro) : x3-velocity component 26c 27c 28c Zdenek Johan, Summer 1990. (Modified from e2ivar.f) 29c Zdenek Johan, Winter 1991. (Fortran 90) 30c Kenneth Jansen, Winter 1997. Primitive Variables 31c---------------------------------------------------------------------- 32c 33 include "common.h" 34c 35c passed arrays 36c 37 dimension yl(npro,nshl,ndof), 38 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 39 & g1yi(npro,nflow), g2yi(npro,nflow), 40 & g3yi(npro,nflow), shg(npro,nshl,nsd), 41 & dxidx(npro,nsd,nsd), WdetJ(npro) 42c 43c local arrays 44c 45 dimension tmp(npro), dxdxi(npro,nsd,nsd) 46 47c 48c.... compute the deformation gradient 49c 50 dxdxi = zero 51c 52 do n = 1, nenl 53 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(:,1,n) 54 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(:,2,n) 55 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(:,3,n) 56 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(:,1,n) 57 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(:,2,n) 58 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(:,3,n) 59 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(:,1,n) 60 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(:,2,n) 61 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(:,3,n) 62 enddo 63c 64c.... compute the inverse of deformation gradient 65c 66 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 67 & - dxdxi(:,3,2) * dxdxi(:,2,3) 68 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 69 & - dxdxi(:,1,2) * dxdxi(:,3,3) 70 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 71 & - dxdxi(:,1,3) * dxdxi(:,2,2) 72 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 73 & + dxidx(:,1,2) * dxdxi(:,2,1) 74 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 75 dxidx(:,1,1) = dxidx(:,1,1) * tmp 76 dxidx(:,1,2) = dxidx(:,1,2) * tmp 77 dxidx(:,1,3) = dxidx(:,1,3) * tmp 78 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 79 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 80 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 81 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 82 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 83 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 84 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 85 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 86 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 87 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 88 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 89 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 90c 91 WdetJ = Qwt(lcsyst,intp)/ tmp 92 93c 94c.... ---------------------> Global Gradients <----------------------- 95c 96 g1yi = zero 97 g2yi = zero 98 g3yi = zero 99c 100c 101 do n = 1, nshl 102c 103c.... compute the global gradient of shape-function 104c 105c ! N_{a,x_i}= N_{a,xi_i} xi_{i,x_j} 106c 107 shg(:,n,1) = shgl(:,1,n) * dxidx(:,1,1) + 108 & shgl(:,2,n) * dxidx(:,2,1) + 109 & shgl(:,3,n) * dxidx(:,3,1) 110 shg(:,n,2) = shgl(:,1,n) * dxidx(:,1,2) + 111 & shgl(:,2,n) * dxidx(:,2,2) + 112 & shgl(:,3,n) * dxidx(:,3,2) 113 shg(:,n,3) = shgl(:,1,n) * dxidx(:,1,3) + 114 & shgl(:,2,n) * dxidx(:,2,3) + 115 & shgl(:,3,n) * dxidx(:,3,3) 116c 117c.... compute the global gradient of Y-variables 118c 119c 120c Y_{,x_i}=SUM_{a=1}^nenl (N_{a,x_i}(int) Ya) 121c 122 g1yi(:,2) = g1yi(:,2) + shg(:,n,1) * yl(:,n,2) 123 g1yi(:,3) = g1yi(:,3) + shg(:,n,1) * yl(:,n,3) 124 g1yi(:,4) = g1yi(:,4) + shg(:,n,1) * yl(:,n,4) 125c 126 g2yi(:,2) = g2yi(:,2) + shg(:,n,2) * yl(:,n,2) 127 g2yi(:,3) = g2yi(:,3) + shg(:,n,2) * yl(:,n,3) 128 g2yi(:,4) = g2yi(:,4) + shg(:,n,2) * yl(:,n,4) 129c 130 g3yi(:,2) = g3yi(:,2) + shg(:,n,3) * yl(:,n,2) 131 g3yi(:,3) = g3yi(:,3) + shg(:,n,3) * yl(:,n,3) 132 g3yi(:,4) = g3yi(:,4) + shg(:,n,3) * yl(:,n,4) 133 134 enddo 135 136c 137c.... return 138c 139 140 return 141 end 142 143c----------------------------------------------------------------------- 144c 145c compute the variables for the local scalar diffusion 146c 147c----------------------------------------------------------------------- 148 subroutine e3qvarSclr (yl, shgl, xl, 149 & gradT, dxidx, WdetJ ) 150 151 include "common.h" 152c 153c passed arrays 154c 155 real*8 yl(npro,nshl,ndof), shp(npro,nshl), 156 & shgl(npro,nsd,nshl), xl(npro,nenl,nsd), 157 & dxidx(npro,nsd,nsd), WdetJ(npro), 158 & gradT(npro,nsd) 159c 160c local arrays 161c 162 real*8 shg(npro,nshl,nsd) 163 164 165 call e3metric( xl, shgl, dxidx, 166 & shg, WdetJ) 167 168 gradT = zero 169 id=5+isclr 170c 171c later, when there are more models than SA we will need a 172c more general function to calculate evisc at a quadrature point 173c 174 do n = 1, nshl 175 gradT(:,1) = gradT(:,1) + shg(:,n,1) * yl(:,n,id) 176 gradT(:,2) = gradT(:,2) + shg(:,n,2) * yl(:,n,id) 177 gradT(:,3) = gradT(:,3) + shg(:,n,3) * yl(:,n,id) 178 enddo 179c 180c.... return 181c 182 183 return 184 end 185 186 187 188 189