1!$Id:$ 2 subroutine resid3d(xsj,shp,sig,d,vl,al,p,ndf,l) 3 4! Purpose: 3-D residual routine 5 6 implicit none 7 8 include 'eldata.h' 9 include 'elplot.h' 10 include 'eltran.h' 11 include 'prld1.h' 12 13 integer ndf,l, j,k 14 real*8 xsj 15 16 real*8 b1,b2,b3,rr 17 real*8 aj1,aj2,aj3,aj0,lfac,cfac 18 19 real*8 d(*),vl(ndf,*),al(ndf,*),p(ndf,*) 20 real*8 shp(4,*),sig(*),ac(3),vc(3) 21 22 save 23 24! Compute stress-divergence vector (p) 25 26 if(int(d(74)).gt.0) then 27 b1 = d(11) + prldv(int(d(74)))*d(71) 28 else 29 b1 = d(11)*dm 30 endif 31 32 if(int(d(75)).gt.0) then 33 b2 = d(12) + prldv(int(d(75)))*d(72) 34 else 35 b2 = d(12)*dm 36 endif 37 38 if(int(d(76)).gt.0) then 39 b3 = d(13) + prldv(int(d(76)))*d(73) 40 else 41 b3 = d(13)*dm 42 endif 43 44 rr = d(4) 45 if(d(7).ge.0.0d0) then 46 cfac = d(7) 47 lfac = 1.d0 - cfac 48 else 49 cfac = 0.0d0 50 lfac = 0.0d0 51 endif 52 53! Store time history plot data for element 54 55 k = 6*(l-1) 56 do j = 1,6 57 tt(j+k) = sig(j) 58 end do ! j 59 60! Compute accelerations 61 62 ac(1) = 0.0d0 63 ac(2) = 0.0d0 64 ac(3) = 0.0d0 65 do j = 1,nel 66 ac(1) = ac(1) + shp(4,j)*al(1,j) 67 ac(2) = ac(2) + shp(4,j)*al(2,j) 68 ac(3) = ac(3) + shp(4,j)*al(3,j) 69 end do ! j 70 ac(1) = rr*ac(1)*cfac 71 ac(2) = rr*ac(2)*cfac 72 ac(3) = rr*ac(3)*cfac 73 74! For Rayleigh Mass Damping: Compute velocity 75 76 if(d(77).ne.0.0d0) then 77 vc(1) = 0.0d0 78 vc(2) = 0.0d0 79 vc(3) = 0.0d0 80 do j = 1,nel 81 vc(1) = vc(1) + shp(4,j)*vl(1,j) 82 vc(2) = vc(2) + shp(4,j)*vl(2,j) 83 vc(3) = vc(3) + shp(4,j)*vl(3,j) 84 end do ! j 85 vc(1) = vc(1)*cfac 86 vc(2) = vc(2)*cfac 87 vc(3) = vc(3)*cfac 88 89 do j = 1,nel 90 aj0 = shp(4,j)*xsj*rr*d(77) 91 p(1,j) = p(1,j) - (vc(1) + lfac*vl(1,j))*aj0 92 p(2,j) = p(2,j) - (vc(2) + lfac*vl(2,j))*aj0 93 p(3,j) = p(3,j) - (vc(3) + lfac*vl(3,j))*aj0 94 end do ! j 95 96 endif 97 98! Loop over rows 99 100 do j = 1,nel 101 aj1 = shp(1,j)*xsj 102 aj2 = shp(2,j)*xsj 103 aj3 = shp(3,j)*xsj 104 aj0 = lfac*rr 105 106! Compute gravity, thermal, inertia, and stress contributions 107 108 p(1,j) = p(1,j) + (b1 - ac(1) - aj0*al(1,j))*shp(4,j)*xsj 109 & - aj1*sig(1) - aj2*sig(4) - aj3*sig(6) 110 p(2,j) = p(2,j) + (b2 - ac(2) - aj0*al(2,j))*shp(4,j)*xsj 111 & - aj1*sig(4) - aj2*sig(2) - aj3*sig(5) 112 p(3,j) = p(3,j) + (b3 - ac(3) - aj0*al(3,j))*shp(4,j)*xsj 113 & - aj1*sig(6) - aj2*sig(5) - aj3*sig(3) 114 115 end do ! j 116 117 end 118