Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d10_nl_matous.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_matous(coor,lmcstet,R_in,d,&
54  s11,s22,s33,s12,s23,s13, &
55  numnp,nstart,nend,numcstet,&
56  p1, p2, yin, softparam, km, k2, gm, g2, vm, nu2, a_eta, a_zeta,&
57  cm, c2, cd, lo, lm, l2, mm, m2, l_bar, strainold)
58 
59 !________________________________________________________________________
60 !
61 ! V3D10 - Performs displacement based computations for Volumetric 3D,
62 ! 10-node Quadratic Linear Elastic Tetrahedron with quadratic
63 ! interpolation functions. (linear strain tetrahedra).
64 ! Large Deformation. Returns the internal force vector R_in.
65 !
66 ! DATE: 10.2001 AUTHOR: SCOT BREITENFELD
67 !________________________________________________________________________
68 
69  IMPLICIT NONE
70 
71 
72 !-----Global variables
73  INTEGER :: numnp ! number of nodes
74  INTEGER :: numcstet ! number of LSTets
75 
76  REAL(KIND=8), DIMENSION(1:6,1:6) :: lo, l_bar
77  REAL(KIND=8) :: km, k2, gm, g2, vm, nu2, a_eta, a_zeta, cm, cb
78 
79  REAL(KIND=8), DIMENSION(1:4,1:numcstet) :: cd, softparam
80 
81  REAL(KIND=8) :: cd_gauss
82 
83 ! added
84 
85  REAL(KIND=8) :: p1, p2, yin, c2, y
86 
87  REAL(KIND=8) :: gy, tot
88 
89  REAL(KIND=8), DIMENSION(1:6,1:6) :: lm, l2, mm, m2
90 
91 
92  real*8, dimension(1:numcstet) :: straintrace
93  real*8 :: straintrace_gauss
94 
95  integer :: ii
96 !-- coordinate array
97  REAL*8, DIMENSION(1:3,1:numnp) :: coor
98 !-- internal force
99  REAL*8, DIMENSION(1:3*numnp) :: r_in
100 !-- displacement vector
101  REAL*8, DIMENSION(1:3*numnp) :: d
102 !-- CSTet stress
103  REAL*8, DIMENSION(1:4,1:numcstet) :: s11, s22, s33, s12, s23, s13
104 !-- connectivity table for CSTet
105  INTEGER, DIMENSION(1:10,1:numcstet) :: lmcstet
106 !-- mat number for CSTet element
107  INTEGER, DIMENSION(1:numcstet) :: matcstet
108 !---- Local variables
109 !-- node numbers
110  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
111 !-- x, y, and z displacements of nodes
112  REAL*8 :: u1,u2,u3,u4,u5,u6,u7,u8,u9,u10
113  REAL*8 :: v1,v2,v3,v4,v5,v6,v7,v8,v9,v10
114  REAL*8 :: w1,w2,w3,w4,w5,w6,w7,w8,w9,w10
115 !-- 6*volume and the volume
116  REAL*8 :: vx6,vx6inv
117 !-- spacial derivatives
118  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10
119  REAL*8 :: b11,b12,b13,b14,b15,b16,b17,b18,b19,b20
120  REAL*8 :: b21,b22,b23,b24,b25,b26,b27,b28,b29,b30
121 !-- partial derivatives of the displacement
122  REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx
123 !-- strains
124  REAL*8 :: e11,e22,e33,e12,e23,e13
125 !-- coordinate holding variable
126  REAL*8 :: x1,x2,x3,x4,x5,x6,x7,x8,x9,x10
127  REAL*8 :: y1,y2,y3,y4,y5,y6,y7,y8,y9,y10
128  REAL*8 :: z1,z2,z3,z4,z5,z6,z7,z8,z9,z10
129 !-- dummy and counters
130  INTEGER :: i,j,nstart,nend
131  REAL*8 :: aux1,aux2,aux3,aux4,aux5,aux6,aux7,aux8,aux9,aux10,aux11,aux12
132 !-- partial internal force
133  REAL*8 :: r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13,r14,r15,r16,r17,r18
134  REAL*8 :: r19,r20,r21,r22,r23,r24,r25,r26,r27,r28,r29,r30
135  REAL*8 :: gaus1, gaus2, gaus3, gaus4
136  REAL*8 :: xn1, xn2, xn3, xn4
137 !-- Coordinate subtractions
138  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
139 !--
140  REAL*8 :: c11, c21, c31
141  INTEGER :: k1n1,k1n2,k1n3,k1n4,k1n5,k1n6,k1n7,k1n8,k1n9,k1n10
142  INTEGER :: k2n1,k2n2,k2n3,k2n4,k2n5,k2n6,k2n7,k2n8,k2n9,k2n10
143  INTEGER :: k3n1,k3n2,k3n3,k3n4,k3n5,k3n6,k3n7,k3n8,k3n9,k3n10
144 !--
145  REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
146  REAL*8, dimension(1:3,1:3) :: f
147  REAL*8, dimension(1:3,1:3) :: s
148 
149  REAL*8 :: strainpt
150 
151  REAL*8,DIMENSION(1:4,1:4) :: gaussintpt = reshape( &
152  (/0.58541020d0,0.13819660d0,0.13819660d0,0.13819660d0, &
153  0.13819660d0,0.58541020d0,0.13819660d0,0.13819660d0, &
154  0.13819660d0,0.13819660d0,0.58541020d0,0.13819660d0, &
155  0.13819660d0,0.13819660d0,0.13819660d0,0.58541020d0/),(/4,4/) )
156 
157  INTEGER :: igpt
158 
159  REAL*8, DIMENSION(1:3,1:3) :: fij, cij
160  INTEGER :: kk,ll
161 
162 
163  REAL*8, DIMENSION(1:6) :: strain, delstrain
164 
165  REAL(KIND=8), DIMENSION(1:4,1:numcstet*6) :: strainold
166 
167 
168  DO i = nstart, nend
169 
170 
171  n1 = lmcstet(1,i)
172  n2 = lmcstet(2,i)
173  n3 = lmcstet(3,i)
174  n4 = lmcstet(4,i)
175  n5 = lmcstet(5,i)
176  n6 = lmcstet(6,i)
177  n7 = lmcstet(7,i)
178  n8 = lmcstet(8,i)
179  n9 = lmcstet(9,i)
180  n10 = lmcstet(10,i)
181 
182  k3n1 = 3*n1
183  k3n2 = 3*n2
184  k3n3 = 3*n3
185  k3n4 = 3*n4
186  k3n5 = 3*n5
187  k3n6 = 3*n6
188  k3n7 = 3*n7
189  k3n8 = 3*n8
190  k3n9 = 3*n9
191  k3n10 = 3*n10
192 
193  k2n1 = k3n1 - 1
194  k2n2 = k3n2 - 1
195  k2n3 = k3n3 - 1
196  k2n4 = k3n4 - 1
197  k2n5 = k3n5 - 1
198  k2n6 = k3n6 - 1
199  k2n7 = k3n7 - 1
200  k2n8 = k3n8 - 1
201  k2n9 = k3n9 - 1
202  k2n10 = k3n10 - 1
203 
204  k1n1 = k3n1 - 2
205  k1n2 = k3n2 - 2
206  k1n3 = k3n3 - 2
207  k1n4 = k3n4 - 2
208  k1n5 = k3n5 - 2
209  k1n6 = k3n6 - 2
210  k1n7 = k3n7 - 2
211  k1n8 = k3n8 - 2
212  k1n9 = k3n9 - 2
213  k1n10 = k3n10 - 2
214 
215  ! k#n# dummy variables replaces:
216  u1 = d(k1n1) ! (3*n1 -2)
217  u2 = d(k1n2) ! (3*n2 -2)
218  u3 = d(k1n3) ! (3*n3 -2)
219  u4 = d(k1n4) ! (3*n4 -2)
220  u5 = d(k1n5) ! (3*n5 -2)
221  u6 = d(k1n6) ! (3*n6 -2)
222  u7 = d(k1n7) ! (3*n7 -2)
223  u8 = d(k1n8) ! (3*n8 -2)
224  u9 = d(k1n9) ! (3*n9 -2)
225  u10 = d(k1n10) ! (3*n10-2)
226  v1 = d(k2n1) ! (3*n1 -1)
227  v2 = d(k2n2) ! (3*n2 -1)
228  v3 = d(k2n3) ! (3*n3 -1)
229  v4 = d(k2n4) ! (3*n4 -1)
230  v5 = d(k2n5) ! (3*n5 -1)
231  v6 = d(k2n6) ! (3*n6 -1)
232  v7 = d(k2n7) ! (3*n7 -1)
233  v8 = d(k2n8) ! (3*n8 -1)
234  v9 = d(k2n9) ! (3*n9 -1)
235  v10 = d(k2n10) ! (3*n10-1)
236  w1 = d(k3n1) ! (3*n1)
237  w2 = d(k3n2) ! (3*n2)
238  w3 = d(k3n3) ! (3*n3)
239  w4 = d(k3n4) ! (3*n4)
240  w5 = d(k3n5) ! (3*n5)
241  w6 = d(k3n6) ! (3*n6)
242  w7 = d(k3n7) ! (3*n7)
243  w8 = d(k3n8) ! (3*n8)
244  w9 = d(k3n9) ! (3*n9)
245  w10 = d(k3n10) ! (3*n10)
246 
247  x1 = coor(1,n1)
248  x2 = coor(1,n2)
249  x3 = coor(1,n3)
250  x4 = coor(1,n4)
251  y1 = coor(2,n1)
252  y2 = coor(2,n2)
253  y3 = coor(2,n3)
254  y4 = coor(2,n4)
255  z1 = coor(3,n1)
256  z2 = coor(3,n2)
257  z3 = coor(3,n3)
258  z4 = coor(3,n4)
259 
260 
261  vx6 = x2*y3*z4 - x2*y4*z3 - y2*x3*z4 + y2*x4*z3 + z2*x3*y4 &
262  - z2*x4*y3 - x1*y3*z4 + x1*y4*z3 + x1*y2*z4 - x1*y2*z3 &
263  - x1*z2*y4 + x1*z2*y3 + y1*x3*z4 - y1*x4*z3 - y1*x2*z4 &
264  + y1*x2*z3 + y1*z2*x4 - y1*z2*x3 - z1*x3*y4 + z1*x4*y3 &
265  + z1*x2*y4 - z1*x2*y3 - z1*y2*x4 + z1*y2*x3
266 
267 
268  vx6inv = 1.d0 / vx6
269 
270  aux1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3)
271  aux2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3)
272  aux3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3)
273  aux4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3)
274  aux5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3)
275  aux6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3)
276  aux7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2)
277  aux8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2)
278  aux9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2)
279  aux10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2)
280  aux11 =-(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2)
281  aux12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2)
282 
283  r1 = 0.d0
284  r2 = 0.d0
285  r3 = 0.d0
286  r4 = 0.d0
287  r5 = 0.d0
288  r6 = 0.d0
289  r7 = 0.d0
290  r8 = 0.d0
291  r9 = 0.d0
292  r10 = 0.d0
293  r11 = 0.d0
294  r12 = 0.d0
295  r13 = 0.d0
296  r14 = 0.d0
297  r15 = 0.d0
298  r16 = 0.d0
299  r17 = 0.d0
300  r18 = 0.d0
301  r19 = 0.d0
302  r20 = 0.d0
303  r21 = 0.d0
304  r22 = 0.d0
305  r23 = 0.d0
306  r24 = 0.d0
307  r25 = 0.d0
308  r26 = 0.d0
309  r27 = 0.d0
310  r28 = 0.d0
311  r29 = 0.d0
312  r30 = 0.d0
313 
314  straintrace_gauss = 0.d0
315 
316 
317  DO igpt = 1, 4
318 ! PRINT*,i,nstart,nend
319  gaus1 = gaussintpt(igpt,1)
320  gaus2 = gaussintpt(igpt,2)
321  gaus3 = gaussintpt(igpt,3)
322  gaus4 = gaussintpt(igpt,4)
323 ! PRINT*,i,nstart,nend
324  xn1 = (4.d0*gaus1-1.d0) ! derivative of shape function
325  xn2 = (4.d0*gaus2-1.d0) ! dN_i/dzeta_i
326  xn3 = (4.d0*gaus3-1.d0)
327  xn4 = (4.d0*gaus4-1.d0)
328 ! xN5 = 4.d0*gaus1*gaus2
329 ! xN6 = 4.d0*gaus2*gaus3
330 ! xN7 = 4.d0*gaus3*gaus1
331 ! xN8 = 4.d0*gaus1*gaus4
332 ! xN9 = 4.d0*gaus2*gaus4
333 ! xN10= 4.d0*gaus3*gaus4
334 
335  b1 = aux1*xn1
336  b2 = aux2*xn1
337  b3 = aux3*xn1
338  b4 = aux4*xn2
339  b5 = aux5*xn2
340  b6 = aux6*xn2
341  b7 = aux7*xn3
342  b8 = aux8*xn3
343  b9 = aux9*xn3
344  b10 = aux10*xn4
345  b11 = aux11*xn4
346  b12 = aux12*xn4
347 
348  b13 = 4.d0*(gaus2*aux1 + gaus1*aux4)
349  b14 = 4.d0*(gaus2*aux2 + gaus1*aux5)
350  b15 = 4.d0*(gaus2*aux3 + gaus1*aux6)
351 
352  b16 = 4.d0*(gaus3*aux4 + gaus2*aux7)
353  b17 = 4.d0*(gaus3*aux5 + gaus2*aux8)
354  b18 = 4.d0*(gaus3*aux6 + gaus2*aux9)
355 
356  b19 = 4.d0*(gaus1*aux7 + gaus3*aux1)
357  b20 = 4.d0*(gaus1*aux8 + gaus3*aux2)
358  b21 = 4.d0*(gaus1*aux9 + gaus3*aux3)
359 
360  b22 = 4.d0*(gaus4*aux1 + gaus1*aux10)
361  b23 = 4.d0*(gaus4*aux2 + gaus1*aux11)
362  b24 = 4.d0*(gaus4*aux3 + gaus1*aux12)
363 
364  b25 = 4.d0*(gaus4*aux4 + gaus2*aux10)
365  b26 = 4.d0*(gaus4*aux5 + gaus2*aux11)
366  b27 = 4.d0*(gaus4*aux6 + gaus2*aux12)
367 
368  b28 = 4.d0*(gaus4*aux7 + gaus3*aux10)
369  b29 = 4.d0*(gaus4*aux8 + gaus3*aux11)
370  b30 = 4.d0*(gaus4*aux9 + gaus3*aux12)
371 !
372 !-----Calculate displacement gradient (H)
373 !
374  dudx = (b1*u1 + b4*u2 + b7*u3 + b10*u4 + b13*u5 + b16*u6 + b19*u7 + b22*u8 + b25*u9 + b28*u10)*vx6inv
375  dvdy = (b2*v1 + b5*v2 + b8*v3 + b11*v4 + b14*v5 + b17*v6 + b20*v7 + b23*v8 + b26*v9 + b29*v10)*vx6inv
376  dwdz = (b3*w1 + b6*w2 + b9*w3 + b12*w4 + b15*w5 + b18*w6 + b21*w7 + b24*w8 + b27*w9 + b30*w10)*vx6inv
377  dudy = (b2*u1 + b5*u2 + b8*u3 + b11*u4 + b14*u5 + b17*u6 + b20*u7 + b23*u8 + b26*u9 + b29*u10)*vx6inv
378  dvdx = (b1*v1 + b4*v2 + b7*v3 + b10*v4 + b13*v5 + b16*v6 + b19*v7 + b22*v8 + b25*v9 + b28*v10)*vx6inv
379  dvdz = (b3*v1 + b6*v2 + b9*v3 + b12*v4 + b15*v5 + b18*v6 + b21*v7 + b24*v8 + b27*v9 + b30*v10)*vx6inv
380  dwdy = (b2*w1 + b5*w2 + b8*w3 + b11*w4 + b14*w5 + b17*w6 + b20*w7 + b23*w8 + b26*w9 + b29*w10)*vx6inv
381  dudz = (b3*u1 + b6*u2 + b9*u3 + b12*u4 + b15*u5 + b18*u6 + b21*u7 + b24*u8 + b27*u9 + b30*u10)*vx6inv
382  dwdx = (b1*w1 + b4*w2 + b7*w3 + b10*w4 + b13*w5 + b16*w6 + b19*w7 + b22*w8 + b25*w9 + b28*w10)*vx6inv
383 
384 
385 ! PRINT*,i,nstart,nend
386 ! calculate the Green-Lagrange Strain (E) strain
387 ! [E] = [B]{d} = E11, E22, E33, E23, E13, E12
388 !
389 
390  strain(1) = dudx + 0.5d0*(dudx*dudx + dvdx*dvdx + dwdx*dwdx) ! E11
391  strain(2) = dvdy + 0.5d0*(dudy*dudy + dvdy*dvdy + dwdy*dwdy) ! E22
392  strain(3) = dwdz + 0.5d0*(dudz*dudz + dvdz*dvdz + dwdz*dwdz) ! E33
393  strain(6) = dudy + dvdx + ( dudx*dudy + dvdx*dvdy + dwdx*dwdy ) ! E12
394  strain(4) = dvdz + dwdy + ( dudy*dudz + dvdy*dvdz + dwdy*dwdz ) ! E23
395  strain(5) = dwdx + dudz + ( dudz*dudx + dvdz*dvdx + dwdz*dwdx ) ! E13
396 
397 ! calculate the strain increment
398 
399 ! PRINT*,i,nstart,nend
400  delstrain(1) = strain(1) - strainold(igpt,(i-1)*6 + 1)
401  delstrain(2) = strain(2) - strainold(igpt,(i-1)*6 + 2)
402  delstrain(3) = strain(3) - strainold(igpt,(i-1)*6 + 3)
403  delstrain(4) = strain(4) - strainold(igpt,(i-1)*6 + 4)
404  delstrain(5) = strain(5) - strainold(igpt,(i-1)*6 + 5)
405  delstrain(6) = strain(6) - strainold(igpt,(i-1)*6 + 6)
406 
407 
408 ! PRINT*,i,igpt,nstart,nend
409  CALL matous_const_model(cd(igpt,i),cb, &
410  p1, p2, yin, softparam(igpt,i), km, k2, gm, g2, vm, nu2, a_eta, a_zeta,&
411  cm, c2, lo, lm, l2, mm, m2, l_bar, strain, delstrain, &
412  s11(igpt,i), s22(igpt,i), s33(igpt,i), s12(igpt,i), s13(igpt,i), s23(igpt,i),i)
413 
414 ! PRINT*,i,igpt,nstart,nend
415  straintrace_gauss = straintrace_gauss + strain(1)+ strain(2)+ strain(3)
416 
417 
418  r1 = r1 + &
419  ( s11(igpt,i)*b1*(1.d0+dudx) + s22(igpt,i)*b2*dudy + s33(igpt,i)*b3*dudz &
420  + s12(igpt,i)*( b2*(1.d0+dudx) + b1*dudy ) &
421  + s23(igpt,i)*( b3*dudy + b2*dudz ) &
422  + s13(igpt,i)*( b3*(1.d0+dudx) + b1*dudz ) )
423  r2 = r2 + &
424  ( s11(igpt,i)*b1*dvdx + s22(igpt,i)*b2*(1.d0+dvdy) + s33(igpt,i)*b3*dvdz &
425  + s12(igpt,i)*( b1*(1.d0+dvdy) + b2*dvdx ) &
426  + s23(igpt,i)*( b3*(1.d0+dvdy) + b2*dvdz ) &
427  + s13(igpt,i)*( b3*dvdx + b1*dvdz ) )
428  r3 = r3 + &
429  ( s11(igpt,i)*b1*dwdx + s22(igpt,i)*b2*dwdy + s33(igpt,i)*b3*(1.d0+dwdz) &
430  + s12(igpt,i)*( b2*dwdx + b1*dwdy ) &
431  + s23(igpt,i)*( b3*dwdy + b2*(1.d0 + dwdz) ) &
432  + s13(igpt,i)*( b3*dwdx + b1*(1.d0 + dwdz) ) )
433 
434  r4 = r4 + &
435  ( s11(igpt,i)*b4*(1.d0+dudx) + s22(igpt,i)*b5*dudy + s33(igpt,i)*b6*dudz &
436  + s12(igpt,i)*( b5*(1.d0+dudx) + b4*dudy ) &
437  + s23(igpt,i)*( b6*dudy + b5*dudz ) &
438  + s13(igpt,i)*( b6*(1.d0+dudx) + b4*dudz ) )
439  r5 = r5 + &
440  ( s11(igpt,i)*b4*dvdx + s22(igpt,i)*b5*(1.d0+dvdy) + s33(igpt,i)*b6*dvdz &
441  + s12(igpt,i)*( b4*(1.d0+dvdy) + b5*dvdx ) &
442  + s23(igpt,i)*( b6*(1.d0+dvdy) + b5*dvdz ) &
443  + s13(igpt,i)*( b6*dvdx + b4*dvdz ) )
444  r6 = r6 + &
445  ( s11(igpt,i)*b4*dwdx + s22(igpt,i)*b5*dwdy + s33(igpt,i)*b6*(1.d0+dwdz) &
446  + s12(igpt,i)*( b5*dwdx + b4*dwdy ) &
447  + s23(igpt,i)*( b6*dwdy + b5*(1.d0 + dwdz) ) &
448  + s13(igpt,i)*( b6*dwdx + b4*(1.d0 + dwdz) ) )
449 
450  r7 = r7 + &
451  ( s11(igpt,i)*b7*(1.d0+dudx) + s22(igpt,i)*b8*dudy + s33(igpt,i)*b9*dudz &
452  + s12(igpt,i)*( b8*(1.d0+dudx) + b7*dudy ) &
453  + s23(igpt,i)*( b9*dudy + b8*dudz ) &
454  + s13(igpt,i)*( b9*(1.d0+dudx) + b7*dudz ) )
455  r8 = r8 + &
456  ( s11(igpt,i)*b7*dvdx + s22(igpt,i)*b8*(1.d0+dvdy) + s33(igpt,i)*b9*dvdz &
457  + s12(igpt,i)*( b7*(1.d0+dvdy) + b8*dvdx ) &
458  + s23(igpt,i)*( b9*(1.d0+dvdy) + b8*dvdz ) &
459  + s13(igpt,i)*( b9*dvdx + b7*dvdz ) )
460  r9 = r9 + &
461  ( s11(igpt,i)*b7*dwdx + s22(igpt,i)*b8*dwdy + s33(igpt,i)*b9*(1.d0+dwdz) &
462  + s12(igpt,i)*( b8*dwdx + b7*dwdy ) &
463  + s23(igpt,i)*( b9*dwdy + b8*(1.d0 + dwdz) ) &
464  + s13(igpt,i)*( b9*dwdx + b7*(1.d0 + dwdz) ) )
465 
466  r10 = r10 + &
467  ( s11(igpt,i)*b10*(1.d0+dudx) + s22(igpt,i)*b11*dudy+s33(igpt,i)*b12*dudz &
468  + s12(igpt,i)*( b11*(1.d0+dudx) + b10*dudy ) &
469  + s23(igpt,i)*( b12*dudy + b11*dudz ) &
470  + s13(igpt,i)*( b12*(1.d0+dudx) + b10*dudz ) )
471  r11 = r11 + &
472  ( s11(igpt,i)*b10*dvdx + s22(igpt,i)*b11*(1.d0+dvdy)+s33(igpt,i)*b12*dvdz &
473  + s12(igpt,i)*( b10*(1.d0+dvdy) + b11*dvdx ) &
474  + s23(igpt,i)*( b12*(1.d0+dvdy) + b11*dvdz ) &
475  + s13(igpt,i)*( b12*dvdx + b10*dvdz ) )
476  r12 = r12 + &
477  ( s11(igpt,i)*b10*dwdx + s22(igpt,i)*b11*dwdy+s33(igpt,i)*b12*(1.d0+dwdz) &
478  + s12(igpt,i)*( b11*dwdx + b10*dwdy ) &
479  + s23(igpt,i)*( b12*dwdy + b11*(1.d0 + dwdz) ) &
480  + s13(igpt,i)*( b12*dwdx + b10*(1.d0 + dwdz) ) )
481 
482  r13 = r13 + &
483  ( s11(igpt,i)*b13*(1.d0+dudx) + s22(igpt,i)*b14*dudy + s33(igpt,i)*b15*dudz &
484  + s12(igpt,i)*( b14*(1.d0+dudx) + b13*dudy ) &
485  + s23(igpt,i)*( b15*dudy + b14*dudz ) &
486  + s13(igpt,i)*( b15*(1.d0+dudx) + b13*dudz ) )
487  r14 = r14 + &
488  ( s11(igpt,i)*b13*dvdx + s22(igpt,i)*b14*(1.d0+dvdy) + s33(igpt,i)*b15*dvdz &
489  + s12(igpt,i)*( b13*(1.d0+dvdy) + b14*dvdx ) &
490  + s23(igpt,i)*( b15*(1.d0+dvdy) + b14*dvdz ) &
491  + s13(igpt,i)*( b15*dvdx + b13*dvdz ) )
492  r15 = r15 + &
493  ( s11(igpt,i)*b13*dwdx + s22(igpt,i)*b14*dwdy + s33(igpt,i)*b15*(1.d0+dwdz) &
494  + s12(igpt,i)*( b14*dwdx + b13*dwdy ) &
495  + s23(igpt,i)*( b15*dwdy + b14*(1.d0 + dwdz) ) &
496  + s13(igpt,i)*( b15*dwdx + b13*(1.d0 + dwdz) ) )
497 
498  r16 = r16 + &
499  ( s11(igpt,i)*b16*(1.d0+dudx) + s22(igpt,i)*b17*dudy + s33(igpt,i)*b18*dudz &
500  + s12(igpt,i)*( b17*(1.d0+dudx) + b16*dudy ) &
501  + s23(igpt,i)*( b18*dudy + b17*dudz ) &
502  + s13(igpt,i)*( b18*(1.d0+dudx) + b16*dudz ) )
503  r17 = r17 + &
504  ( s11(igpt,i)*b16*dvdx + s22(igpt,i)*b17*(1.d0+dvdy) + s33(igpt,i)*b18*dvdz &
505  + s12(igpt,i)*( b16*(1.d0+dvdy) + b17*dvdx ) &
506  + s23(igpt,i)*( b18*(1.d0+dvdy) + b17*dvdz ) &
507  + s13(igpt,i)*( b18*dvdx + b16*dvdz ) )
508  r18 = r18 + &
509  ( s11(igpt,i)*b16*dwdx + s22(igpt,i)*b17*dwdy + s33(igpt,i)*b18*(1.d0+dwdz) &
510  + s12(igpt,i)*( b17*dwdx + b16*dwdy ) &
511  + s23(igpt,i)*( b18*dwdy + b17*(1.d0 + dwdz) ) &
512  + s13(igpt,i)*( b18*dwdx + b16*(1.d0 + dwdz) ) )
513 
514  r19 = r19 + &
515  ( s11(igpt,i)*b19*(1.d0+dudx) + s22(igpt,i)*b20*dudy + s33(igpt,i)*b21*dudz &
516  + s12(igpt,i)*( b20*(1.d0+dudx) + b19*dudy ) &
517  + s23(igpt,i)*( b21*dudy + b20*dudz ) &
518  + s13(igpt,i)*( b21*(1.d0+dudx) + b19*dudz ) )
519  r20 = r20 + &
520  ( s11(igpt,i)*b19*dvdx + s22(igpt,i)*b20*(1.d0+dvdy) + s33(igpt,i)*b21*dvdz &
521  + s12(igpt,i)*( b19*(1.d0+dvdy) + b20*dvdx ) &
522  + s23(igpt,i)*( b21*(1.d0+dvdy) + b20*dvdz ) &
523  + s13(igpt,i)*( b21*dvdx + b19*dvdz ) )
524  r21 = r21 + &
525  ( s11(igpt,i)*b19*dwdx + s22(igpt,i)*b20*dwdy + s33(igpt,i)*b21*(1.d0+dwdz) &
526  + s12(igpt,i)*( b20*dwdx + b19*dwdy ) &
527  + s23(igpt,i)*( b21*dwdy + b20*(1.d0 + dwdz) ) &
528  + s13(igpt,i)*( b21*dwdx + b19*(1.d0 + dwdz) ) )
529 
530  r22 = r22 + &
531  ( s11(igpt,i)*b22*(1.d0+dudx) + s22(igpt,i)*b23*dudy+s33(igpt,i)*b24*dudz &
532  + s12(igpt,i)*( b23*(1.d0+dudx) + b22*dudy ) &
533  + s23(igpt,i)*( b24*dudy + b23*dudz ) &
534  + s13(igpt,i)*( b24*(1.d0+dudx) + b22*dudz ) )
535  r23 = r23 + &
536  ( s11(igpt,i)*b22*dvdx + s22(igpt,i)*b23*(1.d0+dvdy)+s33(igpt,i)*b24*dvdz &
537  + s12(igpt,i)*( b22*(1.d0+dvdy) + b23*dvdx ) &
538  + s23(igpt,i)*( b24*(1.d0+dvdy) + b23*dvdz ) &
539  + s13(igpt,i)*( b24*dvdx + b22*dvdz ) )
540  r24 = r24 + &
541  ( s11(igpt,i)*b22*dwdx + s22(igpt,i)*b23*dwdy+s33(igpt,i)*b24*(1.d0+dwdz) &
542  + s12(igpt,i)*( b23*dwdx + b22*dwdy ) &
543  + s23(igpt,i)*( b24*dwdy + b23*(1.d0 + dwdz) ) &
544  + s13(igpt,i)*( b24*dwdx + b22*(1.d0 + dwdz) ) )
545 
546  r25 = r25 + &
547  ( s11(igpt,i)*b25*(1.d0+dudx) + s22(igpt,i)*b26*dudy+s33(igpt,i)*b27*dudz &
548  + s12(igpt,i)*( b26*(1.d0+dudx) + b25*dudy ) &
549  + s23(igpt,i)*( b27*dudy + b26*dudz ) &
550  + s13(igpt,i)*( b27*(1.d0+dudx) + b25*dudz ) )
551  r26 = r26 + &
552  ( s11(igpt,i)*b25*dvdx + s22(igpt,i)*b26*(1.d0+dvdy)+s33(igpt,i)*b27*dvdz &
553  + s12(igpt,i)*( b25*(1.d0+dvdy) + b26*dvdx ) &
554  + s23(igpt,i)*( b27*(1.d0+dvdy) + b26*dvdz ) &
555  + s13(igpt,i)*( b27*dvdx + b25*dvdz ) )
556  r27 = r27 + &
557  ( s11(igpt,i)*b25*dwdx + s22(igpt,i)*b26*dwdy+s33(igpt,i)*b27*(1.d0+dwdz) &
558  + s12(igpt,i)*( b26*dwdx + b25*dwdy ) &
559  + s23(igpt,i)*( b27*dwdy + b26*(1.d0 + dwdz) ) &
560  + s13(igpt,i)*( b27*dwdx + b25*(1.d0 + dwdz) ) )
561 
562  r28 = r28 + &
563  ( s11(igpt,i)*b28*(1.d0+dudx) + s22(igpt,i)*b29*dudy+s33(igpt,i)*b30*dudz &
564  + s12(igpt,i)*( b29*(1.d0+dudx) + b28*dudy ) &
565  + s23(igpt,i)*( b30*dudy + b29*dudz ) &
566  + s13(igpt,i)*( b30*(1.d0+dudx) + b28*dudz ) )
567  r29 = r29 + &
568  ( s11(igpt,i)*b28*dvdx + s22(igpt,i)*b29*(1.d0+dvdy)+s33(igpt,i)*b30*dvdz &
569  + s12(igpt,i)*( b28*(1.d0+dvdy) + b29*dvdx ) &
570  + s23(igpt,i)*( b30*(1.d0+dvdy) + b29*dvdz ) &
571  + s13(igpt,i)*( b30*dvdx + b28*dvdz ) )
572  r30 = r30 + &
573  ( s11(igpt,i)*b28*dwdx + s22(igpt,i)*b29*dwdy+s33(igpt,i)*b30*(1.d0+dwdz) &
574  + s12(igpt,i)*( b29*dwdx + b28*dwdy ) &
575  + s23(igpt,i)*( b30*dwdy + b29*(1.d0 + dwdz) ) &
576  + s13(igpt,i)*( b30*dwdx + b28*(1.d0 + dwdz) ) )
577 ! PRINT*,i,igpt,nstart,nend
578 
579  ENDDO
580 
581 ! PRINT*,i,nstart,nend
582 
583  straintrace(i) = straintrace_gauss*0.25d0
584 
585 
586 ! Wi (i.e. weight) for 4 guass point integration is 1/4
587 
588 ! Wi * 1/6 because the volume of a reference tetrahedra in
589 ! volume coordinates is 1/6
590 
591 ! ASSEMBLE THE INTERNAL FORCE VECTOR
592 !
593 ! local node 1
594  r_in(k1n1) = r_in(k1n1) - r1*0.04166666666666667d0
595  r_in(k2n1) = r_in(k2n1) - r2*0.04166666666666667d0
596  r_in(k3n1) = r_in(k3n1) - r3*0.04166666666666667d0
597 ! local node 2
598  r_in(k1n2) = r_in(k1n2) - r4*0.04166666666666667d0
599  r_in(k2n2) = r_in(k2n2) - r5*0.04166666666666667d0
600  r_in(k3n2) = r_in(k3n2) - r6*0.04166666666666667d0
601 ! local node 3
602  r_in(k1n3) = r_in(k1n3) - r7*0.04166666666666667d0
603  r_in(k2n3) = r_in(k2n3) - r8*0.04166666666666667d0
604  r_in(k3n3) = r_in(k3n3) - r9*0.04166666666666667d0
605 ! local node 4
606  r_in(k1n4) = r_in(k1n4) - r10*0.04166666666666667d0
607  r_in(k2n4) = r_in(k2n4) - r11*0.04166666666666667d0
608  r_in(k3n4) = r_in(k3n4) - r12*0.04166666666666667d0
609 ! local node 5
610  r_in(k1n5) = r_in(k1n5) - r13*0.04166666666666667d0
611  r_in(k2n5) = r_in(k2n5) - r14*0.04166666666666667d0
612  r_in(k3n5) = r_in(k3n5) - r15*0.04166666666666667d0
613 ! local node 6
614  r_in(k1n6) = r_in(k1n6) - r16*0.04166666666666667d0
615  r_in(k2n6) = r_in(k2n6) - r17*0.04166666666666667d0
616  r_in(k3n6) = r_in(k3n6) - r18*0.04166666666666667d0
617 ! local node 7
618  r_in(k1n7) = r_in(k1n7) - r19*0.04166666666666667d0
619  r_in(k2n7) = r_in(k2n7) - r20*0.04166666666666667d0
620  r_in(k3n7) = r_in(k3n7) - r21*0.04166666666666667d0
621 ! local node 8
622  r_in(k1n8) = r_in(k1n8) - r22*0.04166666666666667d0
623  r_in(k2n8) = r_in(k2n8) - r23*0.04166666666666667d0
624  r_in(k3n8) = r_in(k3n8) - r24*0.04166666666666667d0
625 ! local node 9
626  r_in(k1n9) = r_in(k1n9) - r25*0.04166666666666667d0
627  r_in(k2n9) = r_in(k2n9) - r26*0.04166666666666667d0
628  r_in(k3n9) = r_in(k3n9) - r27*0.04166666666666667d0
629 ! local node 10
630  r_in(k1n10) = r_in(k1n10) - r28*0.04166666666666667d0
631  r_in(k2n10) = r_in(k2n10) - r29*0.04166666666666667d0
632  r_in(k3n10) = r_in(k3n10) - r30*0.04166666666666667d0
633 !
634 ! PRINT*,i,nstart,nend
635  ENDDO
636 
637 ! PRINT*,'(((', nstart, nend
638 
639 
640  RETURN
641 END SUBROUTINE v3d10_nl_matous
642 
d
const NT & d
Definition: rational_rotation.h:73
y
void int int REAL REAL * y
Definition: read.cpp:74
s
double s
Definition: blastest.C:80
p1
NT p1
Definition: rational_rotation.h:108
matous_const_model
subroutine matous_const_model(cd, cb, p1, p2, Yin, SoftParam, Km, K2, Gm, G2, vm, nu2, a_eta, a_zeta, cm, c2, Lo, Lm, L2, Mm, M2, L_bar, strain, DelStrain, s11, s22, s33, s12, s13, s23, ielem)
Definition: Matous_Const_Model.f90:53
i
blockLoc i
Definition: read.cpp:79
v3d10_nl_matous
subroutine v3d10_nl_matous(coor, lmcstet, R_in, d, S11, S22, S33, S12, S23, S13, numnp, nstart, nend, numcstet, p1, p2, Yin, SoftParam, Km, K2, Gm, G2, vm, nu2, a_eta, a_zeta, cm, c2, cd, Lo, Lm, L2, Mm, M2, L_bar, StrainOld)
Definition: v3d10_nl_matous.f90:53
j
j indices j
Definition: Indexing.h:6