Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d4n_NeoHookeanInCompress.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_neohookeanincompress(numnp, numel, coor, disp, nodes, Rnet, &
54  s11ev, s22ev, s33ev, s12ev, s23ev, s13ev, &
55  numelneigh,elconn,alpha, ahat,&
56  xmu,xkappa,myid, &
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  USE rocstar_rocfraccomm
62 
63  IMPLICIT NONE
64 
65  INTEGER :: myid
66 
67  INTEGER :: totnumneighprocs, nprocs,totnumndcomm
68  INTEGER, DIMENSION(1:TotNumNeighProcs) :: neighproclist
69  INTEGER, DIMENSION(1:TotNumNeighProcs) :: numndcomm
70  REAL*8, pointer, DIMENSION(:) :: buf
71 
72  TYPE(rcv_buf), pointer, DIMENSION(:) :: recvdatafrm
73  TYPE(send_buf), DIMENSION(1:TotNumNeighProcs) :: neigh_lst
74 
75  INTEGER :: ierr
76 
77 
78  INTEGER :: numnp, numel
79  REAL*8, DIMENSION(1:numnp) :: ahat
80  REAL*8, DIMENSION(1:numnp*3) :: disp, rnet
81  REAL*8, DIMENSION(1:3,1:numnp) :: coor
82  INTEGER, DIMENSION(1:numnp) ::numelneigh
83  INTEGER, DIMENSION(1:numnp,1:40) :: elconn ! fix 40 should not be hard coded
84  INTEGER,DIMENSION(1:4,1:numel) :: nodes
85  REAL*8, DIMENSION(1:4,1:numel) :: alpha
86  REAL*8 :: xmu, xkappa
87  REAL*8, DIMENSION(1:numel) :: s11ev, s22ev, s33ev, s12ev, s23ev, s13ev
88  REAL*8 :: s11e, s22e, s33e, s12e, s23e, s13e, s21e, s31e, s32e
89 
90 
91  REAL*8 :: onethird = 1.d0/3.d0
92 
93  INTEGER :: j,i,k
94  REAL*8 :: aaa
95 !-- coordinate holding variable
96  REAL*8 :: x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
97  REAL*8 :: x1d,x2d,x3d,x4d,y1d,y2d,y3d,y4d,z1d,z2d,z3d,z4d
98  INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
99  INTEGER :: k3n1,k3n2,k3n3,k3n4
100 
101  REAL*8 :: ahatinv
102  INTEGER :: ielnum
103  REAL*8 :: sixinv
104 !-- x, y, and z displacements of nodes
105  REAL*8 :: u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4
106 !-- partial derivatives of the displacement
107  REAL*8 :: dudx,dvdy,dwdz,dudy,dvdx,dvdz,dwdy,dudz,dwdx
108 
109 !-- Coordinate subtractions
110  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
111 !-- Added these to speed up B calculation
112 
113  REAL*8 :: z12, z13,x12,x13,y12,y13
114  REAL*8 :: c11, c12, c13, c21, c22, c23, c31, c32, c33
115  REAL*8 :: sh1, sh2, sh3, sh4, sh5, sh6, sh7, sh8, sh9, sh10, sh11, sh12
116  REAL*8 :: sh13, sh14, sh15
117 ! -- 6*volume and the volume
118  REAL*8 :: vx6
119 !-- spacial derivatives
120  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
121 !-- strains
122  REAL*8 :: e11,e22,e33,e12,e23,e13
123  REAL*8 :: b13, b14, b15
124  REAL*8 :: sums11n,sums22n,sums33n,sums12n,sums23n,sums13n
125  REAL*8 :: sume11n,sume22n,sume33n,sume12n,sume23n,sume13n
126 
127  REAL*8 :: e11n, e22n, e33n, e12n, e23n, e13n
128  INTEGER :: k3i,k2i,k1i
129  INTEGER :: ix
130  INTEGER :: n1, n2, n3, n4
131  REAL*8 :: volel, volel0,vx6inv,vol
132  REAL*8 :: vainv
133 
134  REAL*8 :: coeff1, coeff2
135  REAL*8 :: jacob, detfa
136  REAL*8 :: f11, f12, f13, f21, f22, f23, f31, f32, f33
137  REAL*8,DIMENSION(1:numnp) :: f11v, f12v, f13v, f21v, f22v, f23v, f31v, f32v, f33v, voldefv
138  REAL*8 :: vx6def, vx6invdef, voldef
139  REAL*8 :: vtotal
140  REAL*8, DIMENSION(1:numnp) :: pa
141  REAL*8 :: cnode(1:3,1:3)
142  REAL*8 :: iiic, doubleproduct, pc
143  REAL*8, DIMENSION(1:numnp) :: s11n, s22n, s33n, s12n, s23n, s13n, s21n, s31n,s32n
144 
145  INTEGER :: k1, j1, k2
146 
147  integer :: mpi_status_size,mpi_comm_rocfrac,mpi_double_precision
148 
149 !-- Non-block receive, Non-block send request arrays
150  INTEGER, DIMENSION(1:TotNumNeighProcs) :: reqrcv, reqsnd
151  INTEGER, DIMENSION(1:MPI_STATUS_SIZE,1:TotNumNeighProcs) :: statrcv, statsnd
152 
153 
154  sixinv = 1.d0/6.d0
155 
156  DO i = 1, numnp ! for each node
157 
158 ! /* Undeformed Volume of node */
159 
160  vainv = 1.d0/ahat(i)
161 
162  ! /* Initialize to 0, Va, Va0, Fa */
163 
164  f11 = 0.d0
165  f12 = 0.d0
166  f13 = 0.d0
167  f21 = 0.d0
168  f22 = 0.d0
169  f23 = 0.d0
170  f31 = 0.d0
171  f32 = 0.d0
172  f33 = 0.d0
173 
174  voldef = 0.d0
175 
176  DO j = 1, numelneigh(i) ! for each element associated with node i
177 
178  ielnum = elconn(i,j)
179 
180  n1 = nodes(1,ielnum)
181  n2 = nodes(2,ielnum)
182  n3 = nodes(3,ielnum)
183  n4 = nodes(4,ielnum)
184 
185  IF(i.EQ.n1) THEN
186  ix = 1
187  ELSE IF(i.EQ.n2)THEN
188  ix = 2
189  ELSE IF(i.EQ.n3)THEN
190  ix = 3
191  ELSE IF(i.EQ.n4)THEN
192  ix = 4
193  ENDIF
194 
195  k3n1 = 3*n1
196  k3n2 = 3*n2
197  k3n3 = 3*n3
198  k3n4 = 3*n4
199 
200  k2n1 = k3n1 - 1
201  k2n2 = k3n2 - 1
202  k2n3 = k3n3 - 1
203  k2n4 = k3n4 - 1
204 
205  k1n1 = k3n1 - 2
206  k1n2 = k3n2 - 2
207  k1n3 = k3n3 - 2
208  k1n4 = k3n4 - 2
209 
210  ! k#n# dummy variables replaces:
211  u1 = disp(k1n1) ! (3*n1-2)
212  u2 = disp(k1n2) ! (3*n2-2)
213  u3 = disp(k1n3) ! (3*n3-2)
214  u4 = disp(k1n4) ! (3*n4-2)
215  v1 = disp(k2n1) ! (3*n1-1)
216  v2 = disp(k2n2) ! (3*n2-1)
217  v3 = disp(k2n3) ! (3*n3-1)
218  v4 = disp(k2n4) ! (3*n4-1)
219  w1 = disp(k3n1) ! (3*n1)
220  w2 = disp(k3n2) ! (3*n2)
221  w3 = disp(k3n3) ! (3*n3)
222  w4 = disp(k3n4) ! (3*n4)
223 
224  x1 = coor(1,n1)
225  x2 = coor(1,n2)
226  x3 = coor(1,n3)
227  x4 = coor(1,n4)
228  y1 = coor(2,n1)
229  y2 = coor(2,n2)
230  y3 = coor(2,n3)
231  y4 = coor(2,n4)
232  z1 = coor(3,n1)
233  z2 = coor(3,n2)
234  z3 = coor(3,n3)
235  z4 = coor(3,n4)
236 
237  x12 = x1 - x2 ! not used in vol. calc
238  x13 = x1 - x3 ! not used in vol. calc
239  x14 = x1 - x4
240  x24 = x2 - x4
241  x34 = x3 - x4
242  y12 = y1 - y2 ! not used in vol. calc
243  y13 = y1 - y3 ! not used in vol. calc
244  y14 = y1 - y4
245  y24 = y2 - y4
246  y34 = y3 - y4
247  z12 = z1 - z2 ! not used in vol. calc
248  z13 = z1 - z3 ! not used in vol. calc
249  z14 = z1 - z4
250  z24 = z2 - z4
251  z34 = z3 - z4
252 
253  c11 = y24*z34 - z24*y34
254  c21 = -( x24*z34 - z24*x34 )
255  c31 = x24*y34 - y24*x34
256 
257  vx6 = -( x14*c11 + y14*c21 + z14*c31 )
258 
259  vx6inv = 1.d0/vx6
260 
261  volel0 = vx6/6.d0 ! undeformed volume of element (Ve_O)
262 
263 ! Compute the Shape functions
264 ! NOTE: Factored for a more equivalent/compact form then maple's
265 
266  b1 = (y34*z24 - y24*z34) * vx6inv
267  b2 = (z34*x24 - z24*x34) * vx6inv
268  b3 = (x34*y24 - x24*y34) * vx6inv
269  b4 = (y13*z14 - y14*z13) * vx6inv
270  b5 = (z13*x14 - z14*x13) * vx6inv
271  b6 = (x13*y14 - x14*y13) * vx6inv
272  b7 = (y14*z12 - y12*z14) * vx6inv
273  b8 = (z14*x12 - z12*x14) * vx6inv
274  b9 = (x14*y12 - x12*y14) * vx6inv
275  b10 = (y12*z13 - y13*z12) * vx6inv
276  b11 = (z12*x13 - z13*x12) * vx6inv
277  b12 = (x12*y13 - x13*y12) * vx6inv
278 !
279 ! deformation gradients F
280 !
281  f11 = f11 + alpha(ix,ielnum)*volel0*(1.d0 + ( b1*u1 + b4*u2 + b7*u3 + b10*u4 )) ! 1 + ( dudx )
282  f22 = f22 + alpha(ix,ielnum)*volel0*(1.d0 + ( b2*v1 + b5*v2 + b8*v3 + b11*v4 )) ! 1 + ( dvdy )
283  f33 = f33 + alpha(ix,ielnum)*volel0*(1.d0 + ( b3*w1 + b6*w2 + b9*w3 + b12*w4 )) ! 1 + ( dwdz )
284  f12 = f12 + alpha(ix,ielnum)*volel0*(b2*u1 + b5*u2 + b8*u3 + b11*u4) ! dudy
285  f21 = f21 + alpha(ix,ielnum)*volel0*(b1*v1 + b4*v2 + b7*v3 + b10*v4) ! dvdx
286  f23 = f23 + alpha(ix,ielnum)*volel0*(b3*v1 + b6*v2 + b9*v3 + b12*v4) ! dvdz
287  f32 = f32 + alpha(ix,ielnum)*volel0*(b2*w1 + b5*w2 + b8*w3 + b11*w4) ! dwdy
288  f13 = f13 + alpha(ix,ielnum)*volel0*(b3*u1 + b6*u2 + b9*u3 + b12*u4) ! dudz
289  f31 = f31 + alpha(ix,ielnum)*volel0*(b1*w1 + b4*w2 + b7*w3 + b10*w4) ! dwdx
290 
291  x1 = x1 + u1
292  x2 = x2 + u2
293  x3 = x3 + u3
294  x4 = x4 + u4
295  y1 = y1 + v1
296  y2 = y2 + v2
297  y3 = y3 + v3
298  y4 = y4 + v4
299  z1 = z1 + w1
300  z2 = z2 + w2
301  z3 = z3 + w3
302  z4 = z4 + w4
303 
304  x14 = x1 - x4
305  x24 = x2 - x4
306  x34 = x3 - x4
307  y14 = y1 - y4
308  y24 = y2 - y4
309  y34 = y3 - y4
310  z14 = z1 - z4
311  z24 = z2 - z4
312  z34 = z3 - z4
313 
314  c11 = y24*z34 - z24*y34
315  c21 = -( x24*z34 - z24*x34 )
316  c31 = x24*y34 - y24*x34
317 
318  vx6def = -( x14*c11 + y14*c21 + z14*c31 )
319 
320 ! calculate the volume
321  voldef = voldef + alpha(ix,ielnum)*vx6def/6.d0
322 
323  ENDDO
324 
325  f11v(i) = f11*vainv
326  f22v(i) = f22*vainv
327  f33v(i) = f33*vainv
328  f12v(i) = f12*vainv
329  f13v(i) = f13*vainv
330  f21v(i) = f21*vainv
331  f23v(i) = f23*vainv
332  f31v(i) = f31*vainv
333  f32v(i) = f32*vainv
334 
335  voldefv(i) = voldef
336 
337  ENDDO
338 
339 
340  ALLOCATE(recvdatafrm(0:nprocs-1))
341 
342 !
343 !----- FORM THE BUFFER CONTAINING COMMUNICATED STRESS MATRIX NODAL VALUES
344 
345  ALLOCATE(buf(1:totnumndcomm/3*10)) ! is really (TotNumNdComm*3 * 3 + 1)
346  k1 = 1
347  DO j1 = 1, totnumneighprocs
348  k = neighproclist(j1)
349  ALLOCATE(recvdatafrm(k)%rcvbuf(1:numndcomm(j1)*10))
350  DO j = 1, numndcomm(j1)
351  k2 = neigh_lst(j1)%NdId(j) !NdCommList(k)%NdId(j)
352  ! print*,myid,k2,S23n(k2)
353  buf(k1 ) = f11v(k2)
354  buf(k1+1) = f12v(k2)
355  buf(k1+2) = f13v(k2)
356  buf(k1+3) = f21v(k2)
357  buf(k1+4) = f22v(k2)
358  buf(k1+5) = f23v(k2)
359  buf(k1+6) = f31v(k2)
360  buf(k1+7) = f32v(k2)
361  buf(k1+8) = f33v(k2)
362  buf(k1+9) = voldefv(k2)
363  k1 = k1 + 10
364  ENDDO
365  ENDDO
366 !
367 !-MPI- RECEIVE THE RECIPRICAL MASS MATRIX DIAGONAL FROM THE NEIGHBORS
368 !
369  DO j1 = 1, totnumneighprocs
370  k = neighproclist(j1)
371  CALL mpi_irecv(recvdatafrm(k)%rcvbuf(1),numndcomm(j1)*10, &
372  mpi_double_precision, k, 10,mpi_comm_rocfrac,reqrcv(j1),ierr)
373  ENDDO
374 !
375 !-MPI- SEND THE RECIPRICAL MASS MATRIX DIAGONAL TO THE NEIGHBORS
376 !
377  k2 = 1
378  DO j1 = 1, totnumneighprocs
379  k = neighproclist(j1)
380 
381 ! IF(myid.EQ.0) print*,'sending',buf(1:)
382  CALL mpi_isend(buf(k2),numndcomm(j1)*10,&
383  mpi_double_precision,k,10,mpi_comm_rocfrac,reqsnd(j1),ierr)
384  k2 = k2 + numndcomm(j1)*10
385  ENDDO
386 !
387 !-MPI- WAIT FOR INTERNAL FORCE VECTOR COMMUNICATION TO COMPLETE
388 !
389  IF(totnumneighprocs.GT.0)THEN
390  CALL mpi_waitall(totnumneighprocs,reqrcv,statrcv,ierr)
391  CALL mpi_waitall(totnumneighprocs,reqsnd,statsnd,ierr)
392  ENDIF
393  DEALLOCATE(buf)
394 
395 !
396 !----- ADD NEIGHBOR'S CONTRIBUTION TO THE INTERNAL FORCE VECTOR
397 !
398  DO j1 = 1, totnumneighprocs
399  k = neighproclist(j1)
400  k1 = 1
401 ! IF(myid.EQ.1) print*,'***',RecvDataFrm(0)%rcvbuf(:)
402  DO j = 1, numndcomm(j1)
403  k2 = neigh_lst(j1)%NdId(j)
404  f11v(k2) = f11v(k2) + recvdatafrm(k)%rcvbuf(k1)
405  f12v(k2) = f12v(k2) + recvdatafrm(k)%rcvbuf(k1+1)
406  f13v(k2) = f13v(k2) + recvdatafrm(k)%rcvbuf(k1+2)
407  f21v(k2) = f21v(k2) + recvdatafrm(k)%rcvbuf(k1+3)
408  f22v(k2) = f22v(k2) + recvdatafrm(k)%rcvbuf(k1+4)
409  f23v(k2) = f23v(k2) + recvdatafrm(k)%rcvbuf(k1+5)
410  f31v(k2) = f31v(k2) + recvdatafrm(k)%rcvbuf(k1+6)
411  f32v(k2) = f32v(k2) + recvdatafrm(k)%rcvbuf(k1+7)
412  f33v(k2) = f33v(k2) + recvdatafrm(k)%rcvbuf(k1+8)
413  voldefv(k2) = voldefv(k2) + recvdatafrm(k)%rcvbuf(k1+9)
414  k1 = k1 + 10
415  ENDDO
416  ENDDO
417 
418  DEALLOCATE(recvdatafrm)
419 
420  DO i = 1, numnp ! for each node
421 
422  f11 = f11v(i)
423  f22 = f22v(i)
424  f33 = f33v(i)
425  f12 = f12v(i)
426  f13 = f13v(i)
427  f21 = f21v(i)
428  f23 = f23v(i)
429  f31 = f31v(i)
430  f32 = f32v(i)
431 
432  voldef = voldefv(i)
433 
434 ! /* Undeformed Volume of node */
435 
436  vainv = 1.d0/ahat(i)
437 
438 ! Jacob = F11*(F22*F33-F23*F32)+F12*(F31*F23-F21*F33)+F13*(-F31*F22+F21*F32)
439 
440  jacob = voldef*vainv
441 
442  IF( jacob.LE.0.d0) THEN
443  WRITE(*,100) i
444  stop
445  ENDIF
446 
447  pa(i) = xkappa*(jacob - 1.d0)
448 
449 ! Third invariant of C
450 
451  detfa = f11*(f22*f33-f23*f32)+f12*(f31*f23-f21*f33)+f13*(-f31*f22+f21*f32)
452 
453  iiic = jacob*jacob
454 
455 !
456 ! Second Piola-Kirchoff tensor
457 ! Eq. (5.28), pg. 124
458 
459 
460  coeff1 = xmu*detfa**(-2.d0/3.d0)
461 !
462 ! F:F
463 !
464  doubleproduct = f11**2+f21**2+f31**2+f12**2+f22**2+f32**2+f13**2+f23**2+f33**2
465 
466  cnode(1,1) = (f22*f33-f23*f32)/detfa*doubleproduct
467  cnode(1,2) = -(f21*f33-f23*f31)/detfa*doubleproduct
468  cnode(1,3) = (f32*f21-f31*f22)/detfa*doubleproduct
469  cnode(2,1) = -(f12*f33-f13*f32)/detfa*doubleproduct
470  cnode(2,2) = (f11*f33-f13*f31)/detfa*doubleproduct
471  cnode(2,3) = -(f11*f32-f31*f12)/detfa*doubleproduct
472  cnode(3,1) = (f12*f23-f22*f13)/detfa*doubleproduct
473  cnode(3,2) = -(f11*f23-f21*f13)/detfa*doubleproduct
474  cnode(3,3) = (f11*f22-f12*f21)/detfa*doubleproduct
475 
476  s11n(i) = coeff1*(f11 - onethird*cnode(1,1))
477  s12n(i) = coeff1*(f12 - onethird*cnode(1,2))
478  s13n(i) = coeff1*(f13 - onethird*cnode(1,3))
479  s21n(i) = coeff1*(f21 - onethird*cnode(2,1))
480  s22n(i) = coeff1*(f22 - onethird*cnode(2,2))
481  s23n(i) = coeff1*(f23 - onethird*cnode(2,3))
482  s31n(i) = coeff1*(f31 - onethird*cnode(3,1))
483  s32n(i) = coeff1*(f32 - onethird*cnode(3,2))
484  s33n(i) = coeff1*(f33 - onethird*cnode(3,3))
485 
486  ENDDO
487 
488 
489 ! IF(myid.eq.0) print*,'1procs',myid,S11n(11),S22n(11),S33n(11),S12n(11),S23n(11),S13n(11)
490 ! IF(myid.EQ.1) print*,myid,S11n(11),S22n(11),S33n(11),S12n(11),S23n(11),S13n(11)
491 ! IF(myid.EQ.0) print*,myid,S11n(3),S22n(3),S33n(3),S12n(3),S23n(3),S13n(3)
492 
493  DO ielnum = 1, numel ! For each element
494 
495  n1 = nodes(1,ielnum)
496  n2 = nodes(2,ielnum)
497  n3 = nodes(3,ielnum)
498  n4 = nodes(4,ielnum)
499 
500  k3n1 = 3*n1
501  k3n2 = 3*n2
502  k3n3 = 3*n3
503  k3n4 = 3*n4
504 
505  k2n1 = k3n1 - 1
506  k2n2 = k3n2 - 1
507  k2n3 = k3n3 - 1
508  k2n4 = k3n4 - 1
509 
510  k1n1 = k3n1 - 2
511  k1n2 = k3n2 - 2
512  k1n3 = k3n3 - 2
513  k1n4 = k3n4 - 2
514 
515  ! k#n# dummy variables replaces:
516  u1 = disp(k1n1) ! (3*n1-2)
517  u2 = disp(k1n2) ! (3*n2-2)
518  u3 = disp(k1n3) ! (3*n3-2)
519  u4 = disp(k1n4) ! (3*n4-2)
520  v1 = disp(k2n1) ! (3*n1-1)
521  v2 = disp(k2n2) ! (3*n2-1)
522  v3 = disp(k2n3) ! (3*n3-1)
523  v4 = disp(k2n4) ! (3*n4-1)
524  w1 = disp(k3n1) ! (3*n1)
525  w2 = disp(k3n2) ! (3*n2)
526  w3 = disp(k3n3) ! (3*n3)
527  w4 = disp(k3n4) ! (3*n4)
528 
529  x1 = coor(1,n1)
530  x2 = coor(1,n2)
531  x3 = coor(1,n3)
532  x4 = coor(1,n4)
533  y1 = coor(2,n1)
534  y2 = coor(2,n2)
535  y3 = coor(2,n3)
536  y4 = coor(2,n4)
537  z1 = coor(3,n1)
538  z2 = coor(3,n2)
539  z3 = coor(3,n3)
540  z4 = coor(3,n4)
541 
542  x1 = x1 + u1
543  x2 = x2 + u2
544  x3 = x3 + u3
545  x4 = x4 + u4
546  y1 = y1 + v1
547  y2 = y2 + v2
548  y3 = y3 + v3
549  y4 = y4 + v4
550  z1 = z1 + w1
551  z2 = z2 + w2
552  z3 = z3 + w3
553  z4 = z4 + w4
554 
555  x12 = x1 - x2 ! not used in vol. calc
556  x13 = x1 - x3 ! not used in vol. calc
557  x14 = x1 - x4
558  x24 = x2 - x4
559  x34 = x3 - x4
560  y12 = y1 - y2 ! not used in vol. calc
561  y13 = y1 - y3 ! not used in vol. calc
562  y14 = y1 - y4
563  y24 = y2 - y4
564  y34 = y3 - y4
565  z12 = z1 - z2 ! not used in vol. calc
566  z13 = z1 - z3 ! not used in vol. calc
567  z14 = z1 - z4
568  z24 = z2 - z4
569  z34 = z3 - z4
570 
571  c11 = y24*z34 - z24*y34
572  c21 = -( x24*z34 - z24*x34 )
573  c31 = x24*y34 - y24*x34
574 
575  vx6def = -( x14*c11 + y14*c21 + z14*c31 )
576 
577  vx6invdef = 1.d0 / vx6def
578 
579 ! calculate the volume
580  voldef = vx6def/6.d0
581 
582  IF(voldef.LE.0.d0) THEN
583  WRITE(*,200) i
584  stop
585  ENDIF
586 
587  ! equation (14)
588 
589  s11e = 0.d0
590  s12e = 0.d0
591  s13e = 0.d0
592  s21e = 0.d0
593  s22e = 0.d0
594  s23e = 0.d0
595  s31e = 0.d0
596  s32e = 0.d0
597  s33e = 0.d0
598  pc = 0.d0
599 
600  DO k = 1, 4
601  s11e = s11e + alpha(k,ielnum)*s11n(nodes(k,ielnum))
602  s12e = s12e + alpha(k,ielnum)*s12n(nodes(k,ielnum))
603  s13e = s13e + alpha(k,ielnum)*s13n(nodes(k,ielnum))
604  s21e = s21e + alpha(k,ielnum)*s21n(nodes(k,ielnum))
605  s22e = s22e + alpha(k,ielnum)*s22n(nodes(k,ielnum))
606  s23e = s23e + alpha(k,ielnum)*s23n(nodes(k,ielnum))
607  s31e = s31e + alpha(k,ielnum)*s31n(nodes(k,ielnum))
608  s32e = s32e + alpha(k,ielnum)*s32n(nodes(k,ielnum))
609  s33e = s33e + alpha(k,ielnum)*s33n(nodes(k,ielnum))
610  pc = pc + alpha(k,ielnum)*pa(nodes(k,ielnum))
611  ENDDO
612 
613  sh1 = (y34*z24 - y24*z34) * vx6invdef
614  sh2 = (z34*x24 - z24*x34) * vx6invdef
615  sh3 = (x34*y24 - x24*y34) * vx6invdef
616  sh4 = (y13*z14 - y14*z13) * vx6invdef
617  sh5 = (z13*x14 - z14*x13) * vx6invdef
618  sh6 = (x13*y14 - x14*y13) * vx6invdef
619  sh7 = (y14*z12 - y12*z14) * vx6invdef
620  sh8 = (z14*x12 - z12*x14) * vx6invdef
621  sh9 = (x14*y12 - x12*y14) * vx6invdef
622  sh10 = (y12*z13 - y13*z12) * vx6invdef
623  sh11 = (z12*x13 - z13*x12) * vx6invdef
624  sh12 = (x12*y13 - x13*y12) * vx6invdef
625 
626 ! ASSEMBLE THE INTERNAL FORCE VECTOR
627 !
628 ! local node 1
629  rnet(k1n1) = rnet(k1n1) - voldef*pc*sh1
630  rnet(k2n1) = rnet(k2n1) - voldef*pc*sh2
631  rnet(k3n1) = rnet(k3n1) - voldef*pc*sh3
632 ! local node 2
633  rnet(k1n2) = rnet(k1n2) - voldef*pc*sh4
634  rnet(k2n2) = rnet(k2n2) - voldef*pc*sh5
635  rnet(k3n2) = rnet(k3n2) - voldef*pc*sh6
636 ! local node 3
637  rnet(k1n3) = rnet(k1n3) - voldef*pc*sh7
638  rnet(k2n3) = rnet(k2n3) - voldef*pc*sh8
639  rnet(k3n3) = rnet(k3n3) - voldef*pc*sh9
640 ! local node 4
641  rnet(k1n4) = rnet(k1n4) - voldef*pc*sh10
642  rnet(k2n4) = rnet(k2n4) - voldef*pc*sh11
643  rnet(k3n4) = rnet(k3n4) - voldef*pc*sh12
644 
645  x1 = coor(1,n1)
646  x2 = coor(1,n2)
647  x3 = coor(1,n3)
648  x4 = coor(1,n4)
649  y1 = coor(2,n1)
650  y2 = coor(2,n2)
651  y3 = coor(2,n3)
652  y4 = coor(2,n4)
653  z1 = coor(3,n1)
654  z2 = coor(3,n2)
655  z3 = coor(3,n3)
656  z4 = coor(3,n4)
657 
658  x12 = x1 - x2 ! not used in vol. calc
659  x13 = x1 - x3 ! not used in vol. calc
660  x14 = x1 - x4
661  x24 = x2 - x4
662  x34 = x3 - x4
663  y12 = y1 - y2 ! not used in vol. calc
664  y13 = y1 - y3 ! not used in vol. calc
665  y14 = y1 - y4
666  y24 = y2 - y4
667  y34 = y3 - y4
668  z12 = z1 - z2 ! not used in vol. calc
669  z13 = z1 - z3 ! not used in vol. calc
670  z14 = z1 - z4
671  z24 = z2 - z4
672  z34 = z3 - z4
673 
674  c11 = y24*z34 - z24*y34
675  c21 = -( x24*z34 - z24*x34 )
676  c31 = x24*y34 - y24*x34
677 
678  vx6 = -( x14*c11 + y14*c21 + z14*c31 )
679 
680  vx6inv = 1.d0/vx6
681 
682  volel0 = vx6/6.d0 ! undeformed volume of element (Ve_O)
683 
684 ! Compute the Shape functions
685 ! NOTE: Factored for a more equivalent/compact form then maple's
686 
687  b1 = (y34*z24 - y24*z34) * vx6inv
688  b2 = (z34*x24 - z24*x34) * vx6inv
689  b3 = (x34*y24 - x24*y34) * vx6inv
690  b4 = (y13*z14 - y14*z13) * vx6inv
691  b5 = (z13*x14 - z14*x13) * vx6inv
692  b6 = (x13*y14 - x14*y13) * vx6inv
693  b7 = (y14*z12 - y12*z14) * vx6inv
694  b8 = (z14*x12 - z12*x14) * vx6inv
695  b9 = (x14*y12 - x12*y14) * vx6inv
696  b10 = (y12*z13 - y13*z12) * vx6inv
697  b11 = (z12*x13 - z13*x12) * vx6inv
698  b12 = (x12*y13 - x13*y12) * vx6inv
699 ! local node 1
700  rnet(k1n1) = rnet(k1n1) - volel0* &
701  ( s11e*b1 + s12e*b2 + s13e*b3 )
702  rnet(k2n1) = rnet(k2n1) - volel0* &
703  ( s21e*b1 + s22e*b2 + s23e*b3 )
704  rnet(k3n1) = rnet(k3n1) - volel0* &
705  ( s31e*b1 + s32e*b2 + s33e*b3 )
706 ! local node 2
707  rnet(k1n2) = rnet(k1n2) - volel0* &
708  ( s11e*b4 + s12e*b5 + s13e*b6 )
709  rnet(k2n2) = rnet(k2n2) - volel0* &
710  ( s21e*b4 + s22e*b5 + s23e*b6 )
711  rnet(k3n2) = rnet(k3n2) - volel0* &
712  ( s31e*b4 + s32e*b5 + s33e*b6 )
713 ! local node 3
714  rnet(k1n3) = rnet(k1n3) - volel0* &
715  ( s11e*b7 + s12e*b8 + s13e*b9 )
716  rnet(k2n3) = rnet(k2n3) - volel0* &
717  ( s21e*b7 + s22e*b8 + s23e*b9 )
718  rnet(k3n3) = rnet(k3n3) - volel0* &
719  ( s31e*b7 + s32e*b8 + s33e*b9 )
720 ! local node 4
721  rnet(k1n4) = rnet(k1n4) - volel0* &
722  ( s11e*b10 + s12e*b11 + s13e*b12 )
723  rnet(k2n4) = rnet(k2n4) - volel0* &
724  ( s21e*b10 + s22e*b11 + s23e*b12 )
725  rnet(k3n4) = rnet(k3n4) - volel0* &
726  ( s31e*b10 + s32e*b11 + s33e*b12 )
727 
728  ENDDO
729 
730  RETURN
731 100 FORMAT(' Negative Jacobian for element: ',i10)
732 200 FORMAT(' Negative Jacobian for element (undef): ',i10)
733 
734 END SUBROUTINE v3d4n_neohookeanincompress
735 
736 
rocstar_rocfraccomm
Definition: ROCSTAR_RocFracComm.f90:53
k
j indices k indices k
Definition: Indexing.h:6
rocstar_rocfraccomm::rcv_buf
Definition: ROCSTAR_RocFracComm.f90:71
v3d4n_neohookeanincompress
subroutine v3d4n_neohookeanincompress(numnp, numel, coor, disp, nodes, Rnet, S11ev, S22ev, S33ev, S12ev, S23ev, S13ev, NumElNeigh, ElConn, alpha, Ahat, xmu, xkappa, myid, nprocs, TotNumNdComm, TotNumNeighProcs, NeighProcList, NumNdComm, neigh_lst, MPI_STATUS_SIZE, MPI_COMM_ROCFRAC, MPI_DOUBLE_PRECISION, ReqRcv, ReqSnd, StatRcv, StatSnd)
Definition: v3d4n_NeoHookeanInCompress.f90:53
i
blockLoc i
Definition: read.cpp:79
rocstar_rocfraccomm::send_buf
Definition: ROCSTAR_RocFracComm.f90:75
j
j indices j
Definition: Indexing.h:6
Color::alpha
unsigned char alpha() const
Definition: Color.h:75