1 subroutine e3qGradV (yl, shp, shgl, 2 & xl, ql, 3 & 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,nsd*nsd) : element RHS diffusion residual 18c 19c---------------------------------------------------------------------- 20c 21 include "common.h" 22c 23 dimension yl(npro,nshl,ndof), 24 & shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 25 & xl(npro,nenl,nsd), 26 & ql(npro,nshl,nsdsq) 27c 28c local arrays 29c 30 dimension g1yi(npro,nflow), g2yi(npro,nflow), 31 & g3yi(npro,nflow), shg(npro,nshl,nsd), 32 & dxidx(npro,nsd,nsd), WdetJ(npro) 33c 34c 35 dimension sgn(npro,nshl), shape(npro,nshl), 36 & shdrv(npro,nsd,nshl) 37 38c 39c.... loop through the integration points 40c 41 42 do intp = 1, ngauss 43 if (Qwt(lcsyst,intp) .eq. zero) cycle ! precaution 44c 45 call getshp(shp, shgl, sgn, 46 & shape, shdrv) 47 48c 49c 50c.... calculate the integration variables necessary for the 51c formation of q 52c 53 54 call e3qvar (yl, shdrv, 55 & xl, g1yi, 56 & g2yi, g3yi, shg, 57 & dxidx, WdetJ ) 58c 59c 60c each element node 61c 62 do i=1,nshl 63 ql(:,i,1 ) = ql(:,i,1 )+ shape(:,i)*WdetJ*g1yi(:,2 ) ! du/dx 64 ql(:,i,2 ) = ql(:,i,2 )+ shape(:,i)*WdetJ*g2yi(:,2 ) ! du/dy 65 ql(:,i,3 ) = ql(:,i,3 )+ shape(:,i)*WdetJ*g3yi(:,2 ) ! du/dz 66 67 ql(:,i,4 ) = ql(:,i,4 )+ shape(:,i)*WdetJ*g1yi(:,3 ) ! dv/dx 68 ql(:,i,5 ) = ql(:,i,5 )+ shape(:,i)*WdetJ*g2yi(:,3 ) ! dv/dy 69 ql(:,i,6 ) = ql(:,i,6 )+ shape(:,i)*WdetJ*g3yi(:,3 ) ! dv/dz 70 71 ql(:,i,7 ) = ql(:,i,7 )+ shape(:,i)*WdetJ*g1yi(:,4 ) ! dw/dx 72 ql(:,i,8 ) = ql(:,i,8 )+ shape(:,i)*WdetJ*g2yi(:,4 ) ! dw/dy 73 ql(:,i,9 ) = ql(:,i,9 )+ shape(:,i)*WdetJ*g3yi(:,4 ) ! dw/dz 74 75 enddo 76c 77c 78c.... end of the loop over integration points 79c 80 enddo 81 82 83c 84c.... return 85c 86 return 87 end 88