Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d10r_nl.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 SUBROUTINE v3d10r_nl(coor,matcstet,lmcstet,R_in,d,ci,cj,&
54  s11,s22,s33,s12,s23,s13, &
55  numnp,nstart,nend,numcstet,numat_vol)
56 
57 !________________________________________________________________________
58 !
59 ! V3D10 - Performs displacement based computations for Volumetric 3D,
60 ! 10-node Quadratic Linear Elastic Tetrahedron with quadratic
61 ! interpolation functions. (linear strain tetrahedra).
62 ! Large Deformation. Returns the internal force vector R_in.
63 !
64 ! DATE: 10.2001 AUTHOR: SCOT BREITENFELD
65 !________________________________________________________________________
66 
67  IMPLICIT NONE
68 !-----Global variables
69  INTEGER :: numnp ! number of nodes
70  INTEGER :: numat_vol ! number of volumetric materials
71  INTEGER :: numcstet ! number of LSTets
72 !-- coordinate array
73  REAL*8, DIMENSION(1:3,1:numnp) :: coor
74 !-- elastic stiffness consts
75  REAL*8, DIMENSION(1:9,1:numat_vol) :: ci,cj
76 !-- internal force
77  REAL*8, DIMENSION(1:3*numnp) :: r_in
78 !-- displacement vector
79  REAL*8, DIMENSION(1:3*numnp) :: d
80 !-- CSTet stress
81  REAL*8, DIMENSION(1:4,1:numcstet) :: s11, s22, s33, s12, s23, s13
82 !-- connectivity table for CSTet
83  INTEGER, DIMENSION(1:10,1:numcstet) :: lmcstet
84 !-- mat number for CSTet element
85  INTEGER, DIMENSION(1:numcstet) :: matcstet
86 !---- Local variables
87 !-- node numbers
88  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
89 !-- x, y, and z displacements of nodes
90  REAL*8 :: u1,u2,u3,u4,u5,u6,u7,u8,u9,u10
91  REAL*8 :: v1,v2,v3,v4,v5,v6,v7,v8,v9,v10
92  REAL*8 :: w1,w2,w3,w4,w5,w6,w7,w8,w9,w10
93 !-- 6*volume and the volume
94  REAL*8 :: vx6,vx6inv
95 !-- spacial derivatives
96  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10
97  REAL*8 :: b11,b12,b13,b14,b15,b16,b17,b18,b19,b20
98  REAL*8 :: b21,b22,b23,b24,b25,b26,b27,b28,b29,b30
99 !-- partial derivatives of the displacement
100  REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx
101 !-- strains
102  REAL*8 :: e11,e22,e33,e12,e23,e13
103 !-- coordinate holding variable
104  REAL*8 :: x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
105  REAL*8 :: y1,y2,y3,y4,y5,y6,y7,y8,y9,y10
106  REAL*8 :: z1,z2,z3,z4,z5,z6,z7,z8,z9,z10
107 !-- dummy and counters
108  INTEGER :: i,j,nstart,nend
109  REAL*8 :: aux1,aux2,aux3,aux4,aux5,aux6,aux7,aux8,aux9,aux10,aux11,aux12
110 !-- partial internal force
111  REAL*8 :: r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16,r17,r18
112  REAL*8 :: r19,r20,r21,r22,r23,r24,r25,r26,r27,r28,r29,r30
113  REAL*8 :: g1, g2, g3, g4
114  REAL*8 :: xn1, xn2, xn3, xn4
115 !-- Coordinate subtractions
116  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
117 !--
118  REAL*8 :: c11, c21, c31
119  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8,k1n9,k1n10
120  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8,k2n9,k2n10
121  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8,k3n9,k3n10
122 
123  DO i = nstart, nend
124 
125  j = matcstet(i)
126 
127  n1 = lmcstet(1,i)
128  n2 = lmcstet(2,i)
129  n3 = lmcstet(3,i)
130  n4 = lmcstet(4,i)
131  n5 = lmcstet(5,i)
132  n6 = lmcstet(6,i)
133  n7 = lmcstet(7,i)
134  n8 = lmcstet(8,i)
135  n9 = lmcstet(9,i)
136  n10 = lmcstet(10,i)
137 
138  k3n1 = 3*n1
139  k3n2 = 3*n2
140  k3n3 = 3*n3
141  k3n4 = 3*n4
142  k3n5 = 3*n5
143  k3n6 = 3*n6
144  k3n7 = 3*n7
145  k3n8 = 3*n8
146  k3n9 = 3*n9
147  k3n10 = 3*n10
148 
149  k2n1 = k3n1 - 1
150  k2n2 = k3n2 - 1
151  k2n3 = k3n3 - 1
152  k2n4 = k3n4 - 1
153  k2n5 = k3n5 - 1
154  k2n6 = k3n6 - 1
155  k2n7 = k3n7 - 1
156  k2n8 = k3n8 - 1
157  k2n9 = k3n9 - 1
158  k2n10 = k3n10 - 1
159 
160  k1n1 = k3n1 - 2
161  k1n2 = k3n2 - 2
162  k1n3 = k3n3 - 2
163  k1n4 = k3n4 - 2
164  k1n5 = k3n5 - 2
165  k1n6 = k3n6 - 2
166  k1n7 = k3n7 - 2
167  k1n8 = k3n8 - 2
168  k1n9 = k3n9 - 2
169  k1n10 = k3n10 - 2
170  ! k#n# dummy variables replaces:
171  u1 = d(k1n1) ! (3*n1 -2)
172  u2 = d(k1n2) ! (3*n2 -2)
173  u3 = d(k1n3) ! (3*n3 -2)
174  u4 = d(k1n4) ! (3*n4 -2)
175  u5 = d(k1n5) ! (3*n5 -2)
176  u6 = d(k1n6) ! (3*n6 -2)
177  u7 = d(k1n7) ! (3*n7 -2)
178  u8 = d(k1n8) ! (3*n8 -2)
179  u9 = d(k1n9) ! (3*n9 -2)
180  u10 = d(k1n10) ! (3*n10-2)
181  v1 = d(k2n1) ! (3*n1 -1)
182  v2 = d(k2n2) ! (3*n2 -1)
183  v3 = d(k2n3) ! (3*n3 -1)
184  v4 = d(k2n4) ! (3*n4 -1)
185  v5 = d(k2n5) ! (3*n5 -1)
186  v6 = d(k2n6) ! (3*n6 -1)
187  v7 = d(k2n7) ! (3*n7 -1)
188  v8 = d(k2n8) ! (3*n8 -1)
189  v9 = d(k2n9) ! (3*n9 -1)
190  v10 = d(k2n10) ! (3*n10-1)
191  w1 = d(k3n1) ! (3*n1)
192  w2 = d(k3n2) ! (3*n2)
193  w3 = d(k3n3) ! (3*n3)
194  w4 = d(k3n4) ! (3*n4)
195  w5 = d(k3n5) ! (3*n5)
196  w6 = d(k3n6) ! (3*n6)
197  w7 = d(k3n7) ! (3*n7)
198  w8 = d(k3n8) ! (3*n8)
199  w9 = d(k3n9) ! (3*n9)
200  w10 = d(k3n10) ! (3*n10)
201 
202  x1 = coor(1,n1)
203  x2 = coor(1,n2)
204  x3 = coor(1,n3)
205  x4 = coor(1,n4)
206  y1 = coor(2,n1)
207  y2 = coor(2,n2)
208  y3 = coor(2,n3)
209  y4 = coor(2,n4)
210  z1 = coor(3,n1)
211  z2 = coor(3,n2)
212  z3 = coor(3,n3)
213  z4 = coor(3,n4)
214 
215  vx6 = x2*y3*z4 - x2*y4*z3 - y2*x3*z4 + y2*x4*z3 + z2*x3*y4 &
216  - z2*x4*y3 - x1*y3*z4 + x1*y4*z3 + x1*y2*z4 - x1*y2*z3 &
217  - x1*z2*y4 + x1*z2*y3 + y1*x3*z4 - y1*x4*z3 - y1*x2*z4 &
218  + y1*x2*z3 + y1*z2*x4 - y1*z2*x3 - z1*x3*y4 + z1*x4*y3 &
219  + z1*x2*y4 - z1*x2*y3 - z1*y2*x4 + z1*y2*x3
220 
221  vx6inv = 1.d0 / vx6
222 
223  aux1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3)
224  aux2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3)
225  aux3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3)
226  aux4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3)
227  aux5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3)
228  aux6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3)
229  aux7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2)
230  aux8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2)
231  aux9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2)
232  aux10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2)
233  aux11 =-(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2)
234  aux12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2)
235 
236 !-- Gauss Point 1
237 
238  g1 = 0.58541020d0
239  g2 = 0.13819660d0
240  g3 = 0.13819660d0
241  g4 = 0.13819660d0
242 
243  xn1 = (4.d0*g1-1.d0) ! derivative of shape function
244  xn2 = (4.d0*g2-1.d0) ! dN_i/dzeta_i
245  xn3 = (4.d0*g3-1.d0)
246  xn4 = (4.d0*g4-1.d0)
247 ! xN5 = 4.d0*g1*g2
248 ! xN6 = 4.d0*g2*g3
249 ! xN7 = 4.d0*g3*g1
250 ! xN8 = 4.d0*g1*g4
251 ! xN9 = 4.d0*g2*g4
252 ! xN10= 4.d0*g3*g4
253 
254  b1 = aux1*xn1
255  b2 = aux2*xn1
256  b3 = aux3*xn1
257  b4 = aux4*xn2
258  b5 = aux5*xn2
259  b6 = aux6*xn2
260  b7 = aux7*xn3
261  b8 = aux8*xn3
262  b9 = aux9*xn3
263  b10 = aux10*xn4
264  b11 = aux11*xn4
265  b12 = aux12*xn4
266 
267  b13 = 4.d0*(g2*aux1 + g1*aux4)
268  b14 = 4.d0*(g2*aux2 + g1*aux5)
269  b15 = 4.d0*(g2*aux3 + g1*aux6)
270 
271  b16 = 4.d0*(g3*aux4 + g2*aux7)
272  b17 = 4.d0*(g3*aux5 + g2*aux8)
273  b18 = 4.d0*(g3*aux6 + g2*aux9)
274 
275  b19 = 4.d0*(g1*aux7 + g3*aux1)
276  b20 = 4.d0*(g1*aux8 + g3*aux2)
277  b21 = 4.d0*(g1*aux9 + g3*aux3)
278 
279  b22 = 4.d0*(g4*aux1 + g1*aux10)
280  b23 = 4.d0*(g4*aux2 + g1*aux11)
281  b24 = 4.d0*(g4*aux3 + g1*aux12)
282 
283  b25 = 4.d0*(g4*aux4 + g2*aux10)
284  b26 = 4.d0*(g4*aux5 + g2*aux11)
285  b27 = 4.d0*(g4*aux6 + g2*aux12)
286 
287  b28 = 4.d0*(g4*aux7 + g3*aux10)
288  b29 = 4.d0*(g4*aux8 + g3*aux11)
289  b30 = 4.d0*(g4*aux9 + g3*aux12)
290 
291 
292 !-----Calculate displacement gradient (H)
293  dudx = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10)*vx6inv
294  dvdy = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10)*vx6inv
295  dwdz = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10)*vx6inv
296  dudy = (b2*u1 + b5*u2 + b8*u3 + b11*u4 + b14*u5 + b17*u6 + b20*u7 + b23*u8 + b26*u9 + b29*u10)*vx6inv
297  dvdx = (b1*v1 + b4*v2 + b7*v3 + b10*v4 + b13*v5 + b16*v6 + b19*v7 + b22*v8 + b25*v9 + b28*v10)*vx6inv
298  dvdz = (b3*v1 + b6*v2 + b9*v3 + b12*v4 + b15*v5 + b18*v6 + b21*v7 + b24*v8 + b27*v9 + b30*v10)*vx6inv
299  dwdy = (b2*w1 + b5*w2 + b8*w3 + b11*w4 + b14*w5 + b17*w6 + b20*w7 + b23*w8 + b26*w9 + b29*w10)*vx6inv
300  dudz = (b3*u1 + b6*u2 + b9*u3 + b12*u4 + b15*u5 + b18*u6 + b21*u7 + b24*u8 + b27*u9 + b30*u10)*vx6inv
301  dwdx = (b1*w1 + b4*w2 + b7*w3 + b10*w4 + b13*w5 + b16*w6 + b19*w7 + b22*w8 + b25*w9 + b28*w10)*vx6inv
302 !
303 ! Green-Lagrange Strain (E)
304 !
305  e11 = dudx + 0.5d0*(dudx*dudx + dvdx*dvdx + dwdx*dwdx)
306  e22 = dvdy + 0.5d0*(dudy*dudy + dvdy*dvdy + dwdy*dwdy)
307  e33 = dwdz + 0.5d0*(dudz*dudz + dvdz*dvdz + dwdz*dwdz)
308  e12 = dudy + dvdx + ( dudx*dudy + dvdx*dvdy + dwdx*dwdy )
309  e23 = dvdz + dwdy + ( dudy*dudz + dvdy*dvdz + dwdy*dwdz )
310  e13 = dwdx + dudz + ( dudz*dudx + dvdz*dvdx + dwdz*dwdx )
311 !
312 ! Second Piola-Kirchhoff stress (S)
313 !
314  s11(1,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
315  s22(1,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
316  s33(1,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
317  s12(1,i) = e12*ci(7,j)
318  s23(1,i) = e23*ci(8,j)
319  s13(1,i) = e13*ci(9,j)
320 
321  r1 = &
322  ( s11(1,i)*b1*(1.d0+dudx) + s22(1,i)*b2*dudy + s33(1,i)*b3*dudz &
323  + s12(1,i)*( b2*(1.d0+dudx) + b1*dudy ) &
324  + s23(1,i)*( b3*dudy + b2*dudz ) &
325  + s13(1,i)*( b3*(1.d0+dudx) + b1*dudz ) )
326  r2 = &
327  ( s11(1,i)*b1*dvdx + s22(1,i)*b2*(1.d0+dvdy) + s33(1,i)*b3*dvdz &
328  + s12(1,i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
329  + s23(1,i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
330  + s13(1,i)*( b3*dvdx + b1*dvdz ) )
331  r3 = &
332  ( s11(1,i)*b1*dwdx + s22(1,i)*b2*dwdy + s33(1,i)*b3*(1.d0+dwdz) &
333  + s12(1,i)*( b2*dwdx + b1*dwdy ) &
334  + s23(1,i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
335  + s13(1,i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
336 
337  r4 = &
338  ( s11(1,i)*b4*(1.d0+dudx) + s22(1,i)*b5*dudy + s33(1,i)*b6*dudz &
339  + s12(1,i)*( b5*(1.d0+dudx) + b4*dudy ) &
340  + s23(1,i)*( b6*dudy + b5*dudz ) &
341  + s13(1,i)*( b6*(1.d0+dudx) + b4*dudz ) )
342  r5 = &
343  ( s11(1,i)*b4*dvdx + s22(1,i)*b5*(1.d0+dvdy) + s33(1,i)*b6*dvdz &
344  + s12(1,i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
345  + s23(1,i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
346  + s13(1,i)*( b6*dvdx + b4*dvdz ) )
347  r6 = &
348  ( s11(1,i)*b4*dwdx + s22(1,i)*b5*dwdy + s33(1,i)*b6*(1.d0+dwdz) &
349  + s12(1,i)*( b5*dwdx + b4*dwdy ) &
350  + s23(1,i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
351  + s13(1,i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
352 
353  r7 = &
354  ( s11(1,i)*b7*(1.d0+dudx) + s22(1,i)*b8*dudy + s33(1,i)*b9*dudz &
355  + s12(1,i)*( b8*(1.d0+dudx) + b7*dudy ) &
356  + s23(1,i)*( b9*dudy + b8*dudz ) &
357  + s13(1,i)*( b9*(1.d0+dudx) + b7*dudz ) )
358  r8 = &
359  ( s11(1,i)*b7*dvdx + s22(1,i)*b8*(1.d0+dvdy) + s33(1,i)*b9*dvdz &
360  + s12(1,i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
361  + s23(1,i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
362  + s13(1,i)*( b9*dvdx + b7*dvdz ) )
363  r9 = &
364  ( s11(1,i)*b7*dwdx + s22(1,i)*b8*dwdy + s33(1,i)*b9*(1.d0+dwdz) &
365  + s12(1,i)*( b8*dwdx + b7*dwdy ) &
366  + s23(1,i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
367  + s13(1,i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
368 
369  r10 = &
370  ( s11(1,i)*b10*(1.d0+dudx) + s22(1,i)*b11*dudy+s33(1,i)*b12*dudz &
371  + s12(1,i)*( b11*(1.d0+dudx) + b10*dudy ) &
372  + s23(1,i)*( b12*dudy + b11*dudz ) &
373  + s13(1,i)*( b12*(1.d0+dudx) + b10*dudz ) )
374  r11 = &
375  ( s11(1,i)*b10*dvdx + s22(1,i)*b11*(1.d0+dvdy)+s33(1,i)*b12*dvdz &
376  + s12(1,i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
377  + s23(1,i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
378  + s13(1,i)*( b12*dvdx + b10*dvdz ) )
379  r12 = &
380  ( s11(1,i)*b10*dwdx + s22(1,i)*b11*dwdy+s33(1,i)*b12*(1.d0+dwdz) &
381  + s12(1,i)*( b11*dwdx + b10*dwdy ) &
382  + s23(1,i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
383  + s13(1,i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )
384 
385  r13 = &
386  ( s11(1,i)*b13*(1.d0+dudx) + s22(1,i)*b14*dudy + s33(1,i)*b15*dudz &
387  + s12(1,i)*( b14*(1.d0+dudx) + b13*dudy ) &
388  + s23(1,i)*( b15*dudy + b14*dudz ) &
389  + s13(1,i)*( b15*(1.d0+dudx) + b13*dudz ) )
390  r14 = &
391  ( s11(1,i)*b13*dvdx + s22(1,i)*b14*(1.d0+dvdy) + s33(1,i)*b15*dvdz &
392  + s12(1,i)*( b13*(1.d0+dvdy) + b14*dvdx ) &
393  + s23(1,i)*( b15*(1.d0+dvdy) + b14*dvdz ) &
394  + s13(1,i)*( b15*dvdx + b13*dvdz ) )
395  r15 = &
396  ( s11(1,i)*b13*dwdx + s22(1,i)*b14*dwdy + s33(1,i)*b15*(1.d0+dwdz) &
397  + s12(1,i)*( b14*dwdx + b13*dwdy ) &
398  + s23(1,i)*( b15*dwdy + b14*(1.d0 + dwdz) ) &
399  + s13(1,i)*( b15*dwdx + b13*(1.d0 + dwdz) ) )
400 
401  r16 = &
402  ( s11(1,i)*b16*(1.d0+dudx) + s22(1,i)*b17*dudy + s33(1,i)*b18*dudz &
403  + s12(1,i)*( b17*(1.d0+dudx) + b16*dudy ) &
404  + s23(1,i)*( b18*dudy + b17*dudz ) &
405  + s13(1,i)*( b18*(1.d0+dudx) + b16*dudz ) )
406  r17 = &
407  ( s11(1,i)*b16*dvdx + s22(1,i)*b17*(1.d0+dvdy) + s33(1,i)*b18*dvdz &
408  + s12(1,i)*( b16*(1.d0+dvdy) + b17*dvdx ) &
409  + s23(1,i)*( b18*(1.d0+dvdy) + b17*dvdz ) &
410  + s13(1,i)*( b18*dvdx + b16*dvdz ) )
411  r18 = &
412  ( s11(1,i)*b16*dwdx + s22(1,i)*b17*dwdy + s33(1,i)*b18*(1.d0+dwdz) &
413  + s12(1,i)*( b17*dwdx + b16*dwdy ) &
414  + s23(1,i)*( b18*dwdy + b17*(1.d0 + dwdz) ) &
415  + s13(1,i)*( b18*dwdx + b16*(1.d0 + dwdz) ) )
416 
417  r19 = &
418  ( s11(1,i)*b19*(1.d0+dudx) + s22(1,i)*b20*dudy + s33(1,i)*b21*dudz &
419  + s12(1,i)*( b20*(1.d0+dudx) + b19*dudy ) &
420  + s23(1,i)*( b21*dudy + b20*dudz ) &
421  + s13(1,i)*( b21*(1.d0+dudx) + b19*dudz ) )
422  r20 = &
423  ( s11(1,i)*b19*dvdx + s22(1,i)*b20*(1.d0+dvdy) + s33(1,i)*b21*dvdz &
424  + s12(1,i)*( b19*(1.d0+dvdy) + b20*dvdx ) &
425  + s23(1,i)*( b21*(1.d0+dvdy) + b20*dvdz ) &
426  + s13(1,i)*( b21*dvdx + b19*dvdz ) )
427  r21 = &
428  ( s11(1,i)*b19*dwdx + s22(1,i)*b20*dwdy + s33(1,i)*b21*(1.d0+dwdz) &
429  + s12(1,i)*( b20*dwdx + b19*dwdy ) &
430  + s23(1,i)*( b21*dwdy + b20*(1.d0 + dwdz) ) &
431  + s13(1,i)*( b21*dwdx + b19*(1.d0 + dwdz) ) )
432 
433  r22 = &
434  ( s11(1,i)*b22*(1.d0+dudx) + s22(1,i)*b23*dudy+s33(1,i)*b24*dudz &
435  + s12(1,i)*( b23*(1.d0+dudx) + b22*dudy ) &
436  + s23(1,i)*( b24*dudy + b23*dudz ) &
437  + s13(1,i)*( b24*(1.d0+dudx) + b22*dudz ) )
438  r23 = &
439  ( s11(1,i)*b22*dvdx + s22(1,i)*b23*(1.d0+dvdy)+s33(1,i)*b24*dvdz &
440  + s12(1,i)*( b22*(1.d0+dvdy) + b23*dvdx ) &
441  + s23(1,i)*( b24*(1.d0+dvdy) + b23*dvdz ) &
442  + s13(1,i)*( b24*dvdx + b22*dvdz ) )
443  r24 = &
444  ( s11(1,i)*b22*dwdx + s22(1,i)*b23*dwdy+s33(1,i)*b24*(1.d0+dwdz) &
445  + s12(1,i)*( b23*dwdx + b22*dwdy ) &
446  + s23(1,i)*( b24*dwdy + b23*(1.d0 + dwdz) ) &
447  + s13(1,i)*( b24*dwdx + b22*(1.d0 + dwdz) ) )
448 
449  r25 = &
450  ( s11(1,i)*b25*(1.d0+dudx) + s22(1,i)*b26*dudy+s33(1,i)*b27*dudz &
451  + s12(1,i)*( b26*(1.d0+dudx) + b25*dudy ) &
452  + s23(1,i)*( b27*dudy + b26*dudz ) &
453  + s13(1,i)*( b27*(1.d0+dudx) + b25*dudz ) )
454  r26 = &
455  ( s11(1,i)*b25*dvdx + s22(1,i)*b26*(1.d0+dvdy)+s33(1,i)*b27*dvdz &
456  + s12(1,i)*( b25*(1.d0+dvdy) + b26*dvdx ) &
457  + s23(1,i)*( b27*(1.d0+dvdy) + b26*dvdz ) &
458  + s13(1,i)*( b27*dvdx + b25*dvdz ) )
459  r27 = &
460  ( s11(1,i)*b25*dwdx + s22(1,i)*b26*dwdy+s33(1,i)*b27*(1.d0+dwdz) &
461  + s12(1,i)*( b26*dwdx + b25*dwdy ) &
462  + s23(1,i)*( b27*dwdy + b26*(1.d0 + dwdz) ) &
463  + s13(1,i)*( b27*dwdx + b25*(1.d0 + dwdz) ) )
464 
465  r28 = &
466  ( s11(1,i)*b28*(1.d0+dudx) + s22(1,i)*b29*dudy+s33(1,i)*b30*dudz &
467  + s12(1,i)*( b29*(1.d0+dudx) + b28*dudy ) &
468  + s23(1,i)*( b30*dudy + b29*dudz ) &
469  + s13(1,i)*( b30*(1.d0+dudx) + b28*dudz ) )
470  r29 = &
471  ( s11(1,i)*b28*dvdx + s22(1,i)*b29*(1.d0+dvdy)+s33(1,i)*b30*dvdz &
472  + s12(1,i)*( b28*(1.d0+dvdy) + b29*dvdx ) &
473  + s23(1,i)*( b30*(1.d0+dvdy) + b29*dvdz ) &
474  + s13(1,i)*( b30*dvdx + b28*dvdz ) )
475  r30 = &
476  ( s11(1,i)*b28*dwdx + s22(1,i)*b29*dwdy+s33(1,i)*b30*(1.d0+dwdz) &
477  + s12(1,i)*( b29*dwdx + b28*dwdy ) &
478  + s23(1,i)*( b30*dwdy + b29*(1.d0 + dwdz) ) &
479  + s13(1,i)*( b30*dwdx + b28*(1.d0 + dwdz) ) )
480 
481 !-- Gauss Point 2
482 
483  g1 = 0.13819660d0
484  g2 = 0.58541020d0
485  g3 = 0.13819660d0
486  g4 = 0.13819660d0
487 
488  xn1 = (4.d0*g1-1.d0) ! derivative of shape function
489  xn2 = (4.d0*g2-1.d0) ! dN_i/dzeta_i
490  xn3 = (4.d0*g3-1.d0)
491  xn4 = (4.d0*g4-1.d0)
492 ! xN5 = 4.d0*g1*g2
493 ! xN6 = 4.d0*g2*g3
494 ! xN7 = 4.d0*g3*g1
495 ! xN8 = 4.d0*g1*g4
496 ! xN9 = 4.d0*g2*g4
497 ! xN10= 4.d0*g3*g4
498 
499  b1 = aux1*xn1
500  b2 = aux2*xn1
501  b3 = aux3*xn1
502  b4 = aux4*xn2
503  b5 = aux5*xn2
504  b6 = aux6*xn2
505  b7 = aux7*xn3
506  b8 = aux8*xn3
507  b9 = aux9*xn3
508  b10 = aux10*xn4
509  b11 = aux11*xn4
510  b12 = aux12*xn4
511 
512  b13 = 4.d0*(g2*aux1 + g1*aux4)
513  b14 = 4.d0*(g2*aux2 + g1*aux5)
514  b15 = 4.d0*(g2*aux3 + g1*aux6)
515 
516  b16 = 4.d0*(g3*aux4 + g2*aux7)
517  b17 = 4.d0*(g3*aux5 + g2*aux8)
518  b18 = 4.d0*(g3*aux6 + g2*aux9)
519 
520  b19 = 4.d0*(g1*aux7 + g3*aux1)
521  b20 = 4.d0*(g1*aux8 + g3*aux2)
522  b21 = 4.d0*(g1*aux9 + g3*aux3)
523 
524  b22 = 4.d0*(g4*aux1 + g1*aux10)
525  b23 = 4.d0*(g4*aux2 + g1*aux11)
526  b24 = 4.d0*(g4*aux3 + g1*aux12)
527 
528  b25 = 4.d0*(g4*aux4 + g2*aux10)
529  b26 = 4.d0*(g4*aux5 + g2*aux11)
530  b27 = 4.d0*(g4*aux6 + g2*aux12)
531 
532  b28 = 4.d0*(g4*aux7 + g3*aux10)
533  b29 = 4.d0*(g4*aux8 + g3*aux11)
534  b30 = 4.d0*(g4*aux9 + g3*aux12)
535 !-----Calculate displacement gradient (H)
536  dudx = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10)*vx6inv
537  dvdy = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10)*vx6inv
538  dwdz = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10)*vx6inv
539  dudy = (b2*u1 + b5*u2 + b8*u3 + b11*u4 + b14*u5 + b17*u6 + b20*u7 + b23*u8 + b26*u9 + b29*u10)*vx6inv
540  dvdx = (b1*v1 + b4*v2 + b7*v3 + b10*v4 + b13*v5 + b16*v6 + b19*v7 + b22*v8 + b25*v9 + b28*v10)*vx6inv
541  dvdz = (b3*v1 + b6*v2 + b9*v3 + b12*v4 + b15*v5 + b18*v6 + b21*v7 + b24*v8 + b27*v9 + b30*v10)*vx6inv
542  dwdy = (b2*w1 + b5*w2 + b8*w3 + b11*w4 + b14*w5 + b17*w6 + b20*w7 + b23*w8 + b26*w9 + b29*w10)*vx6inv
543  dudz = (b3*u1 + b6*u2 + b9*u3 + b12*u4 + b15*u5 + b18*u6 + b21*u7 + b24*u8 + b27*u9 + b30*u10)*vx6inv
544  dwdx = (b1*w1 + b4*w2 + b7*w3 + b10*w4 + b13*w5 + b16*w6 + b19*w7 + b22*w8 + b25*w9 + b28*w10)*vx6inv
545 !
546 ! Green-Lagrange Strain (E)
547 !
548  e11 = dudx + 0.5d0*(dudx*dudx + dvdx*dvdx + dwdx*dwdx)
549  e22 = dvdy + 0.5d0*(dudy*dudy + dvdy*dvdy + dwdy*dwdy)
550  e33 = dwdz + 0.5d0*(dudz*dudz + dvdz*dvdz + dwdz*dwdz)
551  e12 = dudy + dvdx + ( dudx*dudy + dvdx*dvdy + dwdx*dwdy )
552  e23 = dvdz + dwdy + ( dudy*dudz + dvdy*dvdz + dwdy*dwdz )
553  e13 = dwdx + dudz + ( dudz*dudx + dvdz*dvdx + dwdz*dwdx )
554 !
555 ! Second Piola-Kirchhoff stress (S)
556 !
557  s11(2,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
558  s22(2,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
559  s33(2,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
560  s12(2,i) = e12*ci(7,j)
561  s23(2,i) = e23*ci(8,j)
562  s13(2,i) = e13*ci(9,j)
563 
564  r1 = r1 + &
565  ( s11(2,i)*b1*(1.d0+dudx) + s22(2,i)*b2*dudy + s33(2,i)*b3*dudz &
566  + s12(2,i)*( b2*(1.d0+dudx) + b1*dudy ) &
567  + s23(2,i)*( b3*dudy + b2*dudz ) &
568  + s13(2,i)*( b3*(1.d0+dudx) + b1*dudz ) )
569  r2 = r2 + &
570  ( s11(2,i)*b1*dvdx + s22(2,i)*b2*(1.d0+dvdy) + s33(2,i)*b3*dvdz &
571  + s12(2,i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
572  + s23(2,i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
573  + s13(2,i)*( b3*dvdx + b1*dvdz ) )
574  r3 = r3 + &
575  ( s11(2,i)*b1*dwdx + s22(2,i)*b2*dwdy + s33(2,i)*b3*(1.d0+dwdz) &
576  + s12(2,i)*( b2*dwdx + b1*dwdy ) &
577  + s23(2,i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
578  + s13(2,i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
579 
580  r4 = r4 + &
581  ( s11(2,i)*b4*(1.d0+dudx) + s22(2,i)*b5*dudy + s33(2,i)*b6*dudz &
582  + s12(2,i)*( b5*(1.d0+dudx) + b4*dudy ) &
583  + s23(2,i)*( b6*dudy + b5*dudz ) &
584  + s13(2,i)*( b6*(1.d0+dudx) + b4*dudz ) )
585  r5 = r5 + &
586  ( s11(2,i)*b4*dvdx + s22(2,i)*b5*(1.d0+dvdy) + s33(2,i)*b6*dvdz &
587  + s12(2,i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
588  + s23(2,i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
589  + s13(2,i)*( b6*dvdx + b4*dvdz ) )
590  r6 = r6 + &
591  ( s11(2,i)*b4*dwdx + s22(2,i)*b5*dwdy + s33(2,i)*b6*(1.d0+dwdz) &
592  + s12(2,i)*( b5*dwdx + b4*dwdy ) &
593  + s23(2,i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
594  + s13(2,i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
595 
596  r7 = r7 + &
597  ( s11(2,i)*b7*(1.d0+dudx) + s22(2,i)*b8*dudy + s33(2,i)*b9*dudz &
598  + s12(2,i)*( b8*(1.d0+dudx) + b7*dudy ) &
599  + s23(2,i)*( b9*dudy + b8*dudz ) &
600  + s13(2,i)*( b9*(1.d0+dudx) + b7*dudz ) )
601  r8 = r8 + &
602  ( s11(2,i)*b7*dvdx + s22(2,i)*b8*(1.d0+dvdy) + s33(2,i)*b9*dvdz &
603  + s12(2,i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
604  + s23(2,i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
605  + s13(2,i)*( b9*dvdx + b7*dvdz ) )
606  r9 = r9 + &
607  ( s11(2,i)*b7*dwdx + s22(2,i)*b8*dwdy + s33(2,i)*b9*(1.d0+dwdz) &
608  + s12(2,i)*( b8*dwdx + b7*dwdy ) &
609  + s23(2,i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
610  + s13(2,i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
611 
612  r10 = r10 + &
613  ( s11(2,i)*b10*(1.d0+dudx) + s22(2,i)*b11*dudy+s33(2,i)*b12*dudz &
614  + s12(2,i)*( b11*(1.d0+dudx) + b10*dudy ) &
615  + s23(2,i)*( b12*dudy + b11*dudz ) &
616  + s13(2,i)*( b12*(1.d0+dudx) + b10*dudz ) )
617  r11 = r11 + &
618  ( s11(2,i)*b10*dvdx + s22(2,i)*b11*(1.d0+dvdy)+s33(2,i)*b12*dvdz &
619  + s12(2,i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
620  + s23(2,i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
621  + s13(2,i)*( b12*dvdx + b10*dvdz ) )
622  r12 = r12 + &
623  ( s11(2,i)*b10*dwdx + s22(2,i)*b11*dwdy+s33(2,i)*b12*(1.d0+dwdz) &
624  + s12(2,i)*( b11*dwdx + b10*dwdy ) &
625  + s23(2,i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
626  + s13(2,i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )
627 
628  r13 = r13 + &
629  ( s11(2,i)*b13*(1.d0+dudx) + s22(2,i)*b14*dudy + s33(2,i)*b15*dudz &
630  + s12(2,i)*( b14*(1.d0+dudx) + b13*dudy ) &
631  + s23(2,i)*( b15*dudy + b14*dudz ) &
632  + s13(2,i)*( b15*(1.d0+dudx) + b13*dudz ) )
633  r14 = r14 + &
634  ( s11(2,i)*b13*dvdx + s22(2,i)*b14*(1.d0+dvdy) + s33(2,i)*b15*dvdz &
635  + s12(2,i)*( b13*(1.d0+dvdy) + b14*dvdx ) &
636  + s23(2,i)*( b15*(1.d0+dvdy) + b14*dvdz ) &
637  + s13(2,i)*( b15*dvdx + b13*dvdz ) )
638  r15 = r15 + &
639  ( s11(2,i)*b13*dwdx + s22(2,i)*b14*dwdy + s33(2,i)*b15*(1.d0+dwdz) &
640  + s12(2,i)*( b14*dwdx + b13*dwdy ) &
641  + s23(2,i)*( b15*dwdy + b14*(1.d0 + dwdz) ) &
642  + s13(2,i)*( b15*dwdx + b13*(1.d0 + dwdz) ) )
643 
644  r16 = r16 + &
645  ( s11(2,i)*b16*(1.d0+dudx) + s22(2,i)*b17*dudy + s33(2,i)*b18*dudz &
646  + s12(2,i)*( b17*(1.d0+dudx) + b16*dudy ) &
647  + s23(2,i)*( b18*dudy + b17*dudz ) &
648  + s13(2,i)*( b18*(1.d0+dudx) + b16*dudz ) )
649  r17 = r17 + &
650  ( s11(2,i)*b16*dvdx + s22(2,i)*b17*(1.d0+dvdy) + s33(2,i)*b18*dvdz &
651  + s12(2,i)*( b16*(1.d0+dvdy) + b17*dvdx ) &
652  + s23(2,i)*( b18*(1.d0+dvdy) + b17*dvdz ) &
653  + s13(2,i)*( b18*dvdx + b16*dvdz ) )
654  r18 = r18 + &
655  ( s11(2,i)*b16*dwdx + s22(2,i)*b17*dwdy + s33(2,i)*b18*(1.d0+dwdz) &
656  + s12(2,i)*( b17*dwdx + b16*dwdy ) &
657  + s23(2,i)*( b18*dwdy + b17*(1.d0 + dwdz) ) &
658  + s13(2,i)*( b18*dwdx + b16*(1.d0 + dwdz) ) )
659 
660  r19 = r19 + &
661  ( s11(2,i)*b19*(1.d0+dudx) + s22(2,i)*b20*dudy + s33(2,i)*b21*dudz &
662  + s12(2,i)*( b20*(1.d0+dudx) + b19*dudy ) &
663  + s23(2,i)*( b21*dudy + b20*dudz ) &
664  + s13(2,i)*( b21*(1.d0+dudx) + b19*dudz ) )
665  r20 = r20 + &
666  ( s11(2,i)*b19*dvdx + s22(2,i)*b20*(1.d0+dvdy) + s33(2,i)*b21*dvdz &
667  + s12(2,i)*( b19*(1.d0+dvdy) + b20*dvdx ) &
668  + s23(2,i)*( b21*(1.d0+dvdy) + b20*dvdz ) &
669  + s13(2,i)*( b21*dvdx + b19*dvdz ) )
670  r21 = r21 + &
671  ( s11(2,i)*b19*dwdx + s22(2,i)*b20*dwdy + s33(2,i)*b21*(1.d0+dwdz) &
672  + s12(2,i)*( b20*dwdx + b19*dwdy ) &
673  + s23(2,i)*( b21*dwdy + b20*(1.d0 + dwdz) ) &
674  + s13(2,i)*( b21*dwdx + b19*(1.d0 + dwdz) ) )
675 
676  r22 = r22 + &
677  ( s11(2,i)*b22*(1.d0+dudx) + s22(2,i)*b23*dudy+s33(2,i)*b24*dudz &
678  + s12(2,i)*( b23*(1.d0+dudx) + b22*dudy ) &
679  + s23(2,i)*( b24*dudy + b23*dudz ) &
680  + s13(2,i)*( b24*(1.d0+dudx) + b22*dudz ) )
681  r23 = r23 + &
682  ( s11(2,i)*b22*dvdx + s22(2,i)*b23*(1.d0+dvdy)+s33(2,i)*b24*dvdz &
683  + s12(2,i)*( b22*(1.d0+dvdy) + b23*dvdx ) &
684  + s23(2,i)*( b24*(1.d0+dvdy) + b23*dvdz ) &
685  + s13(2,i)*( b24*dvdx + b22*dvdz ) )
686  r24 = r24 + &
687  ( s11(2,i)*b22*dwdx + s22(2,i)*b23*dwdy+s33(2,i)*b24*(1.d0+dwdz) &
688  + s12(2,i)*( b23*dwdx + b22*dwdy ) &
689  + s23(2,i)*( b24*dwdy + b23*(1.d0 + dwdz) ) &
690  + s13(2,i)*( b24*dwdx + b22*(1.d0 + dwdz) ) )
691 
692  r25 = r25 + &
693  ( s11(2,i)*b25*(1.d0+dudx) + s22(2,i)*b26*dudy+s33(2,i)*b27*dudz &
694  + s12(2,i)*( b26*(1.d0+dudx) + b25*dudy ) &
695  + s23(2,i)*( b27*dudy + b26*dudz ) &
696  + s13(2,i)*( b27*(1.d0+dudx) + b25*dudz ) )
697  r26 = r26 + &
698  ( s11(2,i)*b25*dvdx + s22(2,i)*b26*(1.d0+dvdy)+s33(2,i)*b27*dvdz &
699  + s12(2,i)*( b25*(1.d0+dvdy) + b26*dvdx ) &
700  + s23(2,i)*( b27*(1.d0+dvdy) + b26*dvdz ) &
701  + s13(2,i)*( b27*dvdx + b25*dvdz ) )
702  r27 = r27 + &
703  ( s11(2,i)*b25*dwdx + s22(2,i)*b26*dwdy+s33(2,i)*b27*(1.d0+dwdz) &
704  + s12(2,i)*( b26*dwdx + b25*dwdy ) &
705  + s23(2,i)*( b27*dwdy + b26*(1.d0 + dwdz) ) &
706  + s13(2,i)*( b27*dwdx + b25*(1.d0 + dwdz) ) )
707 
708  r28 = r28 + &
709  ( s11(2,i)*b28*(1.d0+dudx) + s22(2,i)*b29*dudy+s33(2,i)*b30*dudz &
710  + s12(2,i)*( b29*(1.d0+dudx) + b28*dudy ) &
711  + s23(2,i)*( b30*dudy + b29*dudz ) &
712  + s13(2,i)*( b30*(1.d0+dudx) + b28*dudz ) )
713  r29 = r29 + &
714  ( s11(2,i)*b28*dvdx + s22(2,i)*b29*(1.d0+dvdy)+s33(2,i)*b30*dvdz &
715  + s12(2,i)*( b28*(1.d0+dvdy) + b29*dvdx ) &
716  + s23(2,i)*( b30*(1.d0+dvdy) + b29*dvdz ) &
717  + s13(2,i)*( b30*dvdx + b28*dvdz ) )
718  r30 = r30 + &
719  ( s11(2,i)*b28*dwdx + s22(2,i)*b29*dwdy+s33(2,i)*b30*(1.d0+dwdz) &
720  + s12(2,i)*( b29*dwdx + b28*dwdy ) &
721  + s23(2,i)*( b30*dwdy + b29*(1.d0 + dwdz) ) &
722  + s13(2,i)*( b30*dwdx + b28*(1.d0 + dwdz) ) )
723 !-- Gauss Point 3
724 
725  g1 = 0.13819660d0
726  g2 = 0.13819660d0
727  g3 = 0.58541020d0
728  g4 = 0.13819660d0
729 
730  xn1 = (4.d0*g1-1.d0) ! derivative of shape function
731  xn2 = (4.d0*g2-1.d0) ! dN_i/dzeta_i
732  xn3 = (4.d0*g3-1.d0)
733  xn4 = (4.d0*g4-1.d0)
734  ! xN5 = 4.d0*g1*g2
735  ! xN6 = 4.d0*g2*g3
736  ! xN7 = 4.d0*g3*g1
737 ! xN8 = 4.d0*g1*g4
738 ! xN9 = 4.d0*g2*g4
739 ! xN10= 4.d0*g3*g4
740 
741  b1 = aux1*xn1
742  b2 = aux2*xn1
743  b3 = aux3*xn1
744  b4 = aux4*xn2
745  b5 = aux5*xn2
746  b6 = aux6*xn2
747  b7 = aux7*xn3
748  b8 = aux8*xn3
749  b9 = aux9*xn3
750  b10 = aux10*xn4
751  b11 = aux11*xn4
752  b12 = aux12*xn4
753 
754  b13 = 4.d0*(g2*aux1 + g1*aux4)
755  b14 = 4.d0*(g2*aux2 + g1*aux5)
756  b15 = 4.d0*(g2*aux3 + g1*aux6)
757 
758  b16 = 4.d0*(g3*aux4 + g2*aux7)
759  b17 = 4.d0*(g3*aux5 + g2*aux8)
760  b18 = 4.d0*(g3*aux6 + g2*aux9)
761 
762  b19 = 4.d0*(g1*aux7 + g3*aux1)
763  b20 = 4.d0*(g1*aux8 + g3*aux2)
764  b21 = 4.d0*(g1*aux9 + g3*aux3)
765 
766  b22 = 4.d0*(g4*aux1 + g1*aux10)
767  b23 = 4.d0*(g4*aux2 + g1*aux11)
768  b24 = 4.d0*(g4*aux3 + g1*aux12)
769 
770  b25 = 4.d0*(g4*aux4 + g2*aux10)
771  b26 = 4.d0*(g4*aux5 + g2*aux11)
772  b27 = 4.d0*(g4*aux6 + g2*aux12)
773 
774  b28 = 4.d0*(g4*aux7 + g3*aux10)
775  b29 = 4.d0*(g4*aux8 + g3*aux11)
776  b30 = 4.d0*(g4*aux9 + g3*aux12)
777 !-----Calculate displacement gradient (H)
778  dudx = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10)*vx6inv
779  dvdy = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10)*vx6inv
780  dwdz = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10)*vx6inv
781  dudy = (b2*u1 + b5*u2 + b8*u3 + b11*u4 + b14*u5 + b17*u6 + b20*u7 + b23*u8 + b26*u9 + b29*u10)*vx6inv
782  dvdx = (b1*v1 + b4*v2 + b7*v3 + b10*v4 + b13*v5 + b16*v6 + b19*v7 + b22*v8 + b25*v9 + b28*v10)*vx6inv
783  dvdz = (b3*v1 + b6*v2 + b9*v3 + b12*v4 + b15*v5 + b18*v6 + b21*v7 + b24*v8 + b27*v9 + b30*v10)*vx6inv
784  dwdy = (b2*w1 + b5*w2 + b8*w3 + b11*w4 + b14*w5 + b17*w6 + b20*w7 + b23*w8 + b26*w9 + b29*w10)*vx6inv
785  dudz = (b3*u1 + b6*u2 + b9*u3 + b12*u4 + b15*u5 + b18*u6 + b21*u7 + b24*u8 + b27*u9 + b30*u10)*vx6inv
786  dwdx = (b1*w1 + b4*w2 + b7*w3 + b10*w4 + b13*w5 + b16*w6 + b19*w7 + b22*w8 + b25*w9 + b28*w10)*vx6inv
787 ! Green-Lagrange Strain (E)
788 !
789  e11 = dudx + 0.5d0*(dudx*dudx + dvdx*dvdx + dwdx*dwdx)
790  e22 = dvdy + 0.5d0*(dudy*dudy + dvdy*dvdy + dwdy*dwdy)
791  e33 = dwdz + 0.5d0*(dudz*dudz + dvdz*dvdz + dwdz*dwdz)
792  e12 = dudy + dvdx + ( dudx*dudy + dvdx*dvdy + dwdx*dwdy )
793  e23 = dvdz + dwdy + ( dudy*dudz + dvdy*dvdz + dwdy*dwdz )
794  e13 = dwdx + dudz + ( dudz*dudx + dvdz*dvdx + dwdz*dwdx )
795 !
796 ! Second Piola-Kirchhoff stress (S)
797 !
798  s11(3,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
799  s22(3,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
800  s33(3,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
801  s12(3,i) = e12*ci(7,j)
802  s23(3,i) = e23*ci(8,j)
803  s13(3,i) = e13*ci(9,j)
804 
805  r1 = r1 + &
806  ( s11(3,i)*b1*(1.d0+dudx) + s22(3,i)*b2*dudy + s33(3,i)*b3*dudz &
807  + s12(3,i)*( b2*(1.d0+dudx) + b1*dudy ) &
808  + s23(3,i)*( b3*dudy + b2*dudz ) &
809  + s13(3,i)*( b3*(1.d0+dudx) + b1*dudz ) )
810  r2 = r2 + &
811  ( s11(3,i)*b1*dvdx + s22(3,i)*b2*(1.d0+dvdy) + s33(3,i)*b3*dvdz &
812  + s12(3,i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
813  + s23(3,i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
814  + s13(3,i)*( b3*dvdx + b1*dvdz ) )
815  r3 = r3 + &
816  ( s11(3,i)*b1*dwdx + s22(3,i)*b2*dwdy + s33(3,i)*b3*(1.d0+dwdz) &
817  + s12(3,i)*( b2*dwdx + b1*dwdy ) &
818  + s23(3,i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
819  + s13(3,i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
820 
821  r4 = r4 + &
822  ( s11(3,i)*b4*(1.d0+dudx) + s22(3,i)*b5*dudy + s33(3,i)*b6*dudz &
823  + s12(3,i)*( b5*(1.d0+dudx) + b4*dudy ) &
824  + s23(3,i)*( b6*dudy + b5*dudz ) &
825  + s13(3,i)*( b6*(1.d0+dudx) + b4*dudz ) )
826  r5 = r5 + &
827  ( s11(3,i)*b4*dvdx + s22(3,i)*b5*(1.d0+dvdy) + s33(3,i)*b6*dvdz &
828  + s12(3,i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
829  + s23(3,i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
830  + s13(3,i)*( b6*dvdx + b4*dvdz ) )
831  r6 = r6 + &
832  ( s11(3,i)*b4*dwdx + s22(3,i)*b5*dwdy + s33(3,i)*b6*(1.d0+dwdz) &
833  + s12(3,i)*( b5*dwdx + b4*dwdy ) &
834  + s23(3,i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
835  + s13(3,i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
836 
837  r7 = r7 + &
838  ( s11(3,i)*b7*(1.d0+dudx) + s22(3,i)*b8*dudy + s33(3,i)*b9*dudz &
839  + s12(3,i)*( b8*(1.d0+dudx) + b7*dudy ) &
840  + s23(3,i)*( b9*dudy + b8*dudz ) &
841  + s13(3,i)*( b9*(1.d0+dudx) + b7*dudz ) )
842  r8 = r8 + &
843  ( s11(3,i)*b7*dvdx + s22(3,i)*b8*(1.d0+dvdy) + s33(3,i)*b9*dvdz &
844  + s12(3,i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
845  + s23(3,i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
846  + s13(3,i)*( b9*dvdx + b7*dvdz ) )
847  r9 = r9 + &
848  ( s11(3,i)*b7*dwdx + s22(3,i)*b8*dwdy + s33(3,i)*b9*(1.d0+dwdz) &
849  + s12(3,i)*( b8*dwdx + b7*dwdy ) &
850  + s23(3,i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
851  + s13(3,i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
852 
853  r10 = r10 + &
854  ( s11(3,i)*b10*(1.d0+dudx) + s22(3,i)*b11*dudy+s33(3,i)*b12*dudz &
855  + s12(3,i)*( b11*(1.d0+dudx) + b10*dudy ) &
856  + s23(3,i)*( b12*dudy + b11*dudz ) &
857  + s13(3,i)*( b12*(1.d0+dudx) + b10*dudz ) )
858  r11 = r11 + &
859  ( s11(3,i)*b10*dvdx + s22(3,i)*b11*(1.d0+dvdy)+s33(3,i)*b12*dvdz &
860  + s12(3,i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
861  + s23(3,i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
862  + s13(3,i)*( b12*dvdx + b10*dvdz ) )
863  r12 = r12 + &
864  ( s11(3,i)*b10*dwdx + s22(3,i)*b11*dwdy+s33(3,i)*b12*(1.d0+dwdz) &
865  + s12(3,i)*( b11*dwdx + b10*dwdy ) &
866  + s23(3,i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
867  + s13(3,i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )
868 
869  r13 = r13 + &
870  ( s11(3,i)*b13*(1.d0+dudx) + s22(3,i)*b14*dudy + s33(3,i)*b15*dudz &
871  + s12(3,i)*( b14*(1.d0+dudx) + b13*dudy ) &
872  + s23(3,i)*( b15*dudy + b14*dudz ) &
873  + s13(3,i)*( b15*(1.d0+dudx) + b13*dudz ) )
874  r14 = r14 + &
875  ( s11(3,i)*b13*dvdx + s22(3,i)*b14*(1.d0+dvdy) + s33(3,i)*b15*dvdz &
876  + s12(3,i)*( b13*(1.d0+dvdy) + b14*dvdx ) &
877  + s23(3,i)*( b15*(1.d0+dvdy) + b14*dvdz ) &
878  + s13(3,i)*( b15*dvdx + b13*dvdz ) )
879  r15 = r15 + &
880  ( s11(3,i)*b13*dwdx + s22(3,i)*b14*dwdy + s33(3,i)*b15*(1.d0+dwdz) &
881  + s12(3,i)*( b14*dwdx + b13*dwdy ) &
882  + s23(3,i)*( b15*dwdy + b14*(1.d0 + dwdz) ) &
883  + s13(3,i)*( b15*dwdx + b13*(1.d0 + dwdz) ) )
884 
885  r16 = r16 + &
886  ( s11(3,i)*b16*(1.d0+dudx) + s22(3,i)*b17*dudy + s33(3,i)*b18*dudz &
887  + s12(3,i)*( b17*(1.d0+dudx) + b16*dudy ) &
888  + s23(3,i)*( b18*dudy + b17*dudz ) &
889  + s13(3,i)*( b18*(1.d0+dudx) + b16*dudz ) )
890  r17 = r17 + &
891  ( s11(3,i)*b16*dvdx + s22(3,i)*b17*(1.d0+dvdy) + s33(3,i)*b18*dvdz &
892  + s12(3,i)*( b16*(1.d0+dvdy) + b17*dvdx ) &
893  + s23(3,i)*( b18*(1.d0+dvdy) + b17*dvdz ) &
894  + s13(3,i)*( b18*dvdx + b16*dvdz ) )
895  r18 = r18 + &
896  ( s11(3,i)*b16*dwdx + s22(3,i)*b17*dwdy + s33(3,i)*b18*(1.d0+dwdz) &
897  + s12(3,i)*( b17*dwdx + b16*dwdy ) &
898  + s23(3,i)*( b18*dwdy + b17*(1.d0 + dwdz) ) &
899  + s13(3,i)*( b18*dwdx + b16*(1.d0 + dwdz) ) )
900 
901  r19 = r19 + &
902  ( s11(3,i)*b19*(1.d0+dudx) + s22(3,i)*b20*dudy + s33(3,i)*b21*dudz &
903  + s12(3,i)*( b20*(1.d0+dudx) + b19*dudy ) &
904  + s23(3,i)*( b21*dudy + b20*dudz ) &
905  + s13(3,i)*( b21*(1.d0+dudx) + b19*dudz ) )
906  r20 = r20 + &
907  ( s11(3,i)*b19*dvdx + s22(3,i)*b20*(1.d0+dvdy) + s33(3,i)*b21*dvdz &
908  + s12(3,i)*( b19*(1.d0+dvdy) + b20*dvdx ) &
909  + s23(3,i)*( b21*(1.d0+dvdy) + b20*dvdz ) &
910  + s13(3,i)*( b21*dvdx + b19*dvdz ) )
911  r21 = r21 + &
912  ( s11(3,i)*b19*dwdx + s22(3,i)*b20*dwdy + s33(3,i)*b21*(1.d0+dwdz) &
913  + s12(3,i)*( b20*dwdx + b19*dwdy ) &
914  + s23(3,i)*( b21*dwdy + b20*(1.d0 + dwdz) ) &
915  + s13(3,i)*( b21*dwdx + b19*(1.d0 + dwdz) ) )
916 
917  r22 = r22 + &
918  ( s11(3,i)*b22*(1.d0+dudx) + s22(3,i)*b23*dudy+s33(3,i)*b24*dudz &
919  + s12(3,i)*( b23*(1.d0+dudx) + b22*dudy ) &
920  + s23(3,i)*( b24*dudy + b23*dudz ) &
921  + s13(3,i)*( b24*(1.d0+dudx) + b22*dudz ) )
922  r23 = r23 + &
923  ( s11(3,i)*b22*dvdx + s22(3,i)*b23*(1.d0+dvdy)+s33(3,i)*b24*dvdz &
924  + s12(3,i)*( b22*(1.d0+dvdy) + b23*dvdx ) &
925  + s23(3,i)*( b24*(1.d0+dvdy) + b23*dvdz ) &
926  + s13(3,i)*( b24*dvdx + b22*dvdz ) )
927  r24 = r24 + &
928  ( s11(3,i)*b22*dwdx + s22(3,i)*b23*dwdy+s33(3,i)*b24*(1.d0+dwdz) &
929  + s12(3,i)*( b23*dwdx + b22*dwdy ) &
930  + s23(3,i)*( b24*dwdy + b23*(1.d0 + dwdz) ) &
931  + s13(3,i)*( b24*dwdx + b22*(1.d0 + dwdz) ) )
932 
933  r25 = r25 + &
934  ( s11(3,i)*b25*(1.d0+dudx) + s22(3,i)*b26*dudy+s33(3,i)*b27*dudz &
935  + s12(3,i)*( b26*(1.d0+dudx) + b25*dudy ) &
936  + s23(3,i)*( b27*dudy + b26*dudz ) &
937  + s13(3,i)*( b27*(1.d0+dudx) + b25*dudz ) )
938  r26 = r26 + &
939  ( s11(3,i)*b25*dvdx + s22(3,i)*b26*(1.d0+dvdy)+s33(3,i)*b27*dvdz &
940  + s12(3,i)*( b25*(1.d0+dvdy) + b26*dvdx ) &
941  + s23(3,i)*( b27*(1.d0+dvdy) + b26*dvdz ) &
942  + s13(3,i)*( b27*dvdx + b25*dvdz ) )
943  r27 = r27 + &
944  ( s11(3,i)*b25*dwdx + s22(3,i)*b26*dwdy+s33(3,i)*b27*(1.d0+dwdz) &
945  + s12(3,i)*( b26*dwdx + b25*dwdy ) &
946  + s23(3,i)*( b27*dwdy + b26*(1.d0 + dwdz) ) &
947  + s13(3,i)*( b27*dwdx + b25*(1.d0 + dwdz) ) )
948 
949  r28 = r28 + &
950  ( s11(3,i)*b28*(1.d0+dudx) + s22(3,i)*b29*dudy+s33(3,i)*b30*dudz &
951  + s12(3,i)*( b29*(1.d0+dudx) + b28*dudy ) &
952  + s23(3,i)*( b30*dudy + b29*dudz ) &
953  + s13(3,i)*( b30*(1.d0+dudx) + b28*dudz ) )
954  r29 = r29 + &
955  ( s11(3,i)*b28*dvdx + s22(3,i)*b29*(1.d0+dvdy)+s33(3,i)*b30*dvdz &
956  + s12(3,i)*( b28*(1.d0+dvdy) + b29*dvdx ) &
957  + s23(3,i)*( b30*(1.d0+dvdy) + b29*dvdz ) &
958  + s13(3,i)*( b30*dvdx + b28*dvdz ) )
959  r30 = r30 + &
960  ( s11(3,i)*b28*dwdx + s22(3,i)*b29*dwdy+s33(3,i)*b30*(1.d0+dwdz) &
961  + s12(3,i)*( b29*dwdx + b28*dwdy ) &
962  + s23(3,i)*( b30*dwdy + b29*(1.d0 + dwdz) ) &
963  + s13(3,i)*( b30*dwdx + b28*(1.d0 + dwdz) ) )
964 
965 !-- Gauss Point 4
966 
967  g1 = 0.13819660d0
968  g2 = 0.13819660d0
969  g3 = 0.13819660d0
970  g4 = 0.58541020d0
971  xn1 = (4.d0*g1-1.d0) ! derivative of shape function
972  xn2 = (4.d0*g2-1.d0) ! dN_i/dzeta_i
973  xn3 = (4.d0*g3-1.d0)
974  xn4 = (4.d0*g4-1.d0)
975 ! xN5 = 4.d0*g1*g2
976 ! xN6 = 4.d0*g2*g3
977 ! xN7 = 4.d0*g3*g1
978 ! xN8 = 4.d0*g1*g4
979 ! xN9 = 4.d0*g2*g4
980 ! xN10= 4.d0*g3*g4
981 
982  b1 = aux1*xn1
983  b2 = aux2*xn1
984  b3 = aux3*xn1
985  b4 = aux4*xn2
986  b5 = aux5*xn2
987  b6 = aux6*xn2
988  b7 = aux7*xn3
989  b8 = aux8*xn3
990  b9 = aux9*xn3
991  b10 = aux10*xn4
992  b11 = aux11*xn4
993  b12 = aux12*xn4
994 
995  b13 = 4.d0*(g2*aux1 + g1*aux4)
996  b14 = 4.d0*(g2*aux2 + g1*aux5)
997  b15 = 4.d0*(g2*aux3 + g1*aux6)
998 
999  b16 = 4.d0*(g3*aux4 + g2*aux7)
1000  b17 = 4.d0*(g3*aux5 + g2*aux8)
1001  b18 = 4.d0*(g3*aux6 + g2*aux9)
1002 
1003  b19 = 4.d0*(g1*aux7 + g3*aux1)
1004  b20 = 4.d0*(g1*aux8 + g3*aux2)
1005  b21 = 4.d0*(g1*aux9 + g3*aux3)
1006 
1007  b22 = 4.d0*(g4*aux1 + g1*aux10)
1008  b23 = 4.d0*(g4*aux2 + g1*aux11)
1009  b24 = 4.d0*(g4*aux3 + g1*aux12)
1010 
1011  b25 = 4.d0*(g4*aux4 + g2*aux10)
1012  b26 = 4.d0*(g4*aux5 + g2*aux11)
1013  b27 = 4.d0*(g4*aux6 + g2*aux12)
1014 
1015  b28 = 4.d0*(g4*aux7 + g3*aux10)
1016  b29 = 4.d0*(g4*aux8 + g3*aux11)
1017  b30 = 4.d0*(g4*aux9 + g3*aux12)
1018 !-----Calculate displacement gradient (H)
1019  dudx = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10)*vx6inv
1020  dvdy = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10)*vx6inv
1021  dwdz = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10)*vx6inv
1022  dudy = (b2*u1 + b5*u2 + b8*u3 + b11*u4 + b14*u5 + b17*u6 + b20*u7 + b23*u8 + b26*u9 + b29*u10)*vx6inv
1023  dvdx = (b1*v1 + b4*v2 + b7*v3 + b10*v4 + b13*v5 + b16*v6 + b19*v7 + b22*v8 + b25*v9 + b28*v10)*vx6inv
1024  dvdz = (b3*v1 + b6*v2 + b9*v3 + b12*v4 + b15*v5 + b18*v6 + b21*v7 + b24*v8 + b27*v9 + b30*v10)*vx6inv
1025  dwdy = (b2*w1 + b5*w2 + b8*w3 + b11*w4 + b14*w5 + b17*w6 + b20*w7 + b23*w8 + b26*w9 + b29*w10)*vx6inv
1026  dudz = (b3*u1 + b6*u2 + b9*u3 + b12*u4 + b15*u5 + b18*u6 + b21*u7 + b24*u8 + b27*u9 + b30*u10)*vx6inv
1027  dwdx = (b1*w1 + b4*w2 + b7*w3 + b10*w4 + b13*w5 + b16*w6 + b19*w7 + b22*w8 + b25*w9 + b28*w10)*vx6inv
1028 !
1029 ! Green-Lagrange Strain (E)
1030 !
1031  e11 = dudx + 0.5d0*(dudx*dudx + dvdx*dvdx + dwdx*dwdx)
1032  e22 = dvdy + 0.5d0*(dudy*dudy + dvdy*dvdy + dwdy*dwdy)
1033  e33 = dwdz + 0.5d0*(dudz*dudz + dvdz*dvdz + dwdz*dwdz)
1034  e12 = dudy + dvdx + ( dudx*dudy + dvdx*dvdy + dwdx*dwdy )
1035  e23 = dvdz + dwdy + ( dudy*dudz + dvdy*dvdz + dwdy*dwdz )
1036  e13 = dwdx + dudz + ( dudz*dudx + dvdz*dvdx + dwdz*dwdx )
1037 !
1038 ! Second Piola-Kirchhoff stress (S)
1039 !
1040  s11(4,i) = e11*ci(1,j) + e22*ci(2,j) + e33*ci(4,j)
1041  s22(4,i) = e11*ci(2,j) + e22*ci(3,j) + e33*ci(5,j)
1042  s33(4,i) = e11*ci(4,j) + e22*ci(5,j) + e33*ci(6,j)
1043  s12(4,i) = e12*ci(7,j)
1044  s23(4,i) = e23*ci(8,j)
1045  s13(4,i) = e13*ci(9,j)
1046 
1047  r1 = r1 + &
1048  ( s11(4,i)*b1*(1.d0+dudx) + s22(4,i)*b2*dudy + s33(4,i)*b3*dudz &
1049  + s12(4,i)*( b2*(1.d0+dudx) + b1*dudy ) &
1050  + s23(4,i)*( b3*dudy + b2*dudz ) &
1051  + s13(4,i)*( b3*(1.d0+dudx) + b1*dudz ) )
1052  r2 = r2 + &
1053  ( s11(4,i)*b1*dvdx + s22(4,i)*b2*(1.d0+dvdy) + s33(4,i)*b3*dvdz &
1054  + s12(4,i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
1055  + s23(4,i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
1056  + s13(4,i)*( b3*dvdx + b1*dvdz ) )
1057  r3 = r3 + &
1058  ( s11(4,i)*b1*dwdx + s22(4,i)*b2*dwdy + s33(4,i)*b3*(1.d0+dwdz) &
1059  + s12(4,i)*( b2*dwdx + b1*dwdy ) &
1060  + s23(4,i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
1061  + s13(4,i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
1062 
1063  r4 = r4 + &
1064  ( s11(4,i)*b4*(1.d0+dudx) + s22(4,i)*b5*dudy + s33(4,i)*b6*dudz &
1065  + s12(4,i)*( b5*(1.d0+dudx) + b4*dudy ) &
1066  + s23(4,i)*( b6*dudy + b5*dudz ) &
1067  + s13(4,i)*( b6*(1.d0+dudx) + b4*dudz ) )
1068  r5 = r5 + &
1069  ( s11(4,i)*b4*dvdx + s22(4,i)*b5*(1.d0+dvdy) + s33(4,i)*b6*dvdz &
1070  + s12(4,i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
1071  + s23(4,i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
1072  + s13(4,i)*( b6*dvdx + b4*dvdz ) )
1073  r6 = r6 + &
1074  ( s11(4,i)*b4*dwdx + s22(4,i)*b5*dwdy + s33(4,i)*b6*(1.d0+dwdz) &
1075  + s12(4,i)*( b5*dwdx + b4*dwdy ) &
1076  + s23(4,i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
1077  + s13(4,i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
1078 
1079  r7 = r7 + &
1080  ( s11(4,i)*b7*(1.d0+dudx) + s22(4,i)*b8*dudy + s33(4,i)*b9*dudz &
1081  + s12(4,i)*( b8*(1.d0+dudx) + b7*dudy ) &
1082  + s23(4,i)*( b9*dudy + b8*dudz ) &
1083  + s13(4,i)*( b9*(1.d0+dudx) + b7*dudz ) )
1084  r8 = r8 + &
1085  ( s11(4,i)*b7*dvdx + s22(4,i)*b8*(1.d0+dvdy) + s33(4,i)*b9*dvdz &
1086  + s12(4,i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
1087  + s23(4,i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
1088  + s13(4,i)*( b9*dvdx + b7*dvdz ) )
1089  r9 = r9 + &
1090  ( s11(4,i)*b7*dwdx + s22(4,i)*b8*dwdy + s33(4,i)*b9*(1.d0+dwdz) &
1091  + s12(4,i)*( b8*dwdx + b7*dwdy ) &
1092  + s23(4,i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
1093  + s13(4,i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
1094 
1095  r10 = r10 + &
1096  ( s11(4,i)*b10*(1.d0+dudx) + s22(4,i)*b11*dudy+s33(4,i)*b12*dudz &
1097  + s12(4,i)*( b11*(1.d0+dudx) + b10*dudy ) &
1098  + s23(4,i)*( b12*dudy + b11*dudz ) &
1099  + s13(4,i)*( b12*(1.d0+dudx) + b10*dudz ) )
1100  r11 = r11 + &
1101  ( s11(4,i)*b10*dvdx + s22(4,i)*b11*(1.d0+dvdy)+s33(4,i)*b12*dvdz &
1102  + s12(4,i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
1103  + s23(4,i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
1104  + s13(4,i)*( b12*dvdx + b10*dvdz ) )
1105  r12 = r12 + &
1106  ( s11(4,i)*b10*dwdx + s22(4,i)*b11*dwdy+s33(4,i)*b12*(1.d0+dwdz) &
1107  + s12(4,i)*( b11*dwdx + b10*dwdy ) &
1108  + s23(4,i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
1109  + s13(4,i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )
1110 
1111  r13 = r13 + &
1112  ( s11(4,i)*b13*(1.d0+dudx) + s22(4,i)*b14*dudy + s33(4,i)*b15*dudz &
1113  + s12(4,i)*( b14*(1.d0+dudx) + b13*dudy ) &
1114  + s23(4,i)*( b15*dudy + b14*dudz ) &
1115  + s13(4,i)*( b15*(1.d0+dudx) + b13*dudz ) )
1116  r14 = r14 + &
1117  ( s11(4,i)*b13*dvdx + s22(4,i)*b14*(1.d0+dvdy) + s33(4,i)*b15*dvdz &
1118  + s12(4,i)*( b13*(1.d0+dvdy) + b14*dvdx ) &
1119  + s23(4,i)*( b15*(1.d0+dvdy) + b14*dvdz ) &
1120  + s13(4,i)*( b15*dvdx + b13*dvdz ) )
1121  r15 = r15 + &
1122  ( s11(4,i)*b13*dwdx + s22(4,i)*b14*dwdy + s33(4,i)*b15*(1.d0+dwdz) &
1123  + s12(4,i)*( b14*dwdx + b13*dwdy ) &
1124  + s23(4,i)*( b15*dwdy + b14*(1.d0 + dwdz) ) &
1125  + s13(4,i)*( b15*dwdx + b13*(1.d0 + dwdz) ) )
1126 
1127  r16 = r16 + &
1128  ( s11(4,i)*b16*(1.d0+dudx) + s22(4,i)*b17*dudy + s33(4,i)*b18*dudz &
1129  + s12(4,i)*( b17*(1.d0+dudx) + b16*dudy ) &
1130  + s23(4,i)*( b18*dudy + b17*dudz ) &
1131  + s13(4,i)*( b18*(1.d0+dudx) + b16*dudz ) )
1132  r17 = r17 + &
1133  ( s11(4,i)*b16*dvdx + s22(4,i)*b17*(1.d0+dvdy) + s33(4,i)*b18*dvdz &
1134  + s12(4,i)*( b16*(1.d0+dvdy) + b17*dvdx ) &
1135  + s23(4,i)*( b18*(1.d0+dvdy) + b17*dvdz ) &
1136  + s13(4,i)*( b18*dvdx + b16*dvdz ) )
1137  r18 = r18 + &
1138  ( s11(4,i)*b16*dwdx + s22(4,i)*b17*dwdy + s33(4,i)*b18*(1.d0+dwdz) &
1139  + s12(4,i)*( b17*dwdx + b16*dwdy ) &
1140  + s23(4,i)*( b18*dwdy + b17*(1.d0 + dwdz) ) &
1141  + s13(4,i)*( b18*dwdx + b16*(1.d0 + dwdz) ) )
1142 
1143  r19 = r19 + &
1144  ( s11(4,i)*b19*(1.d0+dudx) + s22(4,i)*b20*dudy + s33(4,i)*b21*dudz &
1145  + s12(4,i)*( b20*(1.d0+dudx) + b19*dudy ) &
1146  + s23(4,i)*( b21*dudy + b20*dudz ) &
1147  + s13(4,i)*( b21*(1.d0+dudx) + b19*dudz ) )
1148  r20 = r20 + &
1149  ( s11(4,i)*b19*dvdx + s22(4,i)*b20*(1.d0+dvdy) + s33(4,i)*b21*dvdz &
1150  + s12(4,i)*( b19*(1.d0+dvdy) + b20*dvdx ) &
1151  + s23(4,i)*( b21*(1.d0+dvdy) + b20*dvdz ) &
1152  + s13(4,i)*( b21*dvdx + b19*dvdz ) )
1153  r21 = r21 + &
1154  ( s11(4,i)*b19*dwdx + s22(4,i)*b20*dwdy + s33(4,i)*b21*(1.d0+dwdz) &
1155  + s12(4,i)*( b20*dwdx + b19*dwdy ) &
1156  + s23(4,i)*( b21*dwdy + b20*(1.d0 + dwdz) ) &
1157  + s13(4,i)*( b21*dwdx + b19*(1.d0 + dwdz) ) )
1158 
1159  r22 = r22 + &
1160  ( s11(4,i)*b22*(1.d0+dudx) + s22(4,i)*b23*dudy+s33(4,i)*b24*dudz &
1161  + s12(4,i)*( b23*(1.d0+dudx) + b22*dudy ) &
1162  + s23(4,i)*( b24*dudy + b23*dudz ) &
1163  + s13(4,i)*( b24*(1.d0+dudx) + b22*dudz ) )
1164  r23 = r23 + &
1165  ( s11(4,i)*b22*dvdx + s22(4,i)*b23*(1.d0+dvdy)+s33(4,i)*b24*dvdz &
1166  + s12(4,i)*( b22*(1.d0+dvdy) + b23*dvdx ) &
1167  + s23(4,i)*( b24*(1.d0+dvdy) + b23*dvdz ) &
1168  + s13(4,i)*( b24*dvdx + b22*dvdz ) )
1169  r24 = r24 + &
1170  ( s11(4,i)*b22*dwdx + s22(4,i)*b23*dwdy+s33(4,i)*b24*(1.d0+dwdz) &
1171  + s12(4,i)*( b23*dwdx + b22*dwdy ) &
1172  + s23(4,i)*( b24*dwdy + b23*(1.d0 + dwdz) ) &
1173  + s13(4,i)*( b24*dwdx + b22*(1.d0 + dwdz) ) )
1174 
1175  r25 = r25 + &
1176  ( s11(4,i)*b25*(1.d0+dudx) + s22(4,i)*b26*dudy+s33(4,i)*b27*dudz &
1177  + s12(4,i)*( b26*(1.d0+dudx) + b25*dudy ) &
1178  + s23(4,i)*( b27*dudy + b26*dudz ) &
1179  + s13(4,i)*( b27*(1.d0+dudx) + b25*dudz ) )
1180  r26 = r26 + &
1181  ( s11(4,i)*b25*dvdx + s22(4,i)*b26*(1.d0+dvdy)+s33(4,i)*b27*dvdz &
1182  + s12(4,i)*( b25*(1.d0+dvdy) + b26*dvdx ) &
1183  + s23(4,i)*( b27*(1.d0+dvdy) + b26*dvdz ) &
1184  + s13(4,i)*( b27*dvdx + b25*dvdz ) )
1185  r27 = r27 + &
1186  ( s11(4,i)*b25*dwdx + s22(4,i)*b26*dwdy+s33(4,i)*b27*(1.d0+dwdz) &
1187  + s12(4,i)*( b26*dwdx + b25*dwdy ) &
1188  + s23(4,i)*( b27*dwdy + b26*(1.d0 + dwdz) ) &
1189  + s13(4,i)*( b27*dwdx + b25*(1.d0 + dwdz) ) )
1190 
1191  r28 = r28 + &
1192  ( s11(4,i)*b28*(1.d0+dudx) + s22(4,i)*b29*dudy+s33(4,i)*b30*dudz &
1193  + s12(4,i)*( b29*(1.d0+dudx) + b28*dudy ) &
1194  + s23(4,i)*( b30*dudy + b29*dudz ) &
1195  + s13(4,i)*( b30*(1.d0+dudx) + b28*dudz ) )
1196  r29 = r29 + &
1197  ( s11(4,i)*b28*dvdx + s22(4,i)*b29*(1.d0+dvdy)+s33(4,i)*b30*dvdz &
1198  + s12(4,i)*( b28*(1.d0+dvdy) + b29*dvdx ) &
1199  + s23(4,i)*( b30*(1.d0+dvdy) + b29*dvdz ) &
1200  + s13(4,i)*( b30*dvdx + b28*dvdz ) )
1201  r30 = r30 + &
1202  ( s11(4,i)*b28*dwdx + s22(4,i)*b29*dwdy+s33(4,i)*b30*(1.d0+dwdz) &
1203  + s12(4,i)*( b29*dwdx + b28*dwdy ) &
1204  + s23(4,i)*( b30*dwdy + b29*(1.d0 + dwdz) ) &
1205  + s13(4,i)*( b30*dwdx + b28*(1.d0 + dwdz) ) )
1206 
1207 ! Wi (i.e. weight) for 4 guass point integration is 1/4
1208 
1209 ! Wi * 1/6 because the volume of a reference tetrahedra in
1210 ! volume coordinates is 1/6
1211 
1212 ! ASSEMBLE THE INTERNAL FORCE VECTOR
1213 !
1214 ! local node 1
1215  r_in(k1n1) = r_in(k1n1) - r1*0.04166666666666667d0
1216  r_in(k2n1) = r_in(k2n1) - r2*0.04166666666666667d0
1217  r_in(k3n1) = r_in(k3n1) - r3*0.04166666666666667d0
1218 ! local node 2
1219  r_in(k1n2) = r_in(k1n2) - r4*0.04166666666666667d0
1220  r_in(k2n2) = r_in(k2n2) - r5*0.04166666666666667d0
1221  r_in(k3n2) = r_in(k3n2) - r6*0.04166666666666667d0
1222 ! local node 3
1223  r_in(k1n3) = r_in(k1n3) - r7*0.04166666666666667d0
1224  r_in(k2n3) = r_in(k2n3) - r8*0.04166666666666667d0
1225  r_in(k3n3) = r_in(k3n3) - r9*0.04166666666666667d0
1226 ! local node 4
1227  r_in(k1n4) = r_in(k1n4) - r10*0.04166666666666667d0
1228  r_in(k2n4) = r_in(k2n4) - r11*0.04166666666666667d0
1229  r_in(k3n4) = r_in(k3n4) - r12*0.04166666666666667d0
1230 ! local node 5
1231  r_in(k1n5) = r_in(k1n5) - r13*0.04166666666666667d0
1232  r_in(k2n5) = r_in(k2n5) - r14*0.04166666666666667d0
1233  r_in(k3n5) = r_in(k3n5) - r15*0.04166666666666667d0
1234 ! local node 6
1235  r_in(k1n6) = r_in(k1n6) - r16*0.04166666666666667d0
1236  r_in(k2n6) = r_in(k2n6) - r17*0.04166666666666667d0
1237  r_in(k3n6) = r_in(k3n6) - r18*0.04166666666666667d0
1238 ! local node 7
1239  r_in(k1n7) = r_in(k1n7) - r19*0.04166666666666667d0
1240  r_in(k2n7) = r_in(k2n7) - r20*0.04166666666666667d0
1241  r_in(k3n7) = r_in(k3n7) - r21*0.04166666666666667d0
1242 ! local node 8
1243  r_in(k1n8) = r_in(k1n8) - r22*0.04166666666666667d0
1244  r_in(k2n8) = r_in(k2n8) - r23*0.04166666666666667d0
1245  r_in(k3n8) = r_in(k3n8) - r24*0.04166666666666667d0
1246 ! local node 9
1247  r_in(k1n9) = r_in(k1n9) - r25*0.04166666666666667d0
1248  r_in(k2n9) = r_in(k2n9) - r26*0.04166666666666667d0
1249  r_in(k3n9) = r_in(k3n9) - r27*0.04166666666666667d0
1250 ! local node 10
1251  r_in(k1n10) = r_in(k1n10) - r28*0.04166666666666667d0
1252  r_in(k2n10) = r_in(k2n10) - r29*0.04166666666666667d0
1253  r_in(k3n10) = r_in(k3n10) - r30*0.04166666666666667d0
1254 
1255 ! Reduced integration
1256 
1257 !-- Gauss Point 1
1258 
1259  b1 = 0.d0
1260  b2 = 0.d0
1261  b3 = 0.d0
1262  b4 = 0.d0
1263  b5 = 0.d0
1264  b6 = 0.d0
1265  b7 = 0.d0
1266  b8 = 0.d0
1267  b9 = 0.d0
1268  b10 = 0.d0
1269  b11 = 0.d0
1270  b12 = 0.d0
1271 
1272  b13 = aux1 + aux4
1273  b14 = aux2 + aux5
1274  b15 = aux3 + aux6
1275 
1276  b16 = aux4 + aux7
1277  b17 = aux5 + aux8
1278  b18 = aux6 + aux9
1279 
1280  b19 = aux7 + aux1
1281  b20 = aux8 + aux2
1282  b21 = aux9 + aux3
1283 
1284  b22 = aux1 + aux10
1285  b23 = aux2 + aux11
1286  b24 = aux3 + aux12
1287 
1288  b25 = aux4 + aux10
1289  b26 = aux5 + aux11
1290  b27 = aux6 + aux12
1291 
1292  b28 = aux7 + aux10
1293  b29 = aux8 + aux11
1294  b30 = aux9 + aux12
1295 
1296 
1297 !-----Calculate displacement gradient (H)
1298  dudx = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10)*vx6inv
1299  dvdy = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10)*vx6inv
1300  dwdz = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10)*vx6inv
1301  dudy = (b2*u1 + b5*u2 + b8*u3 + b11*u4 + b14*u5 + b17*u6 + b20*u7 + b23*u8 + b26*u9 + b29*u10)*vx6inv
1302  dvdx = (b1*v1 + b4*v2 + b7*v3 + b10*v4 + b13*v5 + b16*v6 + b19*v7 + b22*v8 + b25*v9 + b28*v10)*vx6inv
1303  dvdz = (b3*v1 + b6*v2 + b9*v3 + b12*v4 + b15*v5 + b18*v6 + b21*v7 + b24*v8 + b27*v9 + b30*v10)*vx6inv
1304  dwdy = (b2*w1 + b5*w2 + b8*w3 + b11*w4 + b14*w5 + b17*w6 + b20*w7 + b23*w8 + b26*w9 + b29*w10)*vx6inv
1305  dudz = (b3*u1 + b6*u2 + b9*u3 + b12*u4 + b15*u5 + b18*u6 + b21*u7 + b24*u8 + b27*u9 + b30*u10)*vx6inv
1306  dwdx = (b1*w1 + b4*w2 + b7*w3 + b10*w4 + b13*w5 + b16*w6 + b19*w7 + b22*w8 + b25*w9 + b28*w10)*vx6inv
1307 !
1308 ! Green-Lagrange Strain (E)
1309 !
1310  e11 = dudx + 0.5d0*(dudx*dudx + dvdx*dvdx + dwdx*dwdx)
1311  e22 = dvdy + 0.5d0*(dudy*dudy + dvdy*dvdy + dwdy*dwdy)
1312  e33 = dwdz + 0.5d0*(dudz*dudz + dvdz*dvdz + dwdz*dwdz)
1313  e12 = dudy + dvdx + ( dudx*dudy + dvdx*dvdy + dwdx*dwdy )
1314  e23 = dvdz + dwdy + ( dudy*dudz + dvdy*dvdz + dwdy*dwdz )
1315  e13 = dwdx + dudz + ( dudz*dudx + dvdz*dvdx + dwdz*dwdx )
1316 !
1317 ! Second Piola-Kirchhoff stress (S)
1318 !
1319  s11(1,i) = e11*cj(1,j) + e22*cj(2,j) + e33*cj(4,j)
1320  s22(1,i) = e11*cj(2,j) + e22*cj(3,j) + e33*cj(5,j)
1321  s33(1,i) = e11*cj(4,j) + e22*cj(5,j) + e33*cj(6,j)
1322  s12(1,i) = e12*cj(7,j)
1323  s23(1,i) = e23*cj(8,j)
1324  s13(1,i) = e13*cj(9,j)
1325 
1326  r1 = &
1327  ( s11(1,i)*b1*(1.d0+dudx) + s22(1,i)*b2*dudy + s33(1,i)*b3*dudz &
1328  + s12(1,i)*( b2*(1.d0+dudx) + b1*dudy ) &
1329  + s23(1,i)*( b3*dudy + b2*dudz ) &
1330  + s13(1,i)*( b3*(1.d0+dudx) + b1*dudz ) )
1331  r2 = &
1332  ( s11(1,i)*b1*dvdx + s22(1,i)*b2*(1.d0+dvdy) + s33(1,i)*b3*dvdz &
1333  + s12(1,i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
1334  + s23(1,i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
1335  + s13(1,i)*( b3*dvdx + b1*dvdz ) )
1336  r3 = &
1337  ( s11(1,i)*b1*dwdx + s22(1,i)*b2*dwdy + s33(1,i)*b3*(1.d0+dwdz) &
1338  + s12(1,i)*( b2*dwdx + b1*dwdy ) &
1339  + s23(1,i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
1340  + s13(1,i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
1341 
1342  r4 = &
1343  ( s11(1,i)*b4*(1.d0+dudx) + s22(1,i)*b5*dudy + s33(1,i)*b6*dudz &
1344  + s12(1,i)*( b5*(1.d0+dudx) + b4*dudy ) &
1345  + s23(1,i)*( b6*dudy + b5*dudz ) &
1346  + s13(1,i)*( b6*(1.d0+dudx) + b4*dudz ) )
1347  r5 = &
1348  ( s11(1,i)*b4*dvdx + s22(1,i)*b5*(1.d0+dvdy) + s33(1,i)*b6*dvdz &
1349  + s12(1,i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
1350  + s23(1,i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
1351  + s13(1,i)*( b6*dvdx + b4*dvdz ) )
1352  r6 = &
1353  ( s11(1,i)*b4*dwdx + s22(1,i)*b5*dwdy + s33(1,i)*b6*(1.d0+dwdz) &
1354  + s12(1,i)*( b5*dwdx + b4*dwdy ) &
1355  + s23(1,i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
1356  + s13(1,i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
1357 
1358  r7 = &
1359  ( s11(1,i)*b7*(1.d0+dudx) + s22(1,i)*b8*dudy + s33(1,i)*b9*dudz &
1360  + s12(1,i)*( b8*(1.d0+dudx) + b7*dudy ) &
1361  + s23(1,i)*( b9*dudy + b8*dudz ) &
1362  + s13(1,i)*( b9*(1.d0+dudx) + b7*dudz ) )
1363  r8 = &
1364  ( s11(1,i)*b7*dvdx + s22(1,i)*b8*(1.d0+dvdy) + s33(1,i)*b9*dvdz &
1365  + s12(1,i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
1366  + s23(1,i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
1367  + s13(1,i)*( b9*dvdx + b7*dvdz ) )
1368  r9 = &
1369  ( s11(1,i)*b7*dwdx + s22(1,i)*b8*dwdy + s33(1,i)*b9*(1.d0+dwdz) &
1370  + s12(1,i)*( b8*dwdx + b7*dwdy ) &
1371  + s23(1,i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
1372  + s13(1,i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
1373 
1374  r10 = &
1375  ( s11(1,i)*b10*(1.d0+dudx) + s22(1,i)*b11*dudy+s33(1,i)*b12*dudz &
1376  + s12(1,i)*( b11*(1.d0+dudx) + b10*dudy ) &
1377  + s23(1,i)*( b12*dudy + b11*dudz ) &
1378  + s13(1,i)*( b12*(1.d0+dudx) + b10*dudz ) )
1379  r11 = &
1380  ( s11(1,i)*b10*dvdx + s22(1,i)*b11*(1.d0+dvdy)+s33(1,i)*b12*dvdz &
1381  + s12(1,i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
1382  + s23(1,i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
1383  + s13(1,i)*( b12*dvdx + b10*dvdz ) )
1384  r12 = &
1385  ( s11(1,i)*b10*dwdx + s22(1,i)*b11*dwdy+s33(1,i)*b12*(1.d0+dwdz) &
1386  + s12(1,i)*( b11*dwdx + b10*dwdy ) &
1387  + s23(1,i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
1388  + s13(1,i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )
1389 
1390  r13 = &
1391  ( s11(1,i)*b13*(1.d0+dudx) + s22(1,i)*b14*dudy + s33(1,i)*b15*dudz &
1392  + s12(1,i)*( b14*(1.d0+dudx) + b13*dudy ) &
1393  + s23(1,i)*( b15*dudy + b14*dudz ) &
1394  + s13(1,i)*( b15*(1.d0+dudx) + b13*dudz ) )
1395  r14 = &
1396  ( s11(1,i)*b13*dvdx + s22(1,i)*b14*(1.d0+dvdy) + s33(1,i)*b15*dvdz &
1397  + s12(1,i)*( b13*(1.d0+dvdy) + b14*dvdx ) &
1398  + s23(1,i)*( b15*(1.d0+dvdy) + b14*dvdz ) &
1399  + s13(1,i)*( b15*dvdx + b13*dvdz ) )
1400  r15 = &
1401  ( s11(1,i)*b13*dwdx + s22(1,i)*b14*dwdy + s33(1,i)*b15*(1.d0+dwdz) &
1402  + s12(1,i)*( b14*dwdx + b13*dwdy ) &
1403  + s23(1,i)*( b15*dwdy + b14*(1.d0 + dwdz) ) &
1404  + s13(1,i)*( b15*dwdx + b13*(1.d0 + dwdz) ) )
1405 
1406  r16 = &
1407  ( s11(1,i)*b16*(1.d0+dudx) + s22(1,i)*b17*dudy + s33(1,i)*b18*dudz &
1408  + s12(1,i)*( b17*(1.d0+dudx) + b16*dudy ) &
1409  + s23(1,i)*( b18*dudy + b17*dudz ) &
1410  + s13(1,i)*( b18*(1.d0+dudx) + b16*dudz ) )
1411  r17 = &
1412  ( s11(1,i)*b16*dvdx + s22(1,i)*b17*(1.d0+dvdy) + s33(1,i)*b18*dvdz &
1413  + s12(1,i)*( b16*(1.d0+dvdy) + b17*dvdx ) &
1414  + s23(1,i)*( b18*(1.d0+dvdy) + b17*dvdz ) &
1415  + s13(1,i)*( b18*dvdx + b16*dvdz ) )
1416  r18 = &
1417  ( s11(1,i)*b16*dwdx + s22(1,i)*b17*dwdy + s33(1,i)*b18*(1.d0+dwdz) &
1418  + s12(1,i)*( b17*dwdx + b16*dwdy ) &
1419  + s23(1,i)*( b18*dwdy + b17*(1.d0 + dwdz) ) &
1420  + s13(1,i)*( b18*dwdx + b16*(1.d0 + dwdz) ) )
1421 
1422  r19 = &
1423  ( s11(1,i)*b19*(1.d0+dudx) + s22(1,i)*b20*dudy + s33(1,i)*b21*dudz &
1424  + s12(1,i)*( b20*(1.d0+dudx) + b19*dudy ) &
1425  + s23(1,i)*( b21*dudy + b20*dudz ) &
1426  + s13(1,i)*( b21*(1.d0+dudx) + b19*dudz ) )
1427  r20 = &
1428  ( s11(1,i)*b19*dvdx + s22(1,i)*b20*(1.d0+dvdy) + s33(1,i)*b21*dvdz &
1429  + s12(1,i)*( b19*(1.d0+dvdy) + b20*dvdx ) &
1430  + s23(1,i)*( b21*(1.d0+dvdy) + b20*dvdz ) &
1431  + s13(1,i)*( b21*dvdx + b19*dvdz ) )
1432  r21 = &
1433  ( s11(1,i)*b19*dwdx + s22(1,i)*b20*dwdy + s33(1,i)*b21*(1.d0+dwdz) &
1434  + s12(1,i)*( b20*dwdx + b19*dwdy ) &
1435  + s23(1,i)*( b21*dwdy + b20*(1.d0 + dwdz) ) &
1436  + s13(1,i)*( b21*dwdx + b19*(1.d0 + dwdz) ) )
1437 
1438  r22 = &
1439  ( s11(1,i)*b22*(1.d0+dudx) + s22(1,i)*b23*dudy+s33(1,i)*b24*dudz &
1440  + s12(1,i)*( b23*(1.d0+dudx) + b22*dudy ) &
1441  + s23(1,i)*( b24*dudy + b23*dudz ) &
1442  + s13(1,i)*( b24*(1.d0+dudx) + b22*dudz ) )
1443  r23 = &
1444  ( s11(1,i)*b22*dvdx + s22(1,i)*b23*(1.d0+dvdy)+s33(1,i)*b24*dvdz &
1445  + s12(1,i)*( b22*(1.d0+dvdy) + b23*dvdx ) &
1446  + s23(1,i)*( b24*(1.d0+dvdy) + b23*dvdz ) &
1447  + s13(1,i)*( b24*dvdx + b22*dvdz ) )
1448  r24 = &
1449  ( s11(1,i)*b22*dwdx + s22(1,i)*b23*dwdy+s33(1,i)*b24*(1.d0+dwdz) &
1450  + s12(1,i)*( b23*dwdx + b22*dwdy ) &
1451  + s23(1,i)*( b24*dwdy + b23*(1.d0 + dwdz) ) &
1452  + s13(1,i)*( b24*dwdx + b22*(1.d0 + dwdz) ) )
1453 
1454  r25 = &
1455  ( s11(1,i)*b25*(1.d0+dudx) + s22(1,i)*b26*dudy+s33(1,i)*b27*dudz &
1456  + s12(1,i)*( b26*(1.d0+dudx) + b25*dudy ) &
1457  + s23(1,i)*( b27*dudy + b26*dudz ) &
1458  + s13(1,i)*( b27*(1.d0+dudx) + b25*dudz ) )
1459  r26 = &
1460  ( s11(1,i)*b25*dvdx + s22(1,i)*b26*(1.d0+dvdy)+s33(1,i)*b27*dvdz &
1461  + s12(1,i)*( b25*(1.d0+dvdy) + b26*dvdx ) &
1462  + s23(1,i)*( b27*(1.d0+dvdy) + b26*dvdz ) &
1463  + s13(1,i)*( b27*dvdx + b25*dvdz ) )
1464  r27 = &
1465  ( s11(1,i)*b25*dwdx + s22(1,i)*b26*dwdy+s33(1,i)*b27*(1.d0+dwdz) &
1466  + s12(1,i)*( b26*dwdx + b25*dwdy ) &
1467  + s23(1,i)*( b27*dwdy + b26*(1.d0 + dwdz) ) &
1468  + s13(1,i)*( b27*dwdx + b25*(1.d0 + dwdz) ) )
1469 
1470  r28 = &
1471  ( s11(1,i)*b28*(1.d0+dudx) + s22(1,i)*b29*dudy+s33(1,i)*b30*dudz &
1472  + s12(1,i)*( b29*(1.d0+dudx) + b28*dudy ) &
1473  + s23(1,i)*( b30*dudy + b29*dudz ) &
1474  + s13(1,i)*( b30*(1.d0+dudx) + b28*dudz ) )
1475  r29 = &
1476  ( s11(1,i)*b28*dvdx + s22(1,i)*b29*(1.d0+dvdy)+s33(1,i)*b30*dvdz &
1477  + s12(1,i)*( b28*(1.d0+dvdy) + b29*dvdx ) &
1478  + s23(1,i)*( b30*(1.d0+dvdy) + b29*dvdz ) &
1479  + s13(1,i)*( b30*dvdx + b28*dvdz ) )
1480  r30 = &
1481  ( s11(1,i)*b28*dwdx + s22(1,i)*b29*dwdy+s33(1,i)*b30*(1.d0+dwdz) &
1482  + s12(1,i)*( b29*dwdx + b28*dwdy ) &
1483  + s23(1,i)*( b30*dwdy + b29*(1.d0 + dwdz) ) &
1484  + s13(1,i)*( b30*dwdx + b28*(1.d0 + dwdz) ) )
1485 
1486 ! Wi (i.e. weight) for 1 guass point integration is 1
1487 
1488 ! Wi * 1/6 because the volume of a reference tetrahedra in
1489 ! volume coordinates is 1/6
1490 
1491 ! ASSEMBLE THE INTERNAL FORCE VECTOR
1492 !
1493 ! local node 1
1494  r_in(k1n1) = r_in(k1n1) - r1*0.16666666666666667d0
1495  r_in(k2n1) = r_in(k2n1) - r2*0.16666666666666667d0
1496  r_in(k3n1) = r_in(k3n1) - r3*0.16666666666666667d0
1497 ! local node 2
1498  r_in(k1n2) = r_in(k1n2) - r4*0.16666666666666667d0
1499  r_in(k2n2) = r_in(k2n2) - r5*0.16666666666666667d0
1500  r_in(k3n2) = r_in(k3n2) - r6*0.16666666666666667d0
1501 ! local node 3
1502  r_in(k1n3) = r_in(k1n3) - r7*0.16666666666666667d0
1503  r_in(k2n3) = r_in(k2n3) - r8*0.16666666666666667d0
1504  r_in(k3n3) = r_in(k3n3) - r9*0.16666666666666667d0
1505 ! local node 4
1506  r_in(k1n4) = r_in(k1n4) - r10*0.16666666666666667d0
1507  r_in(k2n4) = r_in(k2n4) - r11*0.16666666666666667d0
1508  r_in(k3n4) = r_in(k3n4) - r12*0.16666666666666667d0
1509 ! local node 5
1510  r_in(k1n5) = r_in(k1n5) - r13*0.16666666666666667d0
1511  r_in(k2n5) = r_in(k2n5) - r14*0.16666666666666667d0
1512  r_in(k3n5) = r_in(k3n5) - r15*0.16666666666666667d0
1513 ! local node 6
1514  r_in(k1n6) = r_in(k1n6) - r16*0.16666666666666667d0
1515  r_in(k2n6) = r_in(k2n6) - r17*0.16666666666666667d0
1516  r_in(k3n6) = r_in(k3n6) - r18*0.16666666666666667d0
1517 ! local node 7
1518  r_in(k1n7) = r_in(k1n7) - r19*0.16666666666666667d0
1519  r_in(k2n7) = r_in(k2n7) - r20*0.16666666666666667d0
1520  r_in(k3n7) = r_in(k3n7) - r21*0.16666666666666667d0
1521 ! local node 8
1522  r_in(k1n8) = r_in(k1n8) - r22*0.16666666666666667d0
1523  r_in(k2n8) = r_in(k2n8) - r23*0.16666666666666667d0
1524  r_in(k3n8) = r_in(k3n8) - r24*0.16666666666666667d0
1525 ! local node 9
1526  r_in(k1n9) = r_in(k1n9) - r25*0.16666666666666667d0
1527  r_in(k2n9) = r_in(k2n9) - r26*0.16666666666666667d0
1528  r_in(k3n9) = r_in(k3n9) - r27*0.16666666666666667d0
1529 ! local node 10
1530  r_in(k1n10) = r_in(k1n10) - r28*0.16666666666666667d0
1531  r_in(k2n10) = r_in(k2n10) - r29*0.16666666666666667d0
1532  r_in(k3n10) = r_in(k3n10) - r30*0.16666666666666667d0
1533 
1534  ENDDO
1535  RETURN
1536 END SUBROUTINE v3d10r_nl
1537 
d
const NT & d
Definition: rational_rotation.h:73
v3d10r_nl
subroutine v3d10r_nl(coor, matcstet, lmcstet, R_in, d, ci, cj, S11, S22, S33, S12, S23, S13, numnp, nstart, nend, numcstet, numat_vol)
Definition: v3d10r_nl.f90:53
i
blockLoc i
Definition: read.cpp:79
j
j indices j
Definition: Indexing.h:6