1 subroutine asithf (y, x, strnrm, ien, fres, shgl, shp, Qwtf) 2 3 include "common.h" 4 5 dimension y(nshg,ndof), fres(nshg,24) 6 dimension x(numnp,nsd), xl(npro,nenl,nsd) 7 dimension ien(npro,nshl), ycl(npro,nshl,ndof), 8 & fresl(npro,24), WdetJ(npro), 9 & u1(npro), u2(npro), 10 & u3(npro), dxdxi(npro,nsd,nsd), 11 & strnrm(npro,maxsh), dxidx(npro,nsd,nsd), 12 & shgl(nsd,nshl,maxsh), shg(npro,nshl,nsd), 13 & shp(nshl,maxsh), 14 & fresli(npro,24), Qwtf(ngaussf) 15 16 dimension tmp(npro) 17 18 call localy (y, ycl, ien, 5, 'gather ') 19 call localx (x, xl, ien, 3, 'gather ') 20c 21 22 if(matflg(1,1).eq.0) then ! compressible 23 ycl (:,:,1) = ycl(:,:,1) / (Rgas * ycl(:,:,5)) !get density 24 else 25 ycl(:,:,1) = one ! Even if density non unity, it would cancel out 26 endif 27 28 fresl = zero 29 30 do intp = 1, ngaussf 31 32 33c calculate the metrics 34c 35c 36c.... ---------------------> Element Metrics <----------------------- 37c 38c.... compute the deformation gradient 39c 40 dxdxi = zero 41c 42 do n = 1, nenl 43 dxdxi(:,1,1) = dxdxi(:,1,1) + xl(:,n,1) * shgl(1,n,intp) 44 dxdxi(:,1,2) = dxdxi(:,1,2) + xl(:,n,1) * shgl(2,n,intp) 45 dxdxi(:,1,3) = dxdxi(:,1,3) + xl(:,n,1) * shgl(3,n,intp) 46 dxdxi(:,2,1) = dxdxi(:,2,1) + xl(:,n,2) * shgl(1,n,intp) 47 dxdxi(:,2,2) = dxdxi(:,2,2) + xl(:,n,2) * shgl(2,n,intp) 48 dxdxi(:,2,3) = dxdxi(:,2,3) + xl(:,n,2) * shgl(3,n,intp) 49 dxdxi(:,3,1) = dxdxi(:,3,1) + xl(:,n,3) * shgl(1,n,intp) 50 dxdxi(:,3,2) = dxdxi(:,3,2) + xl(:,n,3) * shgl(2,n,intp) 51 dxdxi(:,3,3) = dxdxi(:,3,3) + xl(:,n,3) * shgl(3,n,intp) 52 enddo 53c 54c.... compute the inverse of deformation gradient 55c 56 dxidx(:,1,1) = dxdxi(:,2,2) * dxdxi(:,3,3) 57 & - dxdxi(:,3,2) * dxdxi(:,2,3) 58 dxidx(:,1,2) = dxdxi(:,3,2) * dxdxi(:,1,3) 59 & - dxdxi(:,1,2) * dxdxi(:,3,3) 60 dxidx(:,1,3) = dxdxi(:,1,2) * dxdxi(:,2,3) 61 & - dxdxi(:,1,3) * dxdxi(:,2,2) 62 tmp = one / ( dxidx(:,1,1) * dxdxi(:,1,1) 63 & + dxidx(:,1,2) * dxdxi(:,2,1) 64 & + dxidx(:,1,3) * dxdxi(:,3,1) ) 65 dxidx(:,1,1) = dxidx(:,1,1) * tmp 66 dxidx(:,1,2) = dxidx(:,1,2) * tmp 67 dxidx(:,1,3) = dxidx(:,1,3) * tmp 68 dxidx(:,2,1) = (dxdxi(:,2,3) * dxdxi(:,3,1) 69 & - dxdxi(:,2,1) * dxdxi(:,3,3)) * tmp 70 dxidx(:,2,2) = (dxdxi(:,1,1) * dxdxi(:,3,3) 71 & - dxdxi(:,3,1) * dxdxi(:,1,3)) * tmp 72 dxidx(:,2,3) = (dxdxi(:,2,1) * dxdxi(:,1,3) 73 & - dxdxi(:,1,1) * dxdxi(:,2,3)) * tmp 74 dxidx(:,3,1) = (dxdxi(:,2,1) * dxdxi(:,3,2) 75 & - dxdxi(:,2,2) * dxdxi(:,3,1)) * tmp 76 dxidx(:,3,2) = (dxdxi(:,3,1) * dxdxi(:,1,2) 77 & - dxdxi(:,1,1) * dxdxi(:,3,2)) * tmp 78 dxidx(:,3,3) = (dxdxi(:,1,1) * dxdxi(:,2,2) 79 & - dxdxi(:,1,2) * dxdxi(:,2,1)) * tmp 80c 81c wght=Qwt(lcsyst,intp) ! may be different now 82 wght=Qwtf(intp) 83 WdetJ = wght / tmp 84c 85 fresli=zero 86c 87 if(matflg(1,1).eq.0) then ! compressible 88 do i=1,nshl 89 fresli(:,22) = fresli(:,22)+shp(i,intp)*ycl(:,i,1) !density at 90 !qpt 91 enddo 92 else ! incompressible, set density 93 fresli(:,22)= one ! reduce comp2incompr regardless of rho datmat(1,1,1) 94 endif 95c 96 do n = 1,nshl 97 shg(:,n,1) = (shgl(1,n,intp) * dxidx(:,1,1) 98 & + shgl(2,n,intp) * dxidx(:,2,1) 99 & + shgl(3,n,intp) * dxidx(:,3,1)) 100 shg(:,n,2) = (shgl(1,n,intp) * dxidx(:,1,2) 101 & + shgl(2,n,intp) * dxidx(:,2,2) 102 & + shgl(3,n,intp) * dxidx(:,3,2)) 103 shg(:,n,3) = (shgl(1,n,intp) * dxidx(:,1,3) 104 & + shgl(2,n,intp) * dxidx(:,2,3) 105 & + shgl(3,n,intp) * dxidx(:,3,3)) 106 enddo 107 108 do j=10,12 ! normal strainrate u_{i,i} no sum on i 109 ig=j-9 110 iv=j-8 111 do i=1,nshl 112 fresli(:,j) = fresli(:,j)+shg(:,i,ig)*ycl(:,i,iv) 113 enddo 114 enddo 115 116c shear stresses NOTE there may be faster ways to do this 117c check agains CM5 code for speed WTP 118 119 do i=1,nshl 120 fresli(:,13) = fresli(:,13)+shg(:,i,2)*ycl(:,i,2) 121 & +shg(:,i,1)*ycl(:,i,3) 122 fresli(:,14) = fresli(:,14)+shg(:,i,3)*ycl(:,i,2) 123 & +shg(:,i,1)*ycl(:,i,4) 124 fresli(:,15) = fresli(:,15)+shg(:,i,3)*ycl(:,i,3) 125 & +shg(:,i,2)*ycl(:,i,4) 126 enddo 127 128 fresli(:,13) = pt5 * fresli(:,13) 129 fresli(:,14) = pt5 * fresli(:,14) 130 fresli(:,15) = pt5 * fresli(:,15) 131 132 strnrm(:,intp) = fresli(:,22) * sqrt( 133 & two * (fresli(:,10)**2 + fresli(:,11)**2 + fresli(:,12)**2) 134 & + four * ( fresli(:,13)**2 + fresli(:,14)**2 + 135 & fresli(:,15)**2 ) ) 136 137c 138c S_ij 139c 140 141 fresli(:,10) = fresli(:,10) * WdetJ ! u_{1,1}*WdetJ 142 fresli(:,11) = fresli(:,11) * WdetJ ! u_{2,2}*WdetJ 143 fresli(:,12) = fresli(:,12) * WdetJ ! u_{3,3}*WdetJ 144 fresli(:,13) = fresli(:,13) * WdetJ ! (1/2)*(u_{1,2}+u_{2,1})*WdetJ 145 fresli(:,14) = fresli(:,14) * WdetJ ! (1/2)*(u_{1,3}+u_{3,1})*WdetJ 146 fresli(:,15) = fresli(:,15) * WdetJ ! (1/2)*(u_{2,3}+u_{3,2})*WdetJ 147 148 fresli(:,22) = fresli(:,22) * WdetJ !rho * WdetJ 149c fresli(:,24) = fresli(:,24) * WdetJ 150 151 u1=zero 152 u2=zero 153 u3=zero 154 do i=1,nshl 155 u1 = u1 + shp(i,intp)*ycl(:,i,2) 156 u2 = u2 + shp(i,intp)*ycl(:,i,3) 157 u3 = u3 + shp(i,intp)*ycl(:,i,4) 158 enddo 159 160 fresli(:,1) = fresli(:,22) * u1 !rho u1 * WdetJ 161 fresli(:,2) = fresli(:,22) * u2 !rho u2 * WdetJ 162 fresli(:,3) = fresli(:,22) * u3 !rho u3 * WdetJ 163 164 fresli(:,4) = fresli(:,1) * u1 !rho u1 u1 *WdetJ 165 fresli(:,5) = fresli(:,2) * u2 !rho u2 u2 *WdetJ 166 fresli(:,6) = fresli(:,3) * u3 !rho u3 u3 *WdetJ 167 fresli(:,7) = fresli(:,1) * u2 !rho u1 u2 *WdetJ 168 fresli(:,8) = fresli(:,1) * u3 !rho u1 u3 *WdetJ 169 fresli(:,9) = fresli(:,2) * u3 !rho u2 u3 *WdetJ 170 171 fresli(:,16) = strnrm(:,intp) * fresli(:,10) ! rho *|Eps| *Eps11 *WdetJ 172 fresli(:,17) = strnrm(:,intp) * fresli(:,11) ! rho *|Eps| *Eps22 *WdetJ 173 fresli(:,18) = strnrm(:,intp) * fresli(:,12) ! rho *|Eps| *Eps33 *WdetJ 174 fresli(:,19) = strnrm(:,intp) * fresli(:,13) ! rho *|Eps| *Eps12 *WdetJ 175 fresli(:,20) = strnrm(:,intp) * fresli(:,14) ! rho *|Eps| *Eps13 *WdetJ 176 fresli(:,21) = strnrm(:,intp) * fresli(:,15) ! rho *|Eps| *Eps23 *WdetJ 177 178 fresli(:,23) = WdetJ ! Integral of 1 over the element 179c 180 do i = 1, 23 181 fresl(:,i) = fresl(:,i) + fresli(:,i) 182 enddo 183 184 enddo !end of loop over integration points 185c 186 do j = 1,nshl 187 do nel = 1,npro 188 fres(ien(nel,j),:) = fres(ien(nel,j),:) + fresl(nel,:) 189 enddo 190 enddo 191 192 return 193 end 194 195 196 197 198 199 200 201 202 203