Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MP/Source/v3d4n_NeoHookeanCompress.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 v3d4n_neohookeancompress(numnp, numel, coor, disp, nodes, Rnet,&
54  numelneigh, elconn, alpha, ahat,nummatvol, &
55  xmu,xlambda)
56 
57 ! nprocs,TotNumNdComm,TotNumNeighProcs,NeighProcList,NumNdComm,neigh_lst,&
58 ! MPI_STATUS_SIZE,MPI_COMM_ROCFRAC,MPI_DOUBLE_PRECISION, &
59 ! ReqRcv,ReqSnd,StatRcv, StatSnd)
60 
61 
62 !________________________________________________________________________
63 !
64 ! V3D4 - Performs displacement based computations for Volumetric 3D
65 ! 4-node tetrahedra large deformation (i.e. nonlinear) elastic
66 ! elements with linear interpolation functions. (constant strain
67 ! tetrahedra). Returns the internal force vector R_in.
68 !
69 ! DATE: 04.2000 AUTHOR: SCOT BREITENFELD
70 !
71 ! Source:
72 ! "Nonlinear continuum mechanics for finite element analysis"
73 ! Javier Bonet and Richard D. Wood
74 ! ISBN 0-521-57272-X
75 
76 !________________________________________________________________________
77 
78  USE rocstar_rocfraccomm
79 
80  IMPLICIT NONE
81 
82  integer :: nummatvol
83  INTEGER :: numnp, numel
84  REAL*8,dimension(1:NumMatVol) :: xmu, xlambda
85  REAL*8, DIMENSION(1:numnp) :: ahat
86  REAL*8, DIMENSION(1:3*numnp) :: disp, rnet
87  REAL*8, DIMENSION(1:3,1:numnp) :: coor
88  INTEGER, DIMENSION(1:numnp) ::numelneigh
89  INTEGER, DIMENSION(1:numnp,1:40) :: elconn ! fix 40 should not be hard coded
90  INTEGER,DIMENSION(1:4,1:numel) :: nodes
91  REAL*8, DIMENSION(1:4,1:numel) :: alpha
92  REAL*8, DIMENSION(1:numnp) :: s11n, s22n, s33n, s12n, s13n,s23n
93 
94  INTEGER :: i,k
95 !-- x, y and z displacements of nodes
96  REAL*8 :: u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4
97 !-- coordinate holding variable
98  REAL*8 :: x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
99 !-- Coordinate subtractions
100 
101  REAL*8 :: z12, z13
102  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
103  REAL*8 :: x12, x13, y12, y13
104 
105 !-- 6*volume, inverse(6*volume), and the volume
106  REAL*8 :: vx6, vx6inv, vol
107 !-- spacial derivatives
108  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
109 !-- partial derivatives of the displacement
110  REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx
111 !-- strains
112  REAL*8 :: volnd,vainv
113 ! -- node numbers
114  INTEGER :: n1,n2,n3,n4
115 ! -- dummy and counters
116  INTEGER :: j,nstart,nend
117  INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
118  INTEGER :: k3n1,k3n2,k3n3,k3n4
119 ! --
120  REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
121  REAL*8 :: j1,j2,j2inv
122  REAL*8 :: c11, c12, c13, c21, c22, c23, c31, c32, c33
123  INTEGER :: jj,ielnum
124  REAL*8 :: volnd0inv,volel0
125  REAL*8, DIMENSION(1:1) :: s11, s22, s33, s12, s23, s13
126  INTEGER :: k3i,k2i,k1i
127  INTEGER :: ix
128 
129 ! /* For each Node */
130 
131  DO i = 1, numnp ! for each node
132 
133 ! /* Undeformed Volume of node */
134 
135  volnd0inv = 1.d0/ahat(i)
136 
137 ! /* Initialize to 0, Va, Va0, Fa */
138 
139  f11 = 0.d0
140  f22 = 0.d0
141  f33 = 0.d0
142  f12 = 0.d0
143  f13 = 0.d0
144  f21 = 0.d0
145  f23 = 0.d0
146  f31 = 0.d0
147  f32 = 0.d0
148  volnd = 0.d0
149 
150  DO j = 1, numelneigh(i)
151 
152  ielnum = elconn(i,j)
153 
154  n1 = nodes(1,ielnum)
155  n2 = nodes(2,ielnum)
156  n3 = nodes(3,ielnum)
157  n4 = nodes(4,ielnum)
158 
159  IF(i.EQ.n1) THEN
160  ix = 1
161  ELSE IF(i.EQ.n2)THEN
162  ix = 2
163  ELSE IF(i.EQ.n3)THEN
164  ix = 3
165  ELSE IF(i.EQ.n4)THEN
166  ix = 4
167  ENDIF
168 
169  k3n1 = 3*n1
170  k3n2 = 3*n2
171  k3n3 = 3*n3
172  k3n4 = 3*n4
173 
174  k2n1 = k3n1 - 1
175  k2n2 = k3n2 - 1
176  k2n3 = k3n3 - 1
177  k2n4 = k3n4 - 1
178 
179  k1n1 = k3n1 - 2
180  k1n2 = k3n2 - 2
181  k1n3 = k3n3 - 2
182  k1n4 = k3n4 - 2
183  ! k#n# dummy variables replaces:
184  u1 = disp(k1n1) ! (3*n1-2)
185  u2 = disp(k1n2) ! (3*n2-2)
186  u3 = disp(k1n3) ! (3*n3-2)
187  u4 = disp(k1n4) ! (3*n4-2)
188  v1 = disp(k2n1) ! (3*n1-1)
189  v2 = disp(k2n2) ! (3*n2-1)
190  v3 = disp(k2n3) ! (3*n3-1)
191  v4 = disp(k2n4) ! (3*n4-1)
192  w1 = disp(k3n1) ! (3*n1)
193  w2 = disp(k3n2) ! (3*n2)
194  w3 = disp(k3n3) ! (3*n3)
195  w4 = disp(k3n4) ! (3*n4)
196 
197  x1 = coor(1,n1)
198  x2 = coor(1,n2)
199  x3 = coor(1,n3)
200  x4 = coor(1,n4)
201  y1 = coor(2,n1)
202  y2 = coor(2,n2)
203  y3 = coor(2,n3)
204  y4 = coor(2,n4)
205  z1 = coor(3,n1)
206  z2 = coor(3,n2)
207  z3 = coor(3,n3)
208  z4 = coor(3,n4)
209 
210  x12 = x1 - x2 ! not used in vol. calc
211  x13 = x1 - x3 ! not used in vol. calc
212  x14 = x1 - x4
213  x24 = x2 - x4
214  x34 = x3 - x4
215  y12 = y1 - y2 ! not used in vol. calc
216  y13 = y1 - y3 ! not used in vol. calc
217  y14 = y1 - y4
218  y24 = y2 - y4
219  y34 = y3 - y4
220  z12 = z1 - z2 ! not used in vol. calc
221  z13 = z1 - z3 ! not used in vol. calc
222  z14 = z1 - z4
223  z24 = z2 - z4
224  z34 = z3 - z4
225 
226  c11 = y24*z34 - z24*y34
227  c21 = -( x24*z34 - z24*x34 )
228  c31 = x24*y34 - y24*x34
229 
230  vx6 = -( x14*c11 + y14*c21 + z14*c31 )
231 
232  vx6inv = 1.d0/vx6
233 
234  volel0 = vx6/6.d0 ! undeformed volume of element (Ve_O)
235 
236 
237 ! See the maple worksheet 'V3D4.mws' for the derivation of Vx6
238 
239 
240 ! See the maple worksheet 'V3D4.mws' for the derivation of [B]
241 ! Compute the Shape functions
242 ! NOTE: Factored for a more equivalent/compact form then maple's
243 
244  b1 = (y34*z24 - y24*z34) * vx6inv
245  b2 = (z34*x24 - z24*x34) * vx6inv
246  b3 = (x34*y24 - x24*y34) * vx6inv
247  b4 = (y13*z14 - y14*z13) * vx6inv
248  b5 = (z13*x14 - z14*x13) * vx6inv
249  b6 = (x13*y14 - x14*y13) * vx6inv
250  b7 = (y14*z12 - y12*z14) * vx6inv
251  b8 = (z14*x12 - z12*x14) * vx6inv
252  b9 = (x14*y12 - x12*y14) * vx6inv
253  b10 = (y12*z13 - y13*z12) * vx6inv
254  b11 = (z12*x13 - z13*x12) * vx6inv
255  b12 = (x12*y13 - x13*y12) * vx6inv
256 !
257 ! Calculate displacement gradient (H)
258 !
259  dudx = b1*u1 + b4*u2 + b7*u3 + b10*u4
260  dvdy = b2*v1 + b5*v2 + b8*v3 + b11*v4
261  dwdz = b3*w1 + b6*w2 + b9*w3 + b12*w4
262  dudy = b2*u1 + b5*u2 + b8*u3 + b11*u4
263  dvdx = b1*v1 + b4*v2 + b7*v3 + b10*v4
264  dvdz = b3*v1 + b6*v2 + b9*v3 + b12*v4
265  dwdy = b2*w1 + b5*w2 + b8*w3 + b11*w4
266  dudz = b3*u1 + b6*u2 + b9*u3 + b12*u4
267  dwdx = b1*w1 + b4*w2 + b7*w3 + b10*w4
268 !
269 ! deformation gradients F
270 !
271  f11 = f11 + alpha(ix,ielnum)*volel0*(1.d0 + dudx)
272  f12 = f12 + alpha(ix,ielnum)*volel0*dudy
273  f13 = f13 + alpha(ix,ielnum)*volel0*dudz
274  f21 = f21 + alpha(ix,ielnum)*volel0*dvdx
275  f22 = f22 + alpha(ix,ielnum)*volel0*(1.d0 + dvdy)
276  f23 = f23 + alpha(ix,ielnum)*volel0*dvdz
277  f31 = f31 + alpha(ix,ielnum)*volel0*dwdx
278  f32 = f32 + alpha(ix,ielnum)*volel0*dwdy
279  f33 = f33 + alpha(ix,ielnum)*volel0*(1.d0 + dwdz)
280 
281  ENDDO
282 
283 ! 2nd "for each node" in box starts here
284 
285 
286  ! first bullet in box
287 
288  f11 = volnd0inv*f11
289  f22 = volnd0inv*f22
290  f33 = volnd0inv*f33
291  f12 = volnd0inv*f12
292  f21 = volnd0inv*f21
293  f23 = volnd0inv*f23
294  f32 = volnd0inv*f32
295  f13 = volnd0inv*f13
296  f31 = volnd0inv*f31
297 
298  j1 = f11*(f22*f33-f23*f32)+f12*(f31*f23-f21*f33)+f13*(-f31*f22+f21*f32)
299  IF(j1.LE.0.d0) THEN
300  WRITE(*,*)'area has become zero for element',i
301  stop
302  ENDIF
303  j2 = j1*j1
304 
305 ! T
306 ! C = F F = right Cauchy-Green deformation tensor
307 !
308 ! Divided by third invariant squared
309 
310  j2inv = 1.d0/j2
311 
312  c11 = (f11*f11+f21*f21+f31*f31)*j2inv
313  c12 = (f11*f12+f21*f22+f31*f32)*j2inv
314  c13 = (f11*f13+f21*f23+f31*f33)*j2inv
315  c21 = c12
316  c22 = (f12*f12+f22*f22+f32*f32)*j2inv
317  c23 = (f12*f13+f22*f23+f32*f33)*j2inv
318  c31 = c13
319  c32 = c23
320  c33 = (f13*f13+f23*f23+f33*f33)*j2inv
321 
322 !
323 ! Second Piola-Kirchoff tensor
324 ! Eq. (5.28), pg. 124
325 
326  s11n(i) = xmu(1)*(1.d0-(c22*c33-c23*c32)) + xlambda(1)*log(j1)*(c22*c33-c23*c32)
327  s22n(i) = xmu(1)*(1.d0-(c11*c33-c31*c13)) + xlambda(1)*log(j1)*(c11*c33-c31*c13)
328  s33n(i) = xmu(1)*(1.d0-(c11*c22-c12*c21)) + xlambda(1)*log(j1)*(c11*c22-c12*c21)
329  s12n(i) = (-c12*c33+c13*c32)*(xmu(1) - xlambda(1)*log(j1))
330  s23n(i) = (-c11*c23+c13*c21)*(xmu(1) - xlambda(1)*log(j1))
331  s13n(i) = (c12*c23-c13*c22)*(xmu(1) - xlambda(1)*log(j1))
332 
333  ENDDO
334 
335 ! For each node
336 
337  DO i = 1, numnp
338 
339  k3i = 3*i
340  k2i = 3*i-1
341  k1i = 3*i-2
342 
343 ! For each element
344 
345  DO j = 1, numelneigh(i)
346 
347  ielnum = elconn(i,j)
348 
349  n1 = nodes(1,ielnum)
350  n2 = nodes(2,ielnum)
351  n3 = nodes(3,ielnum)
352  n4 = nodes(4,ielnum)
353  k3n1 = 3*n1
354  k3n2 = 3*n2
355  k3n3 = 3*n3
356  k3n4 = 3*n4
357 
358  k2n1 = k3n1 - 1
359  k2n2 = k3n2 - 1
360  k2n3 = k3n3 - 1
361  k2n4 = k3n4 - 1
362 
363  k1n1 = k3n1 - 2
364  k1n2 = k3n2 - 2
365  k1n3 = k3n3 - 2
366  k1n4 = k3n4 - 2
367  ! k#n# dummy variables replaces:
368  u1 = disp(k1n1) ! (3*n1-2)
369  u2 = disp(k1n2) ! (3*n2-2)
370  u3 = disp(k1n3) ! (3*n3-2)
371  u4 = disp(k1n4) ! (3*n4-2)
372  v1 = disp(k2n1) ! (3*n1-1)
373  v2 = disp(k2n2) ! (3*n2-1)
374  v3 = disp(k2n3) ! (3*n3-1)
375  v4 = disp(k2n4) ! (3*n4-1)
376  w1 = disp(k3n1) ! (3*n1)
377  w2 = disp(k3n2) ! (3*n2)
378  w3 = disp(k3n3) ! (3*n3)
379  w4 = disp(k3n4) ! (3*n4)
380 
381  x1 = coor(1,n1)
382  x2 = coor(1,n2)
383  x3 = coor(1,n3)
384  x4 = coor(1,n4)
385  y1 = coor(2,n1)
386  y2 = coor(2,n2)
387  y3 = coor(2,n3)
388  y4 = coor(2,n4)
389  z1 = coor(3,n1)
390  z2 = coor(3,n2)
391  z3 = coor(3,n3)
392  z4 = coor(3,n4)
393 
394  x12 = x1 - x2 ! not used in vol. calc
395  x13 = x1 - x3 ! not used in vol. calc
396  x14 = x1 - x4
397  x24 = x2 - x4
398  x34 = x3 - x4
399  y12 = y1 - y2 ! not used in vol. calc
400  y13 = y1 - y3 ! not used in vol. calc
401  y14 = y1 - y4
402  y24 = y2 - y4
403  y34 = y3 - y4
404  z12 = z1 - z2 ! not used in vol. calc
405  z13 = z1 - z3 ! not used in vol. calc
406  z14 = z1 - z4
407  z24 = z2 - z4
408  z34 = z3 - z4
409 
410  c11 = y24*z34 - z24*y34
411  c21 = -( x24*z34 - z24*x34 )
412  c31 = x24*y34 - y24*x34
413 
414  vx6 = -( x14*c11 + y14*c21 + z14*c31 )
415  volel0 = vx6/6.d0
416 
417 ! See the maple worksheet 'V3D4.mws' for the derivation of Vx6
418 
419  vx6inv = 1.d0 / vx6
420 
421 ! See the maple worksheet 'V3D4.mws' for the derivation of [B]! Compute the Shape functions
422 ! NOTE: Factored for a more equivalent/compact form then maple's
423 
424  b1 = (y34*z24 - y24*z34) * vx6inv
425  b2 = (z34*x24 - z24*x34) * vx6inv
426  b3 = (x34*y24 - x24*y34) * vx6inv
427  b4 = (y13*z14 - y14*z13) * vx6inv
428  b5 = (z13*x14 - z14*x13) * vx6inv
429  b6 = (x13*y14 - x14*y13) * vx6inv
430  b7 = (y14*z12 - y12*z14) * vx6inv
431  b8 = (z14*x12 - z12*x14) * vx6inv
432  b9 = (x14*y12 - x12*y14) * vx6inv
433  b10 = (y12*z13 - y13*z12) * vx6inv
434  b11 = (z12*x13 - z13*x12) * vx6inv
435  b12 = (x12*y13 - x13*y12) * vx6inv
436 !
437 ! Calculate displacement gradient (H)
438 !
439  dudx = b1*u1 + b4*u2 + b7*u3 + b10*u4
440  dvdy = b2*v1 + b5*v2 + b8*v3 + b11*v4
441  dwdz = b3*w1 + b6*w2 + b9*w3 + b12*w4
442  dudy = b2*u1 + b5*u2 + b8*u3 + b11*u4
443  dvdx = b1*v1 + b4*v2 + b7*v3 + b10*v4
444  dvdz = b3*v1 + b6*v2 + b9*v3 + b12*v4
445  dwdy = b2*w1 + b5*w2 + b8*w3 + b11*w4
446  dudz = b3*u1 + b6*u2 + b9*u3 + b12*u4
447  dwdx = b1*w1 + b4*w2 + b7*w3 + b10*w4
448 !
449 ! deformation gradients F
450 !
451  f11 = (1.d0 + dudx)
452  f12 = dudy
453  f13 = dudz
454  f21 = dvdx
455  f22 = (1.d0 + dvdy)
456  f23 = dvdz
457  f31 = dwdx
458  f32 = dwdy
459  f33 = (1.d0 + dwdz)
460 
461 ! 4
462 ! ---- ,
463 ! Evaluate \ P
464 ! / a
465 ! ----
466 ! a = 1
467 
468  s11(1) = 0.d0
469  s22(1) = 0.d0
470  s33(1) = 0.d0
471  s12(1) = 0.d0
472  s23(1) = 0.d0
473  s13(1) = 0.d0
474 
475  DO k = 1, 4
476  s11(1) = s11(1) + alpha(k,ielnum)*s11n(nodes(k,ielnum))
477  s22(1) = s22(1) + alpha(k,ielnum)*s22n(nodes(k,ielnum))
478  s33(1) = s33(1) + alpha(k,ielnum)*s33n(nodes(k,ielnum))
479  s12(1) = s12(1) + alpha(k,ielnum)*s12n(nodes(k,ielnum))
480  s23(1) = s23(1) + alpha(k,ielnum)*s23n(nodes(k,ielnum))
481  s13(1) = s13(1) + alpha(k,ielnum)*s13n(nodes(k,ielnum))
482  ENDDO
483 
484 ! calculate the volume
485 
486  vol = vx6 / 6.d0
487 
488 !!$
489 !!$! ASSEMBLE THE INTERNAL FORCE VECTOR
490 !!$!
491 !!$! local node 1
492 !!$
493 !!$ Rnet(k1i) = Rnet(k1i) - VolEl0* &
494 !!$ ( S11e*B1*(1.d0+dudx) + S22e*B2*dudy + S33e*B3*dudz &
495 !!$ + S12e*( B2*(1.d0+dudx) + B1*dudy ) &
496 !!$ + S23e*( B3*dudy + B2*dudz ) &
497 !!$ + S13e*( B3*(1.d0+dudx) + B1*dudz ) )
498 !!$ Rnet(k2i) = Rnet(k2i) - VolEl0* &
499 !!$ ( S11e*B1*dvdx + S22e*B2*(1.d0+dvdy) + S33e*B3*dvdz &
500 !!$ + S12e*( B1*(1.d0+dvdy) + B2*dvdx ) &
501 !!$ + S23e*( B3*(1.d0+dvdy) + B2*dvdz ) &
502 !!$ + S13e*( B3*dvdx + B1*dvdz ) )
503 !!$ Rnet(k3i) = Rnet(k3i) - VolEl0* &
504 !!$ ( S11e*B1*dwdx + S22e*B2*dwdy + S33e*B3*(1.d0+dwdz) &
505 !!$ + S12e*( B2*dwdx + B1*dwdy ) &
506 !!$ + S23e*( B3*dwdy + B2*(1.d0 + dwdz) ) &
507 !!$ + S13e*( B3*dwdx + B1*(1.d0 + dwdz) ) )
508 !!$! local node 2
509 !!$ Rnet(k1i) = Rnet(k1i) - VolEl0* &
510 !!$ ( S11e*B4*(1.d0+dudx) + S22e*B5*dudy + S33e*B6*dudz &
511 !!$ + S12e*( B5*(1.d0+dudx) + B4*dudy ) &
512 !!$ + S23e*( B6*dudy + B5*dudz ) &
513 !!$ + S13e*( B6*(1.d0+dudx) + B4*dudz ) )
514 !!$ Rnet(k2i) = Rnet(k2i) - VolEl0* &
515 !!$ ( S11e*B4*dvdx + S22e*B5*(1.d0+dvdy) + S33e*B6*dvdz &
516 !!$ + S12e*( B4*(1.d0+dvdy) + B5*dvdx ) &
517 !!$ + S23e*( B6*(1.d0+dvdy) + B5*dvdz ) &
518 !!$ + S13e*( B6*dvdx + B4*dvdz ) )
519 !!$ Rnet(k3i) = Rnet(k3i) - VolEl0* &
520 !!$ ( S11e*B4*dwdx + S22e*B5*dwdy + S33e*B6*(1.d0+dwdz) &
521 !!$ + S12e*( B5*dwdx + B4*dwdy ) &
522 !!$ + S23e*( B6*dwdy + B5*(1.d0 + dwdz) ) &
523 !!$ + S13e*( B6*dwdx + B4*(1.d0 + dwdz) ) )
524 !!$! local node 3
525 !!$ Rnet(k1i) = Rnet(k1i) - VolEl0* &
526 !!$ ( S11e*B7*(1.d0+dudx) + S22e*B8*dudy + S33e*B9*dudz &
527 !!$ + S12e*( B8*(1.d0+dudx) + B7*dudy ) &
528 !!$ + S23e*( B9*dudy + B8*dudz ) &
529 !!$ + S13e*( B9*(1.d0+dudx) + B7*dudz ) )
530 !!$ Rnet(k2i) = Rnet(k2i) - VolEl0* &
531 !!$ ( S11e*B7*dvdx + S22e*B8*(1.d0+dvdy) + S33e*B9*dvdz &
532 !!$ + S12e*( B7*(1.d0+dvdy) + B8*dvdx ) &
533 !!$ + S23e*( B9*(1.d0+dvdy) + B8*dvdz ) &
534 !!$ + S13e*( B9*dvdx + B7*dvdz ) )
535 !!$ Rnet(k3i) = Rnet(k3i) - VolEl0* &
536 !!$ ( S11e*B7*dwdx + S22e*B8*dwdy + S33e*B9*(1.d0+dwdz) &
537 !!$ + S12e*( B8*dwdx + B7*dwdy ) &
538 !!$ + S23e*( B9*dwdy + B8*(1.d0 + dwdz) ) &
539 !!$ + S13e*( B9*dwdx + B7*(1.d0 + dwdz) ) )
540 !!$! local node 4
541 !!$ Rnet(k1i) = Rnet(k1i) - VolEl0* &
542 !!$ ( S11e*B10*(1.d0+dudx) + S22e*B11*dudy+S33e*B12*dudz &
543 !!$ + S12e*( B11*(1.d0+dudx) + B10*dudy ) &
544 !!$ + S23e*( B12*dudy + B11*dudz ) &
545 !!$ + S13e*( B12*(1.d0+dudx) + B10*dudz ) )
546 !!$ Rnet(k2i) = Rnet(k2i) - VolEl0* &
547 !!$ ( S11e*B10*dvdx + S22e*B11*(1.d0+dvdy)+S33e*B12*dvdz &
548 !!$ + S12e*( B10*(1.d0+dvdy) + B11*dvdx ) &
549 !!$ + S23e*( B12*(1.d0+dvdy) + B11*dvdz ) &
550 !!$ + S13e*( B12*dvdx + B10*dvdz ) )
551 !!$ Rnet(k3i) = Rnet(k3i) - VolEl0* &
552 !!$ ( S11e*B10*dwdx + S22e*B11*dwdy+S33e*B12*(1.d0+dwdz) &
553 !!$ + S12e*( B11*dwdx + B10*dwdy ) &
554 !!$ + S23e*( B12*dwdy + B11*(1.d0 + dwdz) ) &
555 !!$ + S13e*( B12*dwdx + B10*(1.d0 + dwdz) ) )
556 
557  rnet(k1n1) = rnet(k1n1) - volel0* &
558  ( s11(1)*b1*f11 + s22(1)*b2*f12 + s33(1)*b3*f13 &
559  + s12(1)*( b2*f11 + b1*f12 ) &
560  + s23(1)*( b3*f12 + b2*f13 ) &
561  + s13(1)*( b3*f11 + b1*f13 ) )
562  rnet(k2n1) = rnet(k2n1) - volel0* &
563  ( s11(1)*b1*f21 + s22(1)*b2*f22 + s33(1)*b3*f23 &
564  + s12(1)*( b1*f22 + b2*f21 ) &
565  + s23(1)*( b3*f22 + b2*f23 ) &
566  + s13(1)*( b3*f21 + b1*f23 ) )
567  rnet(k3n1) = rnet(k3n1) - volel0* &
568  ( s11(1)*b1*f31 + s22(1)*b2*f32 + s33(1)*b3*f33 &
569  + s12(1)*( b2*f31 + b1*f32 ) &
570  + s23(1)*( b3*f32 + b2*f33 ) &
571  + s13(1)*( b3*f31 + b1*f33 ) )
572 ! local node 2
573  rnet(k1n2) = rnet(k1n2) - volel0* &
574  ( s11(1)*b4*f11 + s22(1)*b5*f12 + s33(1)*b6*f13 &
575  + s12(1)*( b5*f11 + b4*f12 ) &
576  + s23(1)*( b6*f12 + b5*f13 ) &
577  + s13(1)*( b6*f11 + b4*f13 ) )
578  rnet(k2n2) = rnet(k2n2) - volel0* &
579  ( s11(1)*b4*f21 + s22(1)*b5*f22 + s33(1)*b6*f23 &
580  + s12(1)*( b4*f22 + b5*f21 ) &
581  + s23(1)*( b6*f22 + b5*f23 ) &
582  + s13(1)*( b6*f21 + b4*f23 ) )
583  rnet(k3n2) = rnet(k3n2) - volel0* &
584  ( s11(1)*b4*f31 + s22(1)*b5*f32 + s33(1)*b6*f33 &
585  + s12(1)*( b5*f31 + b4*f32 ) &
586  + s23(1)*( b6*f32 + b5*f33 ) &
587  + s13(1)*( b6*f31 + b4*f33 ) )
588 ! local node 3
589  rnet(k1n3) = rnet(k1n3) - volel0* &
590  ( s11(1)*b7*f11 + s22(1)*b8*f12 + s33(1)*b9*f13 &
591  + s12(1)*( b8*f11 + b7*f12 ) &
592  + s23(1)*( b9*f12 + b8*f13 ) &
593  + s13(1)*( b9*f11 + b7*f13 ) )
594  rnet(k2n3) = rnet(k2n3) - volel0* &
595  ( s11(1)*b7*f21 + s22(1)*b8*f22 + s33(1)*b9*f23 &
596  + s12(1)*( b7*f22 + b8*f21 ) &
597  + s23(1)*( b9*f22 + b8*f23 ) &
598  + s13(1)*( b9*f21 + b7*f23 ) )
599  rnet(k3n3) = rnet(k3n3) - volel0* &
600  ( s11(1)*b7*f31 + s22(1)*b8*f32 + s33(1)*b9*f33 &
601  + s12(1)*( b8*f31 + b7*f32 ) &
602  + s23(1)*( b9*f32 + b8*f33 ) &
603  + s13(1)*( b9*f31 + b7*f33 ) )
604 ! local node 4
605  rnet(k1n4) = rnet(k1n4) - volel0* &
606  ( s11(1)*b10*f11 + s22(1)*b11*f12+s33(1)*b12*f13 &
607  + s12(1)*( b11*f11 + b10*f12 ) &
608  + s23(1)*( b12*f12 + b11*f13 ) &
609  + s13(1)*( b12*f11 + b10*f13 ) )
610  rnet(k2n4) = rnet(k2n4) - volel0* &
611  ( s11(1)*b10*f21 + s22(1)*b11*f22 + s33(1)*b12*f23 &
612  + s12(1)*( b10*f22 + b11*f21 ) &
613  + s23(1)*( b12*f22 + b11*f23 ) &
614  + s13(1)*( b12*f21 + b10*f23 ) )
615  rnet(k3n4) = rnet(k3n4) - volel0* &
616  ( s11(1)*b10*f31 + s22(1)*b11*f32 + s33(1)*b12*f33 &
617  + s12(1)*( b11*f31 + b10*f32 ) &
618  + s23(1)*( b12*f32 + b11*f33 ) &
619  + s13(1)*( b12*f31 + b10*f33 ) )
620 
621  ENDDO
622  ENDDO
623 
624  RETURN
625 END SUBROUTINE v3d4n_neohookeancompress
626 
rocstar_rocfraccomm
Definition: ROCSTAR_RocFracComm.f90:53
k
j indices k indices k
Definition: Indexing.h:6
i
blockLoc i
Definition: read.cpp:79
j
j indices j
Definition: Indexing.h:6
Color::alpha
unsigned char alpha() const
Definition: Color.h:75
v3d4n_neohookeancompress
subroutine v3d4n_neohookeancompress(numnp, numel, coor, disp, nodes, Rnet, NumElNeigh, ElConn, alpha, Ahat, NumMatVol, xmu, xlambda)
Definition: v3d4n_NeoHookeanCompress.f90:53