Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MP/Source/feminp.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 SUBROUTINE feminp(glb, myid)
54 
55  USE rocstar_rocfrac
56 
57 !!****f* Rocfrac/Source/feminp.f90
58 !!
59 !! NAME
60 !! feminp
61 !!
62 !! FUNCTION
63 !!
64 !! READ INPUT INFORMATION (i.e. Analysis Deck File)
65 !!
66 !! INPUTS
67 !! glb -- global array
68 !! myid -- processor id (starting at 0)
69 !!
70 !!****
71 
72  IMPLICIT NONE
73  include 'mpif.h'
74  TYPE(rocfrac_global) :: glb
75 
76 ! -- local
77  CHARACTER*200 :: keywd ! keyword parameter
78 
79  INTEGER :: ios ! io error
80  INTEGER :: ierr ! mpi error
81 
82  INTEGER :: i ! loop counter
83 ! -- global
84  INTEGER :: myid ! processor id
85  INTEGER :: lll ! length of keyword, with no trailing blanks
86 
87  REAL*8, ALLOCATABLE, DIMENSION(:) :: tmp_e, tmp_xnu, tmp_rho, tmp_alpha
88  INTEGER, ALLOCATABLE,DIMENSION(:) :: tmp_isolntype, tmp_matortho
89  REAL*8, ALLOCATABLE, DIMENSION(:) :: tmp_cp, tmp_kappaht
90  INTEGER, ALLOCATABLE,DIMENSION(:) :: tmp_isolntypeht
91 
92  ALLOCATE(tmp_e(1:10), tmp_xnu(1:10), tmp_rho(1:10), tmp_alpha(1:10), tmp_isolntype(1:10) )
93  ALLOCATE(tmp_cp(1:10), tmp_kappaht(1:10),tmp_isolntypeht(1:10) ,tmp_matortho(1:10))
94 !
95 ! -- Set Default values
96 
97  IF(myid.EQ.0) print*,'ROCFRAC:: INSIDE feminp.f90'
98 
99  glb%DummyTractVal = 0.d0
100  glb%DummyBurnRate = 0.d0
101 
102  glb%NumNodeIO = 0
103  glb%ALEenabled = .false. ! ALE disabled by default
104  glb%ipstatic = .false. ! don't subract the intial pressure
105  glb%ReStart = .false.
106  glb%iElType = 4 ! default to 4 node tets
107  glb%IONEWER = .false.
108  glb%DampEnabled = .false.
109  glb%iElIntgratn = 0
110  glb%NdMassLump = 0
111  glb%EnforceTractionS = .false.! To enforce traction where no fluids pressure
112  glb%NumEntries = 0
113  glb%NumMatVol = 0
114  glb%NumMatVolHT = 0
115  glb%NumMatCoh = 0
116  glb%cd_fastest = 0.d0
117  glb%NdBasedEl = .false.
118  glb%UnDefConfig = .false.
119  glb%HeatTransSoln = .false.
120  glb%Temperature0 = 0.
121  glb%ArtificialDamping = .false.
122  glb%EnforceTractionS = .false.
123  glb%EnforceTractionSF = .false.
124  glb%DebondPart = .false.
125  glb%DebondPart_MATOUS = .false.
126  glb%ThermalExpansion = .false.
127  glb%AmplitudeTable = .false.
128  glb%debug_state = .false.
129  glb%NumProbesEl = 0
130  glb%NumProbesNd = 0
131  glb%OverlayExist = .false.
132  glb%DummyFlux = 0.d0
133 
134  glb%NumElOverlay = 0
135  glb%NumNpOverlay = 0
136 !
137 ! -- Open Analysis Deck File
138 
139  OPEN(glb%io_input,file='./Rocfrac/RocfracControl.txt',status='old',iostat=ios)
140  IF(ios.NE.0)THEN
141  IF(myid.EQ.0) print*, 'ROCFRAC:: Unable to find RocfracControl.txt - STOPPING'
142  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
143  stop
144  ENDIF
145 !
146 ! Input deck summary file
147 
148  IF(myid.EQ.0)THEN
149  OPEN(glb%io_sum,file='Rocfrac/Modin/InputSummary.res',status='unknown',iostat=ios)
150  IF(ios.NE.0)THEN ! Try without Modin
151  OPEN(glb%io_sum,file='Rocfrac/InputSummary.res',status='unknown',iostat=ios)
152  END IF
153 
154  IF(ios.NE.0)THEN
155  IF(myid.EQ.0) print*, 'ROCFRAC:: Unable to find InputSummary.res under Rocfrac/Modin/ or Rocfrac/ - STOPPING'
156  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
157  stop
158  ENDIF
159  ENDIF
160 
161 !
162 ! -- Read Analysis Deck File
163 
164  rewind glb%io_input
165 1 READ(glb%io_input,'(A)',iostat=ios) keywd
166 ! print*,keywd
167  IF(ios.LT.0) THEN ! Negative ios means end-of-file
168  print*,' *END parameter not found - STOPPING'
169  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
170  stop
171  ENDIF
172 
173 2 CONTINUE
174 
175 !
176 ! Comment field
177 
178  IF(keywd(1:1).NE.'*')THEN
179  goto 1
180  ENDIF
181 
182  IF(keywd(1:2).EQ.'**')THEN
183  goto 1
184  ENDIF
185 
186  lll = len_trim(keywd)
187 
188  IF(myid.EQ.0) WRITE(glb%io_sum,'(2X,A,A)') keywd(1:lll)
189 
190  IF(keywd(1:4).EQ.'*END') THEN
191  goto 3
192 
193 ! Amplitude
194 
195  ELSE IF(keywd(1:10).EQ.'*AMPLITUDE') THEN
196  CALL amplitude_sub(glb,keywd)
197  goto 1
198 
199  ELSE IF(keywd(1:13).EQ.'*HYPERELASTIC')THEN
200  CALL matmodel_hyperelastic(glb, keywd, &
201  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype)
202  goto 1
203 
204  ELSE IF(keywd(1:14).EQ.'*HEAT TRANSFER')THEN
205  CALL heat_transfer_sub(glb, keywd, tmp_kappaht, tmp_cp, tmp_isolntypeht)
206  goto 1
207 
208  ELSE IF(keywd(1:8).EQ.'*ELASTIC')THEN
209  CALL matmodel_elastic(glb, keywd, &
210  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype , tmp_matortho)
211  goto 1
212 
213  ELSE IF(keywd(1:8).EQ.'*ELEMENT')THEN
214  CALL element_sub(glb,keywd)
215  goto 1
216  ELSE IF(keywd(1:8).EQ.'*DAMPING') THEN
217  READ(glb%io_input,*) glb%KappaDamp
218  glb%DampEnabled = .true.
219  goto 1
220 !
221 ! Prefix for input/output files
222 !
223  ELSE IF(keywd(1:7).EQ.'*PREFIX') THEN
224  CALL prefix_sub(glb)
225  goto 1
226  ELSE IF(keywd(1:4).EQ.'*ALE') THEN
227  CALL ale_sub(glb)
228  goto 1
229  ELSE IF(keywd(1:5).EQ.'*NRUN') THEN
230  CALL nrun_sub(glb)
231  goto 1
232  ELSE IF(keywd(1:8).EQ.'*DYNAMIC') THEN
233  CALL dynamic_sub(glb,keywd)
234  goto 1
235 !*** OBSOLETE REMOVE
236  ELSE IF(keywd(1:7).EQ.'*MATVOL') THEN ! Read volumetric material props.
237  CALL matvol_sub(glb, &
238  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype)
239  goto 1
240  ELSE IF(keywd(1:7).EQ.'*MATCOH') THEN ! Read cohesive material props
241  CALL matcoh_sub(glb)
242  goto 1
243  ELSE IF(keywd(1:7).EQ.'*PLOAD1') THEN
244  CALL pload1_sub(glb)
245  goto 1
246  ELSE IF(keywd(1:5).EQ.'*NODE') THEN
247  CALL nodeio_sub(myid)
248  goto 1
249  ELSE IF(keywd(1:11).EQ.'*DUMMYTRACT') THEN
250  CALL dummytract_sub(glb,keywd)
251  goto 1
252  ELSE IF(keywd(1:10).EQ.'*DUMMYBURN') THEN
253  READ(glb%io_input,*) glb%DummyBurnRate
254  goto 1
255  ELSE IF(keywd(1:10).EQ.'*DUMMYFLUX') THEN
256  READ(glb%io_input,*) glb%DummyFlux
257  goto 1
258  ELSE IF(keywd(1:11).EQ.'*BOUNDARYMM') THEN
259  CALL boundarymm_sub(glb)
260  goto 1
261  ELSE IF(keywd(1:11).EQ.'*BOUNDARYHT') THEN
262  CALL boundaryht_sub(glb)
263  goto 1
264  ELSE IF(keywd(1:9).EQ.'*BOUNDARY') THEN
265  CALL boundary_sub(glb)
266  goto 1
267  ELSE IF(keywd(1:9).EQ.'*IPSTATIC')THEN ! to subract out the initial pressure
268  glb%ipstatic = .true.
269  goto 1
270  ELSE IF(keywd(1:8).EQ.'*IONEWER')THEN ! new gen 2.5 input format
271  glb%IONEWER = .true.
272  goto 1
273  ELSE IF(keywd(1:10).EQ.'*DEFCONFIG')THEN
274  glb%UnDefConfig = .true.
275  goto 1
276 
277  ELSE IF(keywd(1:16).EQ.'*MICROMECHANICAL') THEN ! micromechanical model
278  CALL micromechanical_sub(glb, keywd, &
279  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype)
280  goto 1
281 
282  ELSE IF(keywd(1:11).EQ.'*ARTDAMPING')THEN
283  glb%ArtificialDamping = .true.
284  goto 1
285  ELSE IF(keywd(1:6).EQ.'*PROBE')THEN
286  CALL probe_sub(glb)
287  goto 1
288  ELSE IF(keywd(1:18).EQ.'*INITIAL CONDITION')THEN
289  CALL initialcondition_sub(glb,keywd)
290  goto 1
291 
292  ELSE IF(keywd(1:6).EQ.'*DEBUG')THEN
293  glb%debug_state = .true.
294  goto 1
295 
296  ELSE IF(keywd(1:1).EQ.'*')THEN
297  print*,'ROCFRAC: ERROR'
298  print*,'ROCFRAC: CONTROL DECK OPTION ',trim(keywd),' NOT SUPPORTED'
299  print*,'ROCFRAC: STOPPING'
300  CALL mpi_barrier(glb%MPI_COMM_ROCFRAC,ierr)
301  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
302  stop
303  ELSE
304  goto 1
305  ENDIF
306 3 CONTINUE
307 
308  CLOSE(glb%io_input)
309 !
310 !---- Write analysis summary file
311 
312  ALLOCATE(glb%E (1:glb%NumMatVol) )
313  ALLOCATE(glb%xnu (1:glb%NumMatVol) )
314  ALLOCATE(glb%rho (1:glb%NumMatVol) )
315  ALLOCATE(glb%alpha (1:glb%NumMatVol) )
316  ALLOCATE(glb%xmu (1:glb%NumMatVol) )
317  ALLOCATE(glb%xlambda (1:glb%NumMatVol) )
318  ALLOCATE(glb%xkappa (1:glb%NumMatVol) )
319  ALLOCATE(glb%MatOrtho(1:glb%NumMatVol) )
320 
321  ALLOCATE(glb%iSolnType(1:glb%NumMatVol) )
322 
323 
324  IF(glb%DebondPart_Matous) glb%rho(1) = tmp_rho(1)
325 
326  DO i = 1, glb%NumMatVol
327  glb%E(i) = tmp_e(i)
328  glb%xnu(i) = tmp_xnu(i)
329  glb%rho(i) = tmp_rho(i)
330  glb%alpha(i) = tmp_alpha(i)
331  glb%iSolnType(i) = tmp_isolntype(i)
332 
333  glb%xmu(i) = glb%E(i)/(2.d0*(1.d0+glb%xnu(i)))
334 
335  glb%xlambda(i) = 2.d0*glb%xmu(i)*glb%xnu(i)/(1.d0-2.d0*glb%xnu(i))
336 
337 ! bulk modulus
338  glb%xkappa(i) = glb%xlambda(i) + 2.d0/3.d0*glb%xmu(i)
339 
340  glb%MatOrtho(i) = tmp_matortho(i)
341 
342  ENDDO
343 
344  IF(glb%HeatTransSoln)THEN
345  glb%ThermalDiffusivity = -1.
346 
347  ALLOCATE(glb%iSolnTypeHT(1:glb%NumMatVolHT) )
348  ALLOCATE(glb%KappaHT (1:glb%NumMatVolHT) )
349  ALLOCATE(glb%Cp (1:glb%NumMatVolHT) )
350 
351  DO i = 1, glb%NumMatVolHT
352  glb%KappaHT(i) = tmp_kappaht(i)
353  glb%Cp(i) = tmp_cp(i)
354  glb%iSolnTypeHT(i) = tmp_isolntypeht(i)
355 
356 
357 
358  glb%ThermalDiffusivity = max( glb%ThermalDiffusivity, glb%KappaHT(i)/(glb%rho(i)*glb%Cp(i)) )
359 
360 
361  ENDDO
362  ENDIF
363 
364  DEALLOCATE(tmp_kappaht, tmp_cp, tmp_isolntypeht)
365  DEALLOCATE(tmp_e,tmp_xnu,tmp_rho,tmp_alpha,tmp_isolntype)
366 
367  IF(myid.EQ.0)THEN
368  WRITE(glb%io_sum,50) glb%prefx
369  WRITE(glb%io_sum,60)
370  IF (glb%restart) WRITE(glb%io_sum,65)
371  WRITE(glb%io_sum,80) glb%NumMatVol,glb%NumMatCoh
372  WRITE(glb%io_sum,85) glb%CourantRatio
373  WRITE(glb%io_sum,86)
374  DO i = 1,glb%NumMatCoh
375  WRITE(glb%io_sum,88) glb%deltan(i),glb%deltat(i),glb%SigmaMax(i),glb%TauMax(i),glb%Sinit(i)
376 
377 
378  ENDDO
379  WRITE(glb%io_sum,90)
380  WRITE(glb%io_sum,91)
381  DO i = 1,glb%NumMatVol
382  WRITE(glb%io_sum,92) i,glb%E(i),glb%xnu(i),glb%rho(i),glb%alpha(i)
383  IF(glb%iSolnType(i).EQ.0)THEN
384  WRITE(glb%io_sum,*) ' Material Model = Arruda-Boyce '
385  ELSE IF(glb%iSolnType(i).EQ.-1)THEN
386  WRITE(glb%io_sum,*) ' Material Model = NeoHookean Incompressible '
387  ELSE IF(glb%iSolnType(i).EQ.1)THEN
388  WRITE(glb%io_sum,*) ' Material Model = Elastic, Large Deformation '
389  ELSE IF(glb%iSolnType(i).EQ.2)THEN
390  WRITE(glb%io_sum,*) ' Material Model = Elastic, Small Deformation '
391  ELSE
392  WRITE(glb%io_sum,*) ' Error, Not a valid Material model, = ', glb%iSolnType(i)
393  ENDIF
394 
395  ENDDO
396 
397  IF(glb%HeatTransSoln) WRITE(glb%io_sum,*) 'Heat Transfer Solution'
398  CLOSE(glb%io_sum)
399  ENDIF
400 
401  RETURN
402 
403 !--------------------------------FORMATS--------------------------------
404 !
405 
406 50 FORMAT(//,'DYNAMIC 3D LINEAR ELASTIC FEA',///,'Job id: ',a20,//)
407 60 FORMAT('*** ISOTROPIC ANALYSIS ***',/)
408 65 FORMAT(' *** RESTART ***',/)
409 80 FORMAT(1x,'Number of material types (NUMAT_VOL) =',i12 &
410  /,1x,'Number of material types (NUMAT_COH) =',i12)
411 85 FORMAT(/,1x,'Steps per characteristic length (STEPS) =',e12.4)
412 86 FORMAT(///,1x,'COHESIVE ELEMENT DATA')
413 88 FORMAT(/,1x,'Characteristic lengths for the cohesive law:', &
414  /,1x,' normal (DELTAN) =',e12.4, &
415  /,1x,' tangential (DELTAT) =',e12.4, &
416  /,1x,'Maximum normal stress (GLB%SIGMAMAX) =',e12.4, &
417  /,1x,'Maximum shearing stress (TAUMAX) =',e12.4, &
418  /,1x,'The initial Sthreshold S(init) =',e12.4)
419 90 FORMAT(///,1x,'VOLUMETRIC ELEMENT DATA')
420 91 FORMAT(/,4x,'MATERIAL SETS')
421 92 FORMAT(/,9x,'SET',4x,i4 &
422  /,12x,'E =',10x,e13.5, &
423  /,12x,'Nu =',10x,e13.5, &
424  /,12x,'rho =',10x,e13.5, &
425  /,12x,'alpha =',10x,e13.5)
426 
427 END SUBROUTINE feminp
428 
429 
430 
431 SUBROUTINE prefix_sub(glb)
432 
433 !!****f* Rocfrac/Source/feminp/PREFIX_SUB
434 !!
435 !! NAME
436 !! PREFIX_SUB
437 !!
438 !! FUNCTION
439 !!
440 !! Reads prefix keyword (i.e. Analysis Deck File)
441 !!
442 !! INPUTS
443 !! glb -- global array
444 !!
445 !!****
446 
447  USE rocstar_rocfrac
448 
449  IMPLICIT NONE
450 
451  TYPE(rocfrac_global) :: glb
452 
453  READ(glb%io_input,'(A)') glb%prefx
454  glb%prefx_lngth = len_trim(glb%prefx)
455 
456  RETURN
457 END SUBROUTINE prefix_sub
458 
459 SUBROUTINE ale_sub(glb)
460 
461 !!****f* Rocfrac/Source/feminp/ALE_SUB
462 !!
463 !! NAME
464 !! ALE_SUB
465 !!
466 !! FUNCTION
467 !!
468 !! ALE keyword turns on ALE routines
469 !!
470 !! INPUTS
471 !! glb -- global array
472 !!
473 !!****
474 
475  USE rocstar_rocfrac
476 
477  IMPLICIT NONE
478 
479  TYPE(rocfrac_global) :: glb
480 
481 ! kappa parameter for mesh motion
482 
483  READ(glb%io_input,*) glb%kappa
484  glb%ALEenabled = .true.
485 
486  !print*,'Running with ALE'
487 
488  RETURN
489 END SUBROUTINE ale_sub
490 
491 SUBROUTINE nrun_sub(glb)
492 
493  USE rocstar_rocfrac
494 
495  IMPLICIT NONE
496 
497  TYPE(rocfrac_global) :: glb
498 
499  READ(glb%io_input,*) glb%CourantRatio
500 
501 !RAF Make this backward compatible with old input decks --
502 !RAF You would not want a ratio > 1, right?
503 
504  IF (glb%CourantRatio > 1.0d0) THEN
505  glb%CourantRatio = 1.0d0 / glb%CourantRatio
506  ENDIF
507 
508  RETURN
509 END SUBROUTINE nrun_sub
510 
511 SUBROUTINE dynamic_sub(glb,keywd)
512 
513 !!****f* Rocfrac/Source/feminp/ALE_SUB
514 !!
515 !! NAME
516 !! ALE_SUB
517 !!
518 !! FUNCTION
519 !!
520 !! Courant limit multiplier
521 !!
522 !! INPUTS
523 !! glb -- global array
524 !! keywd -- keywd for control deck
525 !!
526 !!****
527 
528  USE rocstar_rocfrac
529 
530  IMPLICIT NONE
531 
532  CHARACTER(len=200) :: keywd
533  CHARACTER(len=16) :: scalefactor
534  INTEGER :: k1, k2
535 
536 
537  TYPE(rocfrac_global) :: glb
538 
539  ! Find the scale factor, unless the run is implicit
540  IF ( keywd(1:18).NE.'*DYNAMIC, IMPLICIT' ) THEN
541 
542  CALL locchr(keywd,'SCALE FACTOR ',12,8,k1,k2)
543 
544  scalefactor = keywd(k1:k2)
545 
546  CALL rchar(scalefactor,glb%CourantRatio)
547 
548  ELSE
549 
550  glb%CourantRatio = 1.0
551 
552  END IF
553 
554  RETURN
555 END SUBROUTINE dynamic_sub
556 
557 SUBROUTINE dummytract_sub(glb,keywd)
558 
559 !!****f* Rocfrac/Source/feminp/DUMMYTRACT_SUB
560 !!
561 !! NAME
562 !! DUMMYTRACT_SUB
563 !!
564 !! FUNCTION
565 !! Marks if applying the dummy traction to what surface
566 !!
567 !!
568 !! INPUTS
569 !! glb -- global array
570 !! keywd -- keywd for control deck
571 !!
572 !!****
573 
574  USE rocstar_rocfrac
575 
576  IMPLICIT NONE
577 
578  CHARACTER(len=200) :: keywd
579  INTEGER :: k1, k2
580  CHARACTER :: option*16
581 
582  TYPE(rocfrac_global) :: glb
583 
584  CALL locchr(keywd,'INTERFACE ',9,8,k1,k2)
585 
586  option = keywd(k1:k2)
587 
588  IF(option.EQ.'S') glb%EnforceTractionS = .true.
589  IF(option.EQ.'SF') glb%EnforceTractionSF = .true.
590  IF(option.EQ.'ALL')THEN
591  glb%EnforceTractionS = .true.
592  glb%EnforceTractionSF = .true.
593  ENDIF
594 
595  READ(glb%io_input,*) glb%DummyTractVal
596 
597  RETURN
598 END SUBROUTINE dummytract_sub
599 
600 SUBROUTINE matmodel_hyperelastic(glb,keywd, &
601  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype)
602 
603 !!****f* Rocfrac/Source/feminp/MATMODEL_HYPERELASTIC
604 !!
605 !! NAME
606 !! MATMODEL_HYPERELASTIC
607 !!
608 !! FUNCTION
609 !! reads material model type, and the
610 !! hyperelastic material parameters
611 !!
612 !!
613 !! INPUTS
614 !! glb -- global array
615 !! keywd -- keywd for control deck
616 !!
617 !! OUTPUT
618 !! tmp_E -- Young's Modulus
619 !! tmp_xnu -- Possion's ratio
620 !! tmp_rho -- Density
621 !! tmp_alpha -- thermal coefficient of expansion
622 !! tmp_iSolnType -- material type model
623 !!
624 !!****
625 
626  USE rocstar_rocfrac
627 
628  IMPLICIT NONE
629 
630  TYPE(rocfrac_global) :: glb
631 
632  CHARACTER :: keywd*200
633  INTEGER :: key ! 0 = no key word found, 1 = key word found)
634  CHARACTER*26 :: mattype
635  INTEGER :: nummattype
636  INTEGER :: i,ii,ierr
637  REAL*8, DIMENSION(1:10) :: tmp_e, tmp_xnu, tmp_rho, tmp_alpha
638  INTEGER, DIMENSION(1:10) :: tmp_isolntype
639 
640  CALL conchr(keywd,'NEOHOOKINC ',10,13,key)
641  IF(key.EQ.1) mattype = 'NEOHOOKINC'
642 
643  CALL conchr(keywd,'ARRUDA-BOYCE ',12,13,key)
644  IF(key.EQ.1) mattype = 'ARRUDA-BOYCE'
645 
646  SELECT CASE (trim(mattype))
647 
648  CASE ('ARRUDA-BOYCE')
649 
650  READ(glb%io_input,*) nummattype
651 
652  DO i = 1, nummattype
653  glb%NumMatVol = glb%NumMatVol + 1
654 
655  ii = glb%NumMatVol
656 
657  IF(ii.GT.10)THEN
658  print*,'ROCFRAC :: ERROR'
659  print*,'Number of materials GREATER then 10'
660  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
661  ENDIF
662 
663  READ(glb%io_input,*) tmp_e(ii), tmp_xnu(ii), tmp_rho(ii), tmp_alpha(ii)
664 
665  glb%cd_fastest = max( glb%cd_fastest, &
666  sqrt(tmp_e(ii)*(1.d0-tmp_xnu(ii))/tmp_rho(ii)/(1.d0+tmp_xnu(ii))/(1.d0-2.d0*tmp_xnu(ii) )) )
667  tmp_isolntype(ii) = 0
668 
669  ENDDO
670 
671  CASE ('NEOHOOKINC')
672 
673  READ(glb%io_input,*) nummattype
674 
675  DO i = 1, nummattype
676  glb%NumMatVol = glb%NumMatVol + 1
677 
678  ii = glb%NumMatVol
679 
680  IF(ii.GT.10)THEN
681  print*,'ROCFRAC :: ERROR'
682  print*,'Number of materials GREATER then 10'
683  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
684  ENDIF
685 
686  READ(glb%io_input,*) tmp_e(ii), tmp_xnu(ii), tmp_rho(ii), tmp_alpha(ii)
687 
688  glb%cd_fastest = max( glb%cd_fastest, &
689  sqrt(tmp_e(ii)*(1.d0-tmp_xnu(ii))/tmp_rho(ii)/(1.d0+tmp_xnu(ii))/(1.d0-2.d0*tmp_xnu(ii) )) )
690  tmp_isolntype(ii) = -1
691 
692  ENDDO
693 
694  CASE default
695 
696  print*,'ROCFRAC :: ERROR'
697  print*,'ROCFRAC :: *HYPERELASTIC KEYWORD ',trim(mattype), ' NOT FOUND'
698  print*,'ROCFRAC :: STOPPING'
699  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
700 
701 
702  END SELECT
703 
704 
705 END SUBROUTINE matmodel_hyperelastic
706 
707 SUBROUTINE matmodel_elastic(glb,keywd, &
708  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype, tmp_matortho)
709 
710 !!****f* Rocfrac/Source/feminp/MATMODEL_ELASTIC
711 !!
712 !! NAME
713 !! MATMODEL_ELASTIC
714 !!
715 !! FUNCTION
716 !! reads material model type, and the
717 !! elastic material parameters
718 !!
719 !!
720 !! INPUTS
721 !! glb -- global array
722 !! keywd -- keywd for control deck
723 !!
724 !! OUTPUT
725 !! tmp_E -- Young's Modulus
726 !! tmp_xnu -- Possion's ratio
727 !! tmp_rho -- Density
728 !! tmp_alpha -- thermal coefficient of expansion
729 !! tmp_iSolnType -- material type model
730 !!
731 !!****
732 
733  USE rocstar_rocfrac
734 
735  IMPLICIT NONE
736 
737  TYPE(rocfrac_global) :: glb
738 
739  CHARACTER :: keywd*200
740  INTEGER :: key ! 0 = no key word found, 1 = key word found)
741  CHARACTER*26 :: mattype
742  INTEGER :: nummattype
743  INTEGER :: i,ii,iii,ierr
744  REAL*8, DIMENSION(1:10) :: tmp_e, tmp_xnu, tmp_rho, tmp_alpha
745  INTEGER, DIMENSION(1:10) :: tmp_isolntype, tmp_matortho
746 
747  INTEGER :: k1, k2
748  CHARACTER :: nlgeom*16
749  INTEGER :: nlgeomtype
750 
751  LOGICAL :: ortho
752 
753  CALL locchr(keywd,'NLGEOM ',6,8,k1,k2)
754 
755  nlgeom = keywd(k1:k2)
756 
757  nlgeomtype = 1 ! defaults to nl
758  IF(nlgeom.EQ.'NO')THEN
759  nlgeomtype = 2
760  ENDIF
761 
762 
763  READ(glb%io_input,*) nummattype
764 
765  IF ( keywd(1:25) == '*ELASTIC TYPE=ORTHOTROPIC' ) THEN
766  ortho = .true.
767  ALLOCATE(glb%E11o(1:nummattype))
768  ALLOCATE(glb%E22o(1:nummattype))
769  ALLOCATE(glb%E33o(1:nummattype))
770  ALLOCATE(glb%xnu12o(1:nummattype))
771  ALLOCATE(glb%xnu13o(1:nummattype))
772  ALLOCATE(glb%xnu23o(1:nummattype))
773  ALLOCATE(glb%G12o(1:nummattype))
774  ALLOCATE(glb%G13o(1:nummattype))
775  ALLOCATE(glb%G23o(1:nummattype))
776  ALLOCATE(glb%vx1o(1:nummattype))
777  ALLOCATE(glb%vy1o(1:nummattype))
778  ALLOCATE(glb%vz1o(1:nummattype))
779  ALLOCATE(glb%vx2o(1:nummattype))
780  ALLOCATE(glb%vy2o(1:nummattype))
781  ALLOCATE(glb%vz2o(1:nummattype))
782  ALLOCATE(glb%vx3o(1:nummattype))
783  ALLOCATE(glb%vy3o(1:nummattype))
784  ALLOCATE(glb%vz3o(1:nummattype))
785  ELSE
786  ortho = .false.
787  END IF
788 
789  DO i = 1, nummattype
790 
791  glb%NumMatVol = glb%NumMatVol + 1
792 
793  ii = glb%NumMatVol
794 
795  IF(ii.GT.10)THEN
796  print*,'ROCFRAC :: ERROR'
797  print*,'Number of ELASTIC materials GREATER then 10'
798  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
799  ENDIF
800 
801  tmp_isolntype(ii) = nlgeomtype
802 
803  IF ( ortho ) THEN
804 
805  glb%NumMatOrtho = glb%NumMatOrtho + 1
806 
807  iii = glb%NumMatOrtho
808 
809  tmp_matortho(glb%NumMatVol) = iii
810 
811  READ(glb%io_input,*) tmp_rho(ii), glb%E11o(iii), glb%E22o(iii), glb%E33o(iii), &
812  glb%xnu12o(iii), glb%xnu13o(iii), glb%xnu23o(iii), &
813  glb%G12o(iii), glb%G13o(iii), glb%G23o(iii), &
814  glb%vx1o(iii), glb%vy1o(iii), glb%vz1o(iii), &
815  glb%vx2o(iii), glb%vy2o(iii), glb%vz2o(iii), &
816  glb%vx3o(iii), glb%vy3o(iii), glb%vz3o(iii)
817 
818  ! Put negative values into the isotropic variables so the code blows up if they're used
819  tmp_e(ii) = -9999.9
820  tmp_xnu(ii) = -9999.9
821  tmp_alpha(ii) = -9999.9
822 
823  ELSE
824 
825  tmp_matortho(glb%NumMatVol) = 0
826 
827  READ(glb%io_input,*) tmp_e(ii), tmp_xnu(ii), tmp_rho(ii), tmp_alpha(ii)
828 
829  IF(tmp_alpha(ii).NE.0.d0) glb%ThermalExpansion = .true.
830 
831  glb%cd_fastest = max( glb%cd_fastest, &
832  sqrt(tmp_e(ii)*(1.d0-tmp_xnu(ii))/tmp_rho(ii)/(1.d0+tmp_xnu(ii))/(1.d0-2.d0*tmp_xnu(ii) )) )
833 
834  END IF ! Ortho
835 
836  ENDDO
837 
838 END SUBROUTINE matmodel_elastic
839 
840 !*** OBSOLETE REMOVE
841 
842 SUBROUTINE matvol_sub(glb, &
843  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype)
844 
845  USE rocstar_rocfrac
846 
847  IMPLICIT NONE
848 
849  TYPE(rocfrac_global) :: glb
850 
851  INTEGER :: i,ii,ierr ! loop counter
852  integer :: nummattype
853  REAL*8, DIMENSION(1:10) :: tmp_e, tmp_xnu, tmp_rho, tmp_alpha
854  INTEGER, DIMENSION(1:10) :: tmp_isolntype
855 
856  READ(glb%io_input,*) nummattype
857 
858  DO i = 1, nummattype
859  glb%NumMatVol = glb%NumMatVol + 1
860 
861  ii = glb%NumMatVol
862 
863  IF(ii.GT.10)THEN
864  print*,'ROCFRAC :: ERROR'
865  print*,'Number of materials GREATER then 10'
866  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
867  stop
868  ENDIF
869 
870  READ(glb%io_input,*) tmp_e(ii), tmp_xnu(ii), tmp_rho(ii), tmp_alpha(ii),tmp_isolntype(ii)
871 
872 
873  glb%cd_fastest = max( glb%cd_fastest, &
874  sqrt(tmp_e(ii)*(1.d0-tmp_xnu(ii))/tmp_rho(ii)/(1.d0+tmp_xnu(ii))/(1.d0-2.d0*tmp_xnu(ii) )) )
875  ENDDO
876  RETURN
877 END SUBROUTINE matvol_sub
878 
879 SUBROUTINE matcoh_sub(glb)
880 
881  USE rocstar_rocfrac
882 
883  IMPLICIT NONE
884 
885  TYPE(rocfrac_global) :: glb
886 
887  INTEGER :: i ! loop counter
888 
889  READ(glb%io_input,*) glb%NumMatCoh
890 
891  ALLOCATE(glb%deltan(1:glb%NumMatCoh))
892  ALLOCATE(glb%deltat(1:glb%NumMatCoh))
893  ALLOCATE(glb%SigmaMax(1:glb%NumMatCoh))
894  ALLOCATE(glb%TauMax(1:glb%NumMatCoh))
895  ALLOCATE(glb%Sinit (1:glb%NumMatCoh))
896 
897  DO i = 1, glb%NumMatCoh
898  READ(glb%io_input,*) glb%deltan(i),glb%deltat(i),glb%SigmaMax(i),glb%TauMax(i), glb%Sinit(i)
899  ENDDO
900 
901  RETURN
902 END SUBROUTINE matcoh_sub
903 
904 SUBROUTINE pload1_sub(glb)
905 
906  USE rocstar_rocfrac
907 
908  IMPLICIT NONE
909 
910  TYPE(rocfrac_global):: glb
911 
912  INTEGER :: pload1_input ! fix not donei ! loop counter
913 
914  READ(glb%io_input,*) pload1_input
915 
916  RETURN
917 END SUBROUTINE pload1_sub
918 
919 SUBROUTINE nodeio_sub(myid)
920 
921  USE rocstar_rocfrac
922 
923  IMPLICIT NONE
924 
925  TYPE(rocfrac_global) :: glb
926 
927  INTEGER :: i ! loop counter
928  INTEGER :: myid
929 
930  READ(glb%io_input,*) glb%NumNodeIO, glb%NumNodeIOpid
931 
932  IF(glb%NumNodeIOpid.EQ.myid)THEN
933 
934  ALLOCATE(glb%NodeIO(1:glb%NumNodeIO))
935 
936  DO i = 1, glb%NumNodeIO
937  READ(glb%io_input,*) glb%NodeIO
938  ENDDO
939 
940  ELSE
941 
942  glb%NumNodeIO = 0
943 
944  ENDIF
945 
946  RETURN
947 END SUBROUTINE nodeio_sub
948 
949 
950 SUBROUTINE boundary_sub(glb)
951 
952 !!****f* Rocfrac/Source/feminp/BOUNDARY_SUB
953 !!
954 !! NAME
955 !! BOUNDARY_SUB
956 !!
957 !! FUNCTION
958 !! This option is used to prescibe boundary conditions
959 !! at nodes.
960 !!
961 !!
962 !! INPUTS
963 !! glb -- global array
964 !!
965 !!****
966 
967  USE rocstar_rocfrac
968 
969  IMPLICIT NONE
970 
971  TYPE(rocfrac_global) :: glb
972 
973  INTEGER :: i,iaux ! loop counter
974  INTEGER :: numbcflags
975 
976  READ(glb%io_input,*) numbcflags
977 
978  ALLOCATE(glb%bcCond(1:numbcflags))
979 
980  DO i = 1, numbcflags
981  READ(glb%io_input,*) iaux,glb%bcCond(i)%BCtypeX,glb%bcCond(i)%BCtypeY,glb%bcCond(i)%BCtypeZ, &
982  glb%bcCond(i)%BCvalueX,glb%bcCond(i)%BCvalueY,glb%bcCond(i)%BCvalueZ
983  ENDDO
984 
985  RETURN
986 END SUBROUTINE boundary_sub
987 
988 SUBROUTINE boundaryht_sub(glb)
989 
990 !!****f* Rocfrac/Source/feminp/BOUNDARY_SUB
991 !!
992 !! NAME
993 !! BOUNDARY_SUB
994 !!
995 !! FUNCTION
996 !! This option is used to prescibe boundary conditions
997 !! at nodes.
998 !!
999 !!
1000 !! INPUTS
1001 !! glb -- global array
1002 !!
1003 !!****
1004 
1005  USE rocstar_rocfrac
1006 
1007  IMPLICIT NONE
1008 
1009  TYPE(rocfrac_global) :: glb
1010 
1011  INTEGER :: i,iaux ! loop counter
1012  INTEGER :: numbcflagsht
1013 
1014  READ(glb%io_input,*) numbcflagsht
1015 
1016 
1017  ALLOCATE(glb%bcCondHT(1:numbcflagsht))
1018 
1019  DO i = 1, numbcflagsht
1020  READ(glb%io_input,*) iaux,glb%bcCondHT(i)%BCtypeX,glb%bcCondHT(i)%BCtypeY,glb%bcCondHT(i)%BCtypeZ, &
1021  glb%bcCondHT(i)%BCvalueX,glb%bcCondHT(i)%BCvalueY,glb%bcCondHT(i)%BCvalueZ
1022  ENDDO
1023 
1024  RETURN
1025 END SUBROUTINE boundaryht_sub
1026 
1027 SUBROUTINE boundarymm_sub(glb)
1028 
1029 !!****f* Rocfrac/Source/feminp/BOUNDARY_SUB
1030 !!
1031 !! NAME
1032 !! BOUNDARY_SUB
1033 !!
1034 !! FUNCTION
1035 !! This option is used to prescibe mesh motion boundary
1036 !! conditions at nodes.
1037 !!
1038 !! INPUTS
1039 !! glb -- global array
1040 !!
1041 !!****
1042 
1043  USE rocstar_rocfrac
1044 
1045  IMPLICIT NONE
1046 
1047  TYPE(rocfrac_global) :: glb
1048 
1049  INTEGER :: i,iaux ! loop counter
1050  INTEGER :: numbcflagsmm
1051 
1052  READ(glb%io_input,*) numbcflagsmm
1053 
1054  ALLOCATE(glb%bcCondmm(1:numbcflagsmm))
1055 
1056  DO i = 1, numbcflagsmm
1057  READ(glb%io_input,*) iaux,glb%bcCondmm(i)%BCtypeX,glb%bcCondmm(i)%BCtypeY,glb%bcCondmm(i)%BCtypeZ, &
1058  glb%bcCondmm(i)%BCvalueX,glb%bcCondmm(i)%BCvalueY,glb%bcCondmm(i)%BCvalueZ
1059  ENDDO
1060 
1061  RETURN
1062 END SUBROUTINE boundarymm_sub
1063 
1064 SUBROUTINE amplitude_sub(glb,keywd)
1065 
1066  USE rocstar_rocfrac
1067 
1068  IMPLICIT NONE
1069 
1070  TYPE(rocfrac_global) :: glb
1071 
1072  CHARACTER :: keywd*200
1073  INTEGER :: k1, k2
1074 
1075  INTEGER :: i ! loop counter
1076 
1077  CHARACTER :: amptype*16
1078  REAL*8 :: der1
1079  REAL*8 :: slp, intcpt, rise, run
1080  REAL*8, POINTER, DIMENSION(:,:) :: tableval
1081  REAL*8 :: x1, y1, x2, y2
1082 
1083  CALL locchr(keywd,'DEFINITION ',10,10,k1,k2)
1084 
1085  amptype = keywd(k1:k2)
1086 
1087  IF(amptype.EQ.'TABULAR')THEN
1088 
1089  glb%AmplitudeTable = .true.
1090 
1091  READ(glb%io_input,*) glb%NumEntries
1092 
1093 ! Time
1094 ! 1st Derivative
1095 ! 2nd Derivative
1096 
1097  ALLOCATE(tableval(1:2,glb%NumEntries))
1098  DO i = 1, glb%NumEntries
1099  READ(glb%io_input,*) tableval(1:2,i)
1100  ENDDO
1101 
1102  glb%NumEntries = glb%NumEntries - 1
1103 
1104  ALLOCATE( glb%AmpTable(1:3,glb%NumEntries))
1105 
1106  DO i = 1, glb%NumEntries
1107 
1108  x1 = tableval(1,i)
1109  y1 = tableval(2,i)
1110  x2 = tableval(1,i+1)
1111  y2 = tableval(2,i+1)
1112 
1113  rise = y2-y1
1114  run = x2-x1
1115 
1116  slp = (y2-y1)/(x2-x1)
1117 
1118  IF(abs(slp*x1).LT.1.0e-6)THEN
1119  intcpt = y1
1120  ELSE
1121  intcpt = y1/(slp*x1)
1122  ENDIF
1123 
1124  glb%AmpTable(1,i) = tableval(1,i)
1125  glb%AmpTable(2,i) = slp
1126  glb%AmpTable(3,i) = intcpt
1127 
1128 
1129  ENDDO
1130 
1131  DEALLOCATE(tableval)
1132 
1133  ENDIF
1134 
1135  RETURN
1136 END SUBROUTINE amplitude_sub
1137 
1138 SUBROUTINE element_sub(glb,keywd)
1139 
1140 !!****f* Rocfrac/Source/feminp/ELEMENT_SUB
1141 !!
1142 !! NAME
1143 !! ELEMENT_SUB
1144 !!
1145 !! FUNCTION
1146 !! Specifies the element type
1147 !!
1148 !! INPUTS
1149 !! glb -- global array
1150 !!
1151 !!****
1152 
1153  USE rocstar_rocfrac
1154 
1155  IMPLICIT NONE
1156 
1157  TYPE(rocfrac_global) :: glb
1158 
1159  CHARACTER(len=200) :: keywd
1160  INTEGER :: k1, k2
1161 
1162  CHARACTER(len=16) :: eltype
1163 
1164  CALL locchr(keywd,'TYPE ',4,8,k1,k2)
1165 
1166  eltype = keywd(k1:k2)
1167 
1168  SELECT CASE (trim(eltype))
1169  CASE ('V3D4')
1170  glb%iElType = 4
1171  CASE ('V3D4NCC')
1172  glb%iElType = 4
1173  glb%NdMassLump = 1
1174  glb%NdBasedEl = .true.
1175  CASE ('V3D4N')
1176  glb%iElType = 4
1177  glb%NdBasedEl = .true.
1178  CASE ('V3D10R')
1179  glb%iElType = 10
1180  glb%iElIntgratn = 1
1181  CASE ('V3D10')
1182  glb%iElType = 10
1183  CASE ('V3D10BBAR')
1184  glb%iElType = 10
1185  glb%iElIntgratn = 1
1186  CASE ('V3D8ME')
1187  glb%iElType = 8
1188  CASE default
1189  print*,' ERROR:'
1190  print*,'*ELEMENT TYPE NOT FOUND'
1191  stop
1192  END SELECT
1193 
1194  RETURN
1195 END SUBROUTINE element_sub
1196 
1197 SUBROUTINE heat_transfer_sub(glb, keywd, tmp_KappaHT, tmp_Cp, tmp_iSolnTypeHT)
1198 
1199  USE rocstar_rocfrac
1200 
1201  IMPLICIT NONE
1202 
1203  TYPE(rocfrac_global) :: glb
1204 
1205  CHARACTER :: keywd*200
1206  INTEGER :: key ! 0 = no key word found, 1 = key word found)
1207  CHARACTER*26 :: mattype
1208  INTEGER :: nummattypeht
1209  INTEGER :: i,ii,ierr, j
1210  REAL*8, DIMENSION(1:10) :: tmp_kappaht, tmp_cp
1211  INTEGER, DIMENSION(1:10) :: tmp_isolntypeht
1212 
1213  READ(glb%io_input,*) nummattypeht
1214 
1215  DO i = 1, nummattypeht
1216  glb%NumMatVolHT = glb%NumMatVolHT + 1
1217 
1218  ii = glb%NumMatVolHT
1219 
1220  IF(ii.GT.10)THEN
1221  print*,'ROCFRAC :: ERROR'
1222  print*,'Number of materials GREATER then 10'
1223  CALL mpi_finalize(glb%MPI_COMM_ROCFRAC,ierr)
1224  ENDIF
1225 
1226  READ(glb%io_input,*) j, tmp_kappaht(ii),tmp_cp(ii)
1227 
1228  tmp_isolntypeht(ii) = j ! > 0 heat transfer solution wanted, 0 = no heat transfer solution
1229 
1230 !!$! courant limit for temperature problem
1231 !!$!
1232 !!$! 2
1233 !!$! t <= Dx
1234 !!$! ------ (alpha = Kappa / (rho*cp) )
1235 !!$! 3 alpha
1236 !!$!
1237 !!$ AlphaHT = 3.* (KappaHT/RhoCp)
1238 
1239 
1240  ENDDO
1241 
1242  glb%HeatTransSoln = .true.
1243 
1244 END SUBROUTINE heat_transfer_sub
1245 
1246 
1247 
1248 SUBROUTINE probe_sub(glb)
1249 
1250 !!****f* Rocfrac/Source/feminp/BOUNDARY_SUB
1251 !!
1252 !! NAME
1253 !! PROBE_SUB
1254 !!
1255 !! FUNCTION
1256 !! This option is used to mark a node near a
1257 !! specified coordinate
1258 !!
1259 !!
1260 !! INPUTS
1261 !! glb -- global array
1262 !!
1263 !!****
1264 
1265  USE rocstar_rocfrac
1266 
1267  IMPLICIT NONE
1268 
1269  TYPE(rocfrac_global) :: glb
1270 
1271  INTEGER :: i ! loop counter
1272 
1273  READ(glb%io_input,*) glb%NumProbesNd
1274 
1275  ALLOCATE(glb%ProbeCoorNd(1:3,1:glb%NumProbesNd))
1276  ALLOCATE(glb%ProbeNd(1:glb%NumProbesNd))
1277 
1278  DO i = 1, glb%NumProbesNd
1279  READ(glb%io_input,*) glb%ProbeCoorNd(1:3,i)
1280  ENDDO
1281 
1282  RETURN
1283 END SUBROUTINE probe_sub
1284 
1285 
1286 SUBROUTINE micromechanical_sub(glb, keywd, &
1287  tmp_e, tmp_xnu, tmp_rho, tmp_alpha, tmp_isolntype)
1288 
1289 !!****f* Rocfrac/Source/feminp/MICROMECHANICAL
1290 !!
1291 !! NAME
1292 !! MICROMECHANICAL
1293 !!
1294 !! FUNCTION
1295 !! micro mechanical material model
1296 !!
1297 !! INPUTS
1298 !! glb -- global array
1299 !! keywd -- keywd for control deck
1300 !! OUTPUT
1301 !! tmp_E -- Young's Modulus
1302 !! tmp_xnu -- Possion's ratio
1303 !! tmp_rho -- Density
1304 !! tmp_alpha -- thermal coefficient of expansion
1305 !! tmp_iSolnType -- material type model
1306 !!
1307 !!****
1308 
1309  USE rocstar_rocfrac
1310 
1311  IMPLICIT NONE
1312 
1313  TYPE(rocfrac_global) :: glb
1314 
1315  CHARACTER :: keywd*200
1316  INTEGER :: key ! 0 = no key word found, 1 = key word found)
1317  INTEGER :: k1, k2, i
1318  CHARACTER :: modeltype*16
1319  REAL*8, DIMENSION(1:10) :: tmp_e, tmp_xnu, tmp_rho, tmp_alpha
1320  INTEGER, DIMENSION(1:10) :: tmp_isolntype
1321 
1322  INTEGER :: nummatvol_loc
1323 
1324  integer :: nlgeomtype
1325 
1326  nlgeomtype = 1 ! defaults to nl
1327 
1328  CALL locchr(keywd,'MODEL ',5,8,k1,k2)
1329 
1330  IF(k1.GT.0.and.k2.GT.k1)THEN
1331 
1332  modeltype = keywd(k1:k2)
1333 
1334  IF(modeltype.EQ.'HUANG')THEN
1335  ! currently not used
1336  glb%NSTATEV = 1
1337 
1338 
1339  glb%NMATRIX = 3
1340  ALLOCATE(glb%MATRIX(1:glb%NMATRIX))
1341 
1342  READ(glb%io_input,*) glb%MATRIX(1:glb%NMATRIX), tmp_rho(1)
1343 
1344  READ(glb%io_input,*) glb%NPARTICLETYPE
1345 
1346  glb%NPARTICLE = 4
1347 
1348  allocate(glb%PARTICLE(1:glb%NPARTICLE,1:glb%NPARTICLETYPE))
1349  DO i = 1, glb%NPARTICLETYPE
1350  READ(glb%io_input,*) glb%PARTICLE(1:4,i)
1351  enddo
1352 
1353 !----- two particle sizes:
1354 ! the code in the following is limited to two particle sizes, with the
1355 ! same elastic moduli.
1356 
1357  IF (glb%particle(1,1)-glb%particle(1,2).NE.0.d0 .OR. &
1358  glb%particle(2,1)-glb%particle(2,2).NE.0.d0) THEN
1359  print*,'ROCFRAC: the particles do not have the same moduli.'
1360  stop
1361  END IF
1362 
1363  IF (abs( glb%particle(3,1)+glb%particle(3,2) +glb%matrix(3)-1.d0).GT.0.001) THEN
1364  print*,' the volume fractions of particles and matrix do not add up to 1'
1365  stop
1366  END IF
1367 
1368 
1369  glb%NINTERFAC = 3
1370  allocate(glb%INTERFAC(1:glb%NINTERFAC))
1371  READ(glb%io_input,*) glb%INTERFAC(1:3)
1372 
1373  glb%NumMatVol = 1
1374  glb%DebondPart = .true.
1375 
1376 
1377  CALL homogenizedmat( tmp_rho(1), glb%Matrix(1), glb%Matrix(2), &
1378  glb%PARTICLE(1,1), glb%PARTICLE(2,1), glb%MATRIX(3), glb%cd_fastest )
1379 
1380 !!$ glb%cd_fastest = 0.d0
1381 !!$ DO i = 1, 1
1382 !!$ glb%cd_fastest = MAX( glb%cd_fastest, &
1383 !!$ SQRT(glb%MATRIX(1)*(1.d0-glb%MATRIX(2))/tmp_rho(1)/(1.d0+glb%MATRIX(2))/(1.d0-2.d0*glb%MATRIX(2))) )
1384 !!$ ENDDO
1385 !!$ DO i = 1, glb%NPARTICLETYPE
1386 !!$ glb%cd_fastest = MAX( glb%cd_fastest, &
1387 !!$ SQRT(glb%PARTICLE(1,i)*(1.d0-glb%PARTICLE(2,i))/tmp_rho(1)/(1.d0+glb%PARTICLE(2,i))/(1.d0-2.d0*glb%PARTICLE(2,i))) )
1388 !!$ ENDDO
1389 
1390 
1391  ELSE IF(modeltype.EQ.'MATOUS')THEN
1392 
1393  glb%DebondPart_Matous = .true.
1394  glb%NumMatVol = glb%NumMatVol + 1
1395 
1396  READ(glb%io_input,*) glb%NumMatVol_Part
1397 
1398  ALLOCATE(glb%E1(1:glb%NumMatVol_Part) )
1399  ALLOCATE(glb%E2(1:glb%NumMatVol_Part) )
1400  ALLOCATE(glb%E3(1:glb%NumMatVol_Part) )
1401  ALLOCATE(glb%nu12(1:glb%NumMatVol_Part) )
1402  ALLOCATE(glb%nu13(1:glb%NumMatVol_Part) )
1403  ALLOCATE(glb%nu23(1:glb%NumMatVol_Part) )
1404  ALLOCATE(glb%G12(1:glb%NumMatVol_Part) )
1405  ALLOCATE(glb%G13(1:glb%NumMatVol_Part) )
1406  ALLOCATE(glb%G23(1:glb%NumMatVol_Part) )
1407 
1408  DO i = 1, glb%NumMatVol_Part
1409  READ(glb%io_input,*) glb%E1(i), glb%E2(i), glb%E3(i), &
1410  glb%nu12(i), glb%nu13(i), glb%nu23(i), &
1411  glb%G12(i), glb%G13(i), glb%G23(i), tmp_rho(i)
1412 
1413  ENDDO
1414 
1415  READ(glb%io_input,*) glb%alpha1, glb%alpha2, glb%c2, glb%p1, glb%p2, glb%Yin, glb%a_eta, glb%a_zeta
1416 
1417  glb%cm = 1.d0 - glb%c2
1418  glb%cb = glb%c2
1419 
1420  ENDIF
1421 
1422  ELSE
1423  print*,'ERROR in *MICROMECHANICAL keywd'
1424  print*,'MODEL type not found'
1425 
1426  ENDIF
1427 
1428 
1429  tmp_isolntype(1) = nlgeomtype ! fix should not be 1
1430 
1431 
1432 END SUBROUTINE micromechanical_sub
1433 
1434 SUBROUTINE initialcondition_sub(glb,keywd)
1435 
1436 !!****f* Rocfrac/Source/feminp/ELEMENT_SUB
1437 !!
1438 !! NAME
1439 !! ELEMENT_SUB
1440 !!
1441 !! FUNCTION
1442 !! Specifies the element type
1443 !!
1444 !! INPUTS
1445 !! glb -- global array
1446 !!
1447 !!****
1448 
1449  USE rocstar_rocfrac
1450 
1451  IMPLICIT NONE
1452 
1453  TYPE(rocfrac_global) :: glb
1454 
1455  CHARACTER(len=200) :: keywd
1456  INTEGER :: k1, k2
1457 
1458  CALL locchr(keywd,'TYPE ',4,16,k1,k2)
1459 
1460 
1461  IF(keywd(k1:k2).EQ.'TEMPERATURE')THEN
1462  READ(glb%io_input,*) glb%Temperature0
1463  ENDIF
1464 
1465  RETURN
1466 END SUBROUTINE initialcondition_sub
1467 
matmodel_elastic
subroutine matmodel_elastic(glb, keywd, tmp_E, tmp_xnu, tmp_rho, tmp_alpha, tmp_iSolnType)
Definition: feminp.f90:698
element_sub
subroutine element_sub(glb, keywd)
Definition: feminp.f90:1076
dummytract_sub
subroutine dummytract_sub(glb, keywd)
Definition: feminp.f90:548
rocstar_rocfrac
Definition: ROCSTAR_RocFrac.f90:53
micromechanical_sub
subroutine micromechanical_sub(glb, keywd, tmp_E, tmp_xnu, tmp_rho, tmp_alpha, tmp_iSolnType)
Definition: feminp.f90:1224
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
conchr
subroutine conchr(text, varna, lvari, kpos0, key)
Definition: Source/CharKeyWrd.f90:149
Attribute::status
int status() const
Obtain the status of the attribute.
Definition: Attribute.h:240
matmodel_hyperelastic
subroutine matmodel_hyperelastic(glb, keywd, tmp_E, tmp_xnu, tmp_rho, tmp_alpha, tmp_iSolnType)
Definition: feminp.f90:591
sqrt
double sqrt(double d)
Definition: double.h:73
feminp
subroutine feminp(glb, myid)
Definition: feminp.f90:53
matvol_sub
subroutine matvol_sub(glb, tmp_E, tmp_xnu, tmp_rho, tmp_alpha, tmp_iSolnType)
Definition: feminp.f90:780
amplitude_sub
subroutine amplitude_sub(glb, keywd)
Definition: feminp.f90:1002
pload1_sub
subroutine pload1_sub(glb)
Definition: feminp.f90:842
Action::run
virtual void run(double t, double dt, double alpha)
Definition: Action.h:51
matcoh_sub
subroutine matcoh_sub(glb)
Definition: feminp.f90:817
i
blockLoc i
Definition: read.cpp:79
x
void int int REAL * x
Definition: read.cpp:74
initialcondition_sub
subroutine initialcondition_sub(glb, keywd)
Definition: feminp.f90:1373
boundarymm_sub
subroutine boundarymm_sub(glb)
Definition: feminp.f90:965
probe_sub
subroutine probe_sub(glb)
Definition: feminp.f90:1186
rocstar_rocfrac::rocfrac_global
Definition: ROCSTAR_RocFrac.f90:80
rchar
subroutine rchar(char, key)
Definition: Source/CharKeyWrd.f90:285
locchr
subroutine locchr(text, varna, lvari, kpos0, kpos1, kpos2)
Definition: Source/CharKeyWrd.f90:53
boundary_sub
subroutine boundary_sub(glb)
Definition: feminp.f90:888
Aff_transformation_repS2::print
virtual std::ostream & print(std::ostream &os) const
Definition: Aff_transformationS2.h:277
heat_transfer_sub
subroutine heat_transfer_sub(glb, keywd, tmp_KappaHT, tmp_Cp, tmp_iSolnTypeHT)
Definition: feminp.f90:1135
j
j indices j
Definition: Indexing.h:6
homogenizedmat
subroutine homogenizedmat(Density, Em, PRm, Ep, PRp, cm, cd_fastest)
Definition: HomogenizedMat.f90:54
prefix_sub
subroutine prefix_sub(glb)
Definition: feminp.f90:431
nodeio_sub
subroutine nodeio_sub(myid)
Definition: feminp.f90:857
key
static T_Key key
Definition: vinci_lass.c:76
dynamic_sub
subroutine dynamic_sub(glb, keywd)
Definition: feminp.f90:511
nrun_sub
subroutine nrun_sub(glb)
Definition: feminp.f90:491
boundaryht_sub
subroutine boundaryht_sub(glb)
Definition: feminp.f90:926
ale_sub
subroutine ale_sub(glb)
Definition: feminp.f90:459