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