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