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