1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5! If new header is supported, change according to following method.    !
6!  1) Increase FSTR_CTRL_HEADER_NUMBER                                 !
7!  2) Add new header name to fstr_ctrl_header_names                    !
8!  3) Describe new function to get parameters from control file        !
9!     in fstr_ctrl.f90                                                 !
10!  4) Describe new subroutine to set values of the parameter           !
11!     in fstr_setup.f90                                                !
12!  5) If initial values are necessary, set the value                   !
13!     in subroutine fstr_setup_init in fstr_setup.f90                  !
14!> This module defined coomon data and basic structures for analysis
15module m_fstr
16  use hecmw
17  use m_common_struct
18  use m_step
19  use m_out
20  use m_timepoint
21  use mMechGauss
22  use mContactDef
23
24  implicit none
25
26  public
27
28  !> CONSTANTS
29  !> general
30  integer(kind=kint), parameter :: kYES = 1
31  integer(kind=kint), parameter :: kNO  = 0
32  integer(kind=kint), parameter :: kON  = 1
33  integer(kind=kint), parameter :: kOFF = 0
34
35  !> solution type (st)
36  integer(kind=kint), parameter :: kstPRECHECK = 0
37  integer(kind=kint), parameter :: kstSTATIC   = 1
38  integer(kind=kint), parameter :: kstEIGEN    = 2
39  integer(kind=kint), parameter :: kstHEAT     = 3
40  integer(kind=kint), parameter :: kstDYNAMIC  = 4
41  !integer(kind=kint), parameter :: kstNLSTATIC  =   5
42  integer(kind=kint), parameter :: kstSTATICEIGEN = 6
43  integer(kind=kint), parameter :: kstNZPROF      = 7
44
45  !> solver method (sm)    !CAUTION : (<=100):indirect, (>100):direct
46  integer(kind=kint), parameter :: ksmCG       = 1
47  integer(kind=kint), parameter :: ksmBiCGSTAB = 2
48  integer(kind=kint), parameter :: ksmGMRES    = 3
49  integer(kind=kint), parameter :: ksmGPBiCG   = 4
50  integer(kind=kint), parameter :: ksmDIRECT   = 101
51
52  !> contact analysis algorithm
53  integer(kind=kint), parameter :: kcaSLagrange = 1
54  integer(kind=kint), parameter :: kcaALagrange = 2
55
56  !> boundary condition file type (bcf)
57  integer(kind=kint), parameter :: kbcfFSTR     =   0  ! BC described in fstr control file (default)
58  integer(kind=kint), parameter :: kbcfNASTRAN  =   1  ! nastran file
59
60  integer(kind=kint), parameter :: kbcInitial   =  1
61  integer(kind=kint), parameter :: kbcTransit   =  2
62
63  !> restart type
64  integer(kind=kint), parameter :: restart_outLast = 1
65  integer(kind=kint), parameter :: restart_outAll  = 2
66
67  !> section control
68  integer(kind=kint), parameter :: kel361FI     =  1
69  integer(kind=kint), parameter :: kel361BBAR   =  2
70  integer(kind=kint), parameter :: kel361IC     =  3
71  integer(kind=kint), parameter :: kel361FBAR   =  4
72
73  integer(kind=kint), parameter :: kFLOADTYPE_NODE = 1
74  integer(kind=kint), parameter :: kFLOADTYPE_SURF = 2
75
76  integer(kind=kint), parameter :: kFLOADCASE_RE = 1
77  integer(kind=kint), parameter :: kFLOADCASE_IM = 2
78
79  !> PARALLEL EXECUTION
80  integer(kind = kint) :: myrank
81  integer(kind = kint) :: nprocs
82
83  !> PARALLEL CONTACT FLAG
84  logical :: paraContactFlag = .false.
85
86  !> FILE NAME
87  character(len=HECMW_FILENAME_LEN) :: cntfilNAME
88  character(len=HECMW_FILENAME_LEN) :: restartfilNAME
89
90  !> FILE HANDLER
91  integer(kind=kint), parameter :: ILOG = 16 ! log
92  integer(kind=kint), parameter :: ISTA = 17 ! status
93  integer(kind=kint), parameter :: IUTB = 18 ! utable
94  integer(kind=kint), parameter :: IMSG = 51 ! message (myrank == 0 only)
95  integer(kind=kint), parameter :: IDBG = 52 ! debug
96  integer(kind=kint), parameter :: IFVS = 53 ! visual.ini file
97  integer(kind=kint), parameter :: INEU = 54 ! neutral file (heat)
98  integer(kind=kint), parameter :: IRESOUT = 100 ! ~110, keeping for result output file
99
100  !> SOLVER CONTROL
101  integer(kind=kint) :: svIarray(100)
102  real(kind=kreal)   :: svRarray(100)
103
104  !> FLAG for ECHO/RESULT/POST
105  integer(kind=kint), pointer :: IECHO
106  integer(kind=kint), pointer :: IRESULT
107  integer(kind=kint), pointer :: IVISUAL
108  integer(kind=kint), pointer :: INEUTRAL  ! flag for femap neutral file
109  integer(kind=kint), pointer :: IRRES     ! flag for restart, read
110  integer(kind=kint), pointer :: IWRES     ! flag for restart, write
111  integer(kind=kint), pointer :: NRRES     ! position of restart read
112  integer(kind=kint), pointer :: NPRINT    ! interval of write
113
114  integer(kind=kint), parameter :: kOPSS_SOLUTION = 1
115  integer(kind=kint), parameter :: kOPSS_MATERIAL = 2
116  integer(kind=kint)            :: OPSSTYPE = kOPSS_SOLUTION ! output stress/strain type
117
118
119  !> REFTEMP
120  real(kind=kreal), pointer :: REF_TEMP
121
122  !> ANALYSIS CONTROL for NLGEOM and HEAT
123  real(kind=kreal)   :: DT    ! /=fstr_param%dtime
124  real(kind=kreal)   :: ETIME ! /=fstr_param%etime
125  integer(kind=kint) :: ITMAX
126  real(kind=kreal)   :: EPS   ! /=fstr_param%eps
127
128  type tInitialCondition
129     character(len=HECMW_FILENAME_LEN)          :: cond_name
130     integer                    :: node_elem                    !< node =0;  element =1
131     integer                    :: grpid
132     integer, pointer           :: intval(:)     => null()      !< if -1, not initialized, otherwise dof number
133     real(kind=kreal), pointer  :: realval(:)    => null()      !< initial value
134  end type
135  type( tInitialCondition ), pointer, save :: g_InitialCnd(:) => null()
136
137  !>  FSTR INNER CONTROL PARAMETERS  (fstrPARAM)
138  type fstr_param
139    integer(kind=kint) :: solution_type !< solution type number
140    integer(kind=kint) :: solver_method !< solver method number
141    logical            :: nlgeom        !< is geometrical nonlinearity considered
142
143    !> STATIC !HEAT
144    integer(kind=kint) :: analysis_n      !< Number of analysis
145    real(kind=kreal), pointer  :: dtime(:) !< (=DT)    STEP_DLTIME
146    real(kind=kreal), pointer  :: etime(:) !< (/=ETIME) STEP_EETIME
147    real(kind=kreal), pointer  :: dtmin(:) !< (=DTMIN) STEP_DELMIN
148    real(kind=kreal), pointer  :: delmax(:)!< (=DTMAX) STEP_DELMAX
149    integer(kind=kint), pointer:: itmax(:) !< (/=ITMAX)
150    real(kind=kreal), pointer  :: eps(:)   !< (/=ESP)
151    real(kind=kreal)           :: ref_temp !< (=REF_TEMP)
152    integer(kind=kint)         :: timepoint_id !< time point ID for heat analysis
153
154    !> output control
155    integer(kind=kint) :: fg_echo       !< output echo   (kYES/kNO) (=IECHO)
156    integer(kind=kint) :: fg_result     !< output result (kYES/kNO) (=IRESULT)
157    integer(kind=kint) :: fg_visual     !< visualization (kYES/kNO) (=IVISUAL)
158
159    !> for heat ...
160    integer(kind=kint) :: fg_neutral    !< write by neutral (=INEUTRAL)
161    integer(kind=kint) :: fg_irres      !< restart read     (=IRRES)
162    integer(kind=kint) :: fg_iwres      !< restart write    (=IWRES)
163    integer(kind=kint) :: nrres         !< NRRES
164    integer(kind=kint) :: nprint        !< NPRINT
165
166    !> index table for global node ID sorting
167    integer(kind=kint) :: n_node        !< hecMESH%n_node
168    integer(kind=kint) :: nn_internal   !< hecMESH%nn_internal
169    integer(kind=kint), pointer :: global_local_ID(:,:)  !> (2:nn_internal) (1,:):global, (2,:):local
170
171    !> for couple analysis
172    integer( kind=kint ) :: fg_couple          !< (default:0)
173    integer( kind=kint ) :: fg_couple_type     !< (default:0)
174    integer( kind=kint ) :: fg_couple_first    !< (default:0)
175    integer( kind=kint ) :: fg_couple_window   !< (default:0)
176
177    !> for restart control
178    integer( kind=kint ) :: restart_out_type   !< output type of restart file
179    integer( kind=kint ) :: restart_version    !< version of restart file
180
181    !> for contact analysis
182    integer( kind=kint ) :: contact_algo       !< contact analysis algorithm number(SLagrange or Alagrange)
183
184    !> for auto increment and cutback
185    type(tParamAutoInc), pointer :: ainc(:)        !< auto increment control
186    type(time_points), pointer :: timepoints(:)  !< time points data
187  end type fstr_param
188
189  !> GLOBAL VARIABLE INITIALIZED IN FSTR_SETUP
190  type( fstr_param ),target :: fstrPR
191
192  !> Data for STATIC ANSLYSIS  (fstrSOLID)
193  type fstr_solid_physic_val
194    real(kind=kreal), pointer :: STRESS(:) => null()   !< nodal stress
195    real(kind=kreal), pointer :: STRAIN(:) => null()   !< nodal strain
196    real(kind=kreal), pointer :: MISES(:) => null()    !< nodal MISES
197
198    real(kind=kreal), pointer :: PSTRESS(:) => null()   !< nodal principal stress
199    real(kind=kreal), pointer :: PSTRAIN(:) => null()   !< nodal principal strain
200    real(kind=kreal), pointer :: PSTRESS_VECT(:,:) => null()  !< nodal principal stress vector
201    real(kind=kreal), pointer :: PSTRAIN_VECT(:,:) => null()  !< nodal principal strain vector
202
203    real(kind=kreal), pointer :: ESTRESS(:) => null()  !< elemental stress
204    real(kind=kreal), pointer :: ESTRAIN(:) => null()  !< elemental strain
205    real(kind=kreal), pointer :: EMISES(:) => null()   !< elemental MISES
206
207    real(kind=kreal), pointer :: EPSTRESS(:) => null()  !< elemental principal stress
208    real(kind=kreal), pointer :: EPSTRAIN(:) => null()  !< elemental principal strain
209    real(kind=kreal), pointer :: EPSTRESS_VECT(:,:) => null()  !< elemental principal stress vector
210    real(kind=kreal), pointer :: EPSTRAIN_VECT(:,:) => null()  !< elemental principal strain vector
211    real(kind=kreal), pointer :: ENQM(:) => null()     !< elemental NQM
212
213
214    type(fstr_solid_physic_val), pointer :: LAYER(:) => null()   !< Laminated Shell's layer (1,2,3,4,5,...)
215    type(fstr_solid_physic_val), pointer :: PLUS => null()   !< for SHELL PLUS
216    type(fstr_solid_physic_val), pointer :: MINUS => null()  !< for SHELL MINUS
217  end type fstr_solid_physic_val
218
219  type fstr_solid
220    integer(kind=kint) :: file_type  ! kbcfFSTR or kbcfNASTRAN
221    integer(kind=kint) :: StaticType ! 1:Total, 2:Updated, 3:Infinitesimal
222    integer(kind=kint) :: nstep_tot
223
224    type(step_info), pointer       :: step_ctrl(:)  =>null()   !< step information
225    type(t_output_ctrl),  pointer  :: output_ctrl(:)=>null()   !< output  information
226
227    !> BOUNDARY
228    integer(kind=kint) :: BOUNDARY_ngrp_tot                    !< Following boundary conditions
229    integer(kind=kint), pointer :: BOUNDARY_ngrp_GRPID  (:)  =>null()
230    integer(kind=kint), pointer :: BOUNDARY_ngrp_ID     (:)  =>null()
231    integer(kind=kint), pointer :: BOUNDARY_ngrp_type   (:)  =>null()
232    integer(kind=kint), pointer :: BOUNDARY_ngrp_amp    (:)  =>null()
233    real(kind=kreal), pointer   :: BOUNDARY_ngrp_val    (:)  =>null()
234    integer(kind=kint) :: BOUNDARY_ngrp_rot                   !< number of rotational boundary conditions
235    integer(kind=kint), pointer :: BOUNDARY_ngrp_rotID     (:) =>null()
236    integer(kind=kint), pointer :: BOUNDARY_ngrp_centerID  (:) =>null()
237
238    !> VELOCITY
239    integer(kind=kint) :: VELOCITY_type
240    integer(kind=kint) :: VELOCITY_ngrp_tot                    !< Following velocity boundary condition
241    integer(kind=kint), pointer :: VELOCITY_ngrp_GRPID  (:)  =>null()
242    integer(kind=kint), pointer :: VELOCITY_ngrp_ID     (:)  =>null()
243    integer(kind=kint), pointer :: VELOCITY_ngrp_type   (:)  =>null()
244    integer(kind=kint), pointer :: VELOCITY_ngrp_amp    (:)  =>null()
245    real(kind=kreal), pointer   :: VELOCITY_ngrp_val    (:)  =>null()
246
247    !> ACCELERATION
248    integer(kind=kint) :: ACCELERATION_type
249    integer(kind=kint) :: ACCELERATION_ngrp_tot                !< Following accelerate boundary condition
250    integer(kind=kint), pointer :: ACCELERATION_ngrp_GRPID  (:)  =>null()
251    integer(kind=kint), pointer :: ACCELERATION_ngrp_ID     (:)  =>null()
252    integer(kind=kint), pointer :: ACCELERATION_ngrp_type   (:)  =>null()
253    integer(kind=kint), pointer :: ACCELERATION_ngrp_amp    (:)  =>null()
254    real(kind=kreal), pointer   :: ACCELERATION_ngrp_val    (:)  =>null()
255
256    !> CLOAD
257    integer(kind=kint) :: CLOAD_ngrp_tot                       !< Following concetrated external load
258    integer(kind=kint), pointer :: CLOAD_ngrp_GRPID     (:) =>null()
259    integer(kind=kint), pointer :: CLOAD_ngrp_ID        (:)
260    integer(kind=kint), pointer :: CLOAD_ngrp_DOF       (:)
261    integer(kind=kint), pointer :: CLOAD_ngrp_amp       (:)
262    real(kind=kreal), pointer   :: CLOAD_ngrp_val       (:)
263    integer(kind=kint) :: CLOAD_ngrp_rot                   !< number of torque load conditions
264    integer(kind=kint), pointer :: CLOAD_ngrp_rotID     (:) =>null()
265    integer(kind=kint), pointer :: CLOAD_ngrp_centerID  (:) =>null()
266
267    !> DLOAD
268    integer(kind=kint) :: DLOAD_ngrp_tot                       !< Following distrubuted external load
269    integer(kind=kint) :: DLOAD_follow
270    integer(kind=kint), pointer :: DLOAD_ngrp_GRPID     (:) =>null()
271    integer(kind=kint), pointer :: DLOAD_ngrp_ID        (:)
272    integer(kind=kint), pointer :: DLOAD_ngrp_LID       (:)
273    integer(kind=kint), pointer :: DLOAD_ngrp_amp       (:)
274    real(kind=kreal), pointer   :: DLOAD_ngrp_params    (:,:)
275
276    !> TEMPERATURE
277    integer(kind=kint) :: TEMP_ngrp_tot                        !< Following tempearture conditions
278    integer(kind=kint) :: TEMP_irres
279    integer(kind=kint) :: TEMP_tstep
280    integer(kind=kint) :: TEMP_interval
281    integer(kind=kint) :: TEMP_rtype      ! type of reading result; 1: step-based; 2: time-based
282    real(kind=kreal)   :: TEMP_FACTOR
283    integer(kind=kint), pointer :: TEMP_ngrp_GRPID     (:) =>null()
284    integer(kind=kint), pointer :: TEMP_ngrp_ID        (:)
285    real(kind=kreal), pointer   :: TEMP_ngrp_val       (:)
286
287    !> SPRING
288    integer(kind=kint) :: SPRING_ngrp_tot                      !< Following spring boundary conditions
289    integer(kind=kint), pointer :: SPRING_ngrp_GRPID    (:) =>null()
290    integer(kind=kint), pointer :: SPRING_ngrp_ID       (:)
291    integer(kind=kint), pointer :: SPRING_ngrp_DOF      (:)
292    integer(kind=kint), pointer :: SPRING_ngrp_amp      (:)
293    real(kind=kreal), pointer   :: SPRING_ngrp_val      (:)
294
295    !> for couple analysis
296    integer( kind=kint ) :: COUPLE_ngrp_tot                   !< Following for coupling analysis
297    integer( kind=kint ),pointer :: COUPLE_ngrp_ID(:)
298
299    !> VALUE
300    integer(kind=kint) :: maxn_gauss
301
302    real(kind=kreal), pointer :: STRESS(:)    !< nodal stress
303    real(kind=kreal), pointer :: STRAIN(:)    !< nodal strain
304    real(kind=kreal), pointer :: MISES(:)     !< nodal MISES
305
306    real(kind=kreal), pointer :: PSTRESS(:)   !< nodal principal stress
307    real(kind=kreal), pointer :: PSTRAIN(:)   !< nodal principal strain
308    real(kind=kreal), pointer :: PSTRESS_VECT(:,:)   !< nodal principal stress vector
309    real(kind=kreal), pointer :: PSTRAIN_VECT(:,:)   !< nodal principal strain vector
310
311    real(kind=kreal), pointer :: ESTRESS(:)   !< elemental stress
312    real(kind=kreal), pointer :: ESTRAIN(:)   !< elemental strain
313    real(kind=kreal), pointer :: EMISES(:)    !< elemental MISES
314
315    real(kind=kreal), pointer :: EPSTRESS(:)   !< elemental principal stress
316    real(kind=kreal), pointer :: EPSTRAIN(:)   !< elemental principal strain
317    real(kind=kreal), pointer :: EPSTRESS_VECT(:,:)   !< elemental principal stress vector
318    real(kind=kreal), pointer :: EPSTRAIN_VECT(:,:)   !< elemental principal strain vector
319
320    real(kind=kreal), pointer :: TNSTRAIN(:)   !< thermal nodal strain
321    real(kind=kreal), pointer :: TESTRAIN(:)   !< thermal elemental strain
322
323    real(kind=kreal), pointer :: YIELD_RATIO(:)    !< yield ratio
324
325    real(kind=kreal), pointer :: ENQM(:)      !< elemental NQM
326    real(kind=kreal), pointer :: REACTION(:)    !< reaction_force
327
328    real(kind=kreal), pointer :: CONT_NFORCE(:)  !< contact normal force for output
329    real(kind=kreal), pointer :: CONT_FRIC(:)    !< contact friction force for output
330    real(kind=kreal), pointer :: CONT_RELVEL(:)  !< contact ralative velocity for output
331    real(kind=kreal), pointer :: CONT_STATE(:)   !< contact state for output
332    integer(kind=kint), pointer :: CONT_SGRP_ID(:) !< contact element surf ids for output
333    real(kind=kreal), pointer :: CONT_AREA(:)    !< contact area
334    real(kind=kreal), pointer :: CONT_NTRAC(:)   !< contact normal traction force for output
335    real(kind=kreal), pointer :: CONT_FTRAC(:)   !< contact friction traction force for output
336
337    type(fstr_solid_physic_val), pointer :: SOLID=>null()     !< for solid physical value stracture
338    type(fstr_solid_physic_val), pointer :: SHELL=>null()     !< for shell physical value stracture
339    type(fstr_solid_physic_val), pointer :: BEAM =>null()     !< for beam physical value stracture
340
341    !> ANALYSIS CONTROL for NLGEOM
342    integer(kind=kint) :: restart_nout   !< output interval of restart file
343    !< (if  .gt.0) restart file write
344    !< (if  .lt.0) restart file read and write
345    integer(kind=kint) :: restart_nin    !< input number of restart
346    type(step_info)    :: step_ctrl_restart  !< step information for restart
347
348    integer(kind=kint) :: max_lyr        !< maximun num of layer
349    integer(kind=kint) :: is_33shell
350    integer(kind=kint) :: is_33beam
351    integer(kind=kint) :: is_heat
352    integer(kind=kint), pointer :: is_rot(:) => null()
353    integer(kind=kint) :: elemopt361
354    real(kind=kreal)   :: FACTOR     (2)      !< factor of incrementation
355    !< 1:time t  2: time t+dt
356    !> for increment control
357    integer(kind=kint) :: NRstat_i(10)     !< statistics of newton iteration (integer)
358    real(kind=kreal)   :: NRstat_r(10)     !< statistics of newton iteration (real)
359    integer(kind=kint) :: AutoINC_stat     !< status of auto-increment control
360    integer(kind=kint) :: CutBack_stat     !< status of cutback control
361
362    real(kind=kreal), pointer :: GL          (:)           !< exnternal force
363    real(kind=kreal), pointer :: EFORCE      (:)           !< exnternal force
364    real(kind=kreal), pointer :: QFORCE      (:)           !< equivalent nodal force
365    real(kind=kreal), pointer :: unode(:)      => null()   !< disp at the beginning of curr step
366    real(kind=kreal), pointer :: dunode(:)     => null()   !< curr total disp
367    real(kind=kreal), pointer :: ddunode(:)    => null()   !< =hecMESH%X, disp icrement
368    real(kind=kreal), pointer :: temperature(:)=> null()   !< =temperature
369    real(kind=kreal), pointer :: temp_bak(:)   => null()
370    real(kind=kreal), pointer :: last_temp(:)  => null()
371
372    type( tElement ), pointer :: elements(:)   =>null()  !< elements information
373    type( tMaterial ),pointer :: materials(:)  =>null()  !< material properties
374    type( tContact ), pointer :: contacts(:)   =>null()  !< contact information
375    integer                   :: n_fix_mpc               !< number mpc conditions user defined
376    real(kind=kreal), pointer :: mpc_const(:)  =>null()  !< bakeup of hecmwST_mpc%mpc_const
377    type(tSection), pointer   :: sections(:)   =>null()  !< definition of setion referred by elements(i)%sectionID
378
379    ! for cutback
380    ! ####################### Notice #######################
381    ! # If you add new variables to store analysis status, #
382    ! # - backup variables with postfix "_bkup" here       #
383    ! # - backup process to module m_fstr_Cutback          #
384    ! # must be added if necessary.                        #
385    ! ######################################################
386    real(kind=kreal), pointer :: unode_bkup(:)     => null() !< disp at the beginning of curr step (backup)
387    real(kind=kreal), pointer :: QFORCE_bkup(:)    => null() !< equivalent nodal force (backup)
388    real(kind=kreal), pointer :: last_temp_bkup(:) => null()
389    type( tElement ), pointer :: elements_bkup(:)  =>null()  !< elements information (backup)
390    type( tContact ), pointer :: contacts_bkup(:)  =>null()  !< contact information (backup)
391  end type fstr_solid
392
393  !> Data for HEAT ANSLYSIS  (fstrHEAT)
394  type fstr_heat
395    !> Crank-Nicolson parameter
396    integer(kind=kint) :: is_steady
397    real(kind=kreal)   :: beta
398    logical :: is_iter_max_limit
399
400    !> TIME CONTROL
401    integer(kind=kint) :: STEPtot
402    integer(kind=kint) :: restart_nout
403    real(kind=kreal), pointer :: STEP_DLTIME(:), STEP_EETIME(:)
404    real(kind=kreal), pointer :: STEP_DELMIN(:), STEP_DELMAX(:)
405    integer(kind=kint) :: timepoint_id
406
407    !> MATERIAL
408    integer(kind=kint) :: MATERIALtot
409    integer(kind=kint), pointer :: RHOtab(:), CPtab(:), CONDtab(:)
410    real(kind=kreal), pointer :: RHO(:,:), RHOtemp (:,:)
411    real(kind=kreal), pointer :: CP(:,:),  CPtemp (:,:)
412    real(kind=kreal), pointer :: COND(:,:),CONDtemp (:,:)
413
414    real(kind=kreal), pointer :: RHOfuncA(:,:),  RHOfuncB(:,:)
415    real(kind=kreal), pointer :: CPfuncA (:,:),  CPfuncB(:,:)
416    real(kind=kreal), pointer :: CONDfuncA (:,:),CONDfuncB(:,:)
417
418    !> AMPLITUDE
419    integer(kind=kint) :: AMPLITUDEtot
420    integer(kind=kint), pointer :: AMPLtab(:)
421    real(kind=kreal), pointer :: AMPL(:,:), AMPLtime (:,:)
422    real(kind=kreal), pointer :: AMPLfuncA(:,:),  AMPLfuncB(:,:)
423
424    !> VALUE
425    real(kind=kreal), pointer :: TEMP0(:)
426    real(kind=kreal), pointer :: TEMPC(:)
427    real(kind=kreal), pointer :: TEMP (:)
428
429    !> FIXTEMP
430    integer(kind=kint) :: T_FIX_tot
431    integer(kind=kint), pointer :: T_FIX_node(:)
432    integer(kind=kint), pointer :: T_FIX_ampl(:)
433    real(kind=kreal), pointer :: T_FIX_val(:)
434
435    !> CFLUX
436    integer(kind=kint) :: Q_NOD_tot
437    integer(kind=kint), pointer :: Q_NOD_node(:)
438    integer(kind=kint), pointer :: Q_NOD_ampl(:)
439    real(kind=kreal), pointer :: Q_NOD_val(:)
440
441    !> DFLUX (not used)
442    integer(kind=kint) :: Q_VOL_tot
443    integer(kind=kint), pointer :: Q_VOL_elem(:)
444    integer(kind=kint), pointer :: Q_VOL_ampl(:)
445    real(kind=kreal), pointer :: Q_VOL_val(:)
446
447    !> DFLUX, !SFLUX
448    integer(kind=kint) :: Q_SUF_tot
449    integer(kind=kint), pointer :: Q_SUF_elem(:)
450    integer(kind=kint), pointer :: Q_SUF_ampl(:)
451    integer(kind=kint), pointer :: Q_SUF_surf(:)
452    real(kind=kreal), pointer :: Q_SUF_val(:)
453
454    !> RADIATE, !SRADIATE
455    integer(kind=kint) :: R_SUF_tot
456    integer(kind=kint), pointer :: R_SUF_elem(:)
457    integer(kind=kint), pointer :: R_SUF_ampl(:,:)
458    integer(kind=kint), pointer :: R_SUF_surf(:)
459    real(kind=kreal), pointer :: R_SUF_val(:,:)
460
461    !> FILM, SFILM
462    integer(kind=kint) :: H_SUF_tot
463    integer(kind=kint), pointer :: H_SUF_elem(:)
464    integer(kind=kint), pointer :: H_SUF_ampl(:,:)
465    integer(kind=kint), pointer :: H_SUF_surf(:)
466    real(kind=kreal), pointer :: H_SUF_val(:,:)
467
468    integer(kind=kint)  :: WL_tot
469    type(tWeldLine), pointer :: weldline(:) => null()
470  end type fstr_heat
471
472  !> Data for DYNAMIC ANSLYSIS  (fstrDYNAMIC)
473  type fstr_dynamic
474    !> ANALYSIS TYPE CONTROL
475    integer(kind=kint) :: idx_eqa       ! implicit or explicit
476    integer(kind=kint) :: idx_resp      ! time history or steady-state harmonic response analysis
477
478    !> TIME CONTROL
479    integer(kind=kint) :: n_step        ! total step number of analysis
480    real(kind=kreal)   :: t_start       ! start time of analysis
481    real(kind=kreal)   :: t_curr        ! current time of analysis
482    real(kind=kreal)   :: t_end         ! end time of analysis
483    real(kind=kreal)   :: t_delta       ! time increment
484    integer(kind=kint) :: restart_nout  ! output interval of restart file
485    ! (if  .gt.0) restart file write
486    ! (if  .lt.0) restart file read and write
487    integer(kind=kint) :: restart_nin  !input number of restart file
488
489    !> Newmark-beta parameter
490    real(kind=kreal)   :: ganma         ! Newmark-beta parameter ganma
491    real(kind=kreal)   :: beta          ! Newmark-beta parameter beta
492
493    !> mass matrix control
494    integer(kind=kint) :: idx_mas       ! mass matrix type
495
496    !> damping control
497    integer(kind=kint) :: idx_dmp      ! damping type
498    real(kind=kreal)   :: ray_m         ! Rayleigh damping parameter Rm
499    real(kind=kreal)   :: ray_k         ! Rayleigh damping parameter Rk
500
501    !> initialization control
502    logical           :: VarInitialize  ! initialization flag
503
504    !> OUTPUT CONTROL
505    integer(kind=kint) :: nout           ! output interval of result
506    integer(kind=kint) :: ngrp_monit     ! node of monitoring result
507    integer(kind=kint) :: nout_monit     ! output interval of result monitoring
508    integer(kind=kint) :: i_step         ! step number
509    integer(kind=kint) :: iout_list(6)   ! 0:not output  1:output
510    ! iout_list(1): displacement
511    ! iout_list(2): velocity
512    ! iout_list(3): acceleration
513    ! iout_list(4): reaction force
514    ! iout_list(5): strain
515    ! iout_list(6): stress
516
517    !> VALUE
518    real(kind=kreal), pointer :: DISP  (:,:)     !> Displacement, U(t+dt), U(t), U(t-dt)
519    real(kind=kreal), pointer :: VEL   (:,:)     !> Velocity
520    real(kind=kreal), pointer :: ACC   (:,:)     !> Acceleration
521
522    real(kind=kreal) :: kineticEnergy
523    real(kind=kreal) :: strainEnergy
524    real(kind=kreal) :: totalEnergy
525
526    !> temporary quantity
527    real(kind=kreal), pointer :: VEC1  (:)
528    real(kind=kreal), pointer :: VEC2  (:)
529    real(kind=kreal), pointer :: VEC3  (:)
530
531    integer(kind=kint) :: dynamic_IW4        =   204
532    integer(kind=kint) :: dynamic_IW5        =   205
533    integer(kind=kint) :: dynamic_IW6        =   206
534    integer(kind=kint) :: dynamic_IW7        =   207
535    integer(kind=kint) :: dynamic_IW8        =   208
536    integer(kind=kint) :: dynamic_IW9        =   209
537    integer(kind=kint) :: dynamic_IW10       =   210
538  end type fstr_dynamic
539
540  type fstr_freqanalysis
541    integer(kind=kint)                :: FLOAD_ngrp_tot
542    integer(kind=kint), pointer       :: FLOAD_ngrp_GRPID(:) => null()
543    integer(kind=kint), pointer       :: FLOAD_ngrp_ID(:)    => null()
544    integer(kind=kint), pointer       :: FLOAD_ngrp_TYPE(:)  => null()
545    integer(kind=kint), pointer       :: FLOAD_ngrp_DOF(:)   => null()
546    real(kind=kreal), pointer         :: FLOAD_ngrp_valre(:) => null()
547    real(kind=kreal), pointer         :: FLOAD_ngrp_valim(:) => null()
548    character(len=HECMW_FILENAME_LEN) :: eigenlog_filename
549    integer(kind=kint)                :: start_mode
550    integer(kind=kint)                :: end_mode
551  end type fstr_freqanalysis
552
553  type fstr_freqanalysis_data
554    integer(kind=kint)        :: numMode
555    integer(kind=kint)        :: numNodeDOF
556    real(kind=kreal), pointer :: eigOmega(:)    => null()
557    real(kind=kreal), pointer :: eigVector(:,:) => null()
558    real(kind=kreal)           :: rayAlpha, rayBeta
559  end type fstr_freqanalysis_data
560
561  !> Package of data used by Lanczos eigenvalue solver
562  type fstr_eigen
563    !> Allocatable array, used or Lanczos eigenvalue analysis
564    integer(kind=kint)  :: nget      ! Solved eigen value number (default:5)
565    integer(kind=kint)  :: maxiter   ! Max. Lcz iterations (default:60)
566    integer(kind=kint)  :: iter      ! Max. Lcz iterations (default:60)
567    real   (kind=kreal) :: sigma     ! 0.0
568    real   (kind=kreal) :: tolerance ! Lcz tolerance (default:1.0e-8)
569    real   (kind=kreal) :: totalmass
570    real   (kind=kreal), pointer :: eigval(:)
571    real   (kind=kreal), pointer :: eigvec(:,:)
572    real   (kind=kreal), pointer :: filter(:)
573    real   (kind=kreal), pointer :: mass(:)
574    real   (kind=kreal), pointer :: effmass(:)
575    real   (kind=kreal), pointer :: partfactor(:)
576    logical :: is_free = .false.
577  end type fstr_eigen
578
579  !> Data for coupling analysis
580  type fstr_couple
581    !> for revocap
582    integer( kind=kint )   :: dof              ! == 3
583    integer( kind=kint )   :: ndof             ! total dof (coupled_node_n*dof)
584    integer( kind=kint )   :: coupled_node_n
585    !> note) following types depend on revocap
586    integer, pointer       :: coupled_node(:)  ! local node id sent from revocap
587    real( kind=8 ),pointer :: trac(:)          ! input  (x,y,z,x,y,z ... )
588    real( kind=8 ),pointer :: disp(:)          ! output (x,y,z,x,y,z ... )
589    real( kind=8 ),pointer :: velo(:)          ! output (x,y,z,x,y,z ... )
590    real( kind=8 ),pointer :: accel(:)         ! output (x,y,z,x,y,z ... )
591    !> for inner use
592    integer( kind=kint ),pointer :: index(:)   ! size:total node num.
593    !>  -1:not relation, >1:index of coupled_node
594  end type fstr_couple
595
596  !> Data for weld line
597  type tWeldLine
598    integer(kind=kint) :: egrpid
599    real( kind=kreal ) :: I
600    real( kind=kreal ) :: U
601    real( kind=kreal ) :: coe
602    real( kind=kreal ) :: v
603    integer(kind=kint) :: xyz
604    real(kind=kreal)   :: n1, n2
605    real(kind=kreal)   :: distol
606    real(kind=kreal)   :: tstart
607  end type tWeldLine
608
609  !> Data for section control
610  type tSection
611    !integer              :: mat_ID
612    !integer              :: iset
613    !integer              :: orien_ID
614    !real(kind=kreal)     :: thickness
615    !integer              :: elemopt341
616    !integer              :: elemopt342
617    !integer              :: elemopt351
618    !integer              :: elemopt352
619    integer              :: elemopt361
620    !integer              :: elemopt362
621  end type tSection
622
623contains
624
625  !> NULL POINTER SETTING TO AVOID RUNTIME ERROR
626  subroutine fstr_nullify_fstr_param( P )
627    implicit none
628    type( fstr_param ) :: P
629
630    nullify( P%dtime )
631    nullify( P%etime )
632    nullify( P%dtmin )
633    nullify( P%delmax )
634    nullify( P%itmax )
635    nullify( P%eps )
636    nullify( P%global_local_ID)
637    nullify( P%timepoints )
638  end subroutine fstr_nullify_fstr_param
639
640  subroutine fstr_nullify_fstr_solid( S )
641    implicit none
642    type( fstr_solid ) :: S
643
644    nullify( S%BOUNDARY_ngrp_ID )
645    nullify( S%BOUNDARY_ngrp_type )
646    nullify( S%BOUNDARY_ngrp_amp )
647    nullify( S%BOUNDARY_ngrp_val)
648    nullify( S%BOUNDARY_ngrp_rotID )
649    nullify( S%BOUNDARY_ngrp_centerID )
650    nullify( S%CLOAD_ngrp_ID )
651    nullify( S%CLOAD_ngrp_DOF )
652    nullify( S%CLOAD_ngrp_amp )
653    nullify( S%CLOAD_ngrp_rotID )
654    nullify( S%CLOAD_ngrp_centerID )
655    nullify( S%CLOAD_ngrp_val )
656    nullify( S%DLOAD_ngrp_ID )
657    nullify( S%DLOAD_ngrp_LID )
658    nullify( S%DLOAD_ngrp_amp )
659    nullify( S%DLOAD_ngrp_params )
660    nullify( S%TEMP_ngrp_ID )
661    nullify( S%TEMP_ngrp_val )
662    nullify( S%SPRING_ngrp_ID )
663    nullify( S%SPRING_ngrp_DOF )
664    nullify( S%SPRING_ngrp_amp )
665    nullify( S%SPRING_ngrp_val )
666    nullify( S%STRESS )
667    nullify( S%STRAIN )
668    nullify( S%MISES )
669    nullify( S%PSTRESS )
670    nullify( S%PSTRAIN )
671    nullify( S%PSTRESS_VECT )
672    nullify( S%PSTRAIN_VECT )
673    nullify( S%REACTION )
674    nullify( S%ESTRESS )
675    nullify( S%ESTRAIN )
676    nullify( S%EMISES )
677    nullify( S%EPSTRESS )
678    nullify( S%EPSTRAIN )
679    nullify( S%EPSTRESS_VECT )
680    nullify( S%EPSTRAIN_VECT )
681    nullify( S%ENQM )
682    nullify( S%GL          )
683    nullify( S%QFORCE      )
684    nullify( S%VELOCITY_ngrp_ID )
685    nullify( S%VELOCITY_ngrp_type )
686    nullify( S%VELOCITY_ngrp_amp )
687    nullify( S%VELOCITY_ngrp_val )
688    nullify( S%ACCELERATION_ngrp_ID )
689    nullify( S%ACCELERATION_ngrp_type )
690    nullify( S%ACCELERATION_ngrp_amp )
691    nullify( S%ACCELERATION_ngrp_val )
692    nullify( S%COUPLE_ngrp_ID )
693  end subroutine fstr_nullify_fstr_solid
694
695  subroutine fstr_nullify_fstr_heat( H )
696    implicit none
697    type( fstr_heat ) :: H
698
699    nullify( H%STEP_DLTIME )
700    nullify( H%STEP_EETIME )
701    nullify( H%STEP_DELMIN )
702    nullify( H%STEP_DELMAX )
703    nullify( H%RHO )
704    nullify( H%RHOtemp  )
705    nullify( H%CP )
706    nullify( H%CPtemp  )
707    nullify( H%COND )
708    nullify( H%CONDtemp  )
709    nullify( H%RHOtab )
710    nullify( H%CPtab )
711    nullify( H%CONDtab )
712    nullify( H%RHOfuncA )
713    nullify( H%RHOfuncB )
714    nullify( H%CPfuncA  )
715    nullify( H%CPfuncB )
716    nullify( H%CONDfuncA  )
717    nullify( H%CONDfuncB )
718    nullify( H%AMPL )
719    nullify( H%AMPLtime  )
720    nullify( H%AMPLtab )
721    nullify( H%AMPLfuncA )
722    nullify( H%AMPLfuncB )
723    nullify( H%TEMP0 )
724    nullify( H%TEMPC )
725    nullify( H%TEMP  )
726    nullify( H%T_FIX_node )
727    nullify( H%T_FIX_ampl )
728    nullify( H%T_FIX_val )
729    nullify( H%Q_NOD_node )
730    nullify( H%Q_NOD_ampl )
731    nullify( H%Q_NOD_val )
732    nullify( H%Q_VOL_elem )
733    nullify( H%Q_VOL_ampl )
734    nullify( H%Q_VOL_val )
735    nullify( H%Q_SUF_elem )
736    nullify( H%Q_SUF_ampl )
737    nullify( H%Q_SUF_surf )
738    nullify( H%Q_SUF_val )
739    nullify( H%R_SUF_elem )
740    nullify( H%R_SUF_ampl )
741    nullify( H%R_SUF_surf )
742    nullify( H%R_SUF_val )
743    nullify( H%H_SUF_elem )
744    nullify( H%H_SUF_ampl )
745    nullify( H%H_SUF_surf )
746    nullify( H%H_SUF_val )
747  end subroutine fstr_nullify_fstr_heat
748
749  subroutine fstr_nullify_fstr_dynamic( DY )
750    implicit none
751    type( fstr_dynamic ) :: DY
752
753    nullify( DY%DISP )
754    nullify( DY%VEL  )
755    nullify( DY%ACC  )
756    nullify( DY%VEC1 )
757    nullify( DY%VEC2 )
758    nullify( DY%VEC3 )
759  end subroutine fstr_nullify_fstr_dynamic
760
761  subroutine fstr_nullify_fstr_freqanalysis( f )
762    implicit none
763    type( fstr_freqanalysis ), intent(inout) :: f
764
765    f%FLOAD_ngrp_tot = 0
766    nullify( f%FLOAD_ngrp_GRPID )
767    nullify( f%FLOAD_ngrp_ID )
768    nullify( f%FLOAD_ngrp_TYPE )
769    nullify( f%FLOAD_ngrp_DOF )
770    nullify( f%FLOAD_ngrp_valre )
771    nullify( f%FLOAD_ngrp_valim )
772  end subroutine fstr_nullify_fstr_freqanalysis
773
774  subroutine fstr_nullify_fstr_eigen( E )
775    implicit none
776    type( fstr_eigen ) :: E
777
778    nullify( E%mass )
779  end subroutine fstr_nullify_fstr_eigen
780
781  subroutine fstr_nullify_fstr_couple( C )
782    implicit none
783    type( fstr_couple ) :: C
784
785    nullify( C%coupled_node )
786    nullify( C%trac )
787    nullify( C%velo )
788    nullify( C%index )
789  end subroutine fstr_nullify_fstr_couple
790
791  !> Initializer of structure hecmwST_matrix
792  subroutine fstr_mat_init( hecMAT )
793    implicit none
794    type(hecmwST_matrix) :: hecMAT
795
796    hecMAT%Iarray(1) =  100    ! = nier
797    hecMAT%Iarray(2) =    1    ! = method
798    hecMAT%Iarray(3) =    1    ! = precond
799    hecMAT%Iarray(4) =    0    ! = nset
800    hecMAT%Iarray(5) =    1    ! = iterpremax
801    hecMAT%Iarray(6) =   10    ! = nrest
802    hecMAT%Iarray(7) =    0    ! = scaling
803    hecMAT%Iarray(21)=  kNO    ! = iterlog
804    hecMAT%Iarray(22)=  kNO    ! = timelog
805    hecMAT%Iarray(31)=    0    ! = dumptype
806    hecMAT%Iarray(32)=    0    ! = dumpexit
807    hecMAT%Iarray(33)=    0    ! = usejad
808    hecMAT%Iarray(34)=   10    ! = ncolor_in
809    hecMAT%Iarray(13)=    0    ! = mpc_method
810    hecMAT%Iarray(14)=    0    ! = estcond
811    hecMAT%Iarray(35)=    3    ! = maxrecycle_precond
812    hecMAT%Iarray(41)=    0    ! = solver_opt1
813    hecMAT%Iarray(42)=    0    ! = solver_opt2
814    hecMAT%Iarray(43)=    0    ! = solver_opt3
815    hecMAT%Iarray(44)=    0    ! = solver_opt4
816    hecMAT%Iarray(45)=    0    ! = solver_opt5
817    hecMAT%Iarray(46)=    0    ! = solver_opt6
818
819    hecMAT%Rarray(1) =  1.0e-8 ! = resid
820    hecMAT%Rarray(2) =  1.0    ! = sigma_diag
821    hecMAT%Rarray(3) =  0.0    ! = sigma
822    hecMAT%Rarray(4) =  0.1    ! = thresh
823    hecMAT%Rarray(5) =  0.1    ! = filter
824    hecMAT%Rarray(11)=  1.0e+4 ! = penalty
825
826    hecMAT%Iarray(96) =   0    ! nrecycle_precond
827    hecMAT%Iarray(97) = kYES   ! flag_numfact
828    hecMAT%Iarray(98) = kYES   ! flag_symbfact
829    hecMAT%Iarray(99) = kYES   ! indirect method
830  end subroutine fstr_mat_init
831
832  subroutine hecMAT_init( hecMAT )
833    implicit none
834    type( hecmwST_matrix ) :: hecMAT
835    integer ::  ndof, nn, ierror
836    ndof = hecMAT%NDOF
837    nn = ndof*ndof
838    allocate (hecMAT%AL(nn*hecMAT%NPL)        ,stat=ierror )
839    if( ierror /= 0 ) then
840      write(*,*) "##ERROR : not enough memory"
841      write(idbg,*) 'stop due to allocation error'
842      call flush(idbg)
843      call hecmw_abort( hecmw_comm_get_comm() )
844    end if
845    allocate (hecMAT%AU(nn*hecMAT%NPU)        ,stat=ierror )
846    if( ierror /= 0 ) then
847      write(*,*) "##ERROR : not enough memory"
848      write(idbg,*) 'stop due to allocation error'
849      call flush(idbg)
850      call hecmw_abort( hecmw_comm_get_comm() )
851    end if
852    allocate (hecMAT%B(ndof*hecMAT%NP)          ,stat=ierror )
853    if( ierror /= 0 ) then
854      write(*,*) "##ERROR : not enough memory"
855      write(idbg,*) 'stop due to allocation error'
856      call flush(idbg)
857      call hecmw_abort( hecmw_comm_get_comm() )
858    end if
859    hecMAT%B(:)=0.d0
860    allocate (hecMAT%D(nn*hecMAT%NP)          ,stat=ierror )
861    if( ierror /= 0 ) then
862      write(*,*) "##ERROR : not enough memory"
863      write(idbg,*) 'stop due to allocation error'
864      call flush(idbg)
865      call hecmw_abort( hecmw_comm_get_comm() )
866    end if
867    allocate (hecMAT%X(ndof*hecMAT%NP)          ,stat=ierror )
868    if( ierror /= 0 ) then
869      write(*,*) "##ERROR : not enough memory"
870      write(idbg,*) 'stop due to allocation error'
871      call flush(idbg)
872      call hecmw_abort( hecmw_comm_get_comm() )
873    end if
874    allocate (hecMAT%ALU(nn*hecMAT%N)         ,stat=ierror )
875    if( ierror /= 0 ) then
876      write(*,*) "##ERROR : not enough memory"
877      write(idbg,*) 'stop due to allocation error'
878      call flush(idbg)
879      call hecmw_abort( hecmw_comm_get_comm() )
880    endif
881    call hecmw_cmat_init( hecMAT%cmat )
882    hecMAT%D  = 0.0d0
883    hecMAT%AL = 0.0d0
884    hecMAT%AU = 0.0d0
885    hecMAT%B  = 0.0d0
886    hecMAT%X  = 0.0d0
887    hecMAT%ALU = 0.0d0
888  end subroutine hecMAT_init
889
890  subroutine hecMAT_finalize( hecMAT )
891    implicit none
892    type( hecmwST_matrix ) :: hecMAT
893    integer ::  ndof, nn, ierror
894    ndof = hecMAT%NDOF
895    nn = ndof*ndof
896    if( associated(hecMAT%AL) ) then
897      deallocate(hecMAT%AL                  ,stat=ierror)
898      if( ierror /= 0 ) then
899        write(idbg,*) 'stop due to deallocation error'
900        call flush(idbg)
901        call hecmw_abort( hecmw_comm_get_comm())
902      end if
903    endif
904    if( associated(hecMAT%AU) ) then
905      deallocate(hecMAT%AU                  ,stat=ierror)
906      if( ierror /= 0 ) then
907        write(idbg,*) 'stop due to deallocation error'
908        call flush(idbg)
909        call hecmw_abort( hecmw_comm_get_comm())
910      end if
911    endif
912    if( associated(hecMAT%B) ) then
913      deallocate(hecMAT%B                   ,stat=ierror)
914      if( ierror /= 0 ) then
915        write(idbg,*) 'stop due to deallocation error'
916        call flush(idbg)
917        call hecmw_abort( hecmw_comm_get_comm())
918      end if
919    endif
920    if( associated(hecMAT%D) ) then
921      deallocate(hecMAT%D                   ,stat=ierror)
922      if( ierror /= 0 ) then
923        write(idbg,*) 'stop due to deallocation error'
924        call flush(idbg)
925        call hecmw_abort( hecmw_comm_get_comm())
926      end if
927    endif
928    if( associated(HECMAT%X) ) then
929      deallocate(hecMAT%X                   ,stat=ierror)
930      if( ierror /= 0 ) then
931        write(idbg,*) 'stop due to deallocation error'
932        call flush(idbg)
933        call hecmw_abort( hecmw_comm_get_comm())
934      end if
935    endif
936    if( associated(hecMAT%ALU) ) then
937      deallocate(hecMAT%ALU                 ,stat=ierror)
938      if( ierror /= 0 ) then
939        write(idbg,*) 'stop due to deallocation error'
940        call flush(idbg)
941        call hecmw_abort( hecmw_comm_get_comm())
942      end if
943    endif
944    call hecmw_cmat_finalize(hecmAT%cmat)
945  end subroutine hecMAT_finalize
946
947  !> Initializer of structure fstr_param
948  subroutine fstr_param_init( fstrPARAM, hecMESH )
949    implicit none
950    type(fstr_param) :: fstrPARAM
951    type(hecmwST_local_mesh) :: hecMESH
952    integer(kind=kint) :: i
953    external fstr_sort_index
954
955    fstrPARAM%solution_type = kstSTATIC
956    fstrPARAM%nlgeom        = .false.
957    fstrPARAM%solver_method = ksmCG
958
959    !!STATIC !HEAT
960    fstrPARAM%analysis_n = 0
961    fstrPARAM%ref_temp   = 0
962
963    ! output control
964    fstrPARAM%fg_echo       = kOFF
965    fstrPARAM%fg_result     = kOFF
966    fstrPARAM%fg_visual     = kOFF
967
968    ! for heat ...
969    fstrPARAM%fg_neutral    = kOFF
970    fstrPARAM%fg_irres      = kNO
971    fstrPARAM%fg_iwres      = kNO
972    fstrPARAM%nrres         = 1
973    fstrPARAM%nprint        = 1
974
975    ! for couple
976    fstrPARAM%fg_couple      = 0
977    fstrPARAM%fg_couple_type = 0
978    fstrPARAM%fg_couple_first= 0
979    fstrPARAM%fg_couple_window= 0
980
981    ! for restart control
982    fstrPARAM%restart_version = 5
983
984    ! index table for global node ID sorting
985
986    fstrPARAM%n_node = hecMESH%n_node;
987    fstrPARAM%nn_internal = hecMESH%nn_internal;
988    allocate( fstrPARAM%global_local_ID(2,hecMESH%nn_internal))
989    do i = 1, hecMESH%nn_internal
990      fstrPARAM%global_local_ID(1,i) = hecMESH%global_node_ID(i)
991      fstrPARAM%global_local_ID(2,i) = i
992    end do
993    call fstr_sort_index( fstrPARAM%global_local_ID, hecMESH%nn_internal )
994  end subroutine fstr_param_init
995
996  logical function fstr_isBoundaryActive( fstrSOLID, nbc, cstep )
997    type(fstr_solid)    :: fstrSOLID
998    integer, intent(in) :: nbc    !< group id of boundary condition
999    integer, intent(in) :: cstep  !< current step number
1000    fstr_isBoundaryActive = .true.
1001    if( .not. associated(fstrSOLID%step_ctrl) ) return
1002    if( cstep>fstrSOLID%nstep_tot ) return
1003    fstr_isBoundaryActive = isBoundaryActive( nbc, fstrSOLID%step_ctrl(cstep) )
1004  end function
1005
1006  logical function fstr_isLoadActive( fstrSOLID, nbc, cstep )
1007    type(fstr_solid)    :: fstrSOLID
1008    integer, intent(in) :: nbc    !< group id of boundary condition
1009    integer, intent(in) :: cstep  !< current step number
1010    fstr_isLoadActive = .true.
1011    if( cstep > 0 ) then
1012      if( .not. associated(fstrSOLID%step_ctrl) ) return
1013      if( cstep>fstrSOLID%nstep_tot ) return
1014      fstr_isLoadActive = isLoadActive( nbc, fstrSOLID%step_ctrl(cstep) )
1015    else
1016      fstr_isLoadActive = isLoadActive( nbc, fstrSOLID%step_ctrl_restart )
1017    endif
1018  end function
1019
1020  logical function fstr_isContactActive( fstrSOLID, nbc, cstep )
1021    type(fstr_solid)     :: fstrSOLID !< type fstr_solid
1022    integer, intent(in) :: nbc       !< group id of boundary condition
1023    integer, intent(in) :: cstep     !< current step number
1024    fstr_isContactActive = .true.
1025    if( .not. associated(fstrSOLID%step_ctrl) ) return
1026    if( cstep>fstrSOLID%nstep_tot ) return
1027    fstr_isContactActive = isContactActive( nbc, fstrSOLID%step_ctrl(cstep) )
1028  end function
1029
1030  !> This subroutine fetch coords defined by local coordinate system
1031  subroutine get_coordsys( cdsys_ID, hecMESH, fstrSOLID, coords )
1032    integer, intent(in)             :: cdsys_ID      !< id of local coordinate
1033    type(hecmwST_local_mesh)       :: hecMESH       !< mesh information
1034    type(fstr_solid), intent(inout) :: fstrSOLID     !< fstr_solid
1035    real(kind=kreal), intent(out)   :: coords(3,3)
1036    integer :: ik
1037
1038    coords = 0.d0
1039    if( cdsys_ID>0 ) then
1040      if( isCoordNeeds(g_LocalCoordSys(cdsys_ID)) ) then
1041        coords=g_LocalCoordSys(cdsys_ID)%CoordSys
1042      else
1043        ik=g_LocalCoordSys(cdsys_ID)%node_ID(1)
1044        coords(1,:)= hecMESH%node(3*ik-2:3*ik)+fstrSOLID%unode(3*ik-2:3*ik)  &
1045          + fstrSOLID%dunode(3*ik-2:3*ik)
1046        ik=g_LocalCoordSys(cdsys_ID)%node_ID(2)
1047        coords(2,:)= hecMESH%node(3*ik-2:3*ik)+fstrSOLID%unode(3*ik-2:3*ik)  &
1048          + fstrSOLID%dunode(3*ik-2:3*ik)
1049        ik=g_LocalCoordSys(cdsys_ID)%node_ID(3)
1050        if(ik>0) coords(3,:)= hecMESH%node(3*ik-2:3*ik)+fstrSOLID%unode(3*ik-2:3*ik)  &
1051          + fstrSOLID%dunode(3*ik-2:3*ik)
1052      endif
1053    endif
1054  end subroutine get_coordsys
1055
1056  subroutine fstr_set_current_config_to_mesh(hecMESH,fstrSOLID,coord)
1057    implicit none
1058    type(hecmwST_local_mesh), intent(inout) :: hecMESH
1059    type (fstr_solid), intent(in) :: fstrSOLID
1060    real(kind=kreal), intent(out) :: coord(:)
1061    integer(kind=kint) :: i
1062    if(hecMESH%n_dof == 4) return
1063    do i = 1, hecMESH%nn_internal*min(hecMESH%n_dof,3)
1064      coord(i) = hecMESH%node(i)
1065      hecMESH%node(i) = coord(i)+fstrSOLID%unode(i)+fstrSOLID%dunode(i)
1066    enddo
1067  end subroutine fstr_set_current_config_to_mesh
1068
1069  subroutine fstr_recover_initial_config_to_mesh(hecMESH,fstrSOLID,coord)
1070    implicit none
1071    type(hecmwST_local_mesh), intent(inout) :: hecMESH
1072    type (fstr_solid), intent(in) :: fstrSOLID
1073    real(kind=kreal), intent(in) :: coord(:)
1074    integer(kind=kint) :: i
1075    if(hecMESH%n_dof == 4) return
1076    do i = 1, hecMESH%nn_internal*min(hecMESH%n_dof,3)
1077      hecMESH%node(i) = coord(i)
1078    enddo
1079  end subroutine fstr_recover_initial_config_to_mesh
1080
1081  subroutine fstr_solid_phys_zero(phys)
1082    implicit none
1083    type(fstr_solid_physic_val), pointer :: phys
1084    phys%STRAIN  = 0.0d0
1085    phys%STRESS  = 0.0d0
1086    phys%MISES   = 0.0d0
1087    phys%ESTRAIN = 0.0d0
1088    phys%ESTRESS = 0.0d0
1089    phys%EMISES  = 0.0d0
1090    phys%ENQM    = 0.0d0
1091  end subroutine fstr_solid_phys_zero
1092
1093  subroutine fstr_solid_phys_clear(fstrSOLID)
1094    implicit none
1095    type (fstr_solid)         :: fstrSOLID
1096    integer(kind=kint) :: i
1097
1098    if (associated(fstrSOLID%SOLID)) then
1099      call fstr_solid_phys_zero(fstrSOLID%SOLID)
1100    end if
1101    if (associated(fstrSOLID%SHELL)) then
1102      call fstr_solid_phys_zero(fstrSOLID%SHELL)
1103      do i=1,fstrSOLID%max_lyr
1104        call fstr_solid_phys_zero(fstrSOLID%SHELL%LAYER(i)%PLUS)
1105        call fstr_solid_phys_zero(fstrSOLID%SHELL%LAYER(i)%MINUS)
1106      end do
1107    end if
1108    if (associated(fstrSOLID%BEAM)) then
1109      call fstr_solid_phys_zero(fstrSOLID%BEAM)
1110    end if
1111  end subroutine fstr_solid_phys_clear
1112
1113end module m_fstr
1114