1 subroutine getstrl( y, x, ien, strnrm, shgl, shp ) 2 3 include "common.h" 4 5 dimension y(nshg,5) 6 dimension x(numnp,3), xl(npro,nenl,3) 7 dimension ien(npro,nshl), yl(npro,nshl,5), 8 & u1(npro), u2(npro), 9 & u3(npro), dxdxi(npro,nsd,nsd), 10 & strnrm(npro,maxsh), dxidx(npro,nsd,nsd), 11 & shgl(nsd,nshl,maxsh), shg(npro,nshl,nsd), 12 & shp(nshl,maxsh) 13 dimension tmp(npro), fresli(npro,24) 14 15 call localy (y, yl, ien, 5, 'gather ') 16 call localx (x, xl, ien, 3, 'gather ') 17c 18 19 if(matflg(1,1).eq.0) then ! compressible 20 yl (:,:,1) = yl(:,:,1) / (Rgas * yl(:,:,5)) 21 else 22 yl (:,:,1) = one 23 endif 24 25 do intp = 1, ngauss 26 27c calculate the metrics 28c 29c 30c.... ---------------------> Element Metrics <----------------------- 31c 32c.... compute the deformation gradient 33c 34 dxdxi = zero 35c 36 do n = 1, nenl 37 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(1,n,intp) 38 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(2,n,intp) 39 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(3,n,intp) 40 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(1,n,intp) 41 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(2,n,intp) 42 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(3,n,intp) 43 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(1,n,intp) 44 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(2,n,intp) 45 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(3,n,intp) 46 enddo 47c 48c.... compute the inverse of deformation gradient 49c 50 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 51 & - dxdxi(:,3,2) * dxdxi(:,2,3) 52 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 53 & - dxdxi(:,1,2) * dxdxi(:,3,3) 54 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 55 & - dxdxi(:,1,3) * dxdxi(:,2,2) 56 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 57 & + dxidx(:,1,2) * dxdxi(:,2,1) 58 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 59 dxidx(:,1,1) = dxidx(:,1,1) * tmp 60 dxidx(:,1,2) = dxidx(:,1,2) * tmp 61 dxidx(:,1,3) = dxidx(:,1,3) * tmp 62 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 63 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 64 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 65 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 66 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 67 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 68 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 69 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 70 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 71 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 72 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 73 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 74c 75 76 fresli=zero 77 do i=1,nshl 78 fresli(:,22) = fresli(:,22)+shp(i,intp)*yl(:,i,1) ! density at qpt 79c fresli(:,24) = fresli(:,24)+shp(i,intp)*yl(:,i,5) !temperature at qpt 80 enddo 81c 82c 83c fresli(:,22)=fresli(:,22)*wght 84c fresli(:,24)=fresli(:,24)*wght 85 86 87 do n = 1,nshl 88 shg(:,n,1) = (shgl(1,n,intp) * dxidx(:,1,1) 89 & + shgl(2,n,intp) * dxidx(:,2,1) 90 & + shgl(3,n,intp) * dxidx(:,3,1)) 91 shg(:,n,2) = (shgl(1,n,intp) * dxidx(:,1,2) 92 & + shgl(2,n,intp) * dxidx(:,2,2) 93 & + shgl(3,n,intp) * dxidx(:,3,2)) 94 shg(:,n,3) = (shgl(1,n,intp) * dxidx(:,1,3) 95 & + shgl(2,n,intp) * dxidx(:,2,3) 96 & + shgl(3,n,intp) * dxidx(:,3,3)) 97 enddo 98 99 do j=10,12 ! normal strainrate u_{i,i} no sum on i 100 ig=j-9 101 iv=j-8 102 do i=1,nshl 103 fresli(:,j) = fresli(:,j)+shg(:,i,ig)*yl(:,i,iv) 104 enddo 105 enddo 106 107c shear stresses NOTE there may be faster ways to do this 108c check agains CM5 code for speed WTP 109 110 do i=1,nshl 111 fresli(:,13) = fresli(:,13)+shg(:,i,2)*yl(:,i,2) 112 & +shg(:,i,1)*yl(:,i,3) 113 fresli(:,14) = fresli(:,14)+shg(:,i,3)*yl(:,i,2) 114 & +shg(:,i,1)*yl(:,i,4) 115 fresli(:,15) = fresli(:,15)+shg(:,i,3)*yl(:,i,3) 116 & +shg(:,i,2)*yl(:,i,4) 117 enddo 118 119 120 fresli(:,13) = pt5 * fresli(:,13) 121 fresli(:,14) = pt5 * fresli(:,14) 122 fresli(:,15) = pt5 * fresli(:,15) 123 124 strnrm(:,intp) = fresli(:,22) * sqrt( 125 & two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2) 126 & + four * ( fresli(:,13)**2 + fresli(:,14)**2 + 127 & fresli(:,15)**2 ) ) 128 129 130 enddo !end of loop over integration points 131 132 133 return 134 end 135