fstr_ctrl_common.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module fstr_ctrl_common
  use m_fstr
  use hecmw
  use m_timepoint
  use fstr_ctrl_util_f
 
  implicit none
 
contains
 
  function fstr_ctrl_get_solution( ctrl, type, nlgeom )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: type
    logical            :: nlgeom
    integer(kind=kint) :: fstr_ctrl_get_solution
 
    integer(kind=kint) :: ipt
    character(len=80) :: s
 
    fstr_ctrl_get_solution = -1
 
    s = 'ELEMCHECK,STATIC,EIGEN,HEAT,DYNAMIC,NLSTATIC,STATICEIGEN,NZPROF '
    if( fstr_ctrl_get_param_ex( ctrl,      'TYPE ',     s,    1,   'P',  type )/= 0) return
    type = type -1
 
    ipt=0
    if( fstr_ctrl_get_param_ex( ctrl,    'NONLINEAR ',  '# ',    0,   'E',   ipt )/= 0) return
    if( ipt/=0 .and. ( type == kststatic .or. type == kstdynamic )) nlgeom = .true.
 
    if( type == 5 ) then !if type == NLSTATIC
      type = kststatic
      nlgeom = .true.
    end if
    if( type == kststaticeigen ) nlgeom = .true.
 
    fstr_ctrl_get_solution = 0
  end function fstr_ctrl_get_solution
 
  function fstr_ctrl_get_nonlinear_solver( ctrl, method )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: method
    integer(kind=kint) :: fstr_ctrl_get_nonlinear_solver
 
    character(len=80) :: s
 
    fstr_ctrl_get_nonlinear_solver = -1
 
    s = 'NEWTON,QUASINEWTON '
    if( fstr_ctrl_get_param_ex( ctrl,      'METHOD ',     s,    1,   'P',  method )/= 0) return
    fstr_ctrl_get_nonlinear_solver = 0
  end function fstr_ctrl_get_nonlinear_solver
 
  function fstr_ctrl_get_solver( ctrl, method, precond, nset, iterlog, timelog, steplog, nier, &
      iterpremax, nrest, nBFGS, scaling, &
      dumptype, dumpexit, usejad, ncolor_in, mpc_method, estcond, method2, recyclepre, &
      solver_opt, contact_elim, &
      resid, sigma_diag, omega, thresh, filter, solver_ropt, loglevel, &
      matvec_impl, precond_impl, recompute_residual )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: method
    integer(kind=kint) :: precond
    integer(kind=kint) :: nset
    integer(kind=kint) :: iterlog
    integer(kind=kint) :: timelog
    integer(kind=kint) :: steplog
    integer(kind=kint) :: nier
    integer(kind=kint) :: iterpremax
    integer(kind=kint) :: nrest
    integer(kind=kint) :: nbfgs
    integer(kind=kint) :: scaling
    integer(kind=kint) :: dumptype
    integer(kind=kint) :: dumpexit
    integer(kind=kint) :: usejad
    integer(kind=kint) :: ncolor_in
    integer(kind=kint) :: mpc_method
    integer(kind=kint) :: estcond
    integer(kind=kint) :: method2
    integer(kind=kint) :: recyclepre
    integer(kind=kint) :: solver_opt(10)
    integer(kind=kint) :: contact_elim
    real(kind=kreal) :: resid
    real(kind=kreal) :: sigma_diag
    real(kind=kreal) :: omega
    real(kind=kreal) :: thresh
    real(kind=kreal) :: filter
    real(kind=kreal) :: solver_ropt(10)
    integer(kind=kint) :: loglevel
    integer(kind=kint) :: matvec_impl
    integer(kind=kint) :: precond_impl
    integer(kind=kint) :: recompute_residual
    integer(kind=kint) :: fstr_ctrl_get_solver
 
    character(120) :: mlist = &
      '1,2,3,4,101,CG,BiCGSTAB,GMRES,GPBiCG,GMRESR,GMRESREN,CR,PipeCG,GroppCG,DIRECT,DIRECTmkl,DIRECTlag,MUMPS,MKL,DIRECTmf '
    !character(92) :: mlist = '1,2,3,4,5,101,CG,BiCGSTAB,GMRES,GPBiCG,DIRECT,DIRECTmkl,DIRECTlag,MUMPS,MKL '
    character(24) :: dlist = '0,1,2,3,NONE,MM,CSR,BSR '
 
    integer(kind=kint) :: number_number = 5
    integer(kind=kint) :: indirect_number = 9
    integer(kind=kint) :: iter, time, sclg, dmpt, dmpx, usjd, step, mtxfmt
 
    fstr_ctrl_get_solver = -1
 
    iter = iterlog+1
    time = timelog+1
    step = steplog+1
    sclg = scaling+1
    dmpt = dumptype+1
    dmpx = dumpexit+1
    usjd = usejad+1
    ! 0 = MATRIXFORMAT absent, leaving USEJAD and the build-time default in effect
    mtxfmt = 0
    !* parameter in header line -----------------------------------------------------------------*!
 
    ! JP-0
    if( fstr_ctrl_get_param_ex( ctrl, 'METHOD ',   mlist,              1,   'P',   method  ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'PRECOND ', '1,2,3,4,5,6,7,8,9,10,11,12,20,21,22,30,31,32 ' ,0, 'I', precond ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'NSET ',    '0,-1,+1 ',          0,   'I',   nset    ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'ITERLOG ', 'NO,YES ',           0,   'P',   iter ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'TIMELOG ', 'NO,YES,VERBOSE ',   0,   'P',   time ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'STEPLOG ', 'NO,YES ',           0,   'P',   step ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'SCALING ', 'NO,YES ',           0,   'P',   sclg ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'DUMPTYPE ', dlist,              0,   'P',   dmpt ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'DUMPEXIT ','NO,YES ',           0,   'P',   dmpx ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'USEJAD '  ,'NO,YES ',           0,   'P',   usjd ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'MATRIXFORMAT ','BSR,JAD,CSR,SBLAS ', 0, 'P', mtxfmt ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'MPCMETHOD ','# ',               0, 'I',mpc_method) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'ESTCOND '  ,'# ',               0,   'I',estcond ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'RECOMPUTE_RESIDUAL ','# ',      0,   'I',recompute_residual ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'METHOD2 ',  mlist,              0,   'P',   method2 ) /= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'CONTACT_ELIM ','# ',            0,   'I',contact_elim ) /= 0) return
    ! diagnostic verbosity, independent of TIMELOG.  -1 = "unset": each consumer
    ! (ML / SA-AMG) then falls back to its previous default, so behavior is
    ! unchanged unless LOGLEVEL is given explicitly.
    loglevel = -1
    if( fstr_ctrl_get_param_ex( ctrl, 'LOGLEVEL ','# ',               0,   'I',loglevel ) /= 0) return
    ! JP-1
    if( method > number_number ) then  ! JP-2
      method = method - number_number
      if( method > indirect_number ) then
        ! JP-3
        method = method - indirect_number + 100
        if( method == 103 ) method = 101 ! DIRECTlag => DIRECT
        if( method == 105 ) method = 102 ! MKL => DIRECTmkl
      end if
    end if
    if( method2 > number_number ) then  ! JP-2
      method2 = method2 - number_number
      if( method2 > indirect_number ) then
        ! JP-3
        method2 = method2 - indirect_number + 100
      end if
    end if
 
    dumptype = dmpt - 1
    if( dumptype >= 4 ) then
      dumptype = dumptype - 4
    end if
 
    !* data --------------------------------------------------------------------------------------- *!
    ! JP-4
    if( fstr_ctrl_get_data_ex( ctrl, 1,   'iiiiii ', nier, iterpremax, nrest, ncolor_in, recyclepre, nbfgs )/= 0) return
    if( fstr_ctrl_get_data_ex( ctrl, 2,   'rrr ', resid, sigma_diag, omega )/= 0) return
 
    ! the direct methods own their data lines and come first: a direct solve never uses
    ! PRECOND, so a leftover PRECOND of an iterative control file must not route the
    ! option lines into a preconditioner branch (issue #930)
    if( method == 101 ) then
      if( fstr_ctrl_get_data_ex( ctrl, 3, 'i ', solver_opt(1) )/= 0) return
    else if( method == 106 ) then
      ! DIRECTmf options.  Two optional data lines; trailing entries may be omitted.
      ! 0 = use the built-in default for each.  Integer slot 1 is the ordering, the
      ! same as the single integer of the DIRECT (method=101) data line.
      !  line 3 (integers):
      !    1 ordering (as DIRECT), 2 mode (0=auto/1=force LDLt/2=force LU),
      !    3 BLR (0=off/1=on), 4 max refinement steps, 5 tile size (DOFs),
      !    6 amalgamation width (columns), 7 BLR skip distance eta (-1=off),
      !    8 BLR rank reuse (0=off/1=on), 9-10 RESERVED
      !  line 4 (reals):
      !    1 BLR epsilon, 2 refinement stopping residual, 3 pivot threshold u,
      !    4 zero pivot factor, 5 BLR gain cap beta, 6-8 RESERVED
      solver_opt(1:10) = 0
      solver_ropt(1:10) = 0.0d0
      if( fstr_ctrl_get_data_ex( ctrl, 3, 'iiiiiiiiii ', &
           solver_opt(1), solver_opt(2), solver_opt(3), solver_opt(4), solver_opt(5), &
           solver_opt(6), solver_opt(7), solver_opt(8), solver_opt(9), solver_opt(10) )/= 0) solver_opt(1:10) = 0
      if( fstr_ctrl_get_data_ex( ctrl, 4, 'rrrrrrrr ', &
           solver_ropt(1), solver_ropt(2), solver_ropt(3), solver_ropt(4), solver_ropt(5), &
           solver_ropt(6), solver_ropt(7), solver_ropt(8) )/= 0) solver_ropt(1:10) = 0.0d0
    else if( precond == 20 .or. precond == 21) then
      if( fstr_ctrl_get_data_ex( ctrl, 3, 'rr ', thresh, filter)/= 0) return
    else if( precond == 5 ) then
      if( fstr_ctrl_get_data_ex( ctrl, 3, 'iiiiiiiiii ', &
           solver_opt(1), solver_opt(2), solver_opt(3), solver_opt(4), solver_opt(5), &
           solver_opt(6), solver_opt(7), solver_opt(8), solver_opt(9), solver_opt(10) )/= 0) return
    else if( precond == 22 ) then
      ! SA-AMG options.  Two optional data lines; trailing entries may be omitted.
      ! 0 = use the built-in default for each.  Integer slots 1-7 mirror the ML
      ! (PRECOND=5) layout so an ML option line carries over (hecmw_ML_wrapper.c
      ! reads opt[0..6] = slots 1-7); slots 8-10 are SA-AMG-specific.  ML does not
      ! read the real line at all.
      ! (verbose has no slot here: enable diagnostics via !SOLVER LOGLEVEL>=1.)
      !  line 3 (integers):
      !    1 coarse_solver (0=auto/1=smoother/2=dense/3=mumps), 2 smoother (0/1=Chebyshev only),
      !    3 cycle (0=default(W)/1=V/2=W), 4 max_level,
      !    5 RESERVED (ML CoarsenScheme: ignored with a warning), 6 cheb_deg, 7 coarse_size,
      !    8 max_size, 9 galerkin_lowmem (0=default(2-stage,faster),
      !               >0=low-memory (fuse finest level only)), 10 RESERVED
      !  line 4 (reals):    theta, cheb_alpha, safety,
      !    taper_k (coarsening taper: 0=default(100), >0=use as K, <0=disable),
      !    agg_order (aggregation seed ordering: 0=default(BFS), <0=natural/legacy,
      !               1..4=explicit mode),
      !    min_size, verify (0=off/>0=on), dump_vtk (0=off/>0=on)
      solver_opt(1:10) = 0
      solver_ropt(1:10) = 0.0d0
      if( fstr_ctrl_get_data_ex( ctrl, 3, 'iiiiiiiiii ', &
           solver_opt(1), solver_opt(2), solver_opt(3), solver_opt(4), solver_opt(5), &
           solver_opt(6), solver_opt(7), solver_opt(8), solver_opt(9), solver_opt(10) )/= 0) solver_opt(1:10) = 0
      if( fstr_ctrl_get_data_ex( ctrl, 4, 'rrrrrrrr ', &
           solver_ropt(1), solver_ropt(2), solver_ropt(3), solver_ropt(4), solver_ropt(5), &
           solver_ropt(6), solver_ropt(7), solver_ropt(8) )/= 0) solver_ropt(1:10) = 0.0d0
    end if
 
    iterlog = iter -1
    timelog = time -1
    steplog = step -1
    scaling = sclg -1
    dumpexit = dmpx -1
    usejad = usjd -1
 
    ! MATRIXFORMAT names one storage format for the whole linear solve, so it sets the
    ! matvec and the preconditioner together.  JAD is one of the formats, which makes
    ! USEJAD the older spelling of MATRIXFORMAT=JAD; giving both lets MATRIXFORMAT win.
    ! JAD replaces the matvec only, so the preconditioner keeps the BSR implementation.
    select case (mtxfmt)
      case (1)
        usejad = kno
        matvec_impl  = hecmw_matvec_impl_bsr
        precond_impl = hecmw_precond_impl_bsr
      case (2)
        usejad = kyes
        matvec_impl  = hecmw_matvec_impl_bsr
        precond_impl = hecmw_precond_impl_bsr
      case (3)
        usejad = kno
        matvec_impl  = hecmw_matvec_impl_csr
        precond_impl = hecmw_precond_impl_csr
      case (4)
        usejad = kno
        matvec_impl  = hecmw_matvec_impl_sblas
        precond_impl = hecmw_precond_impl_sblas
    end select
 
    fstr_ctrl_get_solver = 0
 
  end function fstr_ctrl_get_solver
 
 
  function fstr_ctrl_get_step( ctrl, amp, iproc )
    integer(kind=kint) :: ctrl
    character(len=HECMW_NAME_LEN) :: amp
    integer(kind=kint) :: iproc
    integer(kind=kint) :: fstr_ctrl_get_step
 
    integer(kind=kint) :: ipt = 0
    integer(kind=kint) :: ip = 0
 
    fstr_ctrl_get_step = -1
 
    if( fstr_ctrl_get_param_ex( ctrl, 'AMP ',     '# ',  0, 'S', amp )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ',   'STANDARD,NLGEOM ', 0, 'P',   ipt    )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'NLGEOM ',  '# ',           0,    'E',   ip     )/= 0) return
 
    if( ipt == 2 .or. ip == 1 ) iproc = 1
 
    fstr_ctrl_get_step = 0
 
  end function fstr_ctrl_get_step
 
  logical function fstr_ctrl_get_istep( ctrl, hecMESH, steps, tpname, apname, cpname )
    use fstr_setup_util
    use m_step
    integer(kind=kint), intent(in)        :: ctrl
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type(step_info), intent(inout)        :: steps
    character(len=*), intent(out)         :: tpname
    character(len=*), intent(out)         :: apname
    character(len=*), intent(out)         :: cpname
 
    character(len=HECMW_NAME_LEN) :: data_fmt,ss, data_fmt1
    character(len=HECMW_NAME_LEN) :: amp
    character(len=HECMW_NAME_LEN) :: header_name
    integer(kind=kint) :: bcid
    integer(kind=kint) :: i, n, sn, ierr
    integer(kind=kint) :: bc_n, load_n, contact_n, elemact_n
    real(kind=kreal) :: fn, f1, f2, f3
 
    fstr_ctrl_get_istep = .false.
 
    write(ss,*)  hecmw_name_len
    write( data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)), 'I '
    write( data_fmt1, '(a,a,a)') 'S', trim(adjustl(ss)),'rrr '
 
    steps%solution = stepstatic
    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ',   'STATIC,VISCO ', 0, 'P', steps%solution )/= 0) return
    steps%inc_type = stepfixedinc
    if( fstr_ctrl_get_param_ex( ctrl, 'INC_TYPE ', 'FIXED,AUTO ', 0, 'P', steps%inc_type )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'SUBSTEPS ',  '# ',  0, 'I', steps%num_substep )/= 0) return
    steps%initdt = 1.d0/steps%num_substep
    if( fstr_ctrl_get_param_ex( ctrl, 'ITMAX ',  '# ',  0, 'I', steps%max_iter )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'MAXITER ',  '# ',  0, 'I', steps%max_iter )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'MAXCONTITER ',  '# ',  0, 'I', steps%max_contiter )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'CONVERG ',  '# ',  0, 'R', steps%converg )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'CONVERG_LAG ',  '# ',  0, 'R', steps%converg_lag )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'CONVERG_DDISP ',  '# ',  0, 'R', steps%converg_ddisp )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'MAXRES ',  '# ',  0, 'R', steps%maxres )/= 0) return
    amp = ""
    if( fstr_ctrl_get_param_ex( ctrl, 'AMP ',  '# ',  0, 'S', amp )/= 0) return
    if( len( trim(amp) )>0 ) then
      call amp_name_to_id( hecmesh, '!STEP', amp, steps%amp_id )
    endif
    ! AMPLITUDE = RAMP | STEP : default sub-step application of prescribed displacements.
    ! Default preserves current behavior (static=RAMP, visco=STEP).
    steps%amp_default_type = stepampramp
    if( steps%solution == stepvisco ) steps%amp_default_type = stepampstep
    if( fstr_ctrl_get_param_ex( ctrl, 'AMPLITUDE ', 'RAMP,STEP ', 0, 'P', steps%amp_default_type )/= 0) return
    tpname=""
    if( fstr_ctrl_get_param_ex( ctrl, 'TIMEPOINTS ',  '# ',  0, 'S', tpname )/= 0) return
    apname=""
    if( fstr_ctrl_get_param_ex( ctrl, 'AUTOINCPARAM ',  '# ',  0, 'S', apname )/= 0) return
    cpname=""
    if( fstr_ctrl_get_param_ex( ctrl, 'CONVERGPARAM ',  '# ',  0, 'S', cpname )/= 0) return
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) then
      fstr_ctrl_get_istep = .true.;  return
    endif
 
    f2 = steps%mindt
    f3 = steps%maxdt
    if( fstr_ctrl_get_data_ex( ctrl, 1, data_fmt1, ss, f1, f2, f3  )/= 0) return
    read( ss, * , iostat=ierr ) fn
    sn=1
    if( ierr==0 ) then
      steps%initdt = fn
      steps%elapsetime = f1
      if( steps%inc_type == stepautoinc ) then
        steps%mindt = min(f2,steps%initdt)
        steps%maxdt = f3
      endif
      steps%num_substep = max(int((f1+0.999999999d0*fn)/fn),steps%num_substep)
      !if( mod(f1,fn)/=0 ) steps%num_substep =steps%num_substep+1
      sn = 2
    endif
 
    bc_n = 0
    load_n = 0
    contact_n = 0
    elemact_n = 0
    do i=sn,n
      if( fstr_ctrl_get_data_ex( ctrl, i, data_fmt, header_name, bcid  )/= 0) return
      if( trim(header_name) == 'BOUNDARY' ) then
        bc_n = bc_n + 1
      else if( trim(header_name) == 'LOAD' ) then
        load_n = load_n +1
      else if( trim(header_name) == 'CONTACT' ) then
        contact_n = contact_n+1
      else if( trim(header_name) == 'ELEMACT' ) then
        elemact_n = elemact_n+1
      else if( trim(header_name) == 'TEMPERATURE' ) then
        !   steps%Temperature = .true.
      endif
    end do
 
    if( bc_n>0 ) allocate( steps%Boundary(bc_n) )
    if( load_n>0 ) allocate( steps%Load(load_n) )
    if( contact_n>0 ) allocate( steps%Contact(contact_n) )
    if( elemact_n>0 ) allocate( steps%ElemActivation(elemact_n) )
 
    bc_n = 0
    load_n = 0
    contact_n = 0
    elemact_n = 0
    do i=sn,n
      if( fstr_ctrl_get_data_ex( ctrl, i, data_fmt, header_name, bcid  )/= 0) return
      if( trim(header_name) == 'BOUNDARY' ) then
        bc_n = bc_n + 1
        steps%Boundary(bc_n) = bcid
      else if( trim(header_name) == 'LOAD' ) then
        load_n = load_n +1
        steps%Load(load_n) = bcid
      else if( trim(header_name) == 'CONTACT' ) then
        contact_n = contact_n+1
        steps%Contact(contact_n) = bcid
      else if( trim(header_name) == 'ELEMACT' ) then
        elemact_n = elemact_n+1
        steps%ElemActivation(elemact_n) = bcid
      endif
    end do
 
    fstr_ctrl_get_istep = .true.
  end function fstr_ctrl_get_istep
 
  integer function fstr_ctrl_get_section( ctrl, hecMESH, sections )
    use fstr_setup_util
    integer(kind=kint), intent(in)           :: ctrl
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (tsection), pointer, intent(inout)  :: sections(:)
 
    integer(kind=kint)            :: j, k, sect_id, elemopt
    integer(kind=kint),save       :: cache = 1
    character(len=HECMW_NAME_LEN) :: sect_orien
    character(19) :: form341list = 'FI,SELECTIVE_ESNS '
    character(19) :: form361list = 'FI,BBAR,IC,FBAR,UP '
    character(19) :: form611list = 'EULER,TIMOSHENKO '
 
    fstr_ctrl_get_section = -1
 
    if( fstr_ctrl_get_param_ex( ctrl, 'SECNUM ',  '# ',  1, 'I', sect_id )/= 0) return
    if( sect_id > hecmesh%section%n_sect ) return
 
    elemopt = 0
    if( fstr_ctrl_get_param_ex( ctrl, 'FORM341 ',   form341list, 0, 'P', elemopt )/= 0) return
    if( elemopt > 0 ) sections(sect_id)%elemopt341 = elemopt
 
    elemopt = 0
    if( fstr_ctrl_get_param_ex( ctrl, 'FORM361 ',   form361list, 0, 'P', elemopt )/= 0) return
    if( elemopt > 0 ) sections(sect_id)%elemopt361 = elemopt
 
    elemopt = 0
    if( fstr_ctrl_get_param_ex( ctrl, 'FORM611 ', form611list, 0, 'P', elemopt )/= 0) return
    if( elemopt > 0 ) sections(sect_id)%elemopt611 = elemopt
 
    ! sectional orientation ID
    hecmesh%section%sect_orien_ID(sect_id) = -1
    if( fstr_ctrl_get_param_ex( ctrl, 'ORIENTATION ',  '# ',  0, 'S', sect_orien )/= 0) return
 
    if( associated(g_localcoordsys) ) then
      k = size(g_localcoordsys)
 
      if(cache < k)then
        if( sect_orien == g_localcoordsys(cache)%sys_name ) then
          hecmesh%section%sect_orien_ID(sect_id) = cache
          cache = cache + 1
          fstr_ctrl_get_section = 0
          return
        endif
      endif
 
      do j=1, k
        if( sect_orien == g_localcoordsys(j)%sys_name ) then
          hecmesh%section%sect_orien_ID(sect_id) = j
          cache = j + 1
          exit
        endif
      enddo
    endif
 
    fstr_ctrl_get_section = 0
 
  end function fstr_ctrl_get_section
 
 
  function fstr_ctrl_get_write( ctrl, res, visual, femap )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: res
    integer(kind=kint) :: visual
    integer(kind=kint) :: femap
    integer(kind=kint) :: fstr_ctrl_get_write
 
    fstr_ctrl_get_write = -1
 
    ! JP-6
    if( fstr_ctrl_get_param_ex( ctrl, 'RESULT ',  '# ',    0,   'E',   res    )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'VISUAL ',  '# ',    0,   'E',   visual )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'FEMAP ',   '# ',    0,   'E',   femap  )/= 0) return
 
    fstr_ctrl_get_write = 0
 
  end function fstr_ctrl_get_write
 
  function fstr_ctrl_get_echo( ctrl, echo )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: echo
    integer(kind=kint) :: fstr_ctrl_get_echo
 
    echo = kon;
 
    fstr_ctrl_get_echo = 0
 
  end function fstr_ctrl_get_echo
 
  function fstr_ctrl_get_couple( ctrl, fg_type, fg_first, fg_window, surf_id, surf_id_len )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: fg_type
    integer(kind=kint) :: fg_first
    integer(kind=kint) :: fg_window
    character(len=HECMW_NAME_LEN) :: surf_id(:)
    integer(kind=kint) :: surf_id_len
    integer(kind=kint) :: fstr_ctrl_get_couple
 
    character(len=HECMW_NAME_LEN) :: data_fmt,ss
    write(ss,*)  surf_id_len
    write(data_fmt,'(a,a,a)') 'S',trim(adjustl(ss)),' '
 
    fstr_ctrl_get_couple = -1
    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ', '1,2,3,4,5,6 ', 0, 'I', fg_type )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'ISTEP ', '# ', 0, 'I', fg_first )/= 0) return
    if( fstr_ctrl_get_param_ex( ctrl, 'WINDOW ', '# ', 0, 'I', fg_window )/= 0) return
 
    fstr_ctrl_get_couple = &
      fstr_ctrl_get_data_array_ex( ctrl, data_fmt, surf_id )
 
  end function fstr_ctrl_get_couple
 
  function fstr_ctrl_get_mpc( ctrl, penalty )
    integer(kind=kint), intent(in) :: ctrl
    real(kind=kreal), intent(inout)  :: penalty   
    integer(kind=kint) :: fstr_ctrl_get_mpc
 
    fstr_ctrl_get_mpc = fstr_ctrl_get_data_ex( ctrl, 1,   'r ', penalty )
    if( penalty <= 1.0 ) then
      if (myrank == 0) then
        write(imsg,*) "Warging : !MPC : too small penalty: ", penalty
        write(*,*) "Warging : !MPC : too small penalty: ", penalty
      endif
    endif
 
  end function fstr_ctrl_get_mpc
 
  logical function fstr_ctrl_get_outitem( ctrl, hecMESH, outinfo )
    use hecmw_setup_util, only: hecmw_str2index
    use fstr_setup_util
    use m_out
    integer(kind=kint), intent(in)        :: ctrl
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type( output_info ), intent(inout)    :: outinfo
 
    integer(kind=kint) :: rcode, ipos
    integer(kind=kint) :: n, i, j
    character(len=HECMW_NAME_LEN) :: data_fmt, ss
    character(len=HECMW_NAME_LEN), allocatable :: header_name(:), onoff(:), vtype(:)
 
    write( ss, * )  hecmw_name_len
    write( data_fmt, '(a,a,a,a,a)') 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), ' '
    !  write( data_fmt, '(a,a,a,a,a,a,a)') 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), ' '
 
    fstr_ctrl_get_outitem = .false.
 
    outinfo%grp_id_name = "ALL"
    rcode = fstr_ctrl_get_param_ex( ctrl, 'GROUP ', '# ', 0, 'S', outinfo%grp_id_name )
    ipos = 0
    rcode = fstr_ctrl_get_param_ex( ctrl, 'ACTION ', 'SUM ', 0, 'P', ipos )
    outinfo%actn = ipos
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    if( n == 0 ) return
    allocate( header_name(n), onoff(n), vtype(n) )
    header_name(:) = ""; vtype(:) = "";  onoff(:) = ""
    rcode = fstr_ctrl_get_data_array_ex( ctrl, data_fmt, header_name, onoff )
    !  rcode = fstr_ctrl_get_data_array_ex( ctrl, data_fmt, header_name, onoff, vtype )
 
    do i = 1, n
      do j = 1, outinfo%num_items
        if( trim(header_name(i)) == outinfo%keyWord(j) ) then
          outinfo%on(j) = .true.
          if( trim(onoff(i)) == 'OFF' ) outinfo%on(j) = .false.
          if( len( trim(vtype(i)) )>0 ) then
            if( hecmw_str2index( vtype(i), ipos ) ) then
              outinfo%vtype(j) = ipos
            else if( trim(vtype(i)) == "SCALER" ) then
              outinfo%vtype(j) = -1
            else if( trim(vtype(i)) == "VECTOR" ) then
              outinfo%vtype(j) = -2
            else if( trim(vtype(i)) == "SYMTENSOR" ) then
              outinfo%vtype(j) = -3
            else if( trim(vtype(i)) == "TENSOR" ) then
              outinfo%vtype(j) = -4
            endif
          endif
        endif
      enddo
    enddo
 
    deallocate( header_name, onoff, vtype )
    fstr_ctrl_get_outitem = .true.
 
  end function fstr_ctrl_get_outitem
 
  function fstr_ctrl_get_contactalgo( ctrl, algo, augiter, conefollow, eps_fric_band )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: algo
    integer(kind=kint) :: augiter
    logical            :: conefollow
    real(kind=kreal)   :: eps_fric_band
    integer(kind=kint) :: fstr_ctrl_get_contactalgo
 
    integer(kind=kint) :: rcode, icone
    character(len=80) :: s
    s = 'SLAGRANGE,ALAGRANGE '
    rcode = fstr_ctrl_get_param_ex( ctrl, 'TYPE ', s, 0, 'P', algo )
    if( rcode /= 0 ) then
      fstr_ctrl_get_contactalgo = rcode
      return
    endif
    ! FROZEN keeps the ALagrange friction cone at the multiplier of the last augmentation,
    ! FOLLOW lets it follow the normal force the contact element applies
    icone = 2   ! default value
    s = 'FROZEN,FOLLOW '
    rcode = fstr_ctrl_get_param_ex( ctrl, 'FRICTION_CONE ', s, 0, 'P', icone )
    if( rcode /= 0 ) then
      fstr_ctrl_get_contactalgo = rcode
      return
    endif
    conefollow = ( icone == 2 )
    rcode = fstr_ctrl_get_param_ex( ctrl, 'AUGITER ', '# ', 0, 'I', augiter )
    ! Hysteresis half-band of the stick/slip state switch (0 = no band = legacy behavior).
    ! Read here, at the same level as AUGITER, because it controls the augmentation loop's
    ! state machine rather than the geometric tolerances of !CONTACT_PARAM.
    rcode = fstr_ctrl_get_param_ex( ctrl, 'EPS_FRIC_BAND ', '# ', 0, 'R', eps_fric_band )
    if( eps_fric_band<0.d0 .or. 1.d0<=eps_fric_band ) then
      write(*,*) 'fstr control file error : !CONTACT_ALGO : EPS_FRIC_BAND must be 0 <= EPS_FRIC_BAND < 1.'
      write(ilog,*) 'fstr control file error : !CONTACT_ALGO : EPS_FRIC_BAND must be 0 <= EPS_FRIC_BAND < 1.'
      fstr_ctrl_get_contactalgo = -1
      return
    endif
    fstr_ctrl_get_contactalgo = 0
  end function fstr_ctrl_get_contactalgo
 
  logical function fstr_ctrl_get_contact( ctrl, n, contact, np, tp, ntol, ttol, ctAlgo, cpname, smoothing )
    use fstr_setup_util
    integer(kind=kint), intent(in)    :: ctrl
    integer(kind=kint), intent(in)    :: n
    integer(kind=kint), intent(in)    :: ctalgo
    type(tcontact), intent(out)       :: contact(n)
    real(kind=kreal), intent(out)      :: np             
    real(kind=kreal), intent(out)      :: tp             
    real(kind=kreal), intent(out)      :: ntol           
    real(kind=kreal), intent(out)      :: ttol           
    character(len=*), intent(out)      :: cpname
    integer(kind=kint), intent(out)    :: smoothing
 
    integer           :: rcode, mortar, expansion
    character(len=30) :: s1 = 'TIED,GLUED,SSLID,FSLID '
    character(len=HECMW_NAME_LEN) :: data_fmt,ss
    character(len=HECMW_NAME_LEN) :: cp_name(n)
    real(kind=kreal)  :: fcoeff(n),tpenalty(n)
    real(kind=kreal)  :: damp_alpha, damp_gact
 
    write(ss,*)  hecmw_name_len
 
    fstr_ctrl_get_contact = .false.
    contact(1)%ctype = 1   ! pure slave-master contact; default value
    contact(1)%algtype = contactsslid ! small sliding contact; default value
    rcode = fstr_ctrl_get_param_ex( ctrl, 'INTERACTION ', s1, 0, 'P', contact(1)%algtype )
    if( contact(1)%algtype==contactglued ) contact(1)%algtype=contactfslid  ! not complemented yet
    if( fstr_ctrl_get_param_ex( ctrl, 'GRPID ', '# ', 1, 'I', contact(1)%group )/=0) return
    smoothing = kcsnone + 1
    if( fstr_ctrl_get_param_ex( ctrl, 'SMOOTHING ', 'NONE,NAGATA ', 0, 'P', smoothing ) /= 0 ) return
    smoothing = smoothing - 1
    do rcode=2,n
      contact(rcode)%ctype = contact(1)%ctype
      contact(rcode)%group = contact(1)%group
      contact(rcode)%algtype = contact(1)%algtype
    end do
 
    tpenalty = 0.5d0
 
    if( contact(1)%algtype==contactsslid .or. contact(1)%algtype==contactfslid ) then
      write( data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)),'Rr '
      if( fstr_ctrl_get_data_array_ex( ctrl, data_fmt, cp_name, fcoeff, tpenalty ) /= 0 ) return
      do rcode=1,n
        contact(rcode)%pair_name = cp_name(rcode)
        contact(rcode)%fcoeff = fcoeff(rcode)
        contact(rcode)%nPenalty = 5.0d0
        contact(rcode)%tPenalty = tpenalty(rcode)
        contact(rcode)%refStiff = 1.d0
        contact(rcode)%damp_alpha = 0.0d0
        contact(rcode)%damp_gact = 0.0d0
      enddo
    else if( contact(1)%algtype==contacttied ) then
      write( data_fmt, '(a,a)') 'S', trim(adjustl(ss))
      if(  fstr_ctrl_get_data_array_ex( ctrl, data_fmt, cp_name ) /= 0 ) return
      do rcode=1,n
        contact(rcode)%pair_name = cp_name(rcode)
        contact(rcode)%nPenalty = 5.0d0
        contact(rcode)%fcoeff = 0.d0
        contact(rcode)%tPenalty = 1.d0
        contact(rcode)%damp_alpha = 0.0d0
        contact(rcode)%damp_gact = 0.0d0
      enddo
    endif
 
    np = 0.d0;  tp=0.d0
    ntol = 0.d0;  ttol=0.d0
    damp_alpha = 0.0d0;  damp_gact = 0.0d0
    if( fstr_ctrl_get_param_ex( ctrl, 'NPENALTY ',  '# ',  0, 'R', np ) /= 0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'TPENALTY ', '# ', 0, 'R', tp ) /= 0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'NTOL ',  '# ',  0, 'R', ntol ) /= 0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'TTOL ', '# ', 0, 'R', ttol ) /= 0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'DAMP_ALPHA ', '# ', 0, 'R', damp_alpha ) /= 0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'DAMP_GACT ',  '# ', 0, 'R', damp_gact  ) /= 0 ) return
    cpname=""
    if( fstr_ctrl_get_param_ex( ctrl, 'CONTACTPARAM ',  '# ',  0, 'S', cpname )/= 0) return
 
    ! MORTAR = NO / YES -> 1 / 2 = CONTACTN2S / CONTACTS2S
    mortar = 1
    if( fstr_ctrl_get_param_ex( ctrl, 'MORTAR ','NO,YES ', 0, 'P', mortar ) /= 0) return
    contact%method = mortar
 
    ! EXPANSION = NONE / NEIGHBOR -> 1 / 2 ; mapped to SPARSITY_NONE(0) / SPARSITY_NEIGHBOR(1)
    expansion = 1
    if( fstr_ctrl_get_param_ex( ctrl, 'EXPANSION ','NONE,NEIGHBOR ', 0, 'P', expansion ) /= 0) return
    contact%sparsity_expansion = expansion - 1
 
    ! Set penalty coefficients to contact structure if specified (for ALagrange method)
    if( np > 0.d0 ) then
      do rcode=1,n
        contact(rcode)%nPenalty = np
      enddo
    endif
    if( tp > 0.d0 ) then
      do rcode=1,n
        contact(rcode)%tPenalty = tp
      enddo
    endif
    if( damp_alpha > 0.0d0 ) then
      do rcode=1,n
        contact(rcode)%damp_alpha = damp_alpha
      enddo
    endif
    if( damp_gact > 0.0d0 ) then
      do rcode=1,n
        contact(rcode)%damp_gact = damp_gact
      enddo
    endif
    
    fstr_ctrl_get_contact = .true.
  end function fstr_ctrl_get_contact
 
  logical function fstr_ctrl_get_embed( ctrl, n, embed, cpname, smoothing )
    use fstr_setup_util
    integer(kind=kint), intent(in)    :: ctrl
    integer(kind=kint), intent(in)    :: n
    type(tcontact), intent(out)       :: embed(n)
    character(len=*), intent(out)     :: cpname
    integer(kind=kint), intent(out)   :: smoothing
 
    integer           :: rcode
    character(len=HECMW_NAME_LEN) :: data_fmt,ss
    character(len=HECMW_NAME_LEN) :: cp_name(n)
    real(kind=kreal)  :: tpenalty(n)
 
    tpenalty = 1.0d6
 
    write(ss,*)  hecmw_name_len
 
    fstr_ctrl_get_embed = .false.
    embed(1)%ctype = 1   ! pure slave-master contact; default value
    embed(1)%algtype = contacttied ! small sliding contact; default value
    if( fstr_ctrl_get_param_ex( ctrl, 'GRPID ', '# ', 1, 'I', embed(1)%group )/=0) return
    smoothing = kcsnone + 1
    if( fstr_ctrl_get_param_ex( ctrl, 'SMOOTHING ', 'NONE,NAGATA ', 0, 'P', smoothing ) /= 0 ) return
    smoothing = smoothing - 1
    do rcode=2,n
      embed(rcode)%ctype = embed(1)%ctype
      embed(rcode)%group = embed(1)%group
      embed(rcode)%algtype = embed(1)%algtype
    end do
 
    write( data_fmt, '(a,a)') 'S', trim(adjustl(ss))
    if(  fstr_ctrl_get_data_array_ex( ctrl, data_fmt, cp_name ) /= 0 ) return
    do rcode=1,n
      embed(rcode)%pair_name = cp_name(rcode)
    enddo
 
    cpname=""
    if( fstr_ctrl_get_param_ex( ctrl, 'CONTACTPARAM ',  '# ',  0, 'S', cpname )/= 0) return
    fstr_ctrl_get_embed = .true.
  end function 
 
  function fstr_ctrl_get_contactparam( ctrl, contactparam )
    implicit none
    integer(kind=kint)    :: ctrl
    type( tcontactparam ) :: contactparam
    integer(kind=kint) :: fstr_ctrl_get_contactparam
 
    integer(kind=kint) :: rcode, nline
    character(len=HECMW_NAME_LEN) :: data_fmt
    character(len=128) :: msg
    real(kind=kreal) :: clearance, clr_same_elem, clr_difflpos, clr_cal_norm
    real(kind=kreal) :: distclr_init, distclr_free, distclr_nocheck, tensile_force
    real(kind=kreal) :: distclr_c2f
    real(kind=kreal) :: penclr_free, penclr_nocheck
 
    fstr_ctrl_get_contactparam = -1
 
    !parameters
    contactparam%name = ''
    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', contactparam%name ) /=0 ) return
 
    !read first line
    data_fmt = 'rrrr '
    rcode = fstr_ctrl_get_data_ex( ctrl, 1, data_fmt, &
      &  clearance, clr_same_elem, clr_difflpos, clr_cal_norm )
    if( rcode /= 0 ) return
    contactparam%CLEARANCE = clearance
    contactparam%CLR_SAME_ELEM = clr_same_elem
    contactparam%CLR_DIFFLPOS  = clr_difflpos
    contactparam%CLR_CAL_NORM  = clr_cal_norm
 
    !read second line
    data_fmt = 'rrrr '
    rcode = fstr_ctrl_get_data_ex( ctrl, 2, data_fmt, &
      &  distclr_init, distclr_free, distclr_nocheck, tensile_force )
    if( rcode /= 0 ) return
    contactparam%DISTCLR_INIT = distclr_init
    contactparam%DISTCLR_FREE = distclr_free
    contactparam%DISTCLR_NOCHECK = distclr_nocheck
    contactparam%TENSILE_FORCE = tensile_force
 
    !read third line (optional): DISTCLR_C2F
    !  If the 3rd line does not exist, keep the default set in init_ContactParam.
    distclr_c2f = contactparam%DISTCLR_C2F
    nline = fstr_ctrl_get_data_line_n( ctrl )
    if( nline >= 3 ) then
      data_fmt = 'r '
      rcode = fstr_ctrl_get_data_ex( ctrl, 3, data_fmt, distclr_c2f )
      if( rcode /= 0 ) return
      contactparam%DISTCLR_C2F = distclr_c2f
    endif
 
    !read fourth line (optional): PENCLR_FREE, PENCLR_NOCHECK
    penclr_free = contactparam%PENCLR_FREE
    penclr_nocheck = contactparam%PENCLR_NOCHECK
    if( nline >= 4 ) then
      data_fmt = 'rr '
      rcode = fstr_ctrl_get_data_ex( ctrl, 4, data_fmt, penclr_free, penclr_nocheck )
      if( rcode /= 0 ) return
      contactparam%PENCLR_FREE = penclr_free
      contactparam%PENCLR_NOCHECK = penclr_nocheck
    endif
 
    !input check
    rcode = 1
    if( clearance<0.d0 .OR. 1.d0<clearance ) THEN
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : CLEARANCE must be 0 < CLEARANCE < 1.'
    else if( clr_same_elem<0.d0 .or. 1.d0<clr_same_elem ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : CLR_SAME_ELEM must be 0 < CLR_SAME_ELEM < 1.'
    else if( clr_difflpos<0.d0 .or. 1.d0<clr_difflpos ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : CLR_DIFFLPOS must be 0 < CLR_DIFFLPOS < 1.'
    else if( clr_cal_norm<0.d0 .or. 1.d0<clr_cal_norm ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : CLR_CAL_NORM must be 0 < CLR_CAL_NORM < 1.'
    else if( distclr_init<0.d0 .or. 1.d0<distclr_init ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : DISTCLR_INIT must be 0 < DISTCLR_INIT < 1.'
    else if( distclr_free<-1.d0 .or. 1.d0<distclr_free ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : DISTCLR_FREE must be -1 < DISTCLR_FREE < 1.'
    else if( distclr_nocheck<0.5d0 ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : DISTCLR_NOCHECK must be >= 0.5.'
    else if( tensile_force>=0.d0 ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : TENSILE_FORCE must be < 0.'
    else if( contactparam%DISTCLR_C2F<=0.d0 .or. 1.d0<=contactparam%DISTCLR_C2F ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : DISTCLR_C2F must be 0 < DISTCLR_C2F < 1.'
    else if( penclr_free<=0.d0 ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : PENCLR_FREE must be > 0.'
    else if( penclr_nocheck<penclr_free ) then
      write(msg,*) 'fstr control file error : !CONTACT_PARAM : PENCLR_NOCHECK must be >= PENCLR_FREE.'
    else
      rcode =0
    end if
    if( rcode /= 0 ) then
      write(*,*) trim(msg)
      write(ilog,*) trim(msg)
      return
    endif
 
    fstr_ctrl_get_contactparam = 0
  end function fstr_ctrl_get_contactparam
 
  function fstr_ctrl_get_contact_if( ctrl, n, contact_if )
    use fstr_setup_util
    integer(kind=kint), intent(in)    :: ctrl
    integer(kind=kint), intent(in)    :: n
    !
    type(tcontactinterference), intent(inout) :: contact_if(n)
    
    integer           :: i
    character(len=30) :: s1 = 'SLAVE,MASTER '
    character(len=HECMW_NAME_LEN) :: data_fmt,ss
    character(len=HECMW_NAME_LEN) :: cp_name(n)
    real(kind=kreal)              :: init_pos(n), end_pos(n)
    integer(kind=kint)            :: fstr_ctrl_get_contact_if
 
    fstr_ctrl_get_contact_if = -1
    write(ss,*)  hecmw_name_len
    if( fstr_ctrl_get_param_ex( ctrl, 'TYPE ', s1, 0, 'P', contact_if(1)%if_type ) /= 0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'END ',  '# ',  0, 'R', contact_if(1)%etime ) /= 0 ) return
    write( data_fmt, '(a,a,a)') 'S', trim(adjustl(ss)),'rr '
    init_pos = 0.d0; end_pos = 0.d0
    if( fstr_ctrl_get_data_array_ex( ctrl, data_fmt, cp_name, init_pos, end_pos ) /= 0 ) return
    do i = 1, n
      contact_if(i)%if_type     = contact_if(1)%if_type
      contact_if(i)%etime       = contact_if(1)%etime
 
      contact_if(i)%cp_name     = cp_name(i)
      contact_if(i)%initial_pos = - init_pos(i)
      contact_if(i)%end_pos     = - end_pos(i)
      if(contact_if(i)%if_type == c_if_slave .and. init_pos(i) /= 0.d0) contact_if(i)%initial_pos = 0.d0
    end do
    fstr_ctrl_get_contact_if = 0
 
  end function fstr_ctrl_get_contact_if
 
  function fstr_ctrl_get_elemopt( ctrl, elemopt361 )
    integer(kind=kint) :: ctrl
    integer(kind=kint) :: elemopt361
    integer(kind=kint) :: fstr_ctrl_get_elemopt
 
    character(72) :: o361list = 'IC,Bbar '
 
    integer(kind=kint) :: o361
 
    fstr_ctrl_get_elemopt = -1
 
    o361 = elemopt361 + 1
 
    !* parameter in header line -----------------------------------------------------------------*!
    if( fstr_ctrl_get_param_ex( ctrl, '361 ', o361list, 0, 'P', o361 ) /= 0) return
 
    elemopt361 = o361 - 1
 
    fstr_ctrl_get_elemopt = 0
 
  end function fstr_ctrl_get_elemopt
 
 
  function fstr_get_autoinc( ctrl, aincparam )
    implicit none
    integer(kind=kint)    :: ctrl
    type( tparamautoinc ) :: aincparam
    integer(kind=kint) :: fstr_get_autoinc
 
    integer(kind=kint) :: rcode
    character(len=HECMW_NAME_LEN) :: data_fmt
    character(len=128) :: msg
    integer(kind=kint) :: bound_s(10), bound_l(10)
    real(kind=kreal) :: rs, rl
 
    fstr_get_autoinc = -1
 
    bound_s(:) = 0
    bound_l(:) = 0
 
    !parameters
    aincparam%name = ''
    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', aincparam%name ) /=0 ) return
 
    !read first line ( decrease criteria )
    data_fmt = 'riiii '
    rcode = fstr_ctrl_get_data_ex( ctrl, 1, data_fmt, rs, &
      &  bound_s(1), bound_s(2), bound_s(3), aincparam%NRtimes_s )
    if( rcode /= 0 ) return
    aincparam%ainc_Rs = rs
    aincparam%NRbound_s(knstmaxit) = bound_s(1)
    aincparam%NRbound_s(knstsumit) = bound_s(2)
    aincparam%NRbound_s(knstciter) = bound_s(3)
 
    !read second line ( increase criteria )
    data_fmt = 'riiii '
    rcode = fstr_ctrl_get_data_ex( ctrl, 2, data_fmt, rl, &
      &  bound_l(1), bound_l(2), bound_l(3), aincparam%NRtimes_l )
    if( rcode /= 0 ) return
    aincparam%ainc_Rl = rl
    aincparam%NRbound_l(knstmaxit) = bound_l(1)
    aincparam%NRbound_l(knstsumit) = bound_l(2)
    aincparam%NRbound_l(knstciter) = bound_l(3)
 
    !read third line ( cutback criteria )
    data_fmt = 'ri '
    rcode = fstr_ctrl_get_data_ex( ctrl, 3, data_fmt, &
      &  aincparam%ainc_Rc, aincparam%CBbound )
    if( rcode /= 0 ) return
 
    !input check
    rcode = 1
    if( rs<0.d0 .or. rs>1.d0 ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : decrease ratio Rs must 0 < Rs < 1.'
    else if( any(bound_s<0) ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : decrease NR bound must >= 0.'
    else if( aincparam%NRtimes_s < 1 ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : # of times to decrease must > 0.'
    else if( rl<1.d0 ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : increase ratio Rl must > 1.'
    else if( any(bound_l<0) ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : increase NR bound must >= 0.'
    else if( aincparam%NRtimes_l < 1 ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : # of times to increase must > 0.'
    elseif( aincparam%ainc_Rc<0.d0 .or. aincparam%ainc_Rc>1.d0 ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : cutback decrease ratio Rc must 0 < Rc < 1.'
    else if( aincparam%CBbound < 1 ) then
      write(msg,*) 'fstr control file error : !AUTOINC_PARAM : maximum # of cutback times must > 0.'
    else
      rcode =0
    end if
    if( rcode /= 0 ) then
      write(*,*) trim(msg)
      write(ilog,*) trim(msg)
      return
    endif
 
    fstr_get_autoinc = 0
  end function fstr_get_autoinc
 
  function fstr_ctrl_get_convergparam( ctrl, cnvparam )
    implicit none
    integer(kind=kint)    :: ctrl
    type( tparamconverg ) :: cnvparam
    integer(kind=kint) :: fstr_ctrl_get_convergparam
 
    integer(kind=kint) :: rcode, i, n, iq, ig1, ig2, inorm
    character(len=HECMW_NAME_LEN) :: data_fmt, data_fmt_tol, ss
    character(len=HECMW_NAME_LEN) :: quantity, norm, action
    character(len=256) :: msg
    real(kind=kreal) :: tol
 
    fstr_ctrl_get_convergparam = -1
 
    !parameters
    cnvparam%name = ''
    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', cnvparam%name ) /=0 ) return
 
    write( ss, * ) hecmw_name_len
    write( data_fmt, '(a,a,a,a,a,a,a)' ) 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), 'S', trim(adjustl(ss)), ' '
    write( data_fmt_tol, '(a,a)' ) trim(data_fmt), 'R '
 
    n = fstr_ctrl_get_data_line_n( ctrl )
    do i = 1, n
      if( fstr_ctrl_get_data_ex( ctrl, i, data_fmt, quantity, norm, action ) /= 0 ) return
      ! the threshold is required for CHECK only, so it is read in a second pass over the line
      tol = 0.d0
      if( action == 'CHECK' ) then
        if( fstr_ctrl_get_data_ex( ctrl, i, data_fmt_tol, quantity, norm, action, tol ) /= 0 ) then
          write(msg,*) 'fstr control file error : !CONVERG_PARAM : threshold of CHECK must be given : ', &
            & trim(quantity), ', ', trim(norm)
          write(*,*) trim(msg)
          write(ilog,*) trim(msg)
          return
        endif
      endif
 
      ! a quantity of the input covers the translational and rotational DOF groups together
      iq = 0
      ig1 = 0
      ig2 = 0
      if( quantity == 'RESIDUAL' ) then
        iq = kcnvresidual
        ig1 = kcnvtranslation
        ig2 = kcnvrotation
      else if( quantity == 'CORRECTION' ) then
        iq = kcnvcorrection
        ig1 = kcnvtranslation
        ig2 = kcnvrotation
      else if( quantity == 'LAGRANGE' ) then
        iq = kcnvcorrection
        ig1 = kcnvlagrange
        ig2 = kcnvlagrange
      endif
      inorm = 0
      if( norm == 'L2' ) then
        inorm = kcnvl2
      else if( norm == 'MAX' ) then
        inorm = kcnvmax
      endif
 
      !input check
      rcode = 1
      if( iq == 0 ) then
        write(msg,*) 'fstr control file error : !CONVERG_PARAM : quantity must be RESIDUAL, CORRECTION or LAGRANGE : ', &
          & trim(quantity)
      else if( inorm == 0 ) then
        write(msg,*) 'fstr control file error : !CONVERG_PARAM : norm must be L2 or MAX : ', trim(norm)
      else if( action /= 'CHECK' .and. action /= 'SKIP' ) then
        write(msg,*) 'fstr control file error : !CONVERG_PARAM : CHECK or SKIP must be given : ', trim(action)
      else if( cnvparam%given(iq,ig1,inorm) ) then
        write(msg,*) 'fstr control file error : !CONVERG_PARAM : ', trim(quantity), ', ', trim(norm), &
          & ' is given more than once.'
      else if( action == 'CHECK' .and. tol <= 0.d0 ) then
        write(msg,*) 'fstr control file error : !CONVERG_PARAM : threshold of CHECK must be > 0 : ', &
          & trim(quantity), ', ', trim(norm)
      else
        rcode = 0
      end if
      if( rcode /= 0 ) then
        write(*,*) trim(msg)
        write(ilog,*) trim(msg)
        return
      endif
 
      cnvparam%given(iq,ig1:ig2,inorm) = .true.
      cnvparam%check(iq,ig1:ig2,inorm) = ( action == 'CHECK' )
      if( action == 'CHECK' ) cnvparam%tol(iq,ig1:ig2,inorm) = tol
    end do
 
    fstr_ctrl_get_convergparam = 0
  end function fstr_ctrl_get_convergparam
 
  function fstr_ctrl_get_timepoints( ctrl, tp )
    integer(kind=kint) :: ctrl
    type(time_points)  :: tp
    integer(kind=kint) :: fstr_ctrl_get_timepoints
 
    integer(kind=kint) :: i, n
    logical            :: generate
    real(kind=kreal)   :: stime, etime, interval
 
    fstr_ctrl_get_timepoints = -1
 
    tp%name = ''
    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', tp%name ) /=0 ) return
    tp%range_type = 1
    if( fstr_ctrl_get_param_ex( ctrl, 'TIME ', 'STEP,TOTAL ', 0, 'P', tp%range_type ) /= 0 ) return
    generate = .false.
    if( fstr_ctrl_get_param_ex( ctrl, 'GENERATE ',  '# ', 0, 'E', generate ) /= 0) return
 
    if( generate ) then
      stime = 0.d0;  etime = 0.d0;  interval = 1.d0
      if( fstr_ctrl_get_data_ex( ctrl, 1, 'rrr ', stime, etime, interval ) /= 0) return
      tp%n_points = int((etime-stime)/interval)+1
      allocate(tp%points(tp%n_points))
      do i=1,tp%n_points
        tp%points(i) = stime + dble(i-1)*interval
      end do
    else
      n = fstr_ctrl_get_data_line_n( ctrl )
      if( n == 0 ) return
      tp%n_points = n
      allocate(tp%points(tp%n_points))
      if( fstr_ctrl_get_data_array_ex( ctrl, 'r ', tp%points ) /= 0 ) return
      do i=1,tp%n_points-1
        if( tp%points(i) < tp%points(i+1) ) cycle
        write(*,*) 'Error in reading !TIME_POINT: time points must be given in ascending order.'
        return
      end do
    end if
 
    fstr_ctrl_get_timepoints = 0
  end function fstr_ctrl_get_timepoints
 
  function fstr_ctrl_get_amplitude( ctrl, nline, name, type_def, type_time, type_val, n, val, table )
    implicit none
    integer(kind=kint),            intent(in)  :: ctrl
    integer(kind=kint),            intent(in)  :: nline
    character(len=HECMW_NAME_LEN), intent(out) :: name
    integer(kind=kint),            intent(out) :: type_def
    integer(kind=kint),            intent(out) :: type_time
    integer(kind=kint),            intent(out) :: type_val
    integer(kind=kint),            intent(out) :: n
    real(kind=kreal),              pointer     :: val(:)
    real(kind=kreal),              pointer     :: table(:)
    integer(kind=kint) :: fstr_ctrl_get_amplitude
 
    integer(kind=kint) :: t_def, t_time, t_val
    integer(kind=kint) :: i, j
    real(kind=kreal) :: r(4), t(4)
 
    fstr_ctrl_get_amplitude = -1
 
    name = ''
    t_def = 1
    t_time = 1
    t_val = 1
 
    if( fstr_ctrl_get_param_ex( ctrl, 'NAME ', '# ', 1, 'S', name )/=0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'DEFINITION ', 'TABULAR ', 0, 'P', t_def )/=0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'TIME ', 'STEP ', 0, 'P', t_time )/=0 ) return
    if( fstr_ctrl_get_param_ex( ctrl, 'VALUE ', 'RELATIVE,ABSOLUTE ', 0, 'P', t_val )/=0 ) return
 
    if( t_def==1 ) then
      type_def = hecmw_amp_typedef_tabular
    else
      write(*,*) 'Error in reading !AMPLITUDE: invalid value for parameter DEFINITION.'
    endif
    if( t_time==1 ) then
      type_time = hecmw_amp_typetime_step
    else
      write(*,*) 'Error in reading !AMPLITUDE: invalid value for parameter TIME.'
    endif
    if( t_val==1 ) then
      type_val = hecmw_amp_typeval_relative
    elseif( t_val==2 ) then
      type_val = hecmw_amp_typeval_absolute
    else
      write(*,*) 'Error in reading !AMPLITUDE: invalid value for parameter VALUE.'
    endif
 
    n = 0
    do i = 1, nline
      r(:)=huge(0.0d0); t(:)=huge(0.0d0)
      if( fstr_ctrl_get_data_ex( ctrl, i, 'RRrrrrrr ', r(1), t(1), r(2), t(2), r(3), t(3), r(4), t(4) ) /= 0) return
      n = n+1
      val(n) = r(1)
      table(n) = t(1)
      do j = 2, 4
        if (r(j) < huge(0.0d0) .and. t(j) < huge(0.0d0)) then
          n = n+1
          val(n) = r(j)
          table(n) = t(j)
        else
          exit
        endif
      enddo
    enddo
    fstr_ctrl_get_amplitude = 0
 
  end function fstr_ctrl_get_amplitude
 
  function fstr_ctrl_get_element_activation( ctrl, amp, eps, grp_id_name, mode, measure, state, thlow, thup )
    implicit none
    integer(kind=kint) :: ctrl
    character(len=HECMW_NAME_LEN) :: amp
    real(kind=kreal) :: eps
    character(len=HECMW_NAME_LEN),target :: grp_id_name(:)
    integer(kind=kint) :: mode     ! 1=FIXED, 2=AMPLITUDE, 3=DAMAGE
    integer(kind=kint) :: measure  ! 1=NONE, 2=STRESS, 3=STRAIN
    integer(kind=kint) :: state    ! 0=ACTIVE, 1=INACTIVE
    real(kind=kreal), target :: thlow(:), thup(:)
    integer(kind=kint) :: fstr_ctrl_get_element_activation
 
    character(len=HECMW_NAME_LEN),pointer :: element_id_p
    real(kind=kreal),pointer :: thlow_p(:), thup_p(:)
    integer(kind=kint) :: rcode, n
    character(len=HECMW_NAME_LEN) :: data_fmt, s1
 
    fstr_ctrl_get_element_activation = -1
    
    ! MODE (required)
    if( fstr_ctrl_get_param_ex( ctrl, 'MODE ', 'FIXED,AMPLITUDE,DAMAGE ', 1, 'P', mode ) /= 0 ) return
    
    ! Mode-specific parameters
    if( mode == 1 ) then
      ! FIXED: STATE required
      if( fstr_ctrl_get_param_ex( ctrl, 'STATE ', 'ON,OFF ', 1, 'P', state ) /= 0 ) return
      state = state - 1
      measure = 1
      amp = ''
    elseif( mode == 2 ) then
      ! AMPLITUDE: AMP required
      if( fstr_ctrl_get_param_ex( ctrl, 'AMP ', '# ', 1, 'S', amp ) /= 0 ) return
      state = 0
      measure = 1
    elseif( mode == 3 ) then
      ! DAMAGE: MEASURE required
      if( fstr_ctrl_get_param_ex( ctrl, 'MEASURE ', 'STRESS,STRAIN ', 1, 'P', measure ) /= 0 ) return
      measure = measure + 1
      state = 0
      amp = ''
    endif
    
    ! EPSILON (optional)
    eps = 1.0d-6
    if( fstr_ctrl_get_param_ex( ctrl, 'EPSILON ', '# ', 0, 'R', eps ) /= 0 ) return
 
    write(s1,*) hecmw_name_len
    n = fstr_ctrl_get_data_line_n(ctrl)
    element_id_p => grp_id_name(1)
    thlow_p => thlow
    thup_p => thup
 
    if( mode == 3 ) then
      write( data_fmt, '(a,a,a)') 'S', trim(adjustl(s1)), 'RR'
      rcode = fstr_ctrl_get_data_array_ex( ctrl, data_fmt, element_id_p, thlow_p, thup_p )
    else
      write( data_fmt, '(a,a)') 'S', trim(adjustl(s1))
      rcode = fstr_ctrl_get_data_array_ex( ctrl, data_fmt, element_id_p )
    endif
 
    fstr_ctrl_get_element_activation = 0
 
  end function fstr_ctrl_get_element_activation
 
 
end module fstr_ctrl_common