m_fstr.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
! If new header is supported, change according to following method.    !
!  1) Increase FSTR_CTRL_HEADER_NUMBER                                 !
!  2) Add new header name to fstr_ctrl_header_names                    !
!  3) Describe new function to get parameters from control file        !
!     in fstr_ctrl.f90                                                 !
!  4) Describe new subroutine to set values of the parameter           !
!     in fstr_setup.f90                                                !
!  5) If initial values are necessary, set the value                   !
!     in subroutine fstr_setup_init in fstr_setup.f90                  !
module m_fstr
  use hecmw
  use m_common_struct
  use m_step
  use m_out
  use m_timepoint
  use m_elemact
  use mmechgauss
  use mcontactdef
 
  implicit none
 
  public
 
  integer(kind=kint), parameter :: kyes = 1
  integer(kind=kint), parameter :: kno  = 0
  integer(kind=kint), parameter :: kon  = 1
  integer(kind=kint), parameter :: koff = 0
 
  integer(kind=kint), parameter :: kstprecheck = 0
  integer(kind=kint), parameter :: kststatic   = 1
  integer(kind=kint), parameter :: ksteigen    = 2
  integer(kind=kint), parameter :: kstheat     = 3
  integer(kind=kint), parameter :: kstdynamic  = 4
  !integer(kind=kint), parameter :: kstNLSTATIC  = 5
  integer(kind=kint), parameter :: kststaticeigen = 6
  integer(kind=kint), parameter :: kstnzprof      = 7
 
  integer(kind=kint), parameter :: ksmcg       = 1
  integer(kind=kint), parameter :: ksmbicgstab = 2
  integer(kind=kint), parameter :: ksmgmres    = 3
  integer(kind=kint), parameter :: ksmgpbicg   = 4
  integer(kind=kint), parameter :: ksmgmresr   = 5
  integer(kind=kint), parameter :: ksmgmresren = 6
  integer(kind=kint), parameter :: ksmpipecg   = 8
  integer(kind=kint), parameter :: ksmgroppcg  = 9
  integer(kind=kint), parameter :: ksmdirect   = 101
 
  integer(kind=kint), parameter :: knsmnewton      = 1
  integer(kind=kint), parameter :: knsmquasinewton = 2
 
  integer(kind=kint), parameter :: kcaslagrange = 1
  integer(kind=kint), parameter :: kcaalagrange = 2
 
  integer(kind=kint), parameter :: kbcffstr     =  0  ! BC described in fstr control file (default)
  integer(kind=kint), parameter :: kbcfnastran  =  1  ! nastran file
 
  integer(kind=kint), parameter :: kbcinitial   =  1
  integer(kind=kint), parameter :: kbctransit   =  2
 
  integer(kind=kint), parameter :: kmasslumped     = 1
  integer(kind=kint), parameter :: kmassconsistent = 2
 
  integer(kind=kint), parameter :: restart_outlast = 1
  integer(kind=kint), parameter :: restart_outall  = 2
 
  integer(kind=kint), parameter :: kel341fi     =  1
  integer(kind=kint), parameter :: kel341sesns  =  2
 
  integer(kind=kint), parameter :: kel361fi     =  1
  integer(kind=kint), parameter :: kel361bbar   =  2
  integer(kind=kint), parameter :: kel361ic     =  3
  integer(kind=kint), parameter :: kel361fbar   =  4
  integer(kind=kint), parameter :: kel361up     =  5
 
  integer(kind=kint), parameter :: kel611euler      =  1
  integer(kind=kint), parameter :: kel611timoshenko =  2
 
  integer(kind=kint), parameter :: kfloadtype_node = 1
  integer(kind=kint), parameter :: kfloadtype_surf = 2
 
  integer(kind=kint), parameter :: kfloadcase_re = 1
  integer(kind=kint), parameter :: kfloadcase_im = 2
 
  integer(kind=kint), parameter :: kitrcontinue = 0
  integer(kind=kint), parameter :: kitrconverged = 1
  integer(kind=kint), parameter :: kitrdiverged = 2
  integer(kind=kint), parameter :: kitrfloatingerror = 3
 
 
  integer(kind = kint) :: myrank
  integer(kind = kint) :: nprocs
 
  logical :: paracontactflag = .false.
 
  character(len=HECMW_FILENAME_LEN) :: cntfilname
  character(len=HECMW_FILENAME_LEN) :: restartfilname
 
  integer(kind=kint), parameter :: ilog = 16 ! log
  integer(kind=kint), parameter :: ista = 17 ! status
  integer(kind=kint), parameter :: iutb = 18 ! utable
  integer(kind=kint), parameter :: imsg = 51 ! message (myrank == 0 only)
  integer(kind=kint), parameter :: idbg = 52 ! debug
  integer(kind=kint), parameter :: ifvs = 53 ! visual.ini file
  integer(kind=kint), parameter :: ineu = 54 ! neutral file (heat)
  integer(kind=kint), parameter :: iresout = 100 ! ~110, keeping for result output file
 
  integer(kind=kint), parameter :: kopss_solution = 1
  integer(kind=kint), parameter :: kopss_material = 2
  integer(kind=kint)            :: opsstype = kopss_solution ! output stress/strain type
 
  type tinitialcondition
     character(len=HECMW_FILENAME_LEN)          :: cond_name
     integer                    :: node_elem = 0
     integer                    :: grpid = 0
     integer, pointer           :: intval(:)     => null()      
     real(kind=kreal), pointer  :: realval(:)    => null()      
  end type
  type( tinitialcondition ), pointer, save :: g_initialcnd(:) => null()
 
  type fstr_param
    integer(kind=kint) :: solution_type=0
    integer(kind=kint) :: solver_method=0
    integer(kind=kint) :: nlsolver_method=0
    logical            :: nlgeom=.false.             
 
    integer(kind=kint) :: analysis_n=0
    real(kind=kreal), pointer  :: dtime(:)=>null() 
    real(kind=kreal), pointer  :: etime(:)=>null() 
    real(kind=kreal), pointer  :: dtmin(:)=>null() 
    real(kind=kreal), pointer  :: delmax(:)=>null()
    integer(kind=kint), pointer:: itmax(:)=>null() 
    real(kind=kreal), pointer  :: eps(:)=>null()   
    real(kind=kreal)           :: ref_temp=0.d0    
    integer(kind=kint)         :: timepoint_id=0
 
    integer(kind=kint) :: fg_echo=koff       
    integer(kind=kint) :: fg_result=koff     
    integer(kind=kint) :: fg_visual=koff     
 
    integer(kind=kint) :: fg_neutral=koff   
    integer(kind=kint) :: fg_irres=kno      
    integer(kind=kint) :: fg_iwres=kno      
    integer(kind=kint) :: nrres=0
    integer(kind=kint) :: nprint=0
 
    integer(kind=kint) :: n_node=0
    integer(kind=kint) :: nn_internal=0
    integer(kind=kint), pointer :: global_local_id(:,:)=>null()
 
    integer( kind=kint ) :: fg_couple=0
    integer( kind=kint ) :: fg_couple_type=0
    integer( kind=kint ) :: fg_couple_first=0
    integer( kind=kint ) :: fg_couple_window=0
 
    integer( kind=kint ) :: restart_out_type=0
    integer( kind=kint ) :: restart_version=0
 
    integer( kind=kint ) :: contact_algo=0
    integer( kind=kint ) :: augiter=0
    logical              :: fric_cone_follow=.false.  
    real( kind=kreal )   :: eps_fric_band=0.d0        
    type(tcontactparam), pointer :: contactparam(:)=>null()  
    type(tcontactinterference), pointer :: contact_if(:)=>null()  
 
    type(tparamautoinc), pointer :: ainc(:)=>null()      
    type(time_points), pointer :: timepoints(:)=>null()  
 
    type(tparamconverg), pointer :: cnvparam(:)=>null()    
  end type fstr_param
 
  type fstr_solid_physic_val
    real(kind=kreal), pointer :: stress(:) => null()   
    real(kind=kreal), pointer :: strain(:) => null()   
    real(kind=kreal), pointer :: mises(:) => null()    
 
    real(kind=kreal), pointer :: pstress(:) => null()   
    real(kind=kreal), pointer :: pstrain(:) => null()   
    real(kind=kreal), pointer :: pstress_vect(:,:) => null()  
    real(kind=kreal), pointer :: pstrain_vect(:,:) => null()  
 
    real(kind=kreal), pointer :: estress(:) => null()  
    real(kind=kreal), pointer :: estrain(:) => null()  
    real(kind=kreal), pointer :: emises(:) => null()   
    real(kind=kreal), pointer :: eplstrain(:) => null()   
 
    real(kind=kreal), pointer :: epstress(:) => null()  
    real(kind=kreal), pointer :: epstrain(:) => null()  
    real(kind=kreal), pointer :: epstress_vect(:,:) => null()  
    real(kind=kreal), pointer :: epstrain_vect(:,:) => null()  
    real(kind=kreal), pointer :: enqm(:) => null()     
 
 
    type(fstr_solid_physic_val), pointer :: layer(:) => null()   
    type(fstr_solid_physic_val), pointer :: plus => null()   
    type(fstr_solid_physic_val), pointer :: minus => null()  
  end type fstr_solid_physic_val
 
  type fstr_solid
    integer(kind=kint) :: file_type=0         ! kbcfFSTR or kbcfNASTRAN
    integer(kind=kint) :: statictype=0        ! 1:Total, 2:Updated, 3:Infinitesimal
    integer(kind=kint) :: nstep_tot=0
 
    type(step_info), pointer       :: step_ctrl(:)  =>null()   
    type(t_output_ctrl),  pointer  :: output_ctrl(:)=>null()   
 
    integer(kind=kint) :: boundary_ngrp_tot=0
    integer(kind=kint), pointer :: boundary_ngrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_id     (:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_type   (:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: boundary_ngrp_val(:)  =>null()
    integer(kind=kint), pointer :: boundary_ngrp_istot  (:)  =>null()
    integer(kind=kint) :: boundary_ngrp_rot=0
    integer(kind=kint), pointer :: boundary_ngrp_rotid     (:) =>null()
    integer(kind=kint), pointer :: boundary_ngrp_centerid  (:) =>null()
 
    integer(kind=kint) :: velocity_type=0
    integer(kind=kint) :: velocity_ngrp_tot=0
    integer(kind=kint), pointer :: velocity_ngrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: velocity_ngrp_id     (:)  =>null()
    integer(kind=kint), pointer :: velocity_ngrp_type   (:)  =>null()
    integer(kind=kint), pointer :: velocity_ngrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: velocity_ngrp_val(:)  =>null()
    integer(kind=kint) :: velocity_ngrp_rot=0
    integer(kind=kint), pointer :: velocity_ngrp_rotid    (:) =>null()
    integer(kind=kint), pointer :: velocity_ngrp_centerid (:) =>null()
 
    integer(kind=kint) :: acceleration_type=0
    integer(kind=kint) :: acceleration_ngrp_tot=0
    integer(kind=kint), pointer :: acceleration_ngrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: acceleration_ngrp_id     (:)  =>null()
    integer(kind=kint), pointer :: acceleration_ngrp_type   (:)  =>null()
    integer(kind=kint), pointer :: acceleration_ngrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: acceleration_ngrp_val(:)  =>null()
 
    integer(kind=kint) :: cload_ngrp_tot=0
    integer(kind=kint), pointer :: cload_ngrp_grpid     (:) =>null()
    integer(kind=kint), pointer :: cload_ngrp_id        (:) =>null()
    integer(kind=kint), pointer :: cload_ngrp_dof       (:) =>null()
    integer(kind=kint), pointer :: cload_ngrp_amp       (:) =>null()
    real(kind=kreal), pointer   :: cload_ngrp_val(:) =>null()
    integer(kind=kint) :: cload_ngrp_rot=0
    integer(kind=kint), pointer :: cload_ngrp_rotid     (:) =>null()
    integer(kind=kint), pointer :: cload_ngrp_centerid  (:) =>null()
 
    integer(kind=kint) :: dload_ngrp_tot=0
    integer(kind=kint) :: dload_follow=0
    integer(kind=kint), pointer :: dload_ngrp_grpid     (:) =>null()
    integer(kind=kint), pointer :: dload_ngrp_id        (:) =>null()
    integer(kind=kint), pointer :: dload_ngrp_lid       (:) =>null()
    integer(kind=kint), pointer :: dload_ngrp_amp       (:) =>null()
    real(kind=kreal), pointer   :: dload_ngrp_params(:,:) =>null()
 
    integer(kind=kint) :: temp_ngrp_tot=0
    integer(kind=kint) :: temp_irres=0
    integer(kind=kint) :: temp_tstep=0
    integer(kind=kint) :: temp_interval=0
    integer(kind=kint) :: temp_rtype=0      ! type of reading result; 1: step-based; 2: time-based
    real(kind=kreal)   :: temp_factor=0.d0
    integer(kind=kint), pointer :: temp_ngrp_grpid     (:) =>null()
    integer(kind=kint), pointer :: temp_ngrp_id        (:) =>null()
    real(kind=kreal), pointer   :: temp_ngrp_val(:) =>null()
 
    integer(kind=kint) :: spring_ngrp_tot=0
    integer(kind=kint), pointer :: spring_ngrp_grpid    (:) =>null()
    integer(kind=kint), pointer :: spring_ngrp_id       (:) =>null()
    integer(kind=kint), pointer :: spring_ngrp_dof      (:) =>null()
    integer(kind=kint), pointer :: spring_ngrp_amp      (:) =>null()
    real(kind=kreal), pointer   :: spring_ngrp_val(:) =>null()
    integer(kind=kint), pointer :: spring_incremental   (:) =>null()    
 
    type(telemact) :: elemact
 
    integer( kind=kint ) :: couple_ngrp_tot=0
    integer( kind=kint ),pointer :: couple_ngrp_id(:) =>null()
 
    integer(kind=kint) :: maxn_gauss=0
 
    real(kind=kreal), pointer :: stress(:) =>null()    
    real(kind=kreal), pointer :: strain(:) =>null()    
    real(kind=kreal), pointer :: mises(:) =>null()     
 
    real(kind=kreal), pointer :: pstress(:) =>null()   
    real(kind=kreal), pointer :: pstrain(:) =>null()   
    real(kind=kreal), pointer :: pstress_vect(:,:) =>null()   
    real(kind=kreal), pointer :: pstrain_vect(:,:) =>null()   
 
    real(kind=kreal), pointer :: estress(:) =>null()   
    real(kind=kreal), pointer :: estrain(:) =>null()   
    real(kind=kreal), pointer :: emises(:) =>null()    
    real(kind=kreal), pointer :: eplstrain(:) =>null()    
 
    real(kind=kreal), pointer :: epstress(:) =>null()   
    real(kind=kreal), pointer :: epstrain(:) =>null()   
    real(kind=kreal), pointer :: epstress_vect(:,:) =>null()   
    real(kind=kreal), pointer :: epstrain_vect(:,:) =>null()   
 
    real(kind=kreal), pointer :: tnstrain(:) =>null()   
    real(kind=kreal), pointer :: testrain(:) =>null()   
 
    real(kind=kreal), pointer :: yield_ratio(:) =>null()    
 
    real(kind=kreal), pointer :: enqm(:) =>null()      
    real(kind=kreal), pointer :: reaction(:) =>null()    
 
    real(kind=kreal), pointer :: cont_nforce(:) =>null()  
    real(kind=kreal), pointer :: cont_fric(:) =>null()    
    real(kind=kreal), pointer :: cont_relvel(:) =>null()  
    real(kind=kreal), pointer :: cont_state(:) =>null()   
    integer(kind=kint), pointer :: cont_sgrp_id(:) =>null() 
    real(kind=kreal), pointer :: cont_area(:) =>null()    
    real(kind=kreal), pointer :: cont_ntrac(:) =>null()   
    real(kind=kreal), pointer :: cont_ftrac(:) =>null()   
    real(kind=kreal), pointer :: embed_nforce(:) =>null()  
 
    type(fstr_solid_physic_val), pointer :: solid=>null()     
    type(fstr_solid_physic_val), pointer :: shell=>null()     
    type(fstr_solid_physic_val), pointer :: beam =>null()     
 
    integer(kind=kint) :: restart_nout=0
    integer(kind=kint) :: restart_nin=0
    type(step_info)    :: step_ctrl_restart
 
    integer(kind=kint) :: max_lyr=0
    integer(kind=kint) :: is_33shell=0
    integer(kind=kint) :: is_33beam=0
    integer(kind=kint) :: is_heat=0
    integer(kind=kint) :: max_ncon_stf=0
    integer(kind=kint) :: max_ncon=0
    integer(kind=kint), pointer :: is_rot(:) => null()
    integer(kind=kint) :: elemopt361=0
    logical            :: is_smoothing_active=.false.
    real(kind=kreal)   :: factor(2)=0.d0      
    integer(kind=kint) :: nrstat_i(10)=0
    real(kind=kreal)   :: nrstat_r(10)=0.d0  
    integer(kind=kint) :: autoinc_stat=0
    integer(kind=kint) :: cutback_stat=0
 
    real(kind=kreal), pointer :: gl(:) => null()          
    real(kind=kreal), pointer :: gl0(:) => null()           
    real(kind=kreal), pointer :: eforce(:) => null()           
    real(kind=kreal), pointer :: qforce(:) => null()           
    real(kind=kreal), pointer :: qforce_bak(:) => null()           
    real(kind=kreal), pointer :: dforce(:) => null()           
    real(kind=kreal), pointer :: unode(:)      => null()   
    real(kind=kreal), pointer :: unode_bak(:)  => null()   
    real(kind=kreal), pointer :: dunode(:)     => null()   
    real(kind=kreal), pointer :: ddunode(:)    => null()   
    logical :: has_finite_rotation_kinematics=.false.
    logical :: finite_rotation_state_ready=.false.
    integer(kind=kint), pointer :: shell_node_mode(:) => null() 
    integer(kind=kint), pointer :: shell_rot_state(:) => null() 
    real(kind=kreal), pointer :: shell_ref_triad(:) => null()   
    real(kind=kreal), pointer :: shell_triad(:)     => null()   
    real(kind=kreal), pointer :: shell_triad_bak(:) => null()   
    real(kind=kreal), pointer :: shell_dtriad(:)    => null()   
    real(kind=kreal), pointer :: shell_drill(:)     => null()   
    real(kind=kreal), pointer :: shell_drill_bak(:) => null()   
    real(kind=kreal), pointer :: shell_ddrill(:)    => null()   
    real(kind=kreal), pointer :: temperature(:)=> null()   
    real(kind=kreal), pointer :: temp_bak(:)   => null()
    real(kind=kreal), pointer :: last_temp(:)  => null()
 
    type( telement ), pointer :: elements(:)   =>null()  
    type( tmaterial ),pointer :: materials(:)  =>null()  
    integer                   :: n_contacts=0
    type( tcontact ), pointer :: contacts(:)   =>null()  
    integer                   :: n_embeds=0
    type( tcontact ), pointer :: embeds(:)     =>null()  
    integer                   :: n_fix_mpc=0
    real(kind=kreal), pointer :: mpc_const(:)  =>null()  
    type(tsection), pointer   :: sections(:)   =>null()  
 
    ! for cutback
    ! ####################### Notice #######################
    ! # If you add new variables to store analysis status, #
    ! # - backup variables with postfix "_bkup" here       #
    ! # - backup process to module m_fstr_Cutback          #
    ! # must be added if necessary.                        #
    ! ######################################################
    real(kind=kreal), pointer :: unode_bkup(:)     => null() 
    real(kind=kreal), pointer :: qforce_bkup(:)    => null() 
    real(kind=kreal), pointer :: dforce_bkup(:)    => null() 
    real(kind=kreal), pointer :: last_temp_bkup(:) => null()
    real(kind=kreal), pointer :: shell_triad_bkup(:) => null()
    real(kind=kreal), pointer :: shell_drill_bkup(:) => null()
    type( telement ), pointer :: elements_bkup(:)  =>null()  
    type( tcontact ), pointer :: contacts_bkup(:)  =>null()  
    type( tcontact ), pointer :: embeds_bkup(:)  =>null()  
  end type fstr_solid
 
  type fstr_heat
    integer(kind=kint) :: is_steady=0
    real(kind=kreal)   :: beta=0.d0
    logical :: is_iter_max_limit=.false.
 
    integer(kind=kint) :: steptot=0
    integer(kind=kint) :: restart_nout=0
    real(kind=kreal), pointer :: step_dltime(:)=>null(), step_eetime(:)=>null()
    real(kind=kreal), pointer :: step_delmin(:)=>null(), step_delmax(:)=>null()
    integer(kind=kint) :: timepoint_id=0
 
    integer(kind=kint) :: materialtot=0
    integer(kind=kint), pointer :: rhotab(:)=>null(), cptab(:)=>null(), condtab(:)=>null()
    real(kind=kreal), pointer :: rho(:,:)=>null(), rhotemp(:,:)=>null()
    real(kind=kreal), pointer :: cp(:,:)=>null(),  cptemp(:,:)=>null()
    real(kind=kreal), pointer :: cond(:,:)=>null(), condtemp(:,:)=>null()
 
    real(kind=kreal), pointer :: rhofunca(:,:)=>null(), rhofuncb(:,:)=>null()
    real(kind=kreal), pointer :: cpfunca(:,:)=>null(), cpfuncb(:,:)=>null()
    real(kind=kreal), pointer :: condfunca(:,:)=>null(), condfuncb(:,:)=>null()
 
    integer(kind=kint) :: amplitudetot=0
    integer(kind=kint), pointer :: ampltab(:)=>null()
    real(kind=kreal), pointer :: ampl(:,:)=>null(), ampltime(:,:)=>null()
    real(kind=kreal), pointer :: amplfunca(:,:)=>null(),  amplfuncb(:,:)=>null()
 
    real(kind=kreal), pointer :: temp0(:)=>null()
    real(kind=kreal), pointer :: tempc(:)=>null()
    real(kind=kreal), pointer :: temp(:)=>null()
 
    logical, pointer :: is_dummy_node(:)=>null()
 
    integer(kind=kint) :: t_fix_tot=0
    integer(kind=kint), pointer :: t_fix_node(:)=>null()
    integer(kind=kint), pointer :: t_fix_ampl(:)=>null()
    real(kind=kreal), pointer :: t_fix_val(:)=>null()
 
    integer(kind=kint) :: q_nod_tot=0
    integer(kind=kint), pointer :: q_nod_node(:)=>null()
    integer(kind=kint), pointer :: q_nod_ampl(:)=>null()
    real(kind=kreal), pointer :: q_nod_val(:)=>null()
 
    integer(kind=kint) :: q_vol_tot=0
    integer(kind=kint), pointer :: q_vol_elem(:)=>null()
    integer(kind=kint), pointer :: q_vol_ampl(:)=>null()
    real(kind=kreal), pointer :: q_vol_val(:)=>null()
 
    integer(kind=kint) :: q_suf_tot=0
    integer(kind=kint), pointer :: q_suf_elem(:)=>null()
    integer(kind=kint), pointer :: q_suf_ampl(:)=>null()
    integer(kind=kint), pointer :: q_suf_surf(:)=>null()
    real(kind=kreal), pointer :: q_suf_val(:)=>null()
 
    integer(kind=kint) :: r_suf_tot=0
    integer(kind=kint), pointer :: r_suf_elem(:)=>null()
    integer(kind=kint), pointer :: r_suf_ampl(:,:)=>null()
    integer(kind=kint), pointer :: r_suf_surf(:)=>null()
    real(kind=kreal), pointer :: r_suf_val(:,:)=>null()
 
    integer(kind=kint) :: h_suf_tot=0
    integer(kind=kint), pointer :: h_suf_elem(:)=>null()
    integer(kind=kint), pointer :: h_suf_ampl(:,:)=>null()
    integer(kind=kint), pointer :: h_suf_surf(:)=>null()
    real(kind=kreal), pointer :: h_suf_val(:,:)=>null()
 
    integer(kind=kint)  :: wl_tot=0
    type(tweldline), pointer :: weldline(:) => null()
 
    type(telemact) :: elemact
  end type fstr_heat
 
  type fstr_dynamic
    integer(kind=kint) :: idx_eqa=0     ! implicit or explicit
    integer(kind=kint) :: idx_resp=0    ! time history or steady-state harmonic response analysis
 
    integer(kind=kint) :: n_step=0      ! total step number of analysis
    real(kind=kreal)   :: t_start=0.d0  ! start time of analysis
    real(kind=kreal)   :: t_curr=0.d0   ! current time of analysis
    real(kind=kreal)   :: t_end=0.d0    ! end time of analysis
    real(kind=kreal)   :: t_delta=0.d0  ! time increment
    integer(kind=kint) :: restart_nout=0! output interval of restart file
    ! (if  .gt.0) restart file write
    ! (if  .lt.0) restart file read and write
    integer(kind=kint) :: restart_nin=0 !input number of restart file
 
    real(kind=kreal)   :: gamma=0.d0    ! Newmark-beta parameter gamma
    real(kind=kreal)   :: beta=0.d0     ! Newmark-beta parameter beta
 
    real(kind=kreal)   :: alpha=0.d0    ! HHT parameter alpha
 
    integer(kind=kint) :: idx_mas=0     ! kMassLumped or kMassConsistent
 
    integer(kind=kint) :: idx_dmp=0     ! damping type
    real(kind=kreal)   :: ray_m=0.d0    ! Rayleigh damping parameter Rm
    real(kind=kreal)   :: ray_k=0.d0    ! Rayleigh damping parameter Rk
 
    logical           :: varinitialize=.false. ! initialization flag
 
    integer(kind=kint) :: nout=0         ! output interval of result
    integer(kind=kint) :: ngrp_monit=0   ! node of monitoring result
    integer(kind=kint) :: nout_monit=0   ! output interval of result monitoring
    integer(kind=kint) :: iout_list(6)=0 ! 0:not output  1:output
    ! iout_list(1): displacement
    ! iout_list(2): velocity
    ! iout_list(3): acceleration
    ! iout_list(4): reaction force
    ! iout_list(5): strain
    ! iout_list(6): stress
 
    real(kind=kreal), pointer :: disp(:,:) => null()
    real(kind=kreal), pointer :: vel(:,:) => null()
    real(kind=kreal), pointer :: acc(:,:) => null()
    real(kind=kreal), pointer :: hht_force(:) => null()
 
    real(kind=kreal) :: kineticenergy=0.d0
    real(kind=kreal) :: strainenergy=0.d0
    real(kind=kreal) :: totalenergy=0.d0
 
    real(kind=kreal), pointer :: vec1(:) => null()
    real(kind=kreal), pointer :: vec2(:) => null()
    real(kind=kreal), pointer :: vec3(:) => null()
 
    integer(kind=kint) :: dynamic_iw4        =   204
    integer(kind=kint) :: dynamic_iw5        =   205
    integer(kind=kint) :: dynamic_iw6        =   206
    integer(kind=kint) :: dynamic_iw7        =   207
    integer(kind=kint) :: dynamic_iw8        =   208
    integer(kind=kint) :: dynamic_iw9        =   209
    integer(kind=kint) :: dynamic_iw10       =   210
  end type fstr_dynamic
 
  type fstr_freqanalysis
    integer(kind=kint)                :: fload_ngrp_tot=0
    integer(kind=kint), pointer       :: fload_ngrp_grpid(:) => null()
    integer(kind=kint), pointer       :: fload_ngrp_id(:)    => null()
    integer(kind=kint), pointer       :: fload_ngrp_type(:)  => null()
    integer(kind=kint), pointer       :: fload_ngrp_dof(:)   => null()
    real(kind=kreal), pointer         :: fload_ngrp_valre(:) => null()
    real(kind=kreal), pointer         :: fload_ngrp_valim(:) => null()
    !acceleration input given by !FACCELERATION (one entry per data line)
    integer(kind=kint)                :: facc_ngrp_tot=0
    integer(kind=kint), pointer       :: facc_ngrp_id(:)     => null() 
    integer(kind=kint), pointer       :: facc_ngrp_type(:)   => null() 
    real(kind=kreal), pointer         :: facc_ngrp_valre(:)  => null() 
    real(kind=kreal), pointer         :: facc_ngrp_valim(:)  => null() 
    character(len=HECMW_FILENAME_LEN) :: eigenlog_filename = ""
    integer(kind=kint)                :: start_mode=0
    integer(kind=kint)                :: end_mode=0
  end type fstr_freqanalysis
 
  type fstr_freqanalysis_data
    integer(kind=kint)        :: nummode=0
    integer(kind=kint)        :: numnodedof=0
    real(kind=kreal), pointer :: eigomega(:)    => null()
    real(kind=kreal), pointer :: eigvector(:,:) => null()
    real(kind=kreal)           :: rayalpha=0.d0, raybeta=0.d0
  end type fstr_freqanalysis_data
 
  type fstr_eigen
    integer(kind=kint)  :: nget=0    ! Solved eigen value number (default:5)
    integer(kind=kint)  :: maxiter=0 ! Max. Lcz iterations (default:60)
    integer(kind=kint)  :: iter=0    ! Max. Lcz iterations (default:60)
    real   (kind=kreal) :: sigma=0.d0     ! 0.0
    real   (kind=kreal) :: tolerance=0.d0 ! Lcz tolerance (default:1.0e-8)
    real   (kind=kreal) :: totalmass=0.d0
    real   (kind=kreal), pointer :: eigval(:) => null()
    real   (kind=kreal), pointer :: eigvec(:,:) => null()
    real   (kind=kreal), pointer :: resid(:) => null()
    real   (kind=kreal), pointer :: filter(:) => null()
    real   (kind=kreal), pointer :: mass(:) => null()
    real   (kind=kreal), pointer :: effmass(:) => null()
    real   (kind=kreal), pointer :: partfactor(:) => null()
    logical :: is_free = .false.
  end type fstr_eigen
 
  type fstr_couple
    integer( kind=kint )   :: dof=0            ! == 3
    integer( kind=kint )   :: ndof=0           ! total dof (coupled_node_n*dof)
    integer( kind=kint )   :: coupled_node_n=0
    integer, pointer       :: coupled_node(:) => null()  ! local node id sent from revocap
    real( kind=8 ),pointer :: trac(:) => null()          ! input  (x,y,z,x,y,z ... )
    real( kind=8 ),pointer :: disp(:) => null()          ! output (x,y,z,x,y,z ... )
    real( kind=8 ),pointer :: velo(:) => null()          ! output (x,y,z,x,y,z ... )
    real( kind=8 ),pointer :: accel(:) => null()         ! output (x,y,z,x,y,z ... )
    integer( kind=kint ),pointer :: index(:) => null()   ! size:total node num.
  end type fstr_couple
 
  type tweldline
    integer(kind=kint) :: egrpid=0
    real( kind=kreal ) :: i=0.d0
    real( kind=kreal ) :: u=0.d0
    real( kind=kreal ) :: coe=0.d0
    real( kind=kreal ) :: v=0.d0
    integer(kind=kint) :: xyz=0
    real(kind=kreal)   :: n1=0.d0, n2=0.d0
    real(kind=kreal)   :: distol=0.d0
    real(kind=kreal)   :: tstart=0.d0
  end type tweldline
 
  type tsection
    !integer              :: mat_ID=0
    !integer              :: iset=0
    !integer              :: orien_ID=0
    real(kind=kreal)     :: thickness=0.d0
    integer              :: elemopt341=0
    !integer              :: elemopt342=0
    !integer              :: elemopt351=0
    !integer              :: elemopt352=0
    integer              :: elemopt361=0
    !integer              :: elemopt362=0
    integer              :: elemopt611=0
  end type tsection
 
contains
 
  subroutine fstr_abort( code )
    integer(kind=kint), intent(in) :: code
    integer(kind=kint) :: ios
 
    flush( ilog, iostat=ios )
    if( myrank == 0 ) then
      flush( ista, iostat=ios )
      flush( imsg, iostat=ios )
    endif
    call hecmw_abort( hecmw_comm_get_comm(), code )
  end subroutine fstr_abort
 
  subroutine fstr_nullify_initial_condition( cond )
    implicit none
    type( tinitialcondition ) :: cond
    nullify( cond%intval )
    nullify( cond%realval )
  end subroutine
 
  subroutine fstr_nullify_fstr_param( P )
    implicit none
    type( fstr_param ) :: P
 
    nullify( p%dtime )
    nullify( p%etime )
    nullify( p%dtmin )
    nullify( p%delmax )
    nullify( p%itmax )
    nullify( p%eps )
    nullify( p%global_local_ID)
    nullify( p%contactparam )
    nullify( p%contact_if )
    nullify( p%ainc )
    nullify( p%timepoints )
    nullify( p%cnvparam )
  end subroutine fstr_nullify_fstr_param
 
  subroutine fstr_nullify_solid_physic_val( P )
    implicit none
    type( fstr_solid_physic_val ) :: P
    nullify( p%STRESS )
    nullify( p%STRAIN )
    nullify( p%MISES )
 
    nullify( p%PSTRESS )
    nullify( p%PSTRAIN )
    nullify( p%PSTRESS_VECT )
    nullify( p%PSTRAIN_VECT )
 
    nullify( p%ESTRESS )
    nullify( p%ESTRAIN )
    nullify( p%EMISES )
    nullify( p%EPLSTRAIN )
 
    nullify( p%EPSTRESS )
    nullify( p%EPSTRAIN )
    nullify( p%EPSTRESS_VECT )
    nullify( p%EPSTRAIN_VECT )
    nullify( p%ENQM )
 
    nullify( p%LAYER )
    nullify( p%PLUS )
    nullify( p%MINUS )
  end subroutine fstr_nullify_solid_physic_val
 
  subroutine fstr_nullify_fstr_solid( S )
    implicit none
    type( fstr_solid ) :: S
 
    nullify( s%step_ctrl )
    nullify( s%output_ctrl )
 
    nullify( s%BOUNDARY_ngrp_GRPID )
    nullify( s%BOUNDARY_ngrp_ID )
    nullify( s%BOUNDARY_ngrp_type )
    nullify( s%BOUNDARY_ngrp_amp )
    nullify( s%BOUNDARY_ngrp_val )
    nullify( s%BOUNDARY_ngrp_istot )
    nullify( s%BOUNDARY_ngrp_rotID )
    nullify( s%BOUNDARY_ngrp_centerID )
 
    nullify( s%VELOCITY_ngrp_GRPID )
    nullify( s%VELOCITY_ngrp_ID )
    nullify( s%VELOCITY_ngrp_type )
    nullify( s%VELOCITY_ngrp_amp )
    nullify( s%VELOCITY_ngrp_val )
    nullify( s%VELOCITY_ngrp_rotID )
    nullify( s%VELOCITY_ngrp_centerID )
 
    nullify( s%ACCELERATION_ngrp_GRPID )
    nullify( s%ACCELERATION_ngrp_ID )
    nullify( s%ACCELERATION_ngrp_type )
    nullify( s%ACCELERATION_ngrp_amp )
    nullify( s%ACCELERATION_ngrp_val )
 
    nullify( s%CLOAD_ngrp_GRPID )
    nullify( s%CLOAD_ngrp_ID )
    nullify( s%CLOAD_ngrp_DOF )
    nullify( s%CLOAD_ngrp_amp )
    nullify( s%CLOAD_ngrp_val )
    nullify( s%CLOAD_ngrp_rotID )
    nullify( s%CLOAD_ngrp_centerID )
 
    nullify( s%DLOAD_ngrp_GRPID )
    nullify( s%DLOAD_ngrp_ID )
    nullify( s%DLOAD_ngrp_LID )
    nullify( s%DLOAD_ngrp_amp )
    nullify( s%DLOAD_ngrp_params )
 
    nullify( s%TEMP_ngrp_GRPID )
    nullify( s%TEMP_ngrp_ID )
    nullify( s%TEMP_ngrp_val )
 
    nullify( s%SPRING_ngrp_GRPID )
    nullify( s%SPRING_ngrp_ID )
    nullify( s%SPRING_ngrp_DOF )
    nullify( s%SPRING_ngrp_amp )
    nullify( s%SPRING_ngrp_val )
    nullify( s%SPRING_incremental )
 
    nullify( s%COUPLE_ngrp_ID )
 
    nullify( s%STRESS )
    nullify( s%STRAIN )
    nullify( s%MISES )
    nullify( s%PSTRESS )
    nullify( s%PSTRAIN )
    nullify( s%PSTRESS_VECT )
    nullify( s%PSTRAIN_VECT )
    nullify( s%ESTRESS )
    nullify( s%ESTRAIN )
    nullify( s%EMISES )
    nullify( s%EPLSTRAIN )
    nullify( s%EPSTRESS )
    nullify( s%EPSTRAIN )
    nullify( s%EPSTRESS_VECT )
    nullify( s%EPSTRAIN_VECT )
    nullify( s%TNSTRAIN )
    nullify( s%TESTRAIN )
    nullify( s%YIELD_RATIO )
    nullify( s%ENQM )
    nullify( s%REACTION )
 
    nullify( s%CONT_NFORCE )
    nullify( s%CONT_FRIC )
    nullify( s%CONT_RELVEL )
    nullify( s%CONT_STATE )
    nullify( s%CONT_SGRP_ID )
    nullify( s%CONT_AREA )
    nullify( s%CONT_NTRAC )
    nullify( s%CONT_FTRAC )
    nullify( s%EMBED_NFORCE )
 
    nullify( s%SOLID )
    nullify( s%SHELL )
    nullify( s%BEAM )
 
    nullify( s%is_rot )
 
    nullify( s%GL          )
    nullify( s%GL0         )
    nullify( s%EFORCE      )
    nullify( s%QFORCE      )
    nullify( s%QFORCE_bak  )
    nullify( s%DFORCE      )
    nullify( s%unode       )
    nullify( s%unode_bak   )
    nullify( s%dunode      )
    nullify( s%ddunode     )
 
    nullify( s%shell_node_mode )
    nullify( s%shell_rot_state )
    nullify( s%shell_ref_triad )
    nullify( s%shell_triad )
    nullify( s%shell_triad_bak )
    nullify( s%shell_dtriad )
    nullify( s%shell_drill )
    nullify( s%shell_drill_bak )
    nullify( s%shell_ddrill )
    nullify( s%temperature )
    nullify( s%temp_bak )
    nullify( s%last_temp )
 
    nullify( s%elements )
    nullify( s%materials )
    nullify( s%contacts )
    nullify( s%embeds )
    nullify( s%mpc_const )
    nullify( s%sections )
 
    nullify( s%unode_bkup )
    nullify( s%QFORCE_bkup )
    nullify( s%DFORCE_bkup )
    nullify( s%last_temp_bkup )
    nullify( s%shell_triad_bkup )
    nullify( s%shell_drill_bkup )
    nullify( s%elements_bkup )
    nullify( s%contacts_bkup )
    nullify( s%embeds_bkup )
 
  end subroutine fstr_nullify_fstr_solid
 
  subroutine fstr_nullify_fstr_heat( H )
    implicit none
    type( fstr_heat ) :: H
 
    nullify( h%STEP_DLTIME )
    nullify( h%STEP_EETIME )
    nullify( h%STEP_DELMIN )
    nullify( h%STEP_DELMAX )
    nullify( h%RHOtab )
    nullify( h%CPtab )
    nullify( h%CONDtab )
    nullify( h%RHO )
    nullify( h%RHOtemp  )
    nullify( h%CP )
    nullify( h%CPtemp  )
    nullify( h%COND )
    nullify( h%CONDtemp  )
    nullify( h%RHOfuncA )
    nullify( h%RHOfuncB )
    nullify( h%CPfuncA  )
    nullify( h%CPfuncB )
    nullify( h%CONDfuncA  )
    nullify( h%CONDfuncB )
    nullify( h%AMPLtab )
    nullify( h%AMPL )
    nullify( h%AMPLtime  )
    nullify( h%AMPLfuncA )
    nullify( h%AMPLfuncB )
    nullify( h%TEMP0 )
    nullify( h%TEMPC )
    nullify( h%TEMP  )
    nullify( h%is_dummy_node )
    nullify( h%T_FIX_node )
    nullify( h%T_FIX_ampl )
    nullify( h%T_FIX_val )
    nullify( h%Q_NOD_node )
    nullify( h%Q_NOD_ampl )
    nullify( h%Q_NOD_val )
    nullify( h%Q_VOL_elem )
    nullify( h%Q_VOL_ampl )
    nullify( h%Q_VOL_val )
    nullify( h%Q_SUF_elem )
    nullify( h%Q_SUF_ampl )
    nullify( h%Q_SUF_surf )
    nullify( h%Q_SUF_val )
    nullify( h%R_SUF_elem )
    nullify( h%R_SUF_ampl )
    nullify( h%R_SUF_surf )
    nullify( h%R_SUF_val )
    nullify( h%H_SUF_elem )
    nullify( h%H_SUF_ampl )
    nullify( h%H_SUF_surf )
    nullify( h%H_SUF_val )
    nullify( h%weldline )
  end subroutine fstr_nullify_fstr_heat
 
  subroutine fstr_nullify_fstr_dynamic( DY )
    implicit none
    type( fstr_dynamic ) :: DY
 
    nullify( dy%DISP )
    nullify( dy%VEL  )
    nullify( dy%ACC  )
    nullify( dy%HHT_FORCE  )
    nullify( dy%VEC1 )
    nullify( dy%VEC2 )
    nullify( dy%VEC3 )
  end subroutine fstr_nullify_fstr_dynamic
 
  subroutine fstr_nullify_fstr_freqanalysis( f )
    implicit none
    type( fstr_freqanalysis ), intent(inout) :: f
 
    f%FLOAD_ngrp_tot = 0
    nullify( f%FLOAD_ngrp_GRPID )
    nullify( f%FLOAD_ngrp_ID )
    nullify( f%FLOAD_ngrp_TYPE )
    nullify( f%FLOAD_ngrp_DOF )
    nullify( f%FLOAD_ngrp_valre )
    nullify( f%FLOAD_ngrp_valim )
    !no acceleration input unless !FACCELERATION is given
    f%FACC_ngrp_tot = 0
    nullify( f%FACC_ngrp_ID )
    nullify( f%FACC_ngrp_type )
    nullify( f%FACC_ngrp_valre )
    nullify( f%FACC_ngrp_valim )
  end subroutine fstr_nullify_fstr_freqanalysis
 
  subroutine fstr_nullify_fstr_freqanalysis_data( f )
    implicit none
    type( fstr_freqanalysis_data ), intent(inout) :: f
    nullify( f%eigOmega )
    nullify( f%eigVector )
  end subroutine fstr_nullify_fstr_freqanalysis_data
 
  subroutine fstr_nullify_fstr_eigen( E )
    implicit none
    type( fstr_eigen ) :: E
    nullify( e%eigval )
    nullify( e%eigvec )
    nullify( e%resid )
    nullify( e%filter )
    nullify( e%mass )
    nullify( e%effmass )
    nullify( e%partfactor )
  end subroutine fstr_nullify_fstr_eigen
 
  subroutine fstr_nullify_fstr_couple( C )
    implicit none
    type( fstr_couple ) :: C
 
    nullify( c%coupled_node )
    nullify( c%trac )
    nullify( c%disp )
    nullify( c%velo )
    nullify( c%index )
  end subroutine fstr_nullify_fstr_couple
 
  subroutine fstr_mat_init( hecMAT )
    implicit none
    type(hecmwst_matrix) :: hecMAT
 
    hecmat%Iarray(1) =  100    ! = nier
    hecmat%Iarray(2) =    1    ! = method
    hecmat%Iarray(3) =    1    ! = precond
    hecmat%Iarray(4) =    0    ! = nset
    hecmat%Iarray(5) =    1    ! = iterpremax
    hecmat%Iarray(6) =   10    ! = nrest
    hecmat%Iarray(7) =    0    ! = scaling
    hecmat%Iarray(21)=  kno    ! = iterlog
    hecmat%Iarray(22)=  kno    ! = timelog
    hecmat%Iarray(31)=    0    ! = dumptype
    hecmat%Iarray(32)=    0    ! = dumpexit
    hecmat%Iarray(33)=    0    ! = usejad
    hecmat%Iarray(34)=   10    ! = ncolor_in
    hecmat%Iarray(13)=    0    ! = mpc_method
    hecmat%Iarray(14)=    0    ! = estcond
    hecmat%Iarray(16)=    0    ! = recompute_residual
    hecmat%Iarray(35)=    3    ! = maxrecycle_precond
    hecmat%Iarray(36)= hecmw_matvec_impl_default   ! = matvec_impl
    hecmat%Iarray(37)= hecmw_precond_impl_default  ! = precond_impl
    hecmat%Iarray(41)=    0    ! = solver_opt1
    hecmat%Iarray(42)=    0    ! = solver_opt2
    hecmat%Iarray(43)=    0    ! = solver_opt3
    hecmat%Iarray(44)=    0    ! = solver_opt4
    hecmat%Iarray(45)=    0    ! = solver_opt5
    hecmat%Iarray(46)=    0    ! = solver_opt6
 
    hecmat%Rarray(1) =  1.0e-8 ! = resid
    hecmat%Rarray(2) =  1.0    ! = sigma_diag
    hecmat%Rarray(3) =  1.0    ! = omega
    hecmat%Rarray(4) =  0.1    ! = thresh
    hecmat%Rarray(5) =  0.1    ! = filter
    hecmat%Rarray(11)=  1.0e+4 ! = penalty
 
    hecmat%Iarray(96) =   0    ! nrecycle_precond
    hecmat%Iarray(97) = kyes   ! flag_numfact
    hecmat%Iarray(98) = kyes   ! flag_symbfact
    hecmat%Iarray(99) = kyes   ! indirect method
  end subroutine fstr_mat_init
 
  subroutine hecmat_init( hecMAT )
    implicit none
    type( hecmwst_matrix ) :: hecMAT
    integer ::  ndof, nn, ierror
    ndof = hecmat%NDOF
    nn = ndof*ndof
    allocate (hecmat%AL(nn*hecmat%NPL)        ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    allocate (hecmat%AU(nn*hecmat%NPU)        ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    allocate (hecmat%B(ndof*hecmat%NP)          ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    hecmat%B(:)=0.d0
    allocate (hecmat%D(nn*hecmat%NP)          ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    allocate (hecmat%X(ndof*hecmat%NP)          ,stat=ierror )
    if( ierror /= 0 ) then
      write(*,*) "##ERROR : not enough memory"
      write(idbg,*) 'stop due to allocation error'
      call flush(idbg)
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
    hecmat%D  = 0.0d0
    hecmat%AL = 0.0d0
    hecmat%AU = 0.0d0
    hecmat%B  = 0.0d0
    hecmat%X  = 0.0d0
  end subroutine hecmat_init
 
  subroutine hecmat_finalize( hecMAT )
    implicit none
    type( hecmwst_matrix ) :: hecMAT
    integer ::  ndof, nn, ierror
    ndof = hecmat%NDOF
    nn = ndof*ndof
    if( associated(hecmat%AL) ) then
      deallocate(hecmat%AL                  ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%AU) ) then
      deallocate(hecmat%AU                  ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%B) ) then
      deallocate(hecmat%B                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%D) ) then
      deallocate(hecmat%D                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%X) ) then
      deallocate(hecmat%X                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
    if( associated(hecmat%A) ) then
      deallocate(hecmat%A                   ,stat=ierror)
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to deallocation error'
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
    endif
  end subroutine hecmat_finalize
 
  subroutine fstr_param_init( fstrPARAM, hecMESH )
    implicit none
    type(fstr_param) :: fstrPARAM
    type(hecmwst_local_mesh) :: hecMESH
    integer(kind=kint) :: i
    external fstr_sort_index
 
    fstrparam%solution_type = kststatic
    fstrparam%nlgeom        = .false.
    fstrparam%solver_method = ksmcg
    fstrparam%nlsolver_method = knsmnewton
 
    !!STATIC !HEAT
    fstrparam%analysis_n = 0
    fstrparam%ref_temp   = 0
 
    ! output control
    fstrparam%fg_echo       = koff
    fstrparam%fg_result     = koff
    fstrparam%fg_visual     = koff
 
    ! for heat ...
    fstrparam%fg_neutral    = koff
    fstrparam%fg_irres      = kno
    fstrparam%fg_iwres      = kno
    fstrparam%nrres         = 1
    fstrparam%nprint        = 1
 
    ! for couple
    fstrparam%fg_couple      = 0
    fstrparam%fg_couple_type = 0
    fstrparam%fg_couple_first= 0
    fstrparam%fg_couple_window= 0
 
    ! for restart control
    fstrparam%restart_version = 6
 
    ! for contact analysis
    fstrparam%contact_algo = kcaslagrange  ! default: Standard Lagrange
    fstrparam%augiter = 2                  ! default augmentation iteration for ALagrange
    fstrparam%fric_cone_follow = .true.    ! default: cone radius follows the applied normal force
    fstrparam%eps_fric_band = 0.d0         ! default: no stick/slip hysteresis band
 
    ! index table for global node ID sorting
 
    fstrparam%n_node = hecmesh%n_node;
    fstrparam%nn_internal = hecmesh%nn_internal;
    allocate( fstrparam%global_local_ID(2,hecmesh%nn_internal))
    do i = 1, hecmesh%nn_internal
      fstrparam%global_local_ID(1,i) = hecmesh%global_node_ID(i)
      fstrparam%global_local_ID(2,i) = i
    end do
    call fstr_sort_index( fstrparam%global_local_ID, hecmesh%nn_internal )
  end subroutine fstr_param_init
 
  subroutine fstr_initial_condition_finalize( cond )
    implicit none
    type( tinitialcondition ) :: cond
    if(associated(cond%intval)) deallocate(cond%intval)
    if(associated(cond%realval)) deallocate(cond%realval)
  end subroutine fstr_initial_condition_finalize
 
  subroutine fstr_param_finalize( fstrPARAM )
    implicit none
    type(fstr_param) :: fstrPARAM
    integer :: i
    if(associated(fstrparam%dtime)) deallocate(fstrparam%dtime)
    if(associated(fstrparam%etime)) deallocate(fstrparam%etime)
    if(associated(fstrparam%dtmin)) deallocate(fstrparam%dtmin)
    if(associated(fstrparam%delmax)) deallocate(fstrparam%delmax)
    if(associated(fstrparam%itmax)) deallocate(fstrparam%itmax)
    if(associated(fstrparam%eps)) deallocate(fstrparam%eps)
    if(associated(fstrparam%global_local_ID)) deallocate(fstrparam%global_local_ID)
 
    if(associated(fstrparam%contactparam)) deallocate(fstrparam%contactparam)
    if(associated(fstrparam%contact_if)) deallocate(fstrparam%contact_if)
 
    if(associated(fstrparam%ainc)) deallocate(fstrparam%ainc)
    if(associated(fstrparam%timepoints)) then
      do i=1, size(fstrparam%timepoints)
        if(associated(fstrparam%timepoints(i)%points)) &
          deallocate(fstrparam%timepoints(i)%points)
      end do
      deallocate(fstrparam%timepoints)
    end if
 
    if(associated(fstrparam%cnvparam)) deallocate(fstrparam%cnvparam)
  end subroutine fstr_param_finalize
 
  logical function fstr_isboundaryactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)    :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_isboundaryactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_isboundaryactive = isboundaryactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  logical function fstr_isloadactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)    :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_isloadactive = .true.
    if( cstep > 0 ) then
      if( .not. associated(fstrsolid%step_ctrl) ) return
      if( cstep>fstrsolid%nstep_tot ) return
      fstr_isloadactive = isloadactive( nbc, fstrsolid%step_ctrl(cstep) )
    else
      fstr_isloadactive = isloadactive( nbc, fstrsolid%step_ctrl_restart )
    endif
  end function
 
  logical function fstr_iscontactactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)     :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_iscontactactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_iscontactactive = iscontactactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  logical function fstr_isembedactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)     :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_isembedactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_isembedactive = iscontactactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  logical function fstr_iselemactivationactive( fstrSOLID, nbc, cstep )
    type(fstr_solid)    :: fstrsolid
    integer, intent(in) :: nbc
    integer, intent(in) :: cstep
    fstr_iselemactivationactive = .true.
    if( .not. associated(fstrsolid%step_ctrl) ) return
    if( cstep>fstrsolid%nstep_tot ) return
    fstr_iselemactivationactive = iselemactivationactive( nbc, fstrsolid%step_ctrl(cstep) )
  end function
 
  subroutine get_coordsys( cdsys_ID, hecMESH, fstrSOLID, coords, icel )
    integer, intent(in)             :: cdsys_ID
    integer, intent(in)             :: icel
    type(hecmwst_local_mesh)       :: hecMESH
    type(fstr_solid), intent(inout) :: fstrSOLID
    real(kind=kreal), intent(out)   :: coords(3,3)
    integer :: ik, local_nid, iSS
 
    coords = 0.d0
    if( cdsys_id>0 ) then
      if( iscoordneeds(g_localcoordsys(cdsys_id)) ) then
        coords=g_localcoordsys(cdsys_id)%CoordSys
      elseif( g_localcoordsys(cdsys_id)%sys_type == 12 ) then
        iss = hecmesh%elem_node_index(icel-1)
        !
        local_nid=g_localcoordsys(cdsys_id)%node_ID(1)
        ik = hecmesh%elem_node_item(iss+local_nid)
        coords(1,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
        local_nid=g_localcoordsys(cdsys_id)%node_ID(2)
        ik = hecmesh%elem_node_item(iss+local_nid)
        coords(2,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
        local_nid=g_localcoordsys(cdsys_id)%node_ID(3)
        ik = hecmesh%elem_node_item(iss+local_nid)
        coords(3,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
      else
        ik=g_localcoordsys(cdsys_id)%node_ID(1)
        coords(1,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
        ik=g_localcoordsys(cdsys_id)%node_ID(2)
        coords(2,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
        ik=g_localcoordsys(cdsys_id)%node_ID(3)
        if(ik>0) coords(3,:)= hecmesh%node(3*ik-2:3*ik)+fstrsolid%unode(3*ik-2:3*ik)  &
          + fstrsolid%dunode(3*ik-2:3*ik)
      endif
    endif
  end subroutine get_coordsys
 
  subroutine fstr_set_current_config_to_mesh(hecMESH,fstrSOLID,coord)
    implicit none
    type(hecmwst_local_mesh), intent(inout) :: hecMESH
    type (fstr_solid), intent(in) :: fstrSOLID
    real(kind=kreal), intent(out) :: coord(:)
    integer(kind=kint) :: i
    if(hecmesh%n_dof == 4) return
    do i = 1, hecmesh%nn_internal*min(hecmesh%n_dof,3)
      coord(i) = hecmesh%node(i)
      hecmesh%node(i) = coord(i)+fstrsolid%unode(i)+fstrsolid%dunode(i)
    enddo
  end subroutine fstr_set_current_config_to_mesh
 
  subroutine fstr_recover_initial_config_to_mesh(hecMESH,fstrSOLID,coord)
    implicit none
    type(hecmwst_local_mesh), intent(inout) :: hecMESH
    type (fstr_solid), intent(in) :: fstrSOLID
    real(kind=kreal), intent(in) :: coord(:)
    integer(kind=kint) :: i
    if(hecmesh%n_dof == 4) return
    do i = 1, hecmesh%nn_internal*min(hecmesh%n_dof,3)
      hecmesh%node(i) = coord(i)
    enddo
  end subroutine fstr_recover_initial_config_to_mesh
 
  subroutine fstr_solid_phys_zero(phys)
    implicit none
    type(fstr_solid_physic_val), pointer :: phys
    phys%STRAIN  = 0.0d0
    phys%STRESS  = 0.0d0
    phys%MISES   = 0.0d0
    phys%ESTRAIN = 0.0d0
    phys%ESTRESS = 0.0d0
    phys%EMISES  = 0.0d0
    phys%EPLSTRAIN  = 0.0d0
    phys%ENQM    = 0.0d0
  end subroutine fstr_solid_phys_zero
 
  subroutine fstr_solid_phys_clear(fstrSOLID)
    implicit none
    type (fstr_solid)         :: fstrSOLID
    integer(kind=kint) :: i
 
    if (associated(fstrsolid%SOLID)) then
      call fstr_solid_phys_zero(fstrsolid%SOLID)
    end if
    if (associated(fstrsolid%SHELL)) then
      call fstr_solid_phys_zero(fstrsolid%SHELL)
      do i=1,fstrsolid%max_lyr
        call fstr_solid_phys_zero(fstrsolid%SHELL%LAYER(i)%PLUS)
        call fstr_solid_phys_zero(fstrsolid%SHELL%LAYER(i)%MINUS)
      end do
    end if
    if (associated(fstrsolid%BEAM)) then
      call fstr_solid_phys_zero(fstrsolid%BEAM)
    end if
  end subroutine fstr_solid_phys_clear
 
  subroutine fstr_get_amplitude(hecMESH, fstrSOLID, cstep, jj_n_amp, time, value)
    type ( hecmwST_local_mesh ), intent(in) :: hecMESH
    type ( fstr_solid         ), intent(in) :: fstrSOLID
    integer(kind=kint), intent(in)          :: cstep
    integer(kind=kint), intent(in)          :: jj_n_amp
    real(kind=kreal), intent(in)            :: time       
    real(kind=kreal), intent(inout)         :: value      
 
    real(kind=kreal) :: t_eval
 
    if( jj_n_amp <= 0 ) then  ! Amplitude not defined
      if(myrank == 0) then
        write(imsg,*) 'internal error: amplitude table not defined'
      endif
      call hecmw_abort( hecmw_comm_get_comm() )
    endif
 
    if( hecmesh%amp%amp_type_time(jj_n_amp) == hecmw_amp_typetime_total ) then
      t_eval = time
    else
      t_eval = time - fstrsolid%step_ctrl(cstep)%starttime
    endif
 
    call hecmw_get_amplitude_value(hecmesh%amp, jj_n_amp, t_eval, value)
 
  end subroutine fstr_get_amplitude
 
end module m_fstr