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