Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RocFracMain.f90
Go to the documentation of this file.
1 ! * *******************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Version: Sandia Evaluation Distribution *
4 ! * Licensed To: Sandia National Laboratories *
5 ! * License Type: Evaluation *
6 ! * License Expiration: March 13, 2013 *
7 !*********************************************************************
8 ! * *******************************************************************
9 ! * Rocstar Simulation Suite *
10 ! * Copyright@2012, IllinoisRocstar LLC. All rights reserved. *
11 !* *
12 ! * The Rocstar Simulation Suite is the property of IllinoisRocstar *
13 ! * LLC. No use or distribution of this version of the Rocstar *
14 ! * Simulation Suite beyond the license provided through separate *
15 ! * contract is permitted. *
16 !* *
17 ! * IllinoisRocstar LLC *
18 ! * Champaign, IL *
19 ! * www.illinoisrocstar.com *
20 ! * sales@illinoisrocstar.com *
21 ! *********************************************************************
22 ! * *******************************************************************
23 ! * Initial open source Rocstar software developed by *
24 !* Center for Simulation of Advanced Rockets *
25 ! * University of Illinois at Urbana-Champaign *
26 ! * Urbana, IL *
27 !* www.csar.uiuc.edu *
28 !* *
29 ! * Copyright@2008, University of Illinois. All rights reserved. *
30 !* *
31 !* @ Redistributions of source code must retain the above copyright *
32 !* notice, this list of conditions and the following disclaimers. *
33 !* *
34 !* @ Redistributions in binary form must reproduce the above *
35 !* copyright notice, this list of conditions and the following *
36 !* disclaimers in the documentation and/or other materials *
37 !* provided with the distribution. *
38 !* *
39 !* @ Neither the names of the Center for Simulation of Advanced *
40 !* Rockets, the University of Illinois, nor the names of its *
41 !* contributors may be used to endorse or promote products derived *
42 !* from this Software without specific prior written permission. *
43 ! *********************************************************************
44 ! * *******************************************************************
45 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
46 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
47 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
48 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
49 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
50 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
51 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
52 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
53 !*********************************************************************
54 MODULE rocfracmain
55 
56 !
57 !!****h* Rocfrac/Source/RocFracMain
58 !!
59 !! NAME
60 !! RocFracMain.f90
61 !!
62 !! FUNCTION
63 !! 3D Dynamic Explicit Code with ALE formulation for regressing
64 !! boundaries Finite Element Analysis Code with additional fracture
65 !! simulation using cohesive elements.
66 !!
67 !! USAGE
68 !! Finite Element code to solve the 3-Dimensional TRANSIENT structural problem
69 !!
70 !! USES
71 !! RocFracSubInterface, UpdateStructuralSoln, feminp, VolRatio, vol_elem_mat,
72 !! RocFracInterfaceInitial, RocFracInterfaceBuff, UpdateMassMatrix, V3D4_volume,
73 !! max_dt_solid, UpdateRbar,v3d4_ale,V3D10_ALE,FluidPressLoad, TractPressLoad,
74 !! UpdateStructural,principal_stress,bc_enforce
75 !!
76 !! Global variables stored in modules : ROCSTAR_RocFrac,ROCSTAR_RocFracComm,ROCSTAR_RocFracInterp
77 !!
78 !! COPYRIGHT
79 !!
80 !! University of Illinois, Urbana-Champaign, (C) 2003
81 !!
82 !! AUTHOR
83 !!
84 !! principal : M.S. Breitenfeld, P.H. Geubelle
85 !! - email : brtnfld@uiuc.edu, geubelle@uiuc.edu
86 !!
87 !! contributing : Changyu Huang, Amit Acharya
88 !!
89 !!
90 !! CREATION DATE
91 !! 2001
92 !!
93 !!
94 !!*****
95 
96 
97  USE rocstar_rocfrac
98  USE rocstar_rocfraccomm
99  USE rocstar_rocfracinterp
100  USE rocfracsubinterface
101  USE updatestructuralsoln
102 
103  IMPLICIT NONE
104 
105  PUBLIC :: rocfracfinalize, rocfracsoln, rocfracinitialize
106 
107 CONTAINS
108 !
109 ! ------------------------------------------------------------------ RocFracInitialize
110 
111  SUBROUTINE rocfracinitialize( glb, InitialTime, MPI_COMM_ROCSTAR, MAN_init, surfIn, volIn, obtain_attr)
112 
113  IMPLICIT NONE
114  include 'mpif.h'
115  include 'roccomf90.h'
116 
117  TYPE(rocfrac_global),POINTER :: glb
118  REAL*8, INTENT(IN) :: initialtime
119  INTEGER, INTENT(IN) :: mpi_comm_rocstar
120  INTEGER, INTENT(IN) :: man_init, obtain_attr
121  CHARACTER(*), INTENT(IN) :: surfin, volin
122 
123  INTEGER :: i,j,j1,jj,k,k1,k2,idum,iaux,iaux1
124  REAL*8 :: aux1,aux
125  INTEGER :: ndbcflag
126 
127  INTEGER :: myid, numprocs, ierr
128 
129  INTEGER :: id2d
130  INTEGER :: iunit2
131  INTEGER :: kk
132  INTEGER, DIMENSION(3) :: ndsurf
133  REAL*8, ALLOCATABLE, DIMENSION(:) :: buf
134  INTEGER :: idpacket
135  REAL*8 :: mag
136 
137  REAL*8 :: ioversion
138 
139  INTEGER :: ios
140 
141  CHARACTER(len=2) :: chr2
142 
143  CHARACTER*120 :: frachdffname, meshfile
144 
145  character(LEN=1), POINTER, DIMENSION(:) :: names
146  character(LEN=3) :: ichr03
147  character(LEN=4) :: ichr04
148 
149  integer :: numeltypes2d
150 
151  character(LEN=4) :: chreltype
152  integer :: endpt, startpt, chrlngth
153 
154 
155  INTEGER :: numparcomm
156  INTEGER, pointer, dimension(:) :: parcomm
157 
158  INTEGER, pointer :: arraytmp
159  INTEGER, pointer, dimension(:) :: arraytmp1
160  integer :: iptr
161  integer :: icnt1, icnt2, icnt3, icnt4
162 
163  integer :: iprocs
164 
165  logical :: mesherror, mesherrorall, restartall
166 
167  glb%MPI_COMM_ROCFRAC = mpi_comm_rocstar
168  rocstar_communicator = mpi_comm_rocstar
169 
170  CALL mpi_comm_rank(glb%MPI_COMM_ROCFRAC,myid,ierr)
171  CALL mpi_comm_size(glb%MPI_COMM_ROCFRAC,numprocs,ierr)
172  WRITE(glb%MyIdChr,'(i4.4)') myid
173 
174 
175 
176 !
177 ! -- Read the control deck file
178 
179 ! IF(MyId.EQ.0) PRINT*,'ROCFRAC: Reading input deck'
180 
181  glb%io_input = 10
182  glb%io_sum = 11
183 
184  CALL feminp(glb, myid)
185 
186 
187  glb%InterfaceSFNumNodes = 0
188  glb%InterfaceSFnbNumNodes = 0
189 
190 
191  CALL mpi_barrier(glb%MPI_COMM_ROCFRAC,ierr)
192  IF(myid.eq.0 .AND. glb%debug_state) THEN
193  WRITE(6,'(A)') 'Rocfrac: ....Done Reading input deck.'
194  ENDIF
195 
196 ! -- Read in roc_face data structure for extracted 2D mesh.
197 ! mapnode maps the surface nodes to the volume mesh's nodes
198 
199  IF(glb%iElType.EQ.4)THEN
200  glb%iStrGss = 1
201  glb%LwrBnd = 1
202  glb%UppBnd = 3
203  ELSE IF(glb%iElType.EQ.10)THEN
204  glb%iStrGss = 4
205  glb%LwrBnd = 4
206  glb%UppBnd = 6
207  ELSE IF(glb%iElType.EQ.8)THEN
208  glb%iStrGss = 1
209  glb%LwrBnd = 1
210  glb%UppBnd = 4
211  ENDIF
212 
213 ! glb%InterfaceSVbar = 0.d0
214 
215  IF(myid.EQ.0 .AND. glb%debug_state) THEN
216  WRITE(6,'(A)')'Rocfrac: Reading Mesh'
217  ENDIF
218 
219 ! - Start
220 ! - Process the 3D mesh
221 
222 ! Number of Nodes
223  CALL com_new_window( volwin)
224 
225  CALL com_get_size( volin//".nc", myid+1, glb%NumNP)
226  CALL com_set_size( volwin//'.nc', myid+1, glb%NumNP )
227 
228 ! ALLOCATE(glb%coor(1:3,1:glb%NumNP))
229  ALLOCATE(glb%MeshCoor(1:3,1:glb%NumNP))
230  ALLOCATE(glb%xmass(1:glb%NumNP))
231  IF(glb%HeatTransSoln) ALLOCATE(glb%CapctInv(1:glb%NumNP))
232 
233 !!$
234 !!$ ! Nodal Coordinates
235 !!$ DO i=1,glb%NumNP
236 !!$ ! NodeID, X, Y, Z ! , DummyFlag1
237 !!$ READ(14,*) j,glb%coor(1,i),glb%coor(2,i),glb%coor(3,i) ! ,iaux
238 !!$
239 !!$ glb%Coor(1,j) = glb%Coor(1,j)+ mag
240 !!$ glb%Coor(3,j) = glb%Coor(3,j)+ mag
241 !!$ glb%MeshCoor(1:3,i) = glb%coor(1:3,i)
242 !!$ ENDDO
243 !!$
244 
245  CALL com_get_size( volin//".bcnode", myid+1,glb%NumNdsBCcrypt)
246 
247  ALLOCATE(glb%BCFlagCrypt(1:2,1:glb%NumNdsBCcrypt)) ! should this be half NumNdsBCdrypt
248 
249  ALLOCATE(glb%BCValueGlb(1:glb%NumNdsBCcrypt*6))
250 
251 
252 
253 !!$ ! _________________________________________
254 !!$ ! ___ Read Structural Nodal Boundary Flags
255 !!$ ! _________________________________________
256 !!$
257 !!$ READ(14,*) glb%NumNdsBC ! , iaux
258 !!$ ALLOCATE(glb%BCFlag(1:4,1:glb%NumNdsBC),glb%BCvalue(1:3,1:glb%NumNdsBC))
259 !!$ DO i = 1, glb%NumNdsBC
260 !!$ READ(14,*) glb%BCFlag(1,i),NdBCflag !, iaux
261 !!$ glb%BCFlag(2,i) = glb%bcCond(NdBCflag)%BCtypeX
262 !!$ glb%BCFlag(3,i) = glb%bcCond(NdBCflag)%BCtypeY
263 !!$ glb%BCFlag(4,i) = glb%bcCond(NdBCflag)%BCtypeZ
264 !!$ glb%BCvalue(1,i) = glb%bcCond(NdBCflag)%BCvalueX
265 !!$ glb%BCvalue(2,i) = glb%bcCond(NdBCflag)%BCvalueY
266 !!$ glb%BCvalue(3,i) = glb%bcCond(NdBCflag)%BCvalueZ
267 !!$ ENDDO
268 !!$
269 !!$! DEALLOCATE(glb%bcCond)
270 !!$
271 !!$ CASE(4)
272 !!$
273 !!$ ! Packet 04
274 !!$ ! __________________________________________
275 !!$ ! ___ Read Nodal Mesh Motion Boundary Flags
276 !!$ ! __________________________________________
277 !!$
278 !!$
279 !!$ READ(14,*) glb%NumNdsBCmm !, iaux
280 !!$ IF(glb%NumNdsBCmm.NE.0) ALLOCATE(glb%BCFlagmm(1:4,1:glb%NumNdsBCmm),glb%BCvaluemm(1:3,1:glb%NumNdsBCmm))
281 !!$
282 !!$ IF(glb%ALEenabled)THEN
283 !!$
284 !!$ DO i=1,glb%NumNdsBCmm
285 !!$ READ(14,*) glb%BCFlagmm(1,i) ,NdBCflag !, iaux
286 !!$ glb%BCFlagmm(2,i) = glb%bcCondmm(NdBCflag)%BCtypeX
287 !!$ glb%BCFlagmm(3,i) = glb%bcCondmm(NdBCflag)%BCtypeY
288 !!$ glb%BCFlagmm(4,i) = glb%bcCondmm(NdBCflag)%BCtypeZ
289 !!$ glb%BCvaluemm(1,i) = glb%bcCondmm(NdBCflag)%BCvalueX
290 !!$ glb%BCvaluemm(2,i) = glb%bcCondmm(NdBCflag)%BCvalueY
291 !!$ glb%BCvaluemm(3,i) = glb%bcCondmm(NdBCflag)%BCvalueZ
292 !!$ ENDDO
293 !!$ ELSE
294 !!$
295 !!$ DO i=1,glb%NumNdsBCmm
296 !!$ READ(14,'()')
297 !!$ ENDDO
298 !!$
299 !!$ ENDIF
300 !!$
301 !!$
302 !!$! IF(glb%NumNdsBCmm.NE.0) DEALLOCATE(glb%bcCondmm)
303 !!$
304 
305 !!$ ! __________________________________________
306 !!$ ! ___ Read Element Data
307 !!$ ! __________________________________________
308 
309  CALL com_get_connectivities(volin,myid+1,numeltypes2d,names)
310 
311  startpt = 1
312  DO i = 1, numeltypes2d
313  ! Search for the next attribute name
314  endpt = startpt
315  chrlngth = 0
316  DO WHILE (endpt .LE. ubound(names,1))
317  IF (names(endpt) .NE. ' ') THEN
318  chrlngth = chrlngth + 1
319  chreltype(chrlngth:chrlngth) = names(endpt)
320  endpt = endpt + 1
321  ELSE
322  EXIT
323  END IF
324  END DO
325 
326  startpt = endpt + 1
327 
328  mesherror = .false.
329  mesherrorall = .false.
330 
331  IF(chreltype(1:chrlngth).EQ.':T10')THEN
332  CALL com_get_size( volin//".:T10", myid+1, glb%NumElVol)
333  CALL com_set_size( volwin//'.:T10', myid+1, glb%NumElVol)
334  ELSE IF(chreltype(1:chrlngth).EQ.':T4')THEN
335  CALL com_get_size( volin//".:T4", myid+1, glb%NumElVol)
336  CALL com_set_size( volwin//'.:T4', myid+1, glb%NumElVol)
337  ELSE IF(chreltype(1:chrlngth).EQ.':H8')THEN
338  CALL com_get_size( volin//".:H8", myid+1, glb%NumElVol)
339  CALL com_set_size( volwin//'.:H8', myid+1, glb%NumElVol)
340  ELSE
341  mesherror = .true.
342  ENDIF
343 
344  CALL mpi_reduce(mesherror, mesherrorall, 1, mpi_logical, &
345  mpi_lor,0,glb%MPI_COMM_ROCFRAC,ierr)
346 
347  IF( myid.eq.0 .and. mesherrorall .eqv. .true.) THEN
348  WRITE(0,'(A,A)')'Rocfrac: Error: Volume mesh type', &
349  ' element not supported'
350  WRITE(0,'(A,A)')'Rocfrac: Read in Element Type :: ',&
351  chreltype(1:chrlngth)
352  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
353  END IF
354 
355  END DO
356 
357  CALL com_free_buffer(names)
358 
359 !!$ CALL COM_call_function( obtain_attr, 2, &
360 !!$ COM_get_attribute_handle_const( VolIn//".NumElPartBndry"), &
361 !!$ COM_get_attribute_handle( VolIn//".NumElPartBndry"))
362 ! stop
363  CALL com_get_array_const(volin//".NumElPartBndry",myid+1,arraytmp)
364 
365  glb%NumElPartBndry = arraytmp
366 
367 
368 !!$ CALL COM_call_function( obtain_attr, 2, &
369 !!$ COM_get_attribute_handle_const( VolIn//".NumElVolMat"), &
370 !!$ COM_get_attribute_handle( VolIn//".NumElVolMat"))
371  CALL com_get_array_const(volin//".NumElVolMat",myid+1,arraytmp1)
372 
373  ALLOCATE(glb%NumElVolMat(1:glb%NumMatVol))
374 
375  DO i = 1, glb%NumMatVol
376  glb%NumElVolMat(i) = arraytmp1(i)
377  ENDDO
378  CALL com_get_array_const(volin//".NumElPartBndryMat",myid+1,arraytmp1)
379 
380  ALLOCATE(glb%NumElPartBndryMat(1:glb%NumMatVol))
381  DO i = 1, glb%NumMatVol
382  glb%NumElPartBndryMat(i) = arraytmp1(i)
383  ENDDO
384 
385  IF(.NOT.(glb%DebondPart).AND..NOT.(glb%DebondPart_Matous)) THEN
386  ALLOCATE(glb%ci(1:9,1:glb%NumMatVol))
387 
388  ALLOCATE(glb%cj(1:9,1:glb%NumMatVol))
389  ENDIF
390 
391  i = 0
392 
393  DO kk = 1, glb%NumMatVol
394  IF(kk.EQ.1)THEN
395  ALLOCATE(glb%MatIdVol(1:glb%NumElVol))
396  ALLOCATE(glb%ElConnVol(1:glb%iElType,1:glb%NumElVol))
397  ENDIF
398 !
399 ! -- Read element data.
400 
401  IF(glb%NdBasedEl)THEN
402  IF(kk.EQ.1)THEN
403 
404  ALLOCATE(glb%NumElNeigh(1:glb%NumNP))
405  ALLOCATE(glb%ElConnNd(1:glb%NumNP,1:40))
406  ALLOCATE(glb%AlphaR(1:4,1:glb%NumElVol))
407  ALLOCATE(glb%VolUndfmd(1:glb%NumNP))
408 
409  glb%NumElNeigh(:) = 0
410  glb%ElConnNd(:,:) = 0 ! fix not work for more then one material ??
411  ENDIF
412  ENDIF
413  DO jj = 1,glb%NumElVolMat(kk)
414  i = i + 1
415 !!$ IF ( glb%iElType==4) THEN
416 !!$ READ(14,*) glb%MatIdVol(i),glb%ElConnVol(1,i), &
417 !!$ glb%ElConnVol(2,i), glb%ElConnVol(3,i), &
418 !!$ glb%ElConnVol(4,i) !, iaux, iaux
419 !!$ ELSE IF( glb%iElType==10) THEN
420 !!$ READ(14,*) glb%MatIdVol(i),glb%ElConnVol(1,i), &
421 !!$ glb%ElConnVol(2,i), glb%ElConnVol(3,i), &
422 !!$ glb%ElConnVol(4,i), glb%ElConnVol(5,i), &
423 !!$ glb%ElConnVol(6,i), glb%ElConnVol(7,i), &
424 !!$ glb%ElConnVol(8,i), glb%ElConnVol(9,i), &
425 !!$ glb%ElConnVol(10,i) !, iaux, iaux
426 !!$ ELSE IF( glb%iElType==8) THEN
427 !!$ READ(14,*) glb%MatIdVol(i),glb%ElConnVol(1,i), &
428 !!$ glb%ElConnVol(2,i), glb%ElConnVol(3,i), &
429 !!$ glb%ElConnVol(4,i), glb%ElConnVol(5,i), &
430 !!$ glb%ElConnVol(6,i), glb%ElConnVol(7,i), &
431 !!$ glb%ElConnVol(8,i) !, iaux, iaux
432 !!$ END IF
433 
434  IF(glb%NdBasedEl)THEN
435 
436  CALL volratio(glb%ElConnVol(1,i),glb%ElConnVol(2,i),glb%ElConnVol(3,i),glb%ElConnVol(4,i),glb%AlphaR(1:4,i),&
437  glb%MeshCoor,glb%NumNp,glb%NdMassLump)
438 
439  DO j = 1, 4
440  glb%NumElNeigh(glb%ElConnVol(j,i)) = glb%NumElNeigh(glb%ElConnVol(j,i)) + 1
441  glb%ElConnNd(glb%ElConnVol(j,i),glb%NumElNeigh(glb%ElConnVol(j,i))) = i
442  ENDDO
443  ENDIF
444  ENDDO
445  ENDDO
446 ! __________________________________________
447 ! _______ Parallel Communication
448 ! __________________________________________
449 
450  IF(numprocs.NE.0)THEN
451 
452 
453 
454 ! CALL COM_set_size( volWin//'.pconn', ip, NumNeighProcs*2+MaxNumNodesComm)
455 ! CALL COM_set_array(volWin//'.pconn', ip, Pconn_Comm, 1)
456 
457 
458 !!$ CALL COM_new_attribute( volWin//'.NodesToCommunicate', 'p', COM_INTEGER, 1, '')
459 !!$ CALL COM_set_size( volWin//'.NodesToCommunicate', ip, MaxNumNodesComm)
460 !!$ CALL COM_set_array(volWin//'.NodesToCommunicate', ip, NodesToCommunicate, 1)
461 !!$
462 !!$
463 !!$ CALL COM_new_attribute( volWin//'.ID_sendto_List', 'p', COM_INTEGER, 1, '')
464 !!$ CALL COM_set_size( volWin//'.ID_sendto_List', ip, NumNeighProcs)
465 !!$ CALL COM_set_array(volWin//'.ID_sendto_List', ip, ID_sendto_List, 1)
466 !!$
467 !!$ CALL COM_new_attribute( volWin//'.NumNeighProcs_List', 'p', COM_INTEGER, 1, '')
468 !!$ CALL COM_set_size( volWin//'.NumNeighProcs_List', ip, NumNeighProcs)
469 !!$ CALL COM_set_array(volWin//'.NumNeighProcs_List', ip, NumNeighProcs_List, 1)
470 
471 
472  ELSE
473  glb%TotNumNeighProcs = 1
474  ENDIF
475 
476 
477 
478  CALL com_get_size( volin//".pconn", myid+1, numparcomm)
479  CALL com_set_size( volwin//'.pconn', myid+1, numparcomm)
480  CALL com_resize_array(volwin//'.pconn', myid+1, parcomm, 1)
481 
482 !!$!
483 !!$! -- Read lst nodes that need to be sent to other processors
484 !!$!
485 !!$ READ(14,*) glb%TotNumNeighProcs
486 !!$
487 !!$! IF(glb%TotNumNeighProcs.NE.0)THEN
488 !!$ ALLOCATE(glb%NeighProcList(1:glb%TotNumNeighProcs))
489 !!$! ENDIF
490 !!$
491 !!$ ALLOCATE(glb%NumNdComm(1:glb%TotNumNeighProcs))
492 !!$ ALLOCATE(glb%NdCommList(1:glb%TotNumNeighProcs))
493 !!$ glb%TotNumNdComm = 0
494 !!$
495 !!$! NOTE: NeighProcList is the actual processor starting at id = 0, not 1
496 !!$
497 !!$ DO i = 1, glb%TotNumNeighProcs
498 !!$ READ(14,*) glb%NeighProcList(i),glb%NumNdComm(i)
499 !!$ ! print*,glb%NeighProcList(i),glb%NumNdComm(i)
500 !!$
501 !!$ glb%TotNumNdComm = glb%TotNumNdComm + glb%NumNdComm(i)
502 !!$
503 !!$ ALLOCATE(glb%NdCommList(i)%NdId(1:glb%NumNdComm(i)))
504 !!$ DO j=1,glb%NumNdComm(i)
505 !!$ READ(14,*) glb%NdCommList(i)%NdId(j)
506 !!$ ENDDO
507 !!$ ENDDO
508 !!$ glb%TotNumNdComm = glb%TotNumNdComm*3 ! x3 (x,y,z R_in)
509 !!$
510 !!$ CASE(7)
511 !!$
512 !!$ ! Packet 07
513 !!$ ! __________________________________________
514 !!$ ! ___ Nodal History
515 !!$ ! __________________________________________
516 !!$
517 !!$ READ(14,*) glb%NumNodeIO
518 !!$ ALLOCATE(glb%NodeIO(1:glb%NumNodeIO))
519 !!$ DO i = 1, glb%NumNodeIO
520 !!$ READ(14,*) glb%NodeIO(i)
521 !!$ WRITE(chr2,'(i2.2)') i
522 !!$ OPEN(33+i,FILE='Rocfrac/NdHistory.'//glb%MyIdChr//'.'//chr2)
523 !!$ WRITE(33+i,*) '#', glb%coor(1:3,glb%NodeIO(i))
524 !!$ ENDDO
525 !!$
526 !!$ CASE(8)
527 !!$
528 !!$ ! Packet 08
529 !!$ ! __________________________________________
530 !!$ ! ___ Nodal History
531 !!$ ! __________________________________________
532 !!$
533 !!$ READ(14,*) glb%NumNdsBCHT !, iaux
534 !!$ IF(glb%NumNdsBCHT.NE.0) ALLOCATE(glb%BCFlagHT(1:2,1:glb%NumNdsBCHT),glb%BCvalueHT(1,1:glb%NumNdsBCHT))
535 !!$ DO i=1,glb%NumNdsBCHT
536 !!$ READ(14,*) glb%BCFlagHT(1,i) ,NdBCflag !, iaux
537 !!$ glb%BCFlagHT(2,i) = glb%bcCondHT(NdBCflag)%BCtypeX
538 !!$ glb%BCvalueHT(1,i) = glb%bcCondHT(NdBCflag)%BCvalueX
539 !!$ ENDDO
540 !!$ CASE(99)
541 !!$
542 !!$ ! Packet 99
543 !!$ ! _____________________________
544 !!$ ! ___ End of Input File
545 !!$ ! _____________________________
546 !!$
547 !!$ EXIT
548 !!$
549 !!$ CASE DEFAULT
550 !!$ PRINT*, 'ERROR: ROCFRAC'
551 !!$ PRINT*, 'Input Packet = ',IdPacket, 'Is Not A Valid Option'
552 !!$ PRINT*, 'Check {prefix}.####.inp files and try again'
553 !!$
554 !!$ CALL MPI_BARRIER(glb%MPI_COMM_ROCFRAC,ierr)
555 !!$ CALL MPI_FINALIZE(glb%MPI_COMM_ROCFRAC,ierr)
556 !!$ STOP
557 !!$ END SELECT
558 !!$
559 !!$ ENDDO
560 !!$
561 !!$ CLOSE(14)
562 
563  CALL mpi_barrier(glb%MPI_COMM_ROCFRAC,ierr)
564  IF(myid.EQ.0 .AND. glb%debug_state) THEN
565  WRITE(6,'(A)') 'Rocfrac: Finished Reading Solids Mesh'
566  ENDIF
567 
568  IF((glb%DebondPart.eqv..false.).AND.(glb%DebondPart_Matous.eqv..false.))THEN
569 
570  CALL vol_elem_mat(glb%E,glb%xnu,glb%ci,glb%cj,glb%NumMatVol,glb%iElIntgratn)
571 
572  ELSE IF(glb%DebondPart_Matous.eqv..true.)THEN
573 
574  ALLOCATE( glb%StrainOld(1:4,1:glb%NumElVol*6) )
575  ALLOCATE( glb%SoftParam(1:4,1:glb%NumElVol) )
576  ALLOCATE( glb%cd(1:4,1:glb%NumElVol) )
577 
578  CALL vol_elem_mat_matous(glb)
579 
580  ELSE IF(glb%DebondPart.eqv..true.)THEN
581 
582  ALLOCATE(glb%STATEV_Part1(1:glb%NumElVol))
583  ALLOCATE(glb%STATEV_Part2(1:glb%NumElVol))
584 
585  glb%STATEV_Part1(:) = 1
586  glb%STATEV_Part2(:) = 1
587 
588  ALLOCATE(glb%StrainTrace(1:glb%NumElVol))
589  ENDIF
590 
591 
592  IF(glb%iElType.EQ.8)THEN
593 
594  allocate(glb%mixed_map(1:8,1:9,1:12))
595  allocate(glb%enhanced_map(1:8,1:9,1:9))
596  allocate(glb%Aenh(1:9,1:glb%NumElVol))
597 
598  glb%mixed_map = 0.d0
599  glb%enhanced_map = 0.d0
600  glb%Aenh(:,:) = 0.d0
601 
602  CALL enhanced_elem_maps_hex(glb%mixed_map,glb%enhanced_map)
603  allocate(glb%dmat(1:glb%NumMatVol,1:9,1:9))
604 
605  DO i = 1, glb%NumMatVol
606  CALL get_mat_stiffness(glb%E(i),glb%xnu(i),glb%dmat(i,:,:))
607  ENDDO
608 
609  ENDIF
610 
611 !-----Setting initial conditions
612  ALLOCATE(glb%Disp(1:3*glb%NumNP))
613  glb%Disp(:) = 0.d0
614 
615  ALLOCATE(glb%Accel(1:3*glb%NumNP))
616  glb%Accel(:) = 0.d0
617 
618  ALLOCATE(glb%DispOld(1:3*glb%NumNP))
619  ALLOCATE(glb%S11(1:glb%iStrGss,1:glb%NumElVol),glb%S22(1:glb%iStrGss,1:glb%NumElVol))
620  ALLOCATE(glb%S33(1:glb%iStrGss,1:glb%NumElVol) )
621  ALLOCATE(glb%S12(1:glb%iStrGss,1:glb%NumElVol),glb%S23(1:glb%iStrGss,1:glb%NumElVol))
622  ALLOCATE(glb%S13(1:glb%iStrGss,1:glb%NumElVol) )
623  ALLOCATE(glb%SVonMises(1:glb%NumElVol))
624 
625  IF(glb%HeatTransSoln .eqv. .true.)THEN
626  ALLOCATE(glb%Temperature(1:glb%NumNP))
627  glb%Temperature(1:glb%NumNP) = glb%Temperature0
628  ENDIF
629 
630  IF(glb%ArtificialDamping .eqv. .true.)THEN
631  ALLOCATE(glb%DetF_Old(1:glb%iStrGss,1:glb%NumElVol))
632  glb%DetF_Old(:,:) = 1.d0
633  endif
634 
635 ! new: ale formulation
636  ALLOCATE(glb%DispBar(1:3*glb%NumNP),glb%VeloBar(1:3*glb%NumNP),glb%AccelBar(1:3*glb%NumNP))
637  glb%DispBar(:) = 0.d0
638  glb%VeloBar(:) = 0.d0
639  glb%AccelBar(:) = 0.d0
640 
641 
642  ALLOCATE(glb%VeloHalf(1:3*glb%NumNP),glb%VeloBarOld(1:3*glb%NumNP))
643  glb%VeloHalf(:) = 0.d0
644  glb%VeloBarOld(:) = 0.d0
645 
646  ALLOCATE(glb%DispTotal(1:3*glb%NumNP))
647 
648  glb%S11(1:glb%iStrGss,1:glb%NumElVol) = 0.d0
649  glb%S22(1:glb%iStrGss,1:glb%NumElVol) = 0.d0
650  glb%S33(1:glb%iStrGss,1:glb%NumElVol) = 0.d0
651  glb%S12(1:glb%iStrGss,1:glb%NumElVol) = 0.d0
652  glb%S23(1:glb%iStrGss,1:glb%NumElVol) = 0.d0
653  glb%S13(1:glb%iStrGss,1:glb%NumElVol) = 0.d0
654  glb%SVonMises(1:glb%NumElVol) = 0.d0
655  glb%Accel = 0.d0
656  glb%VeloBar = 0.d0
657 ! glb%AccelBndry = 0.d0 ! not zero for restart.
658 ! glb%VeloBndry = 0.d0 ! not zero for restart.
659 
660  glb%DispBar = 0.d0
661  glb%DispTotal = 0.d0
662 
663 ! Create window for HDF output
664 
665  CALL com_new_attribute( volwin//'.disp', 'n', com_double, 3, 'm')
666  CALL com_new_attribute( volwin//'.disp_burn', 'n',com_double, 3, 'm')
667  CALL com_new_attribute( volwin//'.velo', 'n', com_double, 3, 'm/s')
668  CALL com_new_attribute( volwin//'.stress', 'e', com_double, 1, 'Pa')
669  CALL com_new_attribute( volwin//'.accel', 'n', com_double, 3, 'm/s^2')
670  CALL com_new_attribute( volwin//'.vbar', 'n', com_double, 3, 'm/s')
671  CALL com_new_attribute( volwin//'.S11', 'e', com_double, glb%iStrGss, 'Pa')
672  CALL com_new_attribute( volwin//'.S22', 'e', com_double, glb%iStrGss, 'Pa')
673  CALL com_new_attribute( volwin//'.S33', 'e', com_double, glb%iStrGss, 'Pa')
674  CALL com_new_attribute( volwin//'.S12', 'e', com_double, glb%iStrGss, 'Pa')
675  CALL com_new_attribute( volwin//'.S23', 'e', com_double, glb%iStrGss, 'Pa')
676  CALL com_new_attribute( volwin//'.S13', 'e', com_double, glb%iStrGss, 'Pa')
677 
678  CALL com_new_attribute( volwin//'.NumElPartBndry', 'p', com_integer, 1, '')
679  CALL com_new_attribute( volwin//'.NumElVolMat', 'p', com_integer, 1, '')
680  CALL com_new_attribute( volwin//'.NumElPartBndryMat', 'p', com_integer, 1, '')
681 
682  IF(glb%DebondPart .eqv. .true.)THEN
683  CALL com_new_attribute( volwin//'.StrainTrace', 'e', com_double, 1, ' ')
684  CALL com_new_attribute( volwin//'.DebondLg', 'e', com_double, 1, ' ')
685  CALL com_new_attribute( volwin//'.DebondSm', 'e', com_double, 1, ' ')
686  ENDIF
687 
688  IF(glb%DebondPart_Matous .eqv. .true.)THEN
689  CALL com_new_attribute( volwin//'.StrainOld', 'e', com_double, 4, ' ')
690  CALL com_new_attribute( volwin//'.SoftParam', 'e', com_double, 4, ' ')
691  ENDIF
692 
693  IF(glb%HeatTransSoln .eqv. .true.) CALL com_new_attribute( volwin//'.Temp', 'n', com_double, 1, 'K')
694 
695  IF ( glb%NumNP > 0) THEN
696 !
697 ! Registering Coordinates
698 !
699 
700 !!$ CALL COM_init_mesh( volWin//'.nc', MyId+1, glb%MeshCoor, glb%NumNP)
701 
702  CALL com_set_array(volwin//'.nc', myid+1, glb%MeshCoor,3)
703 !
704 ! Registering Element Connectivity
705 !
706  IF(glb%iElType.EQ.4)THEN
707 
708 !!$ CALL COM_init_mesh( volWin//'.T4', MyId+1, glb%ElConnVol, glb%NumElVol)
709 
710  CALL com_set_array( volwin//'.:T4', myid+1, glb%ElConnVol,4)
711 
712  ELSE IF(glb%iElType.EQ.10)THEN
713 
714 !!$ CALL COM_init_mesh( volWin//'.T10', MyId+1, glb%ElConnVol, glb%NumElVol)
715 
716  CALL com_set_array( volwin//'.:T10', myid+1, glb%ElConnVol,10)
717 
718  ELSE IF(glb%iElType.EQ.8)THEN
719 !!$ CALL COM_init_mesh( volWin//'.H8', MyId+1, glb%ElConnVol, glb%NumElVol)
720 
721  CALL com_set_array( volwin//'.:H8', myid+1, glb%ElConnVol,8)
722 
723  ENDIF
724 
725  CALL com_set_array( volwin//'.disp', myid+1, glb%Disp,3) ! was DispTotal
726  CALL com_set_array( volwin//'.disp_burn', myid+1, glb%DispBar,3)
727  CALL com_set_array( volwin//'.velo', myid+1, glb%VeloHalf,3)
728  CALL com_set_array( volwin//'.stress', myid+1, glb%SVonMises,1)
729  CALL com_set_array( volwin//'.accel', myid+1, glb%Accel,3)
730  CALL com_set_array( volwin//'.vbar', myid+1, glb%VeloBar,3)
731  CALL com_set_array( volwin//'.S11', myid+1, glb%S11)
732  CALL com_set_array( volwin//'.S22', myid+1, glb%S22)
733  CALL com_set_array( volwin//'.S33', myid+1, glb%S33)
734  CALL com_set_array( volwin//'.S12', myid+1, glb%S12)
735  CALL com_set_array( volwin//'.S23', myid+1, glb%S23)
736  CALL com_set_array( volwin//'.S13', myid+1, glb%S13)
737 
738  CALL com_set_size( volwin//'.NumElPartBndry', myid+1, 1)
739  CALL com_set_array( volwin//'.NumElPartBndry', myid+1,glb%NumElPartBndry,1 )
740 
741  CALL com_set_size( volwin//'.NumElVolMat', myid+1, glb%NumMatVol)
742  CALL com_set_array( volwin//'.NumElVolMat', myid+1, glb%NumElVolMat, 1)
743 
744  CALL com_set_size( volwin//'.NumElPartBndryMat', myid+1, glb%NumMatVol)
745  CALL com_set_array( volwin//'.NumElPartBndryMat',myid+1, glb%NumElPartBndryMat, 1 )
746 
747  IF(glb%HeatTransSoln.eqv..true.) CALL com_set_array( volwin//'.Temp', myid+1, glb%Temperature,1)
748  IF(glb%DebondPart.eqv..true.) THEN
749  CALL com_set_array( volwin//'.DebondLg', myid+1, glb%STATEV_Part1,1)
750  CALL com_set_array( volwin//'.DebondSm', myid+1, glb%STATEV_Part2,1)
751  CALL com_set_array( volwin//'.StrainTrace', myid+1, glb%StrainTrace,1)
752  ENDIF
753  IF(glb%DebondPart_Matous.eqv..true.) THEN
754  CALL com_set_array( volwin//'.StrainOld', myid+1, glb%StrainOld,4)
755  CALL com_set_array( volwin//'.SoftParam', myid+1, glb%SoftParam,4)
756  ENDIF
757 
758  ENDIF
759 
760  restartall = .false.
761  IF(initialtime.NE.0.d0)THEN
762  glb%ReStart = .true.
763  ENDIF
764 
765  CALL mpi_reduce(glb%ReStart, restartall, 1, mpi_logical, &
766  mpi_lor,0,glb%MPI_COMM_ROCFRAC,ierr)
767 
768  IF(glb%Verb.gt.0 .and. myid.eq.0 .and. &
769  restartall .eqv. .true.) THEN
770  WRITE(6,'(A)') 'Rocfrac: Restarting, Solids'
771  END IF
772 
773  CALL com_new_attribute(volwin//'.BCValue','p',com_double, 1, '')
774  CALL com_set_size( volwin//'.BCValue', myid+1, glb%NumNdsBCcrypt*6)
775  CALL com_set_array( volwin//'.BCValue', myid+1, glb%BCValueGlb, 1)
776 
777 
778 
779  CALL com_new_attribute(volwin//'.bcnode','p',com_integer, 2, '')
780  CALL com_set_size( volwin//'.bcnode', myid+1, glb%NumNdsBCcrypt)
781  CALL com_set_array( volwin//'.bcnode', myid+1, glb%BCFlagCrypt, 2)
782 
783  CALL com_new_attribute(volwin//'.MatType','e',com_integer, 1, '')
784  CALL com_set_array( volwin//'.MatType', myid+1, glb%MatIdVol, 1)
785 
786  glb%NumNdsBC = 0
787  glb%NumNdsBCmm = 0
788  glb%NumNdsBCHT = 0
789 
790  CALL com_window_init_done( volwin)
791 
792  CALL com_call_function( obtain_attr, 2, &
793  com_get_attribute_handle_const( volin//".all"), &
794  com_get_attribute_handle( volwin//".all"))
795 
796 
797 
798 
799 
800  glb%TotNumNeighProcs = 0
801 
802 !
803 ! -- Read lst nodes that need to be sent to other processors
804 !
805 !!$ iptr = 1
806 !!$ DO
807 !!$ IF(iptr.GT.NumParComm)EXIT
808 !!$ glb%TotNumNeighProcs = glb%TotNumNeighProcs + 1
809 !!$ iptr = iptr + ParComm( iptr + 1) + 2
810 !!$ ENDDO
811 
812  ! crashes on copper if no if statement
813  IF(numprocs.GT.1) glb%TotNumNeighProcs = parcomm(1)
814 
815 
816 ! IF(glb%TotNumNeighProcs.NE.0)THEN
817  ALLOCATE(glb%NeighProcList(1:glb%TotNumNeighProcs))
818 ! ENDIF
819 
820 
821  ALLOCATE(glb%NumNdComm(1:glb%TotNumNeighProcs))
822  ALLOCATE(glb%NdCommList(1:glb%TotNumNeighProcs))
823  glb%TotNumNdComm = 0
824 
825 
826 ! NOTE: NeighProcList is the actual processor starting at id = 0, not 1
827 
828  iptr = 1
829  DO i = 1, glb%TotNumNeighProcs
830  iptr = iptr + 1
831  glb%NeighProcList(i) = parcomm(iptr) -1
832  iptr = iptr + 1
833  glb%NumNdComm(i) = parcomm(iptr)
834  glb%TotNumNdComm = glb%TotNumNdComm + glb%NumNdComm(i)
835 
836  ALLOCATE(glb%NdCommList(i)%NdId(1:glb%NumNdComm(i)))
837 
838  DO j=1,glb%NumNdComm(i)
839  iptr = iptr + 1
840  glb%NdCommList(i)%NdId(j) = parcomm(iptr)
841  ENDDO
842 
843  ENDDO
844 
845 
846  glb%TotNumNdComm = glb%TotNumNdComm*3 ! x3 (x,y,z R_in)
847 
848  IF(glb%TotNumNeighProcs.NE.0)THEN
849  ALLOCATE(glb%RecvDataFrm(1:numprocs) )
850  ALLOCATE(glb%ReqRcv (1:glb%TotNumNeighProcs) )
851  ALLOCATE(glb%ReqSnd (1:glb%TotNumNeighProcs) )
852 
853  ALLOCATE(glb%StatSnd (1:mpi_status_size,1:glb%TotNumNeighProcs) )
854  ALLOCATE(glb%StatRcv (1:mpi_status_size,1:glb%TotNumNeighProcs) )
855  ENDIF
856 
857 !!$ IF(glb%HeatTransSoln)THEN
858 !!$ DO j = 1, glb%TotNumNeighProcs
859 !!$ k = glb%NeighProcList(j) + 1
860 !!$ ALLOCATE(glb%RecvDataFrm(k)%rcvbuf(1:glb%NumNdComm(j)*4))
861 !!$ ENDDO
862 !!$ ELSE
863  DO j = 1, glb%TotNumNeighProcs
864  k = glb%NeighProcList(j) + 1
865  ALLOCATE(glb%RecvDataFrm(k)%rcvbuf(1:glb%NumNdComm(j)*3))
866  ENDDO
867 !!$ ENDIF
868 
869 
870  DO i = 1, glb%NumNdsBCcrypt
871  ndbcflag = mod(glb%BCFlagCrypt(2,i),100)
872 ! structural boundary conditions
873  IF(ndbcflag.GT.0) glb%NumNdsBC = glb%NumNdsBC + 1
874 
875  ndbcflag = glb%BCFlagCrypt(2,i)/10000
876 ! mesh motion boundary conditions
877  IF(ndbcflag.GT.0) glb%NumNdsBCmm = glb%NumNdsBCmm + 1
878 
879  ndbcflag = mod(glb%BCFlagCrypt(2,i),10000)/100
880 ! thermal boundary conditions
881  IF(ndbcflag.GT.0) glb%NumNdsBCHT = glb%NumNdsBCHT + 1
882  ENDDO
883 
884  IF(glb%NumNdsBC.NE.0)THEN
885  ALLOCATE(glb%BCFlag(1:4,1:glb%NumNdsBC),glb%BCvalue(1:3,1:glb%NumNdsBC))
886  ALLOCATE(glb%AccelBndry(1:3*glb%NumNdsBC),glb%VeloBndry(1:3*glb%NumNdsBC))
887  ENDIF
888  IF(glb%NumNdsBCmm.NE.0) ALLOCATE(glb%BCFlagmm(1:4,1:glb%NumNdsBCmm),glb%BCvaluemm(1:3,1:glb%NumNdsBCmm))
889  IF(glb%NumNdsBCHT.NE.0) ALLOCATE(glb%BCFlagHT(1:2,1:glb%NumNdsBCHT),glb%BCvalueHT(1,1:glb%NumNdsBCHT))
890 
891  icnt1 = 0
892  icnt2 = 0
893  icnt3 = 0
894  icnt4 = 1
895 
896 !!$ IF ( glb%NumNdsBC > 0) THEN ! boundary conditions
897 !!$
898 !!$ CALL COM_new_attribute( volWin//'.velobndry', 'p', COM_DOUBLE, 3, 'm/s')
899 !!$ CALL COM_new_attribute( volWin//'.accbndry', 'p', COM_DOUBLE, 3, 'm/s^2')
900 !!$
901 !!$ CALL COM_set_size( volWin//'.velobndry', MyId+1, glb%NumNdsBC )
902 !!$ CALL COM_set_array( volWin//'.velobndry', MyId+1, glb%VeloBndry, 3)
903 !!$
904 !!$ CALL COM_set_size( volWin//'.accbndry', MyId+1, glb%NumNdsBC )
905 !!$ CALL COM_set_array( volWin//'.accbndry', MyId+1, glb%AccelBndry, 3)
906 !!$
907 !!$ ENDIF
908 
909 
910 
911  DO i = 1, glb%NumNdsBCcrypt
912 
913 
914  ndbcflag = mod(glb%BCFlagCrypt(2,i),100)
915 
916 
917 
918 ! structural boundary conditions
919  IF(ndbcflag.GT.0) THEN
920  icnt1 = icnt1 + 1
921  glb%BCFlag(1,icnt1) = glb%BCFlagCrypt(1,i)
922 
923  glb%BCFlag(2,icnt1) = glb%bcCond(ndbcflag)%BCtypeX
924  glb%BCFlag(3,icnt1) = glb%bcCond(ndbcflag)%BCtypeY
925  glb%BCFlag(4,icnt1) = glb%bcCond(ndbcflag)%BCtypeZ
926  glb%BCvalue(1,icnt1) = glb%bcCond(ndbcflag)%BCvalueX
927  glb%BCvalue(2,icnt1) = glb%bcCond(ndbcflag)%BCvalueY
928  glb%BCvalue(3,icnt1) = glb%bcCond(ndbcflag)%BCvalueZ
929 
930  glb%AccelBndry(icnt1*3-2:icnt1*3) = glb%BCValueGlb(icnt4:icnt4+2)
931  icnt4 = icnt4 + 3
932  glb%VeloBndry(icnt1*3-2:icnt1*3) = glb%BCValueGlb(icnt4:icnt4+2)
933  icnt4 = icnt4 + 3
934 
935  ENDIF
936 
937  ndbcflag = glb%BCFlagCrypt(2,i)/10000
938 
939 
940 
941 ! mesh motion boundary conditions
942  IF(ndbcflag.GT.0) THEN
943 
944  icnt2 = icnt2 + 1
945 
946  glb%BCFlagmm(1,icnt2) = glb%BCFlagCrypt(1,i)
947 
948 
949  glb%BCFlagmm(2,icnt2) = glb%bcCondmm(ndbcflag)%BCtypeX
950  glb%BCFlagmm(3,icnt2) = glb%bcCondmm(ndbcflag)%BCtypeY
951  glb%BCFlagmm(4,icnt2) = glb%bcCondmm(ndbcflag)%BCtypeZ
952  glb%BCvaluemm(1,icnt2) = glb%bcCondmm(ndbcflag)%BCvalueX
953  glb%BCvaluemm(2,icnt2) = glb%bcCondmm(ndbcflag)%BCvalueY
954  glb%BCvaluemm(3,icnt2) = glb%bcCondmm(ndbcflag)%BCvalueZ
955  ENDIF
956 
957  ndbcflag = mod(glb%BCFlagCrypt(2,i),10000)/100
958 ! thermal boundary conditions
959 
960  IF(ndbcflag.GT.0) THEN
961 
962  icnt3 = icnt3 + 1
963 
964  glb%BCFlagHT(1,icnt3) = glb%BCFlagCrypt(1,i)
965 
966  glb%BCFlagHT(2,icnt3) = glb%bcCondHT(ndbcflag)%BCtypeX
967  glb%BCvalueHT(1,icnt3) = glb%bcCondHT(ndbcflag)%BCvalueX
968  ENDIF
969  ENDDO
970 
971 
972  CALL rocfracinterfaceinitial( glb,obtain_attr,surfin)
973 
974 ! CALL MPI_BARRIER(glb%MPI_COMM_ROCFRAC,ierr)
975 ! STOP
976 
977 !!$ DO j = 1, glb%NumNP
978 !!$ j1 = j*3
979 !!$ glb%MeshCoor(1,j) = glb%coor(1,j) + glb%DispBar(j1 - 2)
980 !!$ glb%MeshCoor(2,j) = glb%coor(2,j) + glb%DispBar(j1 - 1)
981 !!$ glb%MeshCoor(3,j) = glb%coor(3,j) + glb%DispBar(j1 )
982 !!$ ENDDO
983 
984 
985  CALL rocfracinterfacebuff(glb)
986 
987 
988  CALL com_call_function( man_init, 2, surwin, volwin)
989 
990  IF(glb%ALEenabled.eqv..false.)THEN
991  CALL updatemassmatrix(glb)
992 
993  ENDIF
994 11 FORMAT(a,'_',a,a1)
995 
996 !!$ IF(myid.EQ.0) PRINT*,'ROCFRAC :: Finished RocFracInitialize'
997 !!$ CALL MPI_BARRIER(glb%MPI_COMM_ROCFRAC,ierr)
998 !!$ CALL MPI_FINALIZE(glb%MPI_COMM_ROCFRAC,ierr)
999 
1000 ! -- Mass/Volume Conservation
1001 
1002  glb%TotalMassSolidp = 0.d0
1003  glb%TotalGeomVolp = 0.d0
1004 
1005  IF(glb%iElType.EQ.4.OR.glb%iElType.EQ.10 .AND.(.NOT.(glb%NdBasedEl)))THEN
1006  CALL v3d4_volume(glb%MeshCoor,glb%ElConnVol,glb%MatIdVol,glb%rho, &
1007  glb%NumNP,glb%NumElVol,glb%NumMatVol,glb%Disp,1,glb%NumElVol,&
1008  glb%TotalMassSolidp,glb%TotalGeomVolp,glb%TotalGeomUndefVolp, &
1009  glb%iElType)
1010  ENDIF
1011  IF(glb%NumProbesNd.NE.0)THEN
1012  CALL findprobe(glb,myid)
1013  ENDIF
1014  glb%iAmpCnt = 1
1015 
1016  glb%prop = 0.d0
1017  glb%slope = 0.d0
1018 
1019 
1020 
1021 !!$ PRINT*,'sdfsdflkjsdflkj',myid
1022 !!$ PRINT*,'finished readsdv',myid
1023 !!$ CALL MPI_BARRIER(glb%MPI_COMM_ROCFRAC,i)
1024 !!$ stop
1025 
1026  inquire(file='Rocfrac/Rocin/OverlayMappings.txt',exist=glb%OverlayExist)
1027 
1028  IF(glb%OverlayExist)THEN
1029  CALL readsdv(glb,myid)
1030 
1031  OPEN(456,file ='Rocfrac/Rocin/OverlayMappings.txt')
1032  DO
1033 
1034  READ(456,*,iostat=ios) iprocs
1035  IF(ios.LT.0) EXIT
1036 
1037  IF(iprocs.EQ.myid)THEN
1038 
1039 !!$ READ(456,*) glb%nf1
1040 !!$ ALLOCATE(MapFaceEl2Vol1a(1:600),FaceOfVolEL1a(1:600))
1041 !!$
1042 !!$ DO i = 1, glb%nf1
1043 !!$ READ(456,'()')
1044 !!$ ENDDO
1045 !!$
1046 !!$ READ(456,*) iaux
1047 !!$
1048 !!$ DO i = 1, iaux
1049 !!$ READ(456,'()')
1050 !!$ ENDDO
1051 !!$ PRINT*,';;',myid
1052 
1053 !!$ READ(456,*) glb%nf1
1054 !!$ PRINT*,glb%nf1,myid
1055 
1056  READ(456,*) glb%nf1
1057  ALLOCATE(glb%MapFaceEl2Vol1(1:glb%nf1),glb%FaceOfVolEL1(1:glb%nf1))
1058 
1059  DO i = 1, glb%nf1
1060  READ(456,*) glb%MapFaceEl2Vol1(i), glb%FaceOfVolEL1(i)
1061  ENDDO
1062 
1063  READ(456,*) glb%nf2
1064 
1065  ALLOCATE(glb%MapFaceEl2Vol2(1:glb%nf2),glb%FaceOfVolEL2(1:glb%nf2))
1066 
1067  DO i = 1, glb%nf2
1068  READ(456,*) glb%MapFaceEl2Vol2(i), glb%FaceOfVolEL2(i)
1069  ENDDO
1070 !!$
1071 !!$ READ(456,*) iaux
1072 !!$
1073 !!$ DO i = 1, iaux
1074 !!$ READ(456,'()')
1075 !!$ ENDDO
1076 !!$
1077 !!$
1078 !!$ DO i = 1, glb%nf1
1079 !!$ READ(456,'()')
1080 !!$ ENDDO
1081  EXIT
1082 
1083  ELSE
1084  READ(456,*) iaux
1085 
1086  DO i = 1, iaux
1087  READ(456,'()')
1088  ENDDO
1089 
1090  READ(456,*) iaux
1091 
1092  DO i = 1, iaux
1093  READ(456,'()')
1094  ENDDO
1095  ENDIF
1096  ENDDO
1097  close(456)
1098 
1099 endif
1100 
1101 !!$ IF(myid.EQ.1)THEN
1102 !!$ DO i = 1, glb%nf1
1103 !!$ PRINT*, i,glb%MapFaceEl2Vol1(i), glb%FaceOfVolEL1(i)
1104 !!$ ENDDO
1105 !!$
1106 !!$ DO i = 1, glb%nf2
1107 !!$ PRINT*,i, glb%MapFaceEl2Vol2(i), glb%FaceOfVolEL2(i)
1108 !!$ ENDDO
1109 !!$ ENDIF
1110 
1111 !!$ do i = 1, glb%NumNp
1112 !!$ IF( glb%MeshCoor(1,i).LE.1.0001.AND.glb%MeshCoor(2,i).GT.0.02499.AND.glb%MeshCoor(3,i).GT.0.02499)THEN
1113 !!$ glb%NumNodeIO = i
1114 !!$ ENDIF
1115 !!$ enddo
1116 
1117 
1118  END SUBROUTINE rocfracinitialize
1119 
1120 
1121 !
1122 ! ----------------------------------------------------------------------- RocFracFinalize
1123 
1124  SUBROUTINE rocfracfinalize( glb)
1125 
1126  IMPLICIT NONE
1127  include 'roccomf90.h'
1128 
1129  TYPE(rocfrac_global), POINTER :: glb
1130 
1131  CALL com_delete_window( surwin)
1132  CALL com_delete_window( volwin)
1133 
1134  END SUBROUTINE rocfracfinalize
1135 !
1136 ! ----------------------------------------------------------------------- RocFracSoln
1137 
1138  SUBROUTINE rocfracsoln( glb, CurrentTime, CurrentTimeStep, MAN_update_inbuff)
1139 
1140  IMPLICIT NONE
1141  include 'mpif.h'
1142 
1143  TYPE(rocfrac_global), POINTER :: glb
1144  REAL*8, INTENT(IN) :: currenttime, currenttimestep
1145  INTEGER, INTENT(IN) :: man_update_inbuff
1146 
1147  INTEGER, SAVE:: istep
1148 
1149 
1150  INTEGER j,jj,jjj,k,k1,k2,j1,i,j2,j3,k3,k4 ! loop counters & dummy variables
1151 
1152  INTEGER :: ntime_out ! processor's id ( 0-(nprocs-1) )
1153 !-- First section Rnet = R_bar
1154 !-- Force Sum: cohesive traction + (-internal) + (-damping)+ external
1155 !-- i.e. R_co - R_in - R_damp + R_ex
1156  REAL*8, ALLOCATABLE, DIMENSION(:) :: rnet
1157  REAL*8, ALLOCATABLE, DIMENSION(:) :: rnetht
1158 !-- reciprical of mass matrix diagonal
1159  REAL*8 :: tvol, tvol_com
1160  REAL*8, DIMENSION(2) :: trry_dt,trry_vol,trry_calc,trry_dv_com
1161  REAL*8, DIMENSION(2) :: trry_energy,trry_io,trry_coh,trry_vol_com
1162  REAL*8, ALLOCATABLE, DIMENSION(:) :: buf
1163 
1164 
1165 ! Dummy variable
1166  REAL*8 :: a1
1167  CHARACTER*8 :: ichr8 ! output file number (a character)
1168  CHARACTER*3 :: ai,ai1 ! output file number (a character)
1169  INTEGER :: myid, ierr ! mpi error variable
1170 
1171  INTEGER :: isubstep,nsubstep
1172  REAL*8 :: dt_solid_sub, alpha
1173  REAL*8 :: currenttimestepsolid
1174 
1175  REAL*8 :: prin1, prin2, prin3
1176 
1177  INTEGER, SAVE :: ifirststep
1178 
1179  REAL*8 :: maxvalue
1180  INTEGER :: maxlocat
1181  INTEGER :: jaux
1182  INTEGER :: elemstart, elemend, numprocs
1183  REAL*8 :: dispprev
1184 
1185  INTEGER :: ndbcflag, icnt1, icnt5
1186 
1187 
1188  logical :: debug
1189  CHARACTER*4 :: ichr1, ichr2
1190 
1191 ! _____________________________
1192 ! _____ SOLUTION LOOP
1193 ! _____________________________
1194 
1195  CALL mpi_comm_rank(glb%MPI_COMM_ROCFRAC,myid,ierr)
1196  CALL mpi_comm_size(glb%MPI_COMM_ROCFRAC,numprocs,ierr)
1197 
1198  debug = glb%debug_state
1199 
1200  tvol = 0.d0; tvol_com = 0.d0
1201  nsubstep = 1
1202 
1203  currenttimestepsolid = currenttimestep
1204  ifirststep = 0
1205 
1206  glb%DispOld(:) = glb%Disp(:)
1207  ALLOCATE(rnet(1:3*glb%NumNP))
1208  IF(glb%HeatTransSoln) ALLOCATE(rnetht(1:glb%NumNP))
1209 
1210  IF(myid.EQ.0 .AND. glb%Verb.gt.0)THEN
1211  WRITE(6,'(A)') 'Rocfrac:'
1212  WRITE(6,'(A,A21)')'Rocfrac:','(Current DT'
1213  WRITE(6,'(A,A6,4(A15))')'Rocfrac:','Step',&
1214  '- Solid DT)', 'Global DT',&
1215  'Current DT','Global Time'
1216  WRITE(6,'(A,A)')'Rocfrac: ---------------------------',&
1217  '--------------------------------------'
1218  ENDIF
1219 
1220  glb%CurrTime = currenttime
1221 
1222 ! open(3001,file='ndhistory')
1223 
1224  DO WHILE (nsubstep.EQ.1)
1225 
1226  CALL max_dt_solid(dt_solid_sub,glb)
1227 
1228  IF(dt_solid_sub.GE.currenttimestepsolid)THEN
1229  nsubstep = 0
1230  glb%DT = currenttimestepsolid
1231  currenttimestepsolid = 0
1232  ELSE
1233  nsubstep = 1
1234  glb%DT = dt_solid_sub
1235  currenttimestepsolid = currenttimestepsolid - glb%DT
1236  ENDIF
1237 
1238  glb%DTInv = 1.d0 / glb%DT
1239 
1240  istep = istep + 1
1241 
1242  ! Compute alpha here
1243  alpha = 1.d0 - currenttimestepsolid / currenttimestep
1244  CALL com_call_function( man_update_inbuff, 1, alpha)
1245 
1246  glb%CurrTime = currenttime + (currenttimestep-currenttimestepsolid)
1247 
1248  IF(myid.EQ.0 .AND. glb%Verb.gt.0) WRITE(6,'(a,i6,4e15.4)')'Rocfrac:',&
1249  istep,(currenttimestep-currenttimestepsolid),&
1250  glb%DT,currenttimestep,glb%CurrTime
1251 
1252 
1253 !-- (0) INITIALIZE
1254 
1255 
1256  ntime_out = 0
1257  rnet(:) = 0.d0
1258  IF(glb%HeatTransSoln) rnetht(:) = 0.d0
1259 
1260 
1261  IF( glb%AmplitudeTable )THEN
1262 
1263  IF(glb%iAmpCnt.LE.glb%NumEntries-1)THEN
1264  IF(glb%CurrTime.GT.glb%AmpTable(1,glb%iAmpCnt+1))THEN
1265  glb%iAmpCnt = glb%iAmpCnt + 1
1266  ENDIF
1267  ENDIF
1268  glb%slope = glb%AmpTable(2,glb%iAmpCnt)
1269  glb%prop = glb%CurrTime*glb%AmpTable(2,glb%iAmpCnt) + glb%AmpTable(3,glb%iAmpCnt)
1270 
1271  ENDIF
1272 
1273 
1274  IF(glb%ALEenabled)THEN
1275 
1276 ! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
1277 
1278  DO j = 1, glb%InterfaceSFNumNodes
1279  k3 = 3*glb%MapNodeSF(j)
1280  k2 = k3 - 1
1281  k1 = k3 - 2
1282  glb%VeloBar(k1) = glb%InterfaceSFVbar(1,j)
1283  glb%VeloBar(k2) = glb%InterfaceSFVbar(2,j)
1284  glb%VeloBar(k3) = glb%InterfaceSFVbar(3,j)
1285  ENDDO
1286 
1287 
1288 !-- (2) update mesh displacement vector (ale)
1289 
1290  glb%DispBar(:) = glb%DT*glb%VeloBar(:) + glb%DispBar(:)
1291 
1292 !-- (3) UPDATE MESH POSITION
1293 
1294 
1295 
1296  DO j = 1, glb%NumNP
1297  j1 = j*3
1298  glb%MeshCoor(1,j) = glb%Meshcoor(1,j) + glb%DT*glb%VeloBar(j1 - 2)
1299  glb%MeshCoor(2,j) = glb%Meshcoor(2,j) + glb%DT*glb%VeloBar(j1 - 1)
1300  glb%MeshCoor(3,j) = glb%Meshcoor(3,j) + glb%DT*glb%VeloBar(j1 )
1301 
1302  ENDDO
1303 
1304 
1305 
1306  CALL updatemassmatrix(glb)
1307 
1308 !-- (5) CALCULATE R_BAR
1309 
1310  CALL updaterbar(glb,rnet)
1311 
1312 !-- (6) CALCULATE MESH VELOCITY VECTOR
1313  DO j = 1, glb%NumNP
1314  jaux = 3*j
1315  DO k = jaux-2, jaux
1316  glb%VeloBarOld(k) = glb%VeloBar(k)
1317  glb%VeloBar(k) = glb%xmass(j)*glb%rho(1)*(-rnet(k))*glb%kappa ! move - to subroutine
1318  ENDDO
1319  ENDDO
1320 
1321 !---(7) CALCULATE R_CO
1322 
1323 ! ! fix for parallel to do boundary then inside
1324 ! if(numco .gt. 0) then ! change idpressload bounds to 1:2
1325 ! call CST_COH(glb%MeshCoor,lmc,matc,Rnet,d,deltan,deltat,
1326 ! $ sigmax,taumax,glb%Sinit,Sthresh_cst,istep,glb%NumNP,numco,
1327 ! $ numat_coh,delta,numploadelem,idpressload, pressload,
1328 ! $ glb%NumNdsBCmesh, idmesh, rmesh, nboundtype, num_fail_coh,
1329 ! $ index_fail_coh, npress, nmesh, ielem_coh, iface_coh)
1330 ! endif
1331 
1332 !-- (8) Enforce Fixed BC's of Mesh for other boundaries not in contact with fluids
1333 
1334 !!$ DO j = 1, glb%NumNdsBCmm
1335 !!$ k3 = glb%BCFlagmm(1,j)*3
1336 !!$ k2 = k3-1
1337 !!$ k1 = k3-2
1338 !!$ IF(glb%BCFlagmm(2,j).eq.0) glb%VeloBar(k1) = 0.d0 ! rmesh(1,j)*meshvelo
1339 !!$ IF(glb%BCFlagmm(3,j).eq.0) glb%VeloBar(k2) = 0.d0 ! rmesh(2,j)*meshvelo
1340 !!$ IF(glb%BCFlagmm(4,j).eq.0) glb%VeloBar(k3) = 0.d0 ! rmesh(3,j)*meshvelo
1341 !!$ ENDDO
1342 
1343 ! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
1344 
1345  DO j = 1, glb%InterfaceSNumNodes
1346  k3 = 3*glb%MapNodeS(j)
1347  IF(k3.gt.0) THEN
1348  k2 = k3 - 1
1349  k1 = k3 - 2
1350  glb%VeloBar(k1) = glb%InterfaceSVbar(1,j)
1351  glb%VeloBar(k2) = glb%InterfaceSVbar(2,j)
1352  glb%VeloBar(k3) = glb%InterfaceSVbar(3,j)
1353  IF(glb%InterfaceSVbar(1,j).NE.0.d0) glb%VeloBar(k1) = glb%InterfaceSVbar(1,j)
1354  IF(glb%InterfaceSVbar(2,j).NE.0.d0) glb%VeloBar(k2) = glb%InterfaceSVbar(2,j)
1355  IF(glb%InterfaceSVbar(3,j).NE.0.d0) glb%VeloBar(k3) = glb%InterfaceSVbar(3,j)
1356  ENDIF
1357  ENDDO
1358 
1359 
1360 
1361 ! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
1362 
1363  DO j = 1, glb%InterfaceSFNumNodes
1364  k3 = 3*glb%MapNodeSF(j)
1365  k2 = k3 - 1
1366  k1 = k3 - 2
1367  glb%VeloBar(k1) = glb%InterfaceSFVbar(1,j)
1368  glb%VeloBar(k2) = glb%InterfaceSFVbar(2,j)
1369  glb%VeloBar(k3) = glb%InterfaceSFVbar(3,j)
1370 ! IF(glb%InterfaceSFVbar(1,j).NE.0.d0) glb%VeloBar(k1) = glb%InterfaceSFVbar(1,j)
1371 ! IF(glb%InterfaceSFVbar(2,j).NE.0.d0) glb%VeloBar(k2) = glb%InterfaceSFVbar(2,j)
1372 ! IF(glb%InterfaceSFVbar(3,j).NE.0.d0) glb%VeloBar(k3) = glb%InterfaceSFVbar(3,j)
1373  ENDDO
1374 
1375 ! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
1376 
1377 !!$ DO j = 1, glb%InterfaceSNumNodes
1378 !!$ k3 = 3*glb%MapNodeS(j)
1379 !!$ k2 = k3 - 1
1380 !!$ k1 = k3 - 2
1381 !!$ glb%VeloBar(k1) = glb%InterfaceSVbar(1,j)
1382 !!$ glb%VeloBar(k2) = glb%InterfaceSVbar(2,j)
1383 !!$ glb%VeloBar(k3) = glb%InterfaceSVbar(3,j)
1384 !!$! IF(glb%InterfaceSVbar(1,j).NE.0.d0) glb%VeloBar(k1) = glb%InterfaceSVbar(1,j)
1385 !!$! IF(glb%InterfaceSVbar(2,j).NE.0.d0) glb%VeloBar(k2) = glb%InterfaceSVbar(2,j)
1386 !!$! IF(glb%InterfaceSVbar(3,j).NE.0.d0) glb%VeloBar(k3) = glb%InterfaceSVbar(3,j)
1387 !!$ ENDDO
1388 ! Used for mass conservation
1389 !!$ DO j = 1, glb%InterfaceSFNumNodes
1390 !!$ k3 = 3*glb%MapNodeSF(j)
1391 !!$ k2 = k3 - 1
1392 !!$ k1 = k3 - 2
1393 !!$ glb%VeloBar(k1) = 0.d0
1394 !!$ glb%VeloBar(k2) = 0.d0
1395 !!$ glb%VeloBar(k3) = 1.d0
1396 !!$ ENDDO
1397 
1398  DO j = 1, glb%NumNdsBCmm
1399  k3 = glb%BCFlagmm(1,j)*3
1400  k2 = k3-1
1401  k1 = k3-2
1402  IF(glb%BCFlagmm(2,j).EQ.0) glb%VeloBar(k1) = 0.d0 ! rmesh(1,j)*meshvelo
1403  IF(glb%BCFlagmm(3,j).EQ.0) glb%VeloBar(k2) = 0.d0 ! rmesh(2,j)*meshvelo
1404  IF(glb%BCFlagmm(4,j).EQ.0) glb%VeloBar(k3) = 0.d0 ! rmesh(3,j)*meshvelo
1405  ENDDO
1406 
1407  rnet(:) = 0.d0
1408 
1409 ! -- Calculate mesh acceleration
1410 
1411  glb%AccelBar(:) = ( glb%VeloBar(:) - glb%VeloBarOld(:) ) * glb%DTInv
1412  IF(glb%iElType.EQ.4)THEN
1413 
1414  CALL v3d4_ale(glb%VeloBar,glb%AccelBar,glb%Disp,glb%VeloHalf,rnet, &
1415  glb%E,glb%xnu,glb%rho,glb%NumNP,glb%NumMatVol, &
1416  glb%NumElVol,glb%MatIdVol,glb%ElConnVol,glb%MeshCoor, & !Fix not work for more then one material
1417  1,glb%NumElPartBndry)
1418 
1419  ELSE IF(glb%iElType.EQ.10)THEN
1420  CALL v3d10_ale(glb%VeloBar,glb%AccelBar,glb%Disp,glb%VeloHalf,rnet, &
1421  glb%E,glb%xnu,glb%rho,glb%NumNP,glb%NumMatVol, &
1422  glb%NumElVol,glb%MatIdVol,glb%ElConnVol,glb%MeshCoor, & !Fix not work for more then one material
1423  1,glb%NumElPartBndry)
1424 
1425  ENDIF
1426 
1427  ENDIF ! END ALE OPTION
1428 
1429  IF(myid.Eq.0 .AND. debug) THEN
1430  WRITE(6,'(A)') 'Rocfrac: finished ale'
1431  ENDIF
1432 
1433 
1434 
1435 
1436 ! -- Transfer the Tractions Due to Fluid Pressure to the 3D Solid
1437 
1438 
1439 !!$ IF(glb%DefConfig)THEN
1440 
1441 !!$ CALL TractLoadDef(Rnet,glb%NumNP, &
1442 !!$ glb%InterfaceSFElemTract, &
1443 !!$ glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1444 !!$ glb%InterfaceSFElemConn, &
1445 !!$ glb%MapNodeSF,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp)
1446 
1447 
1448 
1449 !!$ ELSE
1450 !!$
1451 ! CALL TractLoad(Rnet,glb%NumNP, &
1452 ! glb%InterfaceSFElemTract, &
1453 ! glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1454 ! glb%InterfaceSFElemConn, &
1455 ! glb%MapNodeSF,glb%LwrBnd,glb%UppBnd,glb%Meshcoor)
1456 !!$
1457 !!$ ENDIF
1458 
1459  IF(myid.Eq.0 .AND. debug) THEN
1460  WRITE(6,'(A)') 'Rocfrac: Applying Pressure Loading'
1461  ENDIF
1462 
1463  IF(glb%iElType.EQ.8)THEN
1464 
1465 ! Burning, Interacting Surfaces Traction
1466 
1467  CALL tractload_hex(rnet,glb%NumNP, & ! Don't imposed pressure on deformed coordinates
1468  glb%InterfaceSFElemTract, &
1469  glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1470  glb%InterfaceSFElemConn, &
1471  glb%MapNodeSF,glb%LwrBnd,glb%UppBnd,glb%Meshcoor)
1472 
1473 ! Non-Burning, Interacting Surfaces Traction
1474 
1475  CALL tractload_hex(rnet,glb%NumNP, & ! Don't imposed pressure on deformed coordinates
1476  glb%InterfaceSFnbElemTract, &
1477  glb%InterfaceSFnbNumElems, glb%InterfaceSFnbNumNodes, &
1478  glb%InterfaceSFnbElemConn, &
1479  glb%MapNodeSFnb,glb%LwrBnd,glb%UppBnd,glb%Meshcoor)
1480 
1481 
1482  ELSE
1483 
1484  CALL fluidpressload(glb%NumNP,rnet, &
1485  glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1486  glb%InterfaceSFElemConn, &
1487  glb%MapNodeSF,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSFElVolEl,&
1488  glb%ElConnVol,glb%iElType,glb%NumElVol,glb%InterfaceSFElemTract)
1489 
1490 
1491  CALL fluidpressload(glb%NumNP,rnet, &
1492  glb%InterfaceSFnbNumElems, glb%InterfaceSFnbNumNodes, &
1493  glb%InterfaceSFnbElemConn, &
1494  glb%MapNodeSFnb,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSFnbElVolEl,&
1495  glb%ElConnVol,glb%iElType,glb%NumElVol,glb%InterfaceSFnbElemTract)
1496  ENDIF
1497 
1498 
1499  IF(myid.Eq.0 .AND. debug) THEN
1500  WRITE(6,'(A)') 'Rocfrac: End Pressure Loading'
1501  ENDIF
1502 !!$ CALL TractLoad(Rnet,glb%NumNP, &
1503 !!$ glb%InterfaceSFElemTract, &
1504 !!$ glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1505 !!$ glb%InterfaceSFElemConn, &
1506 !!$ glb%MapNodeSF,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSFElVolEl,&
1507 !!$ glb%ElConnVol,glb%iElType,glb%NumElVol)
1508 
1509 ! -- Enforce tractions (if any) where there are no fluid's domains.
1510 ! Uses the non-solid/fluid surface mesh.
1511 
1512 
1513  IF(myid.Eq.0 .AND. debug) THEN
1514  WRITE(6,'(A)') 'Rocfrac: Start traction loading'
1515  ENDIF
1516 
1517  IF(glb%EnforceTractionS.OR.glb%EnforceTractionSF)THEN
1518  ! IF(glb%UnDefConfig)THEN
1519 ! if using cauchy stress
1520 ! CALL TractPressLoadDef(Rnet,glb%NumNP, &
1521 ! glb%InterfaceSNumElems, glb%InterfaceSNumNodes, &
1522 ! glb%InterfaceSElemConn, &
1523 ! glb%MapNodeS,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSElVolEl,&
1524 ! glb%ElConnVol,glb%iElType,glb%NumElVol,glb%DummyTractVal*glb%prop)
1525 !
1526 ! IF(myid.EQ.0) PRINT*,'Pressure =', glb%DummyTractVal*glb%prop
1527 
1528 ! ELSE
1529 
1530  IF(glb%iElType.EQ.8)THEN
1531 
1532  ! fix need other option
1533  IF(glb%EnforceTractionS)THEN
1534  CALL tractloadpress_hex(rnet,glb%NumNP, &
1535  glb%InterfaceSNumElems, glb%InterfaceSNumNodes, &
1536  glb%InterfaceSElemConn, &
1537  glb%MapNodeS,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%DummyTractVal*glb%prop)
1538  ENDIF
1539 
1540  ELSE
1541 
1542  IF(glb%EnforceTractionS)THEN
1543 
1544 
1545  CALL tractpressload(rnet,glb%NumNP, &
1546  glb%InterfaceSNumElems, glb%InterfaceSNumNodes, &
1547  glb%InterfaceSElemConn, &
1548  glb%MapNodeS,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSElVolEl,&
1549  glb%ElConnVol,glb%iElType,glb%NumElVol,glb%DummyTractVal*glb%prop)
1550 
1551  IF(myid.EQ.0 .AND. glb%Verb.gt.1) THEN
1552  WRITE(6,'(A,F12.4)') 'Rocfrac: Pressure Solid =', &
1553  glb%DummyTractVal*glb%prop
1554  ENDIF
1555 
1556  ENDIF
1557 
1558  IF(glb%EnforceTractionSF)THEN
1559 
1560  CALL tractpressload(rnet,glb%NumNP, &
1561  glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1562  glb%InterfaceSFElemConn, &
1563  glb%MapNodeSF,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSFElVolEl,&
1564  glb%ElConnVol,glb%iElType,glb%NumElVol,glb%DummyTractVal*glb%prop)
1565 
1566  IF(myid.EQ.0 .AND. glb%Verb.gt.1) THEN
1567  WRITE(6,'(A,F12.4)') 'Rocfrac: Pressure Solid/Fluid =',&
1568  glb%DummyTractVal*glb%prop
1569  ENDIF
1570 
1571 
1572  CALL tractpressload(rnet,glb%NumNP, &
1573  glb%InterfaceSFnbNumElems, glb%InterfaceSFnbNumNodes, &
1574  glb%InterfaceSFnbElemConn, &
1575  glb%MapNodeSFnb,glb%LwrBnd,glb%UppBnd,glb%Meshcoor,glb%Disp,glb%MapSFnbElVolEl,&
1576  glb%ElConnVol,glb%iElType,glb%NumElVol,glb%DummyTractVal*glb%prop)
1577 
1578  ENDIF
1579 
1580  ENDIF
1581 ! ENDIF
1582  ENDIF
1583 
1584 ! Rnet(:) = 0.d0
1585 
1586 ! -- Calculate R_in, R_damp
1587 
1588  IF(myid.EQ.0 .AND. debug) THEN
1589  WRITE(6,'(A)') 'Rocfrac: start UpdateStructural'
1590  ENDIF
1591 
1592  IF(glb%HeatTransSoln)THEN
1593 
1594  CALL heatfluxload(rnetht,glb%NumNP, &
1595  glb%InterfaceSFNumElems, glb%InterfaceSFNumNodes, &
1596  glb%InterfaceSFElemConn, &
1597  glb%MapNodeSF,glb%LwrBnd,glb%UppBnd, glb%MeshCoor, &
1598  glb%MapSFElVolEl,glb%ElConnVol,glb%iElType,glb%NumElVol, &
1599  glb%InterfaceSFHeatFlux)
1600 
1601 
1602 
1603 
1604  CALL updatestructuralht(glb,numprocs,rnet,rnetht)
1605 
1606 
1607 
1608  IF(myid.EQ.0 .AND.glb%Verb.gt.1) THEN
1609  WRITE(6,'(A)') 'Rocfrac: MaxTemperature',maxval(glb%Temperature)
1610  ENDIF
1611 
1612 ! should not you have updatestrural here
1613  ELSE
1614  CALL updatestructural(glb,numprocs,rnet)
1615  ENDIF
1616 
1617 
1618  IF(myid.EQ.0 .AND. debug) THEN
1619  WRITE(6,'(A)') 'Rocfrac: finished UpdateStructural'
1620  ENDIF
1621 
1622  CALL principal_stress(glb%S11,glb%S22,glb%S33, &
1623  glb%S12,glb%S23,glb%S13, &
1624  glb%iStrGss,glb%NumElVol,glb%SVonMises)
1625 
1626 
1627 ! -- (12) UPDATE THE ACCELERATION AND VELOCITY
1628 
1629 
1630  IF(glb%ALEenabled)THEN
1631  DO j = 1, glb%NumNP
1632  jaux = 3*j
1633  DO k = jaux-2, jaux
1634 
1635 !
1636 ! Accelerations are found by
1637 ! .. -1 ext in
1638 ! { D } = [M] * ( { F } - { F } )
1639 ! i+1
1640 !
1641  a1 = rnet(k) * glb%xmass(j)
1642 
1643 ! Store Accelerations
1644 
1645  glb%Accel(k) = a1
1646  ENDDO
1647  ENDDO
1648 
1649 !-- (2) update mesh displacement vector (ale)
1650 
1651 ! glb%DispBar(:) = glb%DT*glb%VeloBar(:) + glb%DispBar(:)
1652 
1653  ELSE IF(.NOT.(glb%DampEnabled))THEN
1654 
1655 
1656  DO j = 1, glb%NumNP
1657  jaux = 3*j
1658  DO k = jaux-2, jaux
1659 !
1660 ! Accelerations are found by
1661 ! .. -1 ext in
1662 ! { D } = [M] * ( { F } - { F } )
1663 ! i+1
1664 !
1665  a1 = rnet(k) * glb%xmass(j)
1666 
1667 
1668 !
1669 ! Velocities are found by
1670 ! . . .. ..
1671 ! { D } = { D } + 0.5*Dt*( { D } + { D } )
1672 ! i+1 i i i+1
1673 
1674  glb%VeloHalf(k) = glb%VeloHalf(k) + glb%DT * ( glb%Accel(k)+a1 ) * 0.5d0
1675 
1676 ! Store Accelerations
1677 
1678  glb%Accel(k) = a1
1679  ENDDO
1680  ENDDO
1681 
1682  ELSE
1683 
1684 ! Damping Enabled
1685 
1686  DO j = 1, glb%NumNP
1687  jaux = 3*j
1688  DO k = jaux-2, jaux
1689  dispprev = glb%Disp(k)
1690 
1691 ! Displacement is found by
1692 !
1693 ! -1 ext int damp .
1694 ! { D } = Dt**2 * [ M ] ( F - F - F ) + { D } + Dt*{ D }
1695 ! i+1 i i i i-1/2
1696 
1697  glb%Disp(k)= glb%DT*glb%DT*rnet(k)*glb%xmass(j) + glb%Disp(k) + glb%DT*glb%VeloHalf(k)
1698 
1699 ! Update the velocity [ at time (i+1/2)Dt ]
1700 !
1701 ! . -1
1702 ! { D } = Dt ( { D } - { D } )
1703 ! i+1/2 n + 1 n
1704 
1705 
1706  glb%VeloHalf(k) = ( glb%Disp(k)-dispprev )* glb%DTInv
1707 
1708  ENDDO
1709  ENDDO
1710  ENDIF
1711 
1712 !
1713 ! -- (13) APPLY BOUNDARY CONDITIONS
1714 !
1715 
1716  IF(myid.EQ.0 .AND. debug) THEN
1717  WRITE(6,'(A)') 'Rocfrac: start bc'
1718  ENDIF
1719 
1720  IF(glb%NumNdsBC.NE.0)THEN
1721  CALL bc_enforce(glb%NumNdsBC,glb%NumNP,glb%BCFlag,glb%BCvalue,glb%slope,glb%prop,&
1722  glb%VeloBndry,glb%AccelBndry,glb%VeloHalf,glb%Accel,glb%Disp,glb%DT,rnet,&
1723  glb%xmass,glb%DampEnabled,glb%CurrTime)
1724  ENDIF
1725 
1726 
1727  IF(myid.EQ.0 .AND. debug) THEN
1728  WRITE(6,'(A)') 'Rocfrac: finished bc'
1729  ENDIF
1730 
1731 !-- (1) UPDATE VELOCITY AND DISPLACEMENT VECTORS (structural)
1732 
1733  IF(glb%ALEenabled)THEN
1734  glb%VeloHalf(:) = glb%DT*glb%Accel(:) + glb%VeloHalf(:)
1735  glb%Disp(:) = glb%DT*glb%VeloHalf(:) + glb%Disp(:)
1736  ELSE IF(.NOT.(glb%DampEnabled))THEN
1737 ! glb%VeloHalf(k) = glb%VeloHalf(k) + glb%DT * ( glb%Accel(k)+a1 ) * 0.5d0
1738 
1739 ! Displacement is found by
1740 !
1741 ! -1 ext int damp .
1742 ! { D } = Dt**2 * [ M ] ( F - F - F ) + { D } + Dt*{ D }
1743 ! i+1 i i i i-1/2
1744 
1745 
1746 ! PRINT*,glb%DT*glb%DT*glb%Accel(:)*0.5d0,glb%DT*glb%VeloHalf(:)
1747 
1748  glb%Disp(:) = glb%DT*glb%DT*glb%Accel(:)*0.5d0 + glb%DT*glb%VeloHalf(:) + glb%Disp(:)
1749  ENDIF
1750 
1751  IF(myid.EQ.0 .AND. glb%Verb.gt.1) THEN
1752  WRITE(6,'(A,F12.4)') 'Rocfrac: MAX DISPLACEMENT =', maxval(glb%Disp(:))
1753  ENDIF
1754 !!$ DO i = 1, 3*glb%numnp
1755 !!$ IF(glb%Disp(i).NE.0.d0) PRINT*,i,glb%Disp(i)
1756 !!$ ENDDO
1757 
1758  glb%DispTotal(:) = glb%Disp(:) + glb%DispBar(:)
1759 
1760 !!$ DO i = 1, glb%NumNodeIO
1761 !!$ WRITE(33+i,*) glb%CurrTime, glb%DispTotal(3*glb%NodeIO(i)-2:3*glb%NodeIO(i))
1762 !!$ ENDDO
1763 
1764 ! IF(glb%NumNodeIO.NE.0)THEN
1765 ! write(400,*) glb%CurrTime, glb%Disp(glb%NumNodeIO*3-2)
1766  ! print*,myid,glb%Meshcoor(1:3,glb%NumNodeIO)
1767 ! endif
1768 
1769 ! PRINT*,'Max Displacement', MAXval(glb%Disp)
1770 ! PRINT*,'Min Displacement', minval(glb%Disp)
1771 
1772 ! -- Mass/Volume Conservation
1773 
1774 
1775  IF(myid.EQ.0 .AND. debug) THEN
1776  WRITE(6,'(A)') 'Rocfrac: finished mass volume conservation'
1777  ENDIF
1778  glb%TotalMassSolidp = 0.d0
1779  glb%TotalGeomVolp = 0.d0
1780 
1781  IF(glb%iElType.EQ.4 .OR. glb%iElType.EQ.10 .AND.(.NOT.(glb%NdBasedEl)) )THEN
1782  CALL v3d4_volume(glb%MeshCoor,glb%ElConnVol,glb%MatIdVol,glb%rho, &
1783  glb%NumNP,glb%NumElVol,glb%NumMatVol,glb%Disp,1,glb%NumElVol,&
1784  glb%TotalMassSolidp,glb%TotalGeomVolp,glb%TotalGeomUndefVolp, &
1785  glb%iElType)
1786  ENDIF
1787 
1788  DO i = 1, glb%NumProbesNd
1789  IF(glb%PointOnProc(i))THEN
1790  WRITE(ichr1,'(i4.4)') i
1791  WRITE(ichr2,'(I4.4)') myid
1792 
1793  OPEN(440+i,file='Rocfrac/Rocout/Probe.'//ichr1//'.'//ichr2,position='APPEND')
1794  WRITE(440+i,*) glb%CurrTime, glb%Disp( glb%ProbeNd(i)*3-2), &
1795  glb%Disp( glb%ProbeNd(i)*3-1) , glb%Disp( glb%ProbeNd(i)*3) ! glb%Temperature(glb%ProbeNd(i))
1796  close(440+i)
1797  ENDIF
1798  ENDDO
1799 
1800  ENDDO
1801 
1802  DEALLOCATE(rnet)
1803 
1804 
1805  IF(myid.EQ.0 .AND. debug) THEN
1806  WRITE(6,'(A)') 'Rocfrac: structural boundary conditions start'
1807  ENDIF
1808 
1809  icnt1 = 1
1810  icnt5 = 0
1811  DO i = 1, glb%NumNdsBCcrypt
1812 
1813  ndbcflag = mod(glb%BCFlagCrypt(2,i),100)
1814 ! structural boundary conditions
1815  IF(ndbcflag.GT.0) THEN
1816 
1817  icnt5 = icnt5 + 1
1818 
1819  glb%BCValueGlb(icnt1:icnt1+2) = glb%AccelBndry(icnt5*3-2:icnt5*3)
1820  icnt1 = icnt1 + 3
1821  glb%BCValueGlb(icnt1:icnt1+2) = glb%VeloBndry(icnt5*3-2:icnt5*3)
1822  icnt1 = icnt1 + 3
1823 
1824  ENDIF
1825  ENDDO
1826 
1827  IF(myid.EQ.0 .AND. debug) THEN
1828  WRITE(6,'(A)') 'Rocfrac: END SOLID STEP'
1829  ENDIF
1830 
1831  CALL rocfracinterfacebuff(glb)
1832 
1833  RETURN
1834  END SUBROUTINE rocfracsoln
1835 
1836 !
1837 ! ----------------------------------------------------------------------- RocFracInterfaceBuff
1838 
1839  SUBROUTINE rocfracinterfacebuff(glb)
1840 
1841 !!****f* Rocfrac/Rocfrac/Source/RocFracMain
1842 !!
1843 !! NAME
1844 !! RocFracInterfaceBuff
1845 !!
1846 !! FUNCTION
1847 !! Passes the variables to the fluid code. It transfers the new mesh
1848 !! velocity and displacement to the interface mesh arrays that are
1849 !! registered with RocCom
1850 !!
1851 !!
1852 !!***
1853 
1854  IMPLICIT NONE
1855 
1856  TYPE(rocfrac_global) :: glb
1857 
1858  INTEGER :: iinterfacenode
1859  INTEGER :: k,k1,k2,k3
1860 
1861  INTEGER :: myid,ierr
1862 
1863  CALL mpi_comm_rank(glb%MPI_COMM_ROCFRAC,myid,ierr)
1864 ! WRITE(*,*) 'Number of Nodes on Burning Interface: ',glb%InterfaceSFNumNodes,myid
1865 ! WRITE(*,*) 'Number of Nodes on Non-Burning Interface: ',glb%InterfaceSFnbNumNodes,myid
1866 ! WRITE(*,*) 'Number of Nodes on noninteracting structures surface: ',glb%InterfaceSNumNodes,myid
1867 
1868 ! WRITE(*,*) 'AAAA'
1869 
1870  DO iinterfacenode = 1, glb%InterfaceSFNumNodes
1871  k = glb%MapNodeSF(iinterfacenode)
1872  k3 = 3*k
1873  k2 = k3 - 1
1874  k1 = k3 - 2
1875 
1876 ! Move the interface mesh to match the volume mesh
1877  glb%InterfaceSFNodalCoors(1:3,iinterfacenode) = glb%MeshCoor(1:3,k)
1878 
1879 ! Total Displacement due to structural tractions (not due to mesh motion)
1880 
1881  glb%InterfaceSFTotalNodalDisps(1,iinterfacenode) = glb%Disp(k1)
1882  glb%InterfaceSFTotalNodalDisps(2,iinterfacenode) = glb%Disp(k2)
1883  glb%InterfaceSFTotalNodalDisps(3,iinterfacenode) = glb%Disp(k3)
1884 
1885 !
1886 ! Incremental Total Displacement
1887 !
1888  glb%InterfaceSFNodalDisps(1,iinterfacenode) = glb%Disp(k1) + glb%DispBar(k1)
1889  glb%InterfaceSFNodalDisps(2,iinterfacenode) = glb%Disp(k2) + glb%DispBar(k2)
1890  glb%InterfaceSFNodalDisps(3,iinterfacenode) = glb%Disp(k3) + glb%DispBar(k3)
1891 
1892 ! Passing the velocity, Vs
1893 ! Structural Velocity ! same
1894 
1895  glb%InterfaceSFNodalVels(1,iinterfacenode) = glb%VeloHalf(k1) ! fix: need more terms
1896  glb%InterfaceSFNodalVels(2,iinterfacenode) = glb%VeloHalf(k2) ! fix: need more terms
1897  glb%InterfaceSFNodalVels(3,iinterfacenode) = glb%VeloHalf(k3) ! fix: need more terms
1898 
1899  ENDDO
1900 
1901  DO iinterfacenode = 1, glb%InterfaceSFnbNumNodes
1902  k = glb%MapNodeSFnb(iinterfacenode)
1903  k3 = 3*k
1904  k2 = k3 - 1
1905  k1 = k3 - 2
1906 
1907 ! Move the interface mesh to match the volume mesh
1908  glb%InterfaceSFnbNodalCoors(1:3,iinterfacenode) = glb%MeshCoor(1:3,k)
1909 
1910 ! Total Displacement due to structural tractions (not due to mesh motion)
1911 
1912  glb%InterfaceSFnbTotalNodalDisps(1,iinterfacenode) = glb%Disp(k1)
1913  glb%InterfaceSFnbTotalNodalDisps(2,iinterfacenode) = glb%Disp(k2)
1914  glb%InterfaceSFnbTotalNodalDisps(3,iinterfacenode) = glb%Disp(k3)
1915 
1916 !
1917 ! Incremental Total Displacement
1918 !
1919  glb%InterfaceSFnbNodalDisps(1,iinterfacenode) = glb%Disp(k1) + glb%DispBar(k1)
1920  glb%InterfaceSFnbNodalDisps(2,iinterfacenode) = glb%Disp(k2) + glb%DispBar(k2)
1921  glb%InterfaceSFnbNodalDisps(3,iinterfacenode) = glb%Disp(k3) + glb%DispBar(k3)
1922 
1923 ! Passing the velocity, Vs
1924 ! Structural Velocity ! same
1925 
1926  glb%InterfaceSFnbNodalVels(1,iinterfacenode) = glb%VeloHalf(k1) ! fix: need more terms
1927  glb%InterfaceSFnbNodalVels(2,iinterfacenode) = glb%VeloHalf(k2) ! fix: need more terms
1928  glb%InterfaceSFnbNodalVels(3,iinterfacenode) = glb%VeloHalf(k3) ! fix: need more terms
1929 
1930  ENDDO
1931 
1932  DO iinterfacenode = 1, glb%InterfaceSFnbNumNodes
1933  k = glb%MapNodeSFnb(iinterfacenode)
1934  k3 = 3*k
1935  k2 = k3 - 1
1936  k1 = k3 - 2
1937 
1938 ! Move the interface mesh to match the volume mesh
1939  glb%InterfaceSFnbNodalCoors(1:3,iinterfacenode) = glb%MeshCoor(1:3,k)
1940 
1941 ! Total Displacement due to structural tractions (not due to mesh motion)
1942 
1943  glb%InterfaceSFnbTotalNodalDisps(1,iinterfacenode) = glb%Disp(k1)
1944  glb%InterfaceSFnbTotalNodalDisps(2,iinterfacenode) = glb%Disp(k2)
1945  glb%InterfaceSFnbTotalNodalDisps(3,iinterfacenode) = glb%Disp(k3)
1946 
1947 !
1948 ! Incremental Total Displacement
1949 !
1950  glb%InterfaceSFnbNodalDisps(1,iinterfacenode) = glb%Disp(k1) + glb%DispBar(k1)
1951  glb%InterfaceSFnbNodalDisps(2,iinterfacenode) = glb%Disp(k2) + glb%DispBar(k2)
1952  glb%InterfaceSFnbNodalDisps(3,iinterfacenode) = glb%Disp(k3) + glb%DispBar(k3)
1953 
1954 ! Passing the velocity, Vs
1955 ! Structural Velocity ! same
1956 
1957  glb%InterfaceSFnbNodalVels(1,iinterfacenode) = glb%VeloHalf(k1) ! fix: need more terms
1958  glb%InterfaceSFnbNodalVels(2,iinterfacenode) = glb%VeloHalf(k2) ! fix: need more terms
1959  glb%InterfaceSFnbNodalVels(3,iinterfacenode) = glb%VeloHalf(k3) ! fix: need more terms
1960 
1961  ENDDO
1962 
1963 
1964 
1965  DO iinterfacenode = 1, glb%InterfaceSNumNodes
1966 
1967  k = abs(glb%MapNodeS(iinterfacenode))
1968 ! SolidsNodeNum = ABS(glb%MapNodeS(iInterfaceNode))
1969  IF((k .le. 0).or.(k .gt. glb%NumNP)) THEN
1970  ! IF((k .le. glb%NumNP) .and. (k .gt. 1)) THEN
1971 ! WRITE(*,*) 'BAD NODE MAPPING',iInterfaceNode,k,myid
1972  ELSE
1973 ! Move the interface mesh to match the volume mesh
1974  glb%InterfaceSNodalCoors(1,iinterfacenode) = glb%MeshCoor(1,k)
1975  glb%InterfaceSNodalCoors(2,iinterfacenode) = glb%MeshCoor(2,k)
1976  glb%InterfaceSNodalCoors(3,iinterfacenode) = glb%MeshCoor(3,k)
1977  ENDIF
1978 ! glb%InterfaceSNodalCoors(1,iInterfaceNode) = 0
1979 ! glb%InterfaceSNodalCoors(2,iInterfaceNode) = 0
1980 ! glb%InterfaceSNodalCoors(3,iInterfaceNode) = 0
1981 ! glb%InterfaceSNodalCoors(1:3,iInterfaceNode) = glb%MeshCoor(1:3,k)
1982 
1983  END DO
1984 
1985 
1986  ! WRITE(*,*) 'AAA'
1987  ! CALL MPI_BARRIER(glb%MPI_COMM_ROCFRAC,ierr)
1988  ! STOP
1989  IF(glb%HeatTransSoln .eqv. .true.)THEN
1990  IF(myid.EQ.0 .AND. glb%Verb .gt. 1) THEN
1991  WRITE(*,'(A)') 'Rocfrac: Doing heat transfer'
1992  ENDIF
1993  DO iinterfacenode = 1, glb%InterfaceSFNumNodes
1994  k = glb%MapNodeSF(iinterfacenode)
1995 
1996  glb%InterfaceSFNodalTemp(iinterfacenode) = glb%Temperature(k)
1997 
1998  ENDDO
1999 
2000  ENDIF
2001 
2002 
2003  RETURN
2004 
2005  END SUBROUTINE rocfracinterfacebuff
2006 
2007  SUBROUTINE rocfracupdateinbuff( glb, alpha)
2008  TYPE(rocfrac_global), POINTER :: glb
2009  REAL*8, INTENT(IN) :: alpha
2010  INTEGER :: triconn(1:3)
2011  INTEGER :: i
2012  REAL*8 :: radius, signx,signy,signz
2013  REAL*8 :: prop
2014  INTEGER, SAVE :: icnt
2015  REAL*8 :: tractval
2016  REAL*8 :: slope,mag
2017 
2018 !!$ icnt = icnt + 1
2019 !!$
2020 !!$ IF(mod(icnt,26).eq.0) THEN
2021 !!$ prop = float(icnt/26)*glb%DummyTractVal
2022 !!$ print*,'**PRESSURE=', prop
2023 !!$ write(44,*) glb%CurrTime,prop
2024 !!$ endif
2025 
2026 ! Wind Block
2027 
2028 !!$ DO i = 1, glb%InterfaceSFNumElems
2029 !!$
2030 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(glb%LwrBnd:glb%UppBnd,i)
2031 !!$
2032 !!$ IF(glb%coor(1,glb%MapNodeSF(TriConn(1))) .LT. 0.276d0 .AND. &
2033 !!$ glb%coor(1,glb%MapNodeSF(TriConn(2))) .LT. 0.276d0 .AND. &
2034 !!$ glb%coor(1,glb%MapNodeSF(TriConn(3))) .LT. 0.276d0 ) THEN
2035 !!$
2036 !!$ glb%InterfaceSFElemTract(1,i) = glb%DummyTractVal
2037 !!$ glb%InterfaceSFElemTract(2:3,i) = 0.d0
2038 !!$ ELSE
2039 !!$ glb%InterfaceSFElemTract(:,i) = 0.d0
2040 !!$ ENDIF
2041 !!$
2042 !!$ ENDDO
2043 !!$
2044 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
2045 !!$
2046 !!$! glb%InterfaceSVbar(:,1:glb%InterfaceSNumNodes) = 0.d0
2047 !!$
2048 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
2049 !!$
2050 !!$ DO i = 1, glb%InterfaceSFNumNodes
2051 !!$
2052 !!$ IF(glb%coor(1,glb%MapNodeSF(i)) .LT. 0.276d0) THEN
2053 !!$ glb%InterfaceSFVbar(1,i) = glb%DummyBurnRate
2054 !!$ glb%InterfaceSFVbar(2:3,i) = 0.d0
2055 !!$ ELSE
2056 !!$ glb%InterfaceSFVbar(:,i) = 0.d0
2057 !!$ ENDIF
2058 !!$
2059 !!$ ENDDO
2060 
2061 ! Scaleability
2062 
2063 !!$ DO i = 1, glb%InterfaceSFNumElems
2064 !!$
2065 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(glb%LwrBnd:glb%UppBnd,i)
2066 !!$
2067 !!$ signx= SUM( glb%coor(1,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2068 !!$ signy= SUM( glb%coor(2,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2069 !!$
2070 !!$ radius = SQRT( signx**2 + signy**2)
2071 !!$
2072 !!$ glb%InterfaceSFElemTract(1,i) = signx/radius*glb%DummyTractVal
2073 !!$ glb%InterfaceSFElemTract(2,i) = signy/radius*glb%DummyTractVal
2074 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2075 !!$
2076 !!$ ENDDO
2077 !!$
2078 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
2079 !!$
2080 !!$! glb%InterfaceSVbar(:,1:glb%InterfaceSNumNodes) = 0.d0
2081 !!$
2082 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
2083 !!$
2084 !!$ DO i = 1, glb%InterfaceSFNumNodes
2085 !!$
2086 !!$
2087 !!$ radius = SQRT( glb%coor(1,glb%MapNodeSF(i))**2 + &
2088 !!$ glb%coor(2,glb%MapNodeSF(i))**2)
2089 !!$
2090 !!$ glb%InterfaceSFVbar(1,i) = glb%coor(1,glb%MapNodeSF(i))/radius*glb%DummyBurnRate
2091 !!$ glb%InterfaceSFVbar(2,i) = glb%coor(2,glb%MapNodeSF(i))/radius*glb%DummyBurnRate
2092 !!$ glb%InterfaceSFVbar(3,i) = 0.d0
2093 !!$
2094 !!$ ENDDO
2095 
2096 ! Mass Conservation
2097 !!$
2098 !!$ DO i = 1, glb%InterfaceSFNumElems
2099 !!$
2100 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(glb%LwrBnd:glb%UppBnd,i)
2101 !!$
2102 !!$ glb%InterfaceSFElemTract(1:2,i) = 0.d0
2103 !!$ glb%InterfaceSFElemTract(3,i) = glb%DummyTractVal
2104 !!$
2105 !!$ ENDDO
2106 !!$
2107 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
2108 !!$
2109 !!$! glb%InterfaceSVbar(:,1:glb%InterfaceSNumNodes) = 0.d0
2110 !!$
2111 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
2112 !!$
2113 !!$ DO i = 1, glb%InterfaceSFNumNodes
2114 !!$
2115 !!$ glb%InterfaceSFVbar(3,i) = glb%DummyBurnRate
2116 !!$ glb%InterfaceSFVbar(1:2,i) = 0.d0
2117 !!$
2118 !!$ ENDDO
2119 
2120 
2121 ! Hollow Sphere
2122 
2123 !!$ DO i = 1, glb%InterfaceSFNumElems
2124 !!$
2125 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(glb%LwrBnd:glb%UppBnd,i)
2126 !!$
2127 !!$ signx= SUM( glb%coor(1,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2128 !!$ signy= SUM( glb%coor(2,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2129 !!$ signz= SUM( glb%coor(3,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2130 !!$
2131 !!$ radius = SQRT( signx**2 + signy**2 + signz**2)
2132 !!$
2133 !!$ glb%InterfaceSFElemTract(1,i) = signx/radius*glb%DummyTractVal
2134 !!$ glb%InterfaceSFElemTract(2,i) = signy/radius*glb%DummyTractVal
2135 !!$ glb%InterfaceSFElemTract(3,i) = signz/radius*glb%DummyTractVal
2136 !!$
2137 !!$ ENDDO
2138 !!$
2139 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
2140 !!$
2141 !!$! glb%InterfaceSVbar(:,1:glb%InterfaceSNumNodes) = 0.d0
2142 !!$
2143 !!$! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
2144 !!$
2145 !!$ DO i = 1, glb%InterfaceSFNumNodes
2146 !!$
2147 !!$ radius = SQRT( glb%coor(1,glb%MapNodeSF(i))**2 + &
2148 !!$ glb%coor(2,glb%MapNodeSF(i))**2 + glb%coor(3,glb%MapNodeSF(i))**2)
2149 !!$
2150 !!$ glb%InterfaceSFVbar(1,i) = glb%coor(1,glb%MapNodeSF(i))/radius*glb%DummyBurnRate
2151 !!$ glb%InterfaceSFVbar(2,i) = glb%coor(2,glb%MapNodeSF(i))/radius*glb%DummyBurnRate
2152 !!$ glb%InterfaceSFVbar(3,i) = glb%coor(3,glb%MapNodeSF(i))/radius*glb%DummyBurnRate
2153 !!$
2154 !!$ ENDDO
2155 ! -- Enforce Mesh Velocity BC's of Mesh Motion (Not on F/S interface)
2156 
2157 ! glb%InterfaceSVbar(:,1:glb%InterfaceSNumNodes) = 0.d0
2158 
2159 ! -- Enforce Mesh Velocity BC's of Mesh Motion Given by the fluids.
2160 
2161  DO i = 1, glb%InterfaceSFNumNodes
2162 
2163  glb%InterfaceSFVbar(1,i) = 0.00075 ! glb%DummyBurnRate
2164  glb%InterfaceSFVbar(2,i) = 0.
2165  glb%InterfaceSFVbar(3,i) = 0.
2166 
2167  ENDDO
2168 
2169 
2170 ! 2 Material Beam
2171 
2172 
2173 !!$ DO i = 1, glb%InterfaceSFNumElems
2174 !!$
2175 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(glb%LwrBnd:glb%UppBnd,i)
2176 !!$
2177 !!$ glb%InterfaceSFElemTract(1,i) = 0.d0
2178 !!$ glb%InterfaceSFElemTract(2,i) = -glb%DummyTractVal
2179 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2180 !!$
2181 !!$ ENDDO
2182 !!$
2183 
2184 ! 5 sided cube
2185 
2186 !!$ glb%InterfaceSFElemTract(:,:) = 0.d0
2187 !!$ IF(glb%CurrTime.LE.0.1d0)THEN
2188 !!$ slope = glb%DummyTractVal/0.1d0
2189 !!$
2190 !!$ prop = slope*glb%CurrTime
2191 !!$ print*,' PRESSURE =', prop
2192 !!$ else
2193 !!$ prop = 0.d0
2194 !!$ endif
2195 !!$
2196 !!$ mag = 0.d0
2197 !!$
2198 !!$ DO i = 1, glb%InterfaceSFNumElems
2199 !!$
2200 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(1:3,i)
2201 !!$
2202 !!$! y = 0.9 face
2203 !!$
2204 !!$ IF(glb%coor(2,glb%MapNodeSF(TriConn(1))) .GT. 0.89999d0 .AND. &
2205 !!$ glb%coor(2,glb%MapNodeSF(TriConn(2))) .GT. 0.89999d0 .AND. &
2206 !!$ glb%coor(2,glb%MapNodeSF(TriConn(3))) .GT. 0.89999d0 ) THEN
2207 !!$
2208 !!$ glb%InterfaceSFElemTract(1,i) = 0.d0
2209 !!$ glb%InterfaceSFElemTract(2,i) = -prop
2210 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2211 !!$
2212 !!$! z = 0.9 face
2213 !!$ ELSE IF(glb%coor(3,glb%MapNodeSF(TriConn(1))) .GT. 0.89999d0+mag .AND. &
2214 !!$ glb%coor(3,glb%MapNodeSF(TriConn(2))) .GT. 0.89999d0+mag .AND. &
2215 !!$ glb%coor(3,glb%MapNodeSF(TriConn(3))) .GT. 0.89999d0+mag ) THEN
2216 !!$
2217 !!$ glb%InterfaceSFElemTract(1,i) = 0.d0
2218 !!$ glb%InterfaceSFElemTract(2,i) = 0.d0
2219 !!$ glb%InterfaceSFElemTract(3,i) = -prop
2220 !!$
2221 !!$! z = 0.0 face
2222 !!$ ELSE IF(glb%coor(3,glb%MapNodeSF(TriConn(1))) .LT. 0.00001d0+mag .AND. &
2223 !!$ glb%coor(3,glb%MapNodeSF(TriConn(2))) .LT. 0.00001d0 +mag.AND. &
2224 !!$ glb%coor(3,glb%MapNodeSF(TriConn(3))) .LT. 0.00001d0+mag ) THEN
2225 !!$
2226 !!$ glb%InterfaceSFElemTract(1,i) = 0.d0
2227 !!$ glb%InterfaceSFElemTract(2,i) = 0.d0
2228 !!$ glb%InterfaceSFElemTract(3,i) = prop
2229 !!$
2230 !!$! x = 0.9 face
2231 !!$ ELSE IF(glb%coor(1,glb%MapNodeSF(TriConn(1))) .GT. 0.89999d0+mag .AND. &
2232 !!$ glb%coor(1,glb%MapNodeSF(TriConn(2))) .GT. 0.89999d0+mag .AND. &
2233 !!$ glb%coor(1,glb%MapNodeSF(TriConn(3))) .GT. 0.89999d0+mag ) THEN
2234 !!$
2235 !!$ glb%InterfaceSFElemTract(1,i) = -prop
2236 !!$ glb%InterfaceSFElemTract(2,i) = 0.d0
2237 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2238 !!$
2239 !!$! x = 0.0 face
2240 !!$ ELSE IF(glb%coor(1,glb%MapNodeSF(TriConn(1))) .LT. 0.00001d0+mag .AND. &
2241 !!$ glb%coor(1,glb%MapNodeSF(TriConn(2))) .LT. 0.00001d0+mag .AND. &
2242 !!$ glb%coor(1,glb%MapNodeSF(TriConn(3))) .LT. 0.00001d0+mag ) THEN
2243 !!$
2244 !!$ glb%InterfaceSFElemTract(1,i) = prop
2245 !!$ glb%InterfaceSFElemTract(2,i) = 0.d0
2246 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2247 !!$
2248 !!$ endif
2249 !!$
2250 !!$ ENDDO
2251 
2252 ! Cantilever beam
2253 
2254 !!$ glb%InterfaceSFElemTract(:,:) = 0.d0
2255 !!$
2256 !!$ IF(glb%CurrTime.LE.0.05d0)THEN
2257 !!$
2258 !!$ slope = glb%DummyTractVal/0.05d0
2259 !!$
2260 !!$ prop = slope*glb%CurrTime
2261 !!$ print*,' PRESSURE =', prop
2262 !!$ ELSE
2263 !!$ prop = 0.d0
2264 !!$ ENDIF
2265 !!$
2266 !!$ DO i = 1, glb%InterfaceSFNumElems
2267 !!$
2268 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(1:3,i)
2269 !!$
2270 !!$! y = 0.25 face
2271 !!$
2272 !!$ IF(glb%coor(2,glb%MapNodeSF(TriConn(1))) .GT. 0.2499d0 .AND. &
2273 !!$ glb%coor(2,glb%MapNodeSF(TriConn(2))) .GT. 0.2499d0 .AND. &
2274 !!$ glb%coor(2,glb%MapNodeSF(TriConn(3))) .GT. 0.2499d0 ) THEN
2275 !!$
2276 !!$ glb%InterfaceSFElemTract(1,i) = 0.d0
2277 !!$ glb%InterfaceSFElemTract(2,i) = -prop
2278 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2279 !!$
2280 !!$ ENDIF
2281 !!$ ENDDO
2282 
2283 ! Propellent/Case
2284 
2285 !!$ DO i = 1, glb%InterfaceSFNumElems
2286 !!$
2287 !!$ TriConn(1:3) = glb%InterfaceSFElemConn(glb%LwrBnd:glb%UppBnd,i)
2288 !!$
2289 !!$ signx= SUM( glb%coor(1,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2290 !!$ signy= SUM( glb%coor(2,glb%MapNodeSF(TriConn(1:3))) ) / 3.0d0
2291 !!$
2292 !!$ radius = SQRT( signx**2 + signy**2)
2293 !!$
2294 !!$ TractVal = MIN(1.d0, Prop/25.d0 )*glb%DummyTractVal
2295 !!$
2296 !!$ glb%InterfaceSFElemTract(1,i) = signx/radius*TractVal
2297 !!$ glb%InterfaceSFElemTract(2,i) = signy/radius*TractVal
2298 !!$ glb%InterfaceSFElemTract(3,i) = 0.d0
2299 !!$
2300 !!$ ENDDO
2301 
2302  END SUBROUTINE rocfracupdateinbuff
2303 
2304 END MODULE rocfracmain
2305 
rocstar_rocfraccomm
Definition: ROCSTAR_RocFracComm.f90:53
rocstar_rocfrac
Definition: ROCSTAR_RocFrac.f90:53
vol_elem_mat_matous
subroutine vol_elem_mat_matous(glb)
Definition: vol_elem_mat_matous.f90:53
vol_elem_mat
subroutine vol_elem_mat(e, xnu, ci, cj, numat_vol, Integration)
Definition: vol_elem_mat.f90:53
k
j indices k indices k
Definition: Indexing.h:6
updatestructuralsoln::updatestructural
subroutine updatestructural(glb, NumProcs, Rnet)
Definition: UpdateStructuralSoln.f90:61
tractpressload
subroutine tractpressload(R_ex, numnp, InterfaceSFNumElems, InterfaceSFNumNodes, InterfaceSFElemConn, MapNodeSF, LwrBnd, UppBnd, coor, Disp, MapSFElVolEl, ElConnVol, iElType, NumElVol, TractPress)
Definition: TractPressLoad.f90:53
rocfracmain::rocfracinitialize
subroutine, public rocfracinitialize(glb, InitialTime, MPI_COMM_ROCSTAR, MAN_init, surfIn, volIn, obtain_attr)
Definition: RocFracMain.f90:111
rocfracmain::rocfracinterfacebuff
subroutine rocfracinterfacebuff(glb)
Definition: RocFracMain.f90:1839
feminp
subroutine feminp(glb, myid)
Definition: feminp.f90:53
tractloadpress_hex
subroutine tractloadpress_hex(R_ex, numnp, InterfaceNumElems, InterfaceNumNodes, InterfaceElemConn, MapNode, LwrBnd, UppBnd, coor, InterfaceElemTract)
Definition: TractLoadPress_Hex.f90:53
heatfluxload
subroutine heatfluxload(Rnet, NumNP, InterfaceNumElems, InterfaceNumNodes, InterfaceElemConn, MapNode, LwrBnd, UppBnd, Coor, MapSFElVolEl, ElConnVol, iElType, NumElVol, BCqflux)
Definition: HeatFluxLoad.f90:53
updatestructuralsoln::updatestructuralht
subroutine updatestructuralht(glb, NumProcs, Rnet, RnetHT)
Definition: UpdateStructuralSoln.f90:324
fluidpressload
subroutine fluidpressload(NumNP, R_ex, InterfaceSFNumElems, InterfaceSFNumNodes, InterfaceSFElemConn, MapNodeSF, LwrBnd, UppBnd, coor, Disp, MapSFElVolEl, ElConnVol, iElType, NumElVol, TractPress)
Definition: FluidPressLoad.f90:53
rocstar_rocfracinterp
Definition: ROCSTAR_RocFracInterp.f90:53
i
blockLoc i
Definition: read.cpp:79
rocfracmain::rocfracupdateinbuff
subroutine rocfracupdateinbuff(glb, alpha)
Definition: RocFracMain.f90:2007
principal_stress
subroutine principal_stress(s11, s22, s33, s12, s23, s13, istrgss, NumElVol, SVonMises)
Definition: principal_stress.f90:53
bc_enforce
subroutine bc_enforce(numbound, numnp, id, r, slope, prop, vb, ab, v, a, d, delta, Rnet, xm, DampEnabled, CurrTime)
Definition: bc_enforce.f90:53
volratio
subroutine volratio(n1, n2, n3, n4, AlphaR, coor, numnp, NdMassLump)
Definition: VolRatio.f90:53
rocstar_rocfrac::rocfrac_global
Definition: ROCSTAR_RocFrac.f90:80
max_dt_solid
subroutine max_dt_solid(dt_courant, glb)
Definition: max_dt_solid.f90:53
vector3d::mag
real mag(void) const
Definition: vector3d.h:130
updatemassmatrix
subroutine updatemassmatrix(glb)
Definition: UpdateMassMatrix.f90:53
tractload_hex
subroutine tractload_hex(R_ex, numnp, InterfaceElemTract, InterfaceNumElems, InterfaceNumNodes, InterfaceElemConn, MapNode, LwrBnd, UppBnd, coor)
Definition: TractLoad_Hex.f90:53
updaterbar
subroutine updaterbar(glb, Rnet)
Definition: UpdateR_bar.f90:53
GeoPrim::CLine::slope
CVector & slope()
Definition: GeoPrimitives.H:1041
enhanced_elem_maps_hex
subroutine enhanced_elem_maps_hex(mixed_map, enhanced_map)
Definition: enhanced_elem_maps_hex.f90:53
j
j indices j
Definition: Indexing.h:6
GEM_DomainBoundary::debug
bool debug(bool s=true)
Definition: GEM.H:193
get_mat_stiffness
subroutine get_mat_stiffness(e, dnu, dmat)
Definition: get_mat_stiffness.f90:54
rocfracsubinterface::rocfracinterfaceinitial
subroutine, public rocfracinterfaceinitial(glb, obtain_attr, surfIn)
Definition: RocFracSubInterface.f90:63
Color::alpha
unsigned char alpha() const
Definition: Color.h:75
v3d10_ale
subroutine v3d10_ale(v_bar, a_bar, d, vhalf, Rnet, E, xnu, rho, numnp, numat_vol, numlstet, matlstet, lmlstet, meshcoor, nstart, nend)
Definition: v3d10_ale.f90:65
rocfracmain::rocfracsoln
subroutine, public rocfracsoln(glb, CurrentTime, CurrentTimeStep, MAN_update_inbuff)
Definition: RocFracMain.f90:1138
readsdv
program readsdv
Definition: readsdv.f90:23
rocfracmain::rocfracfinalize
subroutine, public rocfracfinalize(glb)
Definition: RocFracMain.f90:1124
findprobe
subroutine findprobe(glb, myid)
Definition: FindProbe.f90:53
Line_2::a
RT a() const
Definition: Line_2.h:140
v3d4_ale
subroutine v3d4_ale(v_bar, a_bar, d, vhalf, Rnet, E, xnu, rho, numnp, numat_vol, numcstet, matcstet, lmcstet, meshcoor, nstart, nend)
Definition: v3d4_ale.f90:53
v3d4_volume
subroutine v3d4_volume(coor, lmcstet, matcstet, rho, numnp, numcstet, numat_vol, Disp, nstart, nend, TotalMass, TotalGeomVolp, TotalGeomUndefVolp, NumVertx)
Definition: v3d4_volume.f90:53