Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
v3d4n.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(coor,nodes, NumElNeigh,ElConn,Rnet,alpha,disp,ci,&
54  numnp,numel,ahat,numat_vol,&
55  nprocs,totnumndcomm,totnumneighprocs,neighproclist,numndcomm,neigh_lst)
56 
57  USE rocstar_rocfraccomm
58 
59  IMPLICIT NONE
60  include 'mpif.h'
61 
62  INTEGER :: k1, k2, j1, k
63 
64 ! INTEGER :: MPI_COMM_WORLD ! Communicator of Rocfrac
65 ! INTEGER :: MPI_DOUBLE_PRECISION
66  INTEGER :: totnumneighprocs, nprocs,totnumndcomm
67  INTEGER, DIMENSION(1:TotNumNeighProcs) :: neighproclist
68  INTEGER, DIMENSION(1:TotNumNeighProcs) ::numndcomm
69  REAL*8, allocatable, dimension(:) :: buf
70 
71  TYPE(rcv_buf), ALLOCATABLE, DIMENSION(:) :: recvdatafrm
72  TYPE(send_buf), DIMENSION(1:TotNumNeighProcs) :: neigh_lst
73 
74  INTEGER :: ierr
75 !-- Non-block receive, Non-block send request arrays
76  INTEGER, ALLOCATABLE, DIMENSION(:) :: req_rcv_lst, req_snd_lst
77  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: stat_rcv_lst,stat_snd_lst
78 
79 !!$
80  INTEGER :: numnp, numel,i
81  REAL*8, DIMENSION(1:numnp) :: ahat
82  REAL*8, DIMENSION(1:numnp*3) :: disp, rnet
83  REAL*8, DIMENSION(1:3,1:numnp) :: coor
84  INTEGER, DIMENSION(1:numnp) :: numelneigh
85  INTEGER, DIMENSION(1:numnp,1:40) :: elconn ! fix 40 should not be hard coded
86  INTEGER,DIMENSION(1:4,1:numel) :: nodes
87  REAL*8, DIMENSION(1:4,1:numel) :: alpha
88  INTEGER :: numat_vol ! number of volumetric materials
89 !-- elastic stiffness consts
90  REAL*8, DIMENSION(1:9,1:numat_vol) :: ci
91 
92  INTEGER :: j
93  REAL*8 :: aaa
94 !-- coordinate holding variable
95  REAL*8 :: x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4
96  INTEGER :: k1n1,k1n2,k1n3,k1n4,k2n1,k2n2,k2n3,k2n4
97  INTEGER :: k3n1,k3n2,k3n3,k3n4
98 
99  REAL*8 :: ahatinv
100  INTEGER :: ielnum
101  REAL*8 :: sixinv
102 !-- x, y, and z displacements of nodes
103  REAL*8 :: u1,u2,u3,u4,v1,v2,v3,v4,w1,w2,w3,w4
104 
105 !-- Coordinate subtractions
106  REAL*8 :: x14, x24, x34, y14, y24, y34, z14, z24, z34
107  REAL*8 :: c11, c21, c31
108 ! -- 6*volume and the volume
109  REAL*8 :: vx6
110 !-- spacial derivatives
111  REAL*8 :: b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12
112 !-- strains
113  REAL*8 :: e11,e22,e33,e12,e23,e13
114  REAL*8 :: b13, b14, b15
115  REAL*8 :: sums11n,sums22n,sums33n,sums12n,sums23n,sums13n
116  INTEGER :: k3i,k2i,k1i
117  INTEGER :: ix
118  INTEGER :: n1, n2, n3, n4
119  REAL*8, DIMENSION(1:numnp) :: s11n, s22n, s33n, s12n, s23n, s13n
120 
121 ! CALL MPI_INIT(ierr)
122 ! CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
123 ! CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
124 
125  INTEGER :: myid
126  IF(totnumneighprocs.NE.0)THEN
127  ALLOCATE(req_rcv_lst(1:totnumneighprocs) )
128  ALLOCATE(req_snd_lst(1:totnumneighprocs) )
129 
130  ALLOCATE(stat_snd_lst(1:mpi_status_size,1:totnumneighprocs) )
131  ALLOCATE(stat_rcv_lst(1:mpi_status_size,1:totnumneighprocs) )
132  ENDIF
133 
134  sixinv = 1.d0/6.d0
135 
136  DO i = 1, numnp
137  ahatinv = 1.d0/ahat(i)
138 
139  sums11n = 0.d0
140  sums22n = 0.d0
141  sums33n = 0.d0
142  sums12n = 0.d0
143  sums23n = 0.d0
144  sums13n = 0.d0
145 
146  DO j = 1, numelneigh(i)
147 
148  ielnum = elconn(i,j)
149 
150  n1 = nodes(1,ielnum)
151  n2 = nodes(2,ielnum)
152  n3 = nodes(3,ielnum)
153  n4 = nodes(4,ielnum)
154 
155  IF(i.EQ.n1) THEN
156  ix = 1
157  ELSE IF(i.EQ.n2)THEN
158  ix = 2
159  ELSE IF(i.EQ.n3)THEN
160  ix = 3
161  ELSE IF(i.EQ.n4)THEN
162  ix = 4
163  ENDIF
164 
165  k3n1 = 3*n1
166  k3n2 = 3*n2
167  k3n3 = 3*n3
168  k3n4 = 3*n4
169 
170  k2n1 = k3n1 - 1
171  k2n2 = k3n2 - 1
172  k2n3 = k3n3 - 1
173  k2n4 = k3n4 - 1
174 
175  k1n1 = k3n1 - 2
176  k1n2 = k3n2 - 2
177  k1n3 = k3n3 - 2
178  k1n4 = k3n4 - 2
179  ! k#n# dummy variables replaces:
180  u1 = disp(k1n1) ! (3*n1-2)
181  u2 = disp(k1n2) ! (3*n2-2)
182  u3 = disp(k1n3) ! (3*n3-2)
183  u4 = disp(k1n4) ! (3*n4-2)
184  v1 = disp(k2n1) ! (3*n1-1)
185  v2 = disp(k2n2) ! (3*n2-1)
186  v3 = disp(k2n3) ! (3*n3-1)
187  v4 = disp(k2n4) ! (3*n4-1)
188  w1 = disp(k3n1) ! (3*n1)
189  w2 = disp(k3n2) ! (3*n2)
190  w3 = disp(k3n3) ! (3*n3)
191  w4 = disp(k3n4) ! (3*n4)
192 
193  x1 = coor(1,n1)
194  x2 = coor(1,n2)
195  x3 = coor(1,n3)
196  x4 = coor(1,n4)
197  y1 = coor(2,n1)
198  y2 = coor(2,n2)
199  y3 = coor(2,n3)
200  y4 = coor(2,n4)
201  z1 = coor(3,n1)
202  z2 = coor(3,n2)
203  z3 = coor(3,n3)
204  z4 = coor(3,n4)
205 
206 ! See the maple worksheet 'V3D4.mws' for the derivation of [B]
207 
208  b1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3) * sixinv
209  b2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3) * sixinv
210  b3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3) * sixinv
211  b4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3) * sixinv
212  b5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3) * sixinv
213  b6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3) * sixinv
214  b7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2) * sixinv
215  b8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2) * sixinv
216  b9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2) * sixinv
217  b10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2) * sixinv
218  b11 = -(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2) * sixinv
219  b12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2) * sixinv
220 
221 ! calculate the strain
222 ! [E] = [B]{d}
223 !
224  e11 = alpha(ix,ielnum)*(b1*u1 + b4*u2 + b7*u3 + b10*u4)
225  e22 = alpha(ix,ielnum)*(b2*v1 + b5*v2 + b8*v3 + b11*v4)
226  e33 = alpha(ix,ielnum)*(b3*w1 + b6*w2 + b9*w3 + b12*w4)
227  e12 = alpha(ix,ielnum)*(b2*u1 + b1*v1 + b5*u2 + b4*v2 + b8*u3 + b7*v3 + b11*u4 + b10*v4)
228  e23 = alpha(ix,ielnum)*(b3*v1 + b2*w1 + b6*v2 + b5*w2 + b9*v3 + b8*w3 + b12*v4 + b11*w4)
229  e13 = alpha(ix,ielnum)*(b3*u1 + b1*w1 + b6*u2 + b4*w2 + b9*u3 + b7*w3 + b12*u4 + b10*w4)
230 ! calculate the stress -1
231 ! [S] = [C] {E}
232 !
233  sums11n = sums11n + e11*ci(1,1) + e22*ci(2,1) + e33*ci(4,1)
234  sums22n = sums22n + e11*ci(2,1) + e22*ci(3,1) + e33*ci(5,1)
235  sums33n = sums33n + e11*ci(4,1) + e22*ci(5,1) + e33*ci(6,1)
236  sums12n = sums12n + e12*ci(7,1)
237  sums23n = sums23n + e23*ci(8,1)
238  sums13n = sums13n + e13*ci(9,1)
239 
240  ENDDO
241  s11n(i) = ahatinv*sums11n
242  s22n(i) = ahatinv*sums22n
243  s33n(i) = ahatinv*sums33n
244  s12n(i) = ahatinv*sums12n
245  s23n(i) = ahatinv*sums23n
246  s13n(i) = ahatinv*sums13n ! stress at node
247 
248 ! print*,myid,i,AhatInv*SumS23n
249 
250  ENDDO
251 
252  ALLOCATE(recvdatafrm(0:nprocs-1))
253 
254 !
255 !----- FORM THE BUFFER CONTAINING COMMUNICATED STRESS MATRIX NODAL VALUES
256 
257  ALLOCATE(buf(1:totnumndcomm*2))
258  k1 = 1
259  DO j1 = 1, totnumneighprocs
260  k = neighproclist(j1)
261  ALLOCATE(recvdatafrm(k)%rcvbuf(1:numndcomm(j1)*6))
262  DO j = 1, numndcomm(j1)
263  k2 = neigh_lst(j1)%NdID(j) !NdCommList(k)%NdId(j)
264 ! print*,myid,k2,S23n(k2)
265  buf(k1 ) = s11n(k2)
266  buf(k1+1) = s22n(k2)
267  buf(k1+2) = s33n(k2)
268  buf(k1+3) = s12n(k2)
269  buf(k1+4) = s23n(k2)
270  buf(k1+5) = s13n(k2)
271  k1 = k1 + 6
272  ENDDO
273  ENDDO
274 !
275 !-MPI- RECEIVE THE RECIPRICAL MASS MATRIX DIAGONAL FROM THE NEIGHBORS
276 !
277  DO j1 = 1, totnumneighprocs
278  k = neighproclist(j1)
279  CALL mpi_irecv(recvdatafrm(k)%rcvbuf(1),numndcomm(j1)*6, &
280  mpi_double_precision, k, 10, rocstar_communicator,req_rcv_lst(j1),ierr)
281  ENDDO
282 !
283 !-MPI- SEND THE RECIPRICAL MASS MATRIX DIAGONAL TO THE NEIGHBORS
284 !
285  k2 = 1
286  DO j1 = 1, totnumneighprocs
287  k = neighproclist(j1)
288 
289 ! IF(myid.EQ.0) print*,'sending',buf(1:)
290  CALL mpi_isend(buf(k2),numndcomm(j1)*6,&
291  mpi_double_precision,k,10,rocstar_communicator,req_snd_lst(j1),ierr)
292  k2 = k2 + numndcomm(j1)
293  ENDDO
294 !
295 !-MPI- WAIT FOR INTERNAL FORCE VECTOR COMMUNICATION TO COMPLETE
296 !
297  IF(totnumneighprocs.GT.0)THEN
298  CALL mpi_waitall(totnumneighprocs,req_rcv_lst,stat_rcv_lst,ierr)
299  CALL mpi_waitall(totnumneighprocs,req_snd_lst,stat_snd_lst,ierr)
300  ENDIF
301  DEALLOCATE(buf)
302 
303 !
304 !----- ADD NEIGHBOR'S CONTRIBUTION TO THE INTERNAL FORCE VECTOR
305 !
306  DO j1 = 1, totnumneighprocs
307  k = neighproclist(j1)
308  k1 = 1
309 ! IF(myid.EQ.1) print*,'***',RecvDataFrm(0)%rcvbuf(:)
310  DO j = 1, numndcomm(j1)
311  k2 = neigh_lst(j1)%NdID(j)
312  s11n(k2) = s11n(k2) + recvdatafrm(k)%rcvbuf(k1)
313  s22n(k2) = s22n(k2) + recvdatafrm(k)%rcvbuf(k1+1)
314  s33n(k2) = s33n(k2) + recvdatafrm(k)%rcvbuf(k1+2)
315  s12n(k2) = s12n(k2) + recvdatafrm(k)%rcvbuf(k1+3)
316  s23n(k2) = s23n(k2) + recvdatafrm(k)%rcvbuf(k1+4)
317  s13n(k2) = s13n(k2) + recvdatafrm(k)%rcvbuf(k1+5)
318  k1 = k1 + 6
319  ENDDO
320  ENDDO
321 
322  DEALLOCATE(recvdatafrm)
323 
324 ! IF(myid.eq.0) print*,'1procs',myid,S11n(11),S22n(11),S33n(11),S12n(11),S23n(11),S13n(11)
325 ! IF(myid.EQ.1) print*,myid,S11n(11),S22n(11),S33n(11),S12n(11),S23n(11),S13n(11)
326 ! IF(myid.EQ.0) print*,myid,S11n(3),S22n(3),S33n(3),S12n(3),S23n(3),S13n(3)
327 
328  DO i = 1, numnp
329 
330  k3i = 3*i
331  k2i = 3*i-1
332  k1i = 3*i-2
333 
334  DO j = 1, numelneigh(i)
335 
336  ielnum = elconn(i,j)
337 
338  n1 = nodes(1,ielnum)
339  n2 = nodes(2,ielnum)
340  n3 = nodes(3,ielnum)
341  n4 = nodes(4,ielnum)
342 
343  k3n1 = 3*n1
344  k3n2 = 3*n2
345  k3n3 = 3*n3
346  k3n4 = 3*n4
347 
348  k2n1 = k3n1 - 1
349  k2n2 = k3n2 - 1
350  k2n3 = k3n3 - 1
351  k2n4 = k3n4 - 1
352 
353  k1n1 = k3n1 - 2
354  k1n2 = k3n2 - 2
355  k1n3 = k3n3 - 2
356  k1n4 = k3n4 - 2
357 
358  x1 = coor(1,n1)
359  x2 = coor(1,n2)
360  x3 = coor(1,n3)
361  x4 = coor(1,n4)
362  y1 = coor(2,n1)
363  y2 = coor(2,n2)
364  y3 = coor(2,n3)
365  y4 = coor(2,n4)
366  z1 = coor(3,n1)
367  z2 = coor(3,n2)
368  z3 = coor(3,n3)
369  z4 = coor(3,n4)
370 
371 ! See the maple worksheet 'V3D4.mws' for the derivation of [B]
372 
373 !!$ B1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3) * SixInv
374 !!$ B2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3) * SixInv
375 !!$ B3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3) * SixInv
376 !!$ B4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3) * SixInv
377 !!$ B5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3) * SixInv
378 !!$ B6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3) * SixInv
379 !!$ B7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2) * SixInv
380 !!$ B8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2) * SixInv
381 !!$ B9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2) * SixInv
382 !!$ B10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2) * SixInv
383 !!$ B11 = -(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2) * SixInv
384 !!$ B12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2) * SixInv
385 
386  b1 = ( (y3-y4)*(z2-z4) - (y2-y4)*(z3-z4) ) * sixinv
387  b2 = ( (z3-z4)*(x2-x4) - (z2-z4)*(x3-x4) ) * sixinv
388  b3 = ( (x3-x4)*(y2-y4) - (x2-x4)*(y3-y4) ) * sixinv
389  b4 = ( (y1-y3)*(z1-z4) - (y1-y4)*(z1-z3) ) * sixinv
390  b5 = ( (z1-z3)*(x1-x4) - (z1-z4)*(x1-x3) ) * sixinv
391  b6 = ( (x1-x3)*(y1-y4) - (x1-x4)*(y1-y3) ) * sixinv
392  b7 = ( (y1-y4)*(z1-z2) - (y1-y2)*(z1-z4) ) * sixinv
393  b8 = ( (z1-z4)*(x1-x2) - (z1-z2)*(x1-x4) ) * sixinv
394  b9 = ( (x1-x4)*(y1-y2) - (x1-x2)*(y1-y4) ) * sixinv
395  b10 = ( (y1-y2)*(z1-z3) - (y1-y3)*(z1-z2) ) * sixinv
396  b11 = ( (z1-z2)*(x1-x3) - (z1-z3)*(x1-x2) ) * sixinv
397  b12 = ( (x1-x2)*(y1-y3) - (x1-x3)*(y1-y2) ) * sixinv
398 
399  IF(n1.EQ.i)THEN
400  b13 = b1
401  b14 = b2
402  b15 = b3
403  ELSE IF (n2.EQ.i)THEN
404  b13 = b4
405  b14 = b5
406  b15 = b6
407  ELSE IF(n3.EQ.i)THEN
408  b13 = b7
409  b14 = b8
410  b15 = b9
411  ELSE IF(n4.EQ.i)THEN
412  b13 = b10
413  b14 = b11
414  b15 = b12
415  ENDIF
416 
417 ! ASSEMBLE THE INTERNAL FORCE VECTOR
418 
419  rnet(k1i) = rnet(k1i) - &
420  (s11n(n1)*b13 + s12n(n1)*b14 + s13n(n1)*b15)*alpha(1,ielnum) - & ! local node 1
421  (s11n(n2)*b13 + s12n(n2)*b14 + s13n(n2)*b15)*alpha(2,ielnum) - & ! local node 2
422  (s11n(n3)*b13 + s12n(n3)*b14 + s13n(n3)*b15)*alpha(3,ielnum) - & ! local node 3
423  (s11n(n4)*b13 + s12n(n4)*b14 + s13n(n4)*b15)*alpha(4,ielnum) ! local node 4
424 
425  rnet(k2i) = rnet(k2i) - &
426  (s22n(n1)*b14 + s12n(n1)*b13 + s23n(n1)*b15)*alpha(1,ielnum) - &
427  (s22n(n2)*b14 + s12n(n2)*b13 + s23n(n2)*b15)*alpha(2,ielnum) - &
428  (s22n(n3)*b14 + s12n(n3)*b13 + s23n(n3)*b15)*alpha(3,ielnum) - &
429  (s22n(n4)*b14 + s12n(n4)*b13 + s23n(n4)*b15)*alpha(4,ielnum)
430 
431  rnet(k3i) = rnet(k3i) - &
432  (s33n(n1)*b15 + s23n(n1)*b14 + s13n(n1)*b13)*alpha(1,ielnum) - &
433  (s33n(n2)*b15 + s23n(n2)*b14 + s13n(n2)*b13)*alpha(2,ielnum) - &
434  (s33n(n3)*b15 + s23n(n3)*b14 + s13n(n3)*b13)*alpha(3,ielnum) - &
435  (s33n(n4)*b15 + s23n(n4)*b14 + s13n(n4)*b13)*alpha(4,ielnum)
436 
437 !!$ IF(i.EQ.11.and.Myid.EQ.0)print*,'1proc',(S22n(n1)*B14 + S12n(n1)*B13 + S23n(n1)*B15)*alpha(1,iElNum) - &
438 !!$ (S22n(n2)*B14 + S12n(n2)*B13 + S23n(n2)*B15)*alpha(2,iElNum) - &
439 !!$ (S22n(n3)*B14 + S12n(n3)*B13 + S23n(n3)*B15)*alpha(3,iElNum) - &
440 !!$ (S22n(n4)*B14 + S12n(n4)*B13 + S23n(n4)*B15)*alpha(4,iElNum)
441 !!$
442 !!$ IF(i.EQ.3.AND.myid.EQ.0)print*,(S22n(n1)*B14 + S12n(n1)*B13 + S23n(n1)*B15)*alpha(1,iElNum) - &
443 !!$ (S22n(n2)*B14 + S12n(n2)*B13 + S23n(n2)*B15)*alpha(2,iElNum) - &
444 !!$ (S22n(n3)*B14 + S12n(n3)*B13 + S23n(n3)*B15)*alpha(3,iElNum) - &
445 !!$ (S22n(n4)*B14 + S12n(n4)*B13 + S23n(n4)*B15)*alpha(4,iElNum)
446 !!$
447 !!$ IF(i.EQ.11.AND.myid.EQ.1)print*,(S22n(n1)*B14 + S12n(n1)*B13 + S23n(n1)*B15)*alpha(1,iElNum) - &
448 !!$ (S22n(n2)*B14 + S12n(n2)*B13 + S23n(n2)*B15)*alpha(2,iElNum) - &
449 !!$ (S22n(n3)*B14 + S12n(n3)*B13 + S23n(n3)*B15)*alpha(3,iElNum) - &
450 !!$ (S22n(n4)*B14 + S12n(n4)*B13 + S23n(n4)*B15)*alpha(4,iElNum)
451 !!$
452 
453 !!$ IF(myid.EQ.0.and.k3i.eq.9) print*,- &
454 !!$ (S33n(n1)*B15 + S23n(n1)*B14 + S13n(n1)*B13)*alpha(1,iElNum) - &
455 !!$ (S33n(n2)*B15 + S23n(n2)*B14 + S13n(n2)*B13)*alpha(2,iElNum) - &
456 !!$ (S33n(n3)*B15 + S23n(n3)*B14 + S13n(n3)*B13)*alpha(3,iElNum) - &
457 !!$ (S33n(n4)*B15 + S23n(n4)*B14 + S13n(n4)*B13)*alpha(4,iElNum)
458 !!$ IF(myid.EQ.0.and.k3i.eq.33) print*,'1procs',- &
459 !!$ (S33n(n1)*B15 + S23n(n1)*B14 + S13n(n1)*B13)*alpha(1,iElNum) - &
460 !!$ (S33n(n2)*B15 + S23n(n2)*B14 + S13n(n2)*B13)*alpha(2,iElNum) - &
461 !!$ (S33n(n3)*B15 + S23n(n3)*B14 + S13n(n3)*B13)*alpha(3,iElNum) - &
462 !!$ (S33n(n4)*B15 + S23n(n4)*B14 + S13n(n4)*B13)*alpha(4,iElNum)
463 !!$ IF(myid.EQ.1.and.k3i.eq.33) print*,'proc2',- &
464 !!$ (S33n(n1)*B15 + S23n(n1)*B14 + S13n(n1)*B13)*alpha(1,iElNum) - &
465 !!$ (S33n(n2)*B15 + S23n(n2)*B14 + S13n(n2)*B13)*alpha(2,iElNum) - &
466 !!$ (S33n(n3)*B15 + S23n(n3)*B14 + S13n(n3)*B13)*alpha(3,iElNum) - &
467 !!$ (S33n(n4)*B15 + S23n(n4)*B14 + S13n(n4)*B13)*alpha(4,iElNum)
468 
469  ENDDO
470  ENDDO
471 
472 ! IF(myid.eq.0) print*,'***1procs',myid,Rnet(3*11)
473 !
474 ! IF(myid.EQ.0) print*,'***ll',myid,Rnet(3*3)
475 ! IF(myid.EQ.1) print*,'***ll',myid,Rnet(3*11)
476 
477 
478 
479 !!$
480 !!$ SixInv = 1.d0/6.d0
481 !!$
482 !!$ DO i = 1, numnp
483 !!$ AhatInv = 1.d0/Ahat(i)
484 !!$
485 !!$ SumS11n = 0.d0
486 !!$ SumS22n = 0.d0
487 !!$ SumS33n = 0.d0
488 !!$ SumS12n = 0.d0
489 !!$ SumS23n = 0.d0
490 !!$ SumS13n = 0.d0
491 !!$
492 !!$ DO j = 1, NumElNeigh(i)
493 !!$
494 !!$ iElNum = ElConn(i,j)
495 !!$
496 !!$ n1 = nodes(1,iElNum)
497 !!$ n2 = nodes(2,iElNum)
498 !!$ n3 = nodes(3,iElNum)
499 !!$ n4 = nodes(4,iElNum)
500 !!$
501 !!$ IF(i.EQ.n1) THEN
502 !!$ ix = 1
503 !!$ ELSE IF(i.EQ.n2)THEN
504 !!$ ix = 2
505 !!$ ELSE IF(i.EQ.n3)THEN
506 !!$ ix = 3
507 !!$ ELSE IF(i.EQ.n4)THEN
508 !!$ ix = 4
509 !!$ ENDIF
510 !!$
511 !!$ k3n1 = 3*n1
512 !!$ k3n2 = 3*n2
513 !!$ k3n3 = 3*n3
514 !!$ k3n4 = 3*n4
515 !!$
516 !!$ k2n1 = k3n1 - 1
517 !!$ k2n2 = k3n2 - 1
518 !!$ k2n3 = k3n3 - 1
519 !!$ k2n4 = k3n4 - 1
520 !!$
521 !!$ k1n1 = k3n1 - 2
522 !!$ k1n2 = k3n2 - 2
523 !!$ k1n3 = k3n3 - 2
524 !!$ k1n4 = k3n4 - 2
525 !!$ ! k#n# dummy variables replaces:
526 !!$ u1 = disp(k1n1) ! (3*n1-2)
527 !!$ u2 = disp(k1n2) ! (3*n2-2)
528 !!$ u3 = disp(k1n3) ! (3*n3-2)
529 !!$ u4 = disp(k1n4) ! (3*n4-2)
530 !!$ v1 = disp(k2n1) ! (3*n1-1)
531 !!$ v2 = disp(k2n2) ! (3*n2-1)
532 !!$ v3 = disp(k2n3) ! (3*n3-1)
533 !!$ v4 = disp(k2n4) ! (3*n4-1)
534 !!$ w1 = disp(k3n1) ! (3*n1)
535 !!$ w2 = disp(k3n2) ! (3*n2)
536 !!$ w3 = disp(k3n3) ! (3*n3)
537 !!$ w4 = disp(k3n4) ! (3*n4)
538 !!$
539 !!$ x1 = coor(1,n1)
540 !!$ x2 = coor(1,n2)
541 !!$ x3 = coor(1,n3)
542 !!$ x4 = coor(1,n4)
543 !!$ y1 = coor(2,n1)
544 !!$ y2 = coor(2,n2)
545 !!$ y3 = coor(2,n3)
546 !!$ y4 = coor(2,n4)
547 !!$ z1 = coor(3,n1)
548 !!$ z2 = coor(3,n2)
549 !!$ z3 = coor(3,n3)
550 !!$ z4 = coor(3,n4)
551 !!$
552 !!$! See the maple worksheet 'V3D4.mws' for the derivation of [B]
553 !!$
554 !!$ B1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3) * SixInv
555 !!$ B2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3) * SixInv
556 !!$ B3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3) * SixInv
557 !!$ B4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3) * SixInv
558 !!$ B5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3) * SixInv
559 !!$ B6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3) * SixInv
560 !!$ B7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2) * SixInv
561 !!$ B8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2) * SixInv
562 !!$ B9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2) * SixInv
563 !!$ B10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2) * SixInv
564 !!$ B11 = -(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2) * SixInv
565 !!$ B12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2) * SixInv
566 !!$
567 !!$! calculate the strain
568 !!$! [E] = [B]{d}
569 !!$!
570 !!$ E11 = alpha(ix,iElNum)*(B1*u1 + B4*u2 + B7*u3 + B10*u4)
571 !!$ E22 = alpha(ix,iElNum)*(B2*v1 + B5*v2 + B8*v3 + B11*v4)
572 !!$ E33 = alpha(ix,iElNum)*(B3*w1 + B6*w2 + B9*w3 + B12*w4)
573 !!$ E12 = alpha(ix,iElNum)*(B2*u1 + B1*v1 + B5*u2 + B4*v2 + B8*u3 + B7*v3 + B11*u4 + B10*v4)
574 !!$ E23 = alpha(ix,iElNum)*(B3*v1 + B2*w1 + B6*v2 + B5*w2 + B9*v3 + B8*w3 + B12*v4 + B11*w4)
575 !!$ E13 = alpha(ix,iElNum)*(B3*u1 + B1*w1 + B6*u2 + B4*w2 + B9*u3 + B7*w3 + B12*u4 + B10*w4)
576 !!$! calculate the stress -1
577 !!$! [S] = [C] {E}
578 !!$!
579 !!$ SumS11n = SumS11n + E11*ci(1,1) + E22*ci(2,1) + E33*ci(4,1)
580 !!$ SumS22n = SumS22n + E11*ci(2,1) + E22*ci(3,1) + E33*ci(5,1)
581 !!$ SumS33n = SumS33n + E11*ci(4,1) + E22*ci(5,1) + E33*ci(6,1)
582 !!$ SumS12n = SumS12n + E12*ci(7,1)
583 !!$ SumS23n = SumS23n + E23*ci(8,1)
584 !!$ SumS13n = SumS13n + E13*ci(9,1)
585 !!$
586 !!$ ENDDO
587 !!$ S11n(i) = AhatInv*SumS11n
588 !!$ S22n(i) = AhatInv*SumS22n
589 !!$ S33n(i) = AhatInv*SumS33n
590 !!$ S12n(i) = AhatInv*SumS12n
591 !!$ S23n(i) = AhatInv*SumS23n
592 !!$ S13n(i) = AhatInv*SumS13n ! stress at node
593 !!$ ENDDO
594 !!$
595 !!$ DO j = 1, NumEl
596 !!$
597 !!$ iElNum = j
598 !!$
599 !!$ n1 = nodes(1,iElNum)
600 !!$ n2 = nodes(2,iElNum)
601 !!$ n3 = nodes(3,iElNum)
602 !!$ n4 = nodes(4,iElNum)
603 !!$
604 !!$ k3n1 = 3*n1
605 !!$ k3n2 = 3*n2
606 !!$ k3n3 = 3*n3
607 !!$ k3n4 = 3*n4
608 !!$
609 !!$ k2n1 = k3n1 - 1
610 !!$ k2n2 = k3n2 - 1
611 !!$ k2n3 = k3n3 - 1
612 !!$ k2n4 = k3n4 - 1
613 !!$
614 !!$ k1n1 = k3n1 - 2
615 !!$ k1n2 = k3n2 - 2
616 !!$ k1n3 = k3n3 - 2
617 !!$ k1n4 = k3n4 - 2
618 !!$
619 !!$ x1 = coor(1,n1)
620 !!$ x2 = coor(1,n2)
621 !!$ x3 = coor(1,n3)
622 !!$ x4 = coor(1,n4)
623 !!$ y1 = coor(2,n1)
624 !!$ y2 = coor(2,n2)
625 !!$ y3 = coor(2,n3)
626 !!$ y4 = coor(2,n4)
627 !!$ z1 = coor(3,n1)
628 !!$ z2 = coor(3,n2)
629 !!$ z3 = coor(3,n3)
630 !!$ z4 = coor(3,n4)
631 !!$
632 !!$! See the maple worksheet 'V3D4.mws' for the derivation of [B]
633 !!$
634 !!$ B1 = -(y3*z4 - y4*z3 - y2*z4 + y2*z3 + z2*y4 - z2*y3) * SixInv
635 !!$ B2 = (x3*z4 - x4*z3 - x2*z4 + x2*z3 + z2*x4 - z2*x3) * SixInv
636 !!$ B3 = -(x3*y4 - x4*y3 - x2*y4 + x2*y3 + y2*x4 - y2*x3) * SixInv
637 !!$ B4 = (y3*z4 - y4*z3 - y1*z4 + y1*z3 + z1*y4 - z1*y3) * SixInv
638 !!$ B5 = -(x3*z4 - x4*z3 - x1*z4 + x1*z3 + z1*x4 - z1*x3) * SixInv
639 !!$ B6 = (x3*y4 - x4*y3 - x1*y4 + x1*y3 + y1*x4 - y1*x3) * SixInv
640 !!$ B7 = -(y2*z4 - z2*y4 - y1*z4 + y1*z2 + z1*y4 - z1*y2) * SixInv
641 !!$ B8 = (x2*z4 - z2*x4 - x1*z4 + x1*z2 + z1*x4 - z1*x2) * SixInv
642 !!$ B9 = -(x2*y4 - y2*x4 - x1*y4 + x1*y2 + y1*x4 - y1*x2) * SixInv
643 !!$ B10 = (y2*z3 - z2*y3 - y1*z3 + y1*z2 + z1*y3 - z1*y2) * SixInv
644 !!$ B11 = -(x2*z3 - z2*x3 - x1*z3 + x1*z2 + z1*x3 - z1*x2) * SixInv
645 !!$ B12 = (x2*y3 - y2*x3 - x1*y3 + x1*y2 + y1*x3 - y1*x2) * SixInv
646 !!$
647 !!$! local node 1
648 !!$ Rnet(k1n1) = Rnet(k1n1) - (S11n(n1)*B1 + S12n(n1)*B2 + S13n(n1)*B3)*alpha(1,iElNum)
649 !!$ Rnet(k2n1) = Rnet(k2n1) - (S22n(n1)*B2 + S12n(n1)*B1 + S23n(n1)*B3)*alpha(1,iElNum)
650 !!$ Rnet(k3n1) = Rnet(k3n1) - (S33n(n1)*B3 + S23n(n1)*B2 + S13n(n1)*B1)*alpha(1,iElNum)
651 !!$! local node 2
652 !!$ Rnet(k1n2) = Rnet(k1n2) - (S11n(n2)*B4 + S12n(n2)*B5 + S13n(n2)*B6)*alpha(2,iElNum)
653 !!$ Rnet(k2n2) = Rnet(k2n2) - (S22n(n2)*B5 + S12n(n2)*B4 + S23n(n2)*B6)*alpha(2,iElNum)
654 !!$ Rnet(k3n2) = Rnet(k3n2) - (S33n(n2)*B6 + S23n(n2)*B5 + S13n(n2)*B4)*alpha(2,iElNum)
655 !!$! local node 3
656 !!$ Rnet(k1n3) = Rnet(k1n3) - (S11n(n3)*B7 + S12n(n3)*B8 + S13n(n3)*B9 )*alpha(3,iElNum)
657 !!$ Rnet(k2n3) = Rnet(k2n3) - (S22n(n3)*B8 + S12n(n3)*B7 + S23n(n3)*B9 )*alpha(3,iElNum)
658 !!$ Rnet(k3n3) = Rnet(k3n3) - (S33n(n3)*B9 + S23n(n3)*B8 + S13n(n3)*B7 )*alpha(3,iElNum)
659 !!$! local node 4
660 !!$ Rnet(k1n4) = Rnet(k1n4) - (S11n(n4)*B10 + S12n(n4)*B11 + S13n(n4)*B12 )*alpha(4,iElNum)
661 !!$ Rnet(k2n4) = Rnet(k2n4) - (S22n(n4)*B11 + S12n(n4)*B10 + S23n(n4)*B12 )*alpha(4,iElNum)
662 !!$ Rnet(k3n4) = Rnet(k3n4) - (S33n(n4)*B12 + S23n(n4)*B11 + S13n(n4)*B10 )*alpha(4,iElNum)
663 !!$
664 
665 !!$! ASSEMBLE THE INTERNAL FORCE VECTOR
666 !!$
667 !!$ Rnet(k1i) = Rnet(k1i) - &
668 !!$ (S11n(n1)*B13 + S12n(n1)*B14 + S13n(n1)*B15)*alpha(1,iElNum) - & ! local node 1
669 !!$ (S11n(n2)*B13 + S12n(n2)*B14 + S13n(n2)*B15)*alpha(2,iElNum) - & ! local node 2
670 !!$ (S11n(n3)*B13 + S12n(n3)*B14 + S13n(n3)*B15)*alpha(3,iElNum) - & ! local node 3
671 !!$ (S11n(n4)*B13 + S12n(n4)*B14 + S13n(n4)*B15)*alpha(4,iElNum) ! local node 4
672 !!$
673 !!$ Rnet(k2i) = Rnet(k2i) - &
674 !!$ (S22n(n1)*B14 + S12n(n1)*B13 + S23n(n1)*B15)*alpha(1,iElNum) - &
675 !!$ (S22n(n2)*B14 + S12n(n2)*B13 + S23n(n2)*B15)*alpha(2,iElNum) - &
676 !!$ (S22n(n3)*B14 + S12n(n3)*B13 + S23n(n3)*B15)*alpha(3,iElNum) - &
677 !!$ (S22n(n4)*B14 + S12n(n4)*B13 + S23n(n4)*B15)*alpha(4,iElNum)
678 !!$
679 !!$ Rnet(k3i) = Rnet(k3i) - &
680 !!$ (S33n(n1)*B15 + S23n(n1)*B14 + S13n(n1)*B13)*alpha(1,iElNum) - &
681 !!$ (S33n(n2)*B15 + S23n(n2)*B14 + S13n(n2)*B13)*alpha(2,iElNum) - &
682 !!$ (S33n(n3)*B15 + S23n(n3)*B14 + S13n(n3)*B13)*alpha(3,iElNum) - &
683 !!$ (S33n(n4)*B15 + S23n(n4)*B14 + S13n(n4)*B13)*alpha(4,iElNum)
684 !!$
685 
686 !!$ ENDDO
687 !!$ENDDO
688 
689  END SUBROUTINE v3d4n
690 
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
i
blockLoc i
Definition: read.cpp:79
rocstar_rocfraccomm::send_buf
Definition: ROCSTAR_RocFracComm.f90:75
j
j indices j
Definition: Indexing.h:6
v3d4n
subroutine v3d4n(coor, nodes, NumElNeigh, ElConn, Rnet, alpha, disp, ci, numnp, numel, Ahat, numat_vol, nprocs, TotNumNdComm, TotNumNeighProcs, NeighProcList, NumNdComm, neigh_lst)
Definition: v3d4n.f90:53
Color::alpha
unsigned char alpha() const
Definition: Color.h:75