fstr_setup_finish.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 m_fstr_setup_finish
  use m_fstr
  use fstr_setup_util, only &
    : fstr_param_pack &
    , fstr_setup_util_err_stop &
    , append_node_grp_from_surf_grp &
    , reallocate_integer &
    , reallocate_real
  use m_fstr_solid_setup, only &
    : fstr_solid_alloc &
    , fstr_element_init
  use mcontactparam
  use m_static_get_prop
  implicit none
 
  type fstr_setup_counts
    integer(kind=kint) :: contact_pairs = 0
    integer(kind=kint) :: embed_pairs = 0
    integer(kind=kint) :: contact_interferences = 0
    integer(kind=kint) :: contact_params = 0
    integer(kind=kint) :: weld_lines = 0
    integer(kind=kint) :: initial_conditions = 0
    integer(kind=kint) :: steps = 0
    integer(kind=kint) :: local_coords = 0
    integer(kind=kint) :: autoinc_params = 0
    integer(kind=kint) :: converg_params = 0
    integer(kind=kint) :: time_points = 0
    integer(kind=kint) :: materials = 0
  end type fstr_setup_counts
 
contains
 
  subroutine fstr_setup_prepare( P, counts )
    use mmaterial
    implicit none
    type(fstr_param_pack) :: P
    type(fstr_setup_counts), intent(in) :: counts
 
    type(hecmwst_local_mesh), pointer :: hecMESH
    type(fstr_param), pointer :: fstrPARAM
    type(fstr_solid), pointer :: fstrSOLID
    type(fstr_heat), pointer :: fstrHEAT
    type(fstr_dynamic), pointer :: fstrDYNAMIC
    integer(kind=kint) :: n, i, cid, n_totlyr
    real(kind=kreal) :: ee, pp, rho, alpha, thick, alpha_over_mu
    real(kind=kreal) :: beam_radius,                          &
      beam_angle1, beam_angle2, beam_angle3,&
      beam_angle4, beam_angle5, beam_angle6
    type(tshellmat),pointer :: shmat(:)
    character(len=HECMW_FILENAME_LEN) :: logfilename
 
    hecmesh => p%MESH
    fstrparam => p%PARAM
    fstrsolid => p%SOLID
    fstrheat => p%HEAT
    fstrdynamic => p%DYN
    write( logfilename, '(i5,''.log'')' ) myrank
 
    fstrsolid%n_contacts = counts%contact_pairs
    if( counts%contact_pairs>0 ) then
      allocate( fstrsolid%contacts( counts%contact_pairs ) )
      ! convert SURF_SURF contact to NODE_SURF contact
      call fstr_convert_contact_type( p%MESH )
    endif
    call fstr_convert_rotation_dof( p )
    fstrsolid%n_embeds = counts%embed_pairs
    if( counts%embed_pairs>0 )  allocate( fstrsolid%embeds( counts%embed_pairs ) )
    if( counts%weld_lines>0 ) allocate( fstrheat%weldline( counts%weld_lines ) )
    if( counts%initial_conditions>0 ) allocate( g_initialcnd( counts%initial_conditions ) )
    if( counts%steps>0 ) then
      allocate( fstrsolid%step_ctrl( counts%steps ) )
      do i=1, counts%steps
        call init_stepinfo( fstrsolid%step_ctrl(i) )
        if( p%PARAM%solution_type==kstdynamic ) then
          fstrsolid%step_ctrl(i)%num_substep = fstrdynamic%n_step
          fstrsolid%step_ctrl(i)%initdt      = fstrdynamic%t_delta
          fstrsolid%step_ctrl(i)%elapsetime  = dble(fstrdynamic%n_step) * fstrdynamic%t_delta
          fstrsolid%step_ctrl(i)%mindt       = fstrdynamic%t_delta
          fstrsolid%step_ctrl(i)%maxdt       = fstrdynamic%t_delta
        endif
      end do
    endif
    if( counts%local_coords>0 ) allocate( g_localcoordsys(counts%local_coords) )
    allocate( fstrparam%ainc(0:counts%autoinc_params) )
    do i=0,counts%autoinc_params
      call init_aincparam( fstrparam%ainc(i) )
    end do
    allocate( fstrparam%cnvparam(0:counts%converg_params) )
    do i=0,counts%converg_params
      call init_convergparam( fstrparam%cnvparam(i) )
    end do
    if( counts%time_points>0 ) allocate( fstrparam%timepoints(counts%time_points) )
    allocate( fstrparam%contactparam(0:counts%contact_params) )
    do i=0,counts%contact_params
      call init_contactparam( fstrparam%contactparam(i) )
    end do
    if( counts%contact_interferences>0 )then
      allocate( fstrparam%contact_if( counts%contact_interferences ) )
      do i=1,counts%contact_interferences
        call init_contact_if( fstrparam%contact_if(i) )
      end do
    end if
 
    p%SOLID%is_33shell = 0
    p%SOLID%is_33beam  = 0
 
    do i=1,hecmesh%n_elem_type
      n =  hecmesh%elem_type_item(i)
      if (n == 781 .or. n == 761)then
        p%SOLID%is_33shell = 1
      elseif (n == 641)then
        p%SOLID%is_33beam  = 1
      endif
    enddo
 
    n = counts%materials
    if( hecmesh%material%n_mat>n ) n= hecmesh%material%n_mat
    if( n==0 ) stop "material property not defined!"
    allocate( fstrsolid%materials( n ) )
    do i = 1, n
      call initmaterial(fstrsolid%materials(i))
    enddo
 
    ! for section control
    allocate( fstrsolid%sections(hecmesh%section%n_sect) )
    do i=1,hecmesh%section%n_sect
      ! set default 361 element formulation
      if( p%PARAM%solution_type==kststatic .or. p%PARAM%solution_type==kstdynamic ) then
        if( p%PARAM%nlgeom ) then
          fstrsolid%sections(i)%elemopt361 = kel361fbar
        else
          fstrsolid%sections(i)%elemopt361 = kel361ic
        end if
      else if( p%PARAM%solution_type==ksteigen ) then
        fstrsolid%sections(i)%elemopt361 = kel361ic
      else if( p%PARAM%solution_type==kststaticeigen ) then
        fstrsolid%sections(i)%elemopt361 = kel361fbar
      else
        fstrsolid%sections(i)%elemopt361 = kel361fi
      end if
      fstrsolid%sections(i)%elemopt341 = kel341fi
      fstrsolid%sections(i)%elemopt611 = kel611euler
      ! INTERFACE sections are skipped in the loop below and keep this thickness
      fstrsolid%sections(i)%thickness = 1.0d0
    enddo
 
    if( hecmesh%section%n_sect >0 ) then
      do i=1,hecmesh%section%n_sect
        if( hecmesh%section%sect_type(i) == 4 ) cycle
        cid = hecmesh%section%sect_mat_ID_item(i)
        if( cid>n ) stop "Error in material property definition!"
        if( fstrparam%nlgeom .or. fstrparam%solution_type==kststaticeigen ) &
          fstrsolid%materials(cid)%nlgeom_flag = 1
        nullify(shmat)
        call fstr_get_prop(hecmesh,shmat,i,ee,pp,rho,alpha,thick,&
          n_totlyr,alpha_over_mu, &
          beam_radius,beam_angle1,beam_angle2,beam_angle3, &
          beam_angle4,beam_angle5,beam_angle6)
        fstrsolid%sections(i)%thickness = thick
        fstrsolid%materials(cid)%name = hecmesh%material%mat_name(cid)
        fstrsolid%materials(cid)%variables(m_youngs)=ee
        fstrsolid%materials(cid)%variables(m_poisson)=pp
        fstrsolid%materials(cid)%variables(m_density)=rho
        fstrsolid%materials(cid)%variables(m_expansion)=alpha
        fstrsolid%materials(cid)%variables(m_alpha_over_mu)= alpha_over_mu
        fstrsolid%materials(cid)%variables(m_beam_radius)=beam_radius
        fstrsolid%materials(cid)%variables(m_beam_angle1)=beam_angle1
        fstrsolid%materials(cid)%variables(m_beam_angle2)=beam_angle2
        fstrsolid%materials(cid)%variables(m_beam_angle3)=beam_angle3
        fstrsolid%materials(cid)%variables(m_beam_angle4)=beam_angle4
        fstrsolid%materials(cid)%variables(m_beam_angle5)=beam_angle5
        fstrsolid%materials(cid)%variables(m_beam_angle6)=beam_angle6
        fstrsolid%materials(cid)%mtype = elastic
        if( hecmesh%section%sect_type(i) == 2 ) then
          fstrsolid%materials(cid)%totallyr = n_totlyr
          fstrsolid%materials(cid)%shell_var => shmat
        endif
      enddo
    endif
 
    allocate( fstrsolid%output_ctrl( 4 ) )
    call fstr_init_outctrl(fstrsolid%output_ctrl(1))
    fstrsolid%output_ctrl( 1 )%filename = trim(logfilename)
    fstrsolid%output_ctrl( 1 )%filenum = ilog
    call fstr_init_outctrl(fstrsolid%output_ctrl(2))
    call fstr_init_outctrl(fstrsolid%output_ctrl(3))
    call fstr_init_outctrl(fstrsolid%output_ctrl(4))
  end subroutine fstr_setup_prepare
 
  subroutine fstr_setup_finish( P, counts )
    implicit none
    type(fstr_param_pack) :: P
    type(fstr_setup_counts), intent(in) :: counts
 
    type(hecmwst_local_mesh), pointer :: hecMESH
    type(fstr_param), pointer :: fstrPARAM
    type(fstr_solid), pointer :: fstrSOLID
    type(fstr_dynamic), pointer :: fstrDYNAMIC
    integer(kind=kint) :: n, i, cid, ierror, n_totlyr
 
    hecmesh => p%MESH
    fstrparam => p%PARAM
    fstrsolid => p%SOLID
    fstrdynamic => p%DYN
 
    ! ----- material type judgement. in case of infinitive analysis, nlgeom_flag=0
    if( .not. p%PARAM%nlgeom ) then
      do i=1, counts%materials
        fstrsolid%materials(i)%nlgeom_flag = 0
      enddo
    endif
 
    if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
      allocate ( fstrsolid%temperature( hecmesh%n_node )      ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, TEMPERATURE>'
        write(idbg,*) '  rank = ', myrank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      fstrsolid%temperature = fstrparam%ref_temp
      allocate ( fstrsolid%last_temp( hecmesh%n_node )      ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, LAST_TEMP>'
        write(idbg,*) '  rank = ', myrank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      fstrsolid%last_temp = 0.d0
      allocate ( fstrsolid%temp_bak( hecmesh%n_node )      ,stat=ierror )
      if( ierror /= 0 ) then
        write(idbg,*) 'stop due to allocation error <FSTR_SOLID, TEMP_BAK>'
        write(idbg,*) '  rank = ', myrank,'  ierror = ',ierror
        call flush(idbg)
        call hecmw_abort( hecmw_comm_get_comm())
      end if
      fstrsolid%temp_bak = 0.d0
    endif
 
    if( associated(fstrsolid%step_ctrl) )  then
      fstrsolid%nstep_tot = size(fstrsolid%step_ctrl)
      call setup_stepinfo_starttime( fstrsolid%step_ctrl )
    else
      if( p%PARAM%solution_type==kststatic .and. p%PARAM%nlgeom ) then
        write( *,* ) " ERROR: STEP not defined!"
        write( idbg,* ) "ERROR: STEP not defined!"
        call flush(idbg)
        call fstr_abort( hecmw_exit_model )
      endif
 
      if( myrank==0 ) write(*,*)"Step control not defined! Using default step=1"
      fstrsolid%nstep_tot = 1
      allocate( fstrsolid%step_ctrl(1) )
      call init_stepinfo( fstrsolid%step_ctrl(1) )
      if( p%PARAM%solution_type==kstdynamic ) then
        fstrsolid%step_ctrl(1)%num_substep = fstrdynamic%n_step
        fstrsolid%step_ctrl(1)%initdt      = fstrdynamic%t_delta
        fstrsolid%step_ctrl(1)%elapsetime  = dble(fstrdynamic%n_step) * fstrdynamic%t_delta
        fstrsolid%step_ctrl(1)%mindt       = fstrdynamic%t_delta
        fstrsolid%step_ctrl(1)%maxdt       = fstrdynamic%t_delta
      endif
      n =  fstrsolid%BOUNDARY_ngrp_tot
      if( n>0 ) allocate( fstrsolid%step_ctrl(1)%Boundary(n) )
      do i = 1, n
        fstrsolid%step_ctrl(1)%Boundary(i) = fstrsolid%BOUNDARY_ngrp_GRPID(i)
      enddo
      n = fstrsolid%CLOAD_ngrp_tot + fstrsolid%DLOAD_ngrp_tot + fstrsolid%TEMP_ngrp_tot + fstrsolid%SPRING_ngrp_tot
      if( n>0 ) allocate( fstrsolid%step_ctrl(1)%Load(n) )
      n = 0
      do i = 1, fstrsolid%CLOAD_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%CLOAD_ngrp_GRPID(i)
      enddo
      do i = 1, fstrsolid%DLOAD_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%DLOAD_ngrp_GRPID(i)
      enddo
      do i = 1, fstrsolid%TEMP_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%TEMP_ngrp_GRPID(i)
      enddo
      do i = 1, fstrsolid%SPRING_ngrp_tot
        n = n + 1
        fstrsolid%step_ctrl(1)%Load(n) = fstrsolid%SPRING_ngrp_GRPID(i)
      enddo
      n = fstrsolid%elemact%ELEMACT_egrp_tot
      if( n>0 ) allocate( fstrsolid%step_ctrl(1)%ElemActivation(n) )
      do i = 1, n
        fstrsolid%step_ctrl(1)%ElemActivation(i) = fstrsolid%elemact%ELEMACT_egrp_GRPID(i)
      enddo
    endif
 
    call fstr_element_init( hecmesh, fstrsolid, p%PARAM%solution_type )
    call fstr_check_shell_j2_inputs(hecmesh, fstrsolid, fstrparam)
    call fstr_check_up_inputs(hecmesh, fstrsolid)
    call fstr_check_shell_gauss_output(hecmesh, fstrsolid)
    if( p%PARAM%solution_type==kststatic .or. p%PARAM%solution_type==kstdynamic .or.   &
      p%PARAM%solution_type==ksteigen  .or. p%PARAM%solution_type==kststaticeigen )  &
      call fstr_solid_alloc( hecmesh, fstrsolid )
 
    if( p%PARAM%solution_type == kstheat) then
      p%PARAM%fg_irres = fstrsolid%output_ctrl(3)%frequency
      p%PARAM%fg_iwres = fstrsolid%output_ctrl(4)%frequency
      p%HEAT%elemact = p%SOLID%elemact
    endif
 
    n_totlyr = 1
    do i=1,hecmesh%section%n_sect
      cid = hecmesh%section%sect_mat_ID_item(i)
      n  = fstrsolid%materials(cid)%totallyr
      if (n > n_totlyr)then
        n_totlyr = n
      endif
    enddo
    p%SOLID%max_lyr = n_totlyr
 
    call fstr_setup_post( p )
  end subroutine fstr_setup_finish
 
  subroutine fstr_check_up_inputs(hecMESH, fstrSOLID)
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(in) :: fstrSOLID
    integer(kind=kint) :: i, isect, mtype
 
    do i = 1, hecmesh%n_elem
      if( hecmesh%elem_type(i) /= 361 ) cycle
      isect = hecmesh%section_ID(i)
      if( fstrsolid%sections(isect)%elemopt361 /= kel361up ) cycle
      if( .not.associated(fstrsolid%elements(i)%gausses) ) cycle
      mtype = fstrsolid%elements(i)%gausses(1)%pMaterial%mtype
      if( .not.iselastoplastic(mtype) .and. .not.iscreep(mtype) ) cycle
      write(*,*) '###ERROR### : Elastoplastic and creep materials are not supported by FORM361=UP'
      write(ilog,*) '###ERROR### : Elastoplastic and creep materials are not supported by FORM361=UP'
      call hecmw_abort(hecmw_comm_get_comm())
    enddo
  end subroutine fstr_check_up_inputs
 
  subroutine fstr_check_shell_gauss_output(hecMESH, fstrSOLID)
    use mmechgauss, only: fstr_shell_output_point_count
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(in) :: fstrSOLID
    integer(kind=kint) :: i
 
    if( hecmesh%n_dof /= 6 ) return
    if( .not.any(fstrsolid%output_ctrl(3)%outinfo%on(9:10)) ) return
    do i=1,hecmesh%n_elem
      if( fstr_shell_output_point_count(fstrsolid%elements(i))>0 ) cycle
      write(*,*) '###ERROR### : ISTRAIN/ISTRESS require shell thickness-point output support; element type = ', &
        hecmesh%elem_type(i)
      call hecmw_abort(hecmw_comm_get_comm())
    enddo
  end subroutine fstr_check_shell_gauss_output
 
  subroutine fstr_check_shell_j2_inputs(hecMESH, fstrSOLID, fstrPARAM)
    use m_elastoplastic, only: checkplanestressj2material, planestressj2errormessage
    use elementinfo, only: numofshellthicknessquadpoints
    use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(in) :: fstrSOLID
    type(fstr_param), intent(in) :: fstrPARAM
    type(tmaterial), pointer :: material
    integer(kind=kint) :: i, nn, ierr, isect, index, ilayer
    real(kind=kreal) :: elastic(2)
    logical :: missing
    character(len=160) :: message
 
    do i = 1, hecmesh%n_elem
      if( hecmesh%elem_type(i) /= 741 ) cycle
      if( .not.associated(fstrsolid%elements(i)%gausses) ) cycle
      material => fstrsolid%elements(i)%gausses(1)%pMaterial
      if( .not.iselastoplastic(material%mtype) ) cycle
      nn = hecmesh%elem_node_index(i)-hecmesh%elem_node_index(i-1)
      message = ''
      call checkplanestressj2material(material, ierr)
      if( ierr /= 0 ) then
        message = planestressj2errormessage(ierr)
      else if( .not.fstr_uses_plane_stress_j2_shell(741, nn, material) ) then
        message = 'Elastoplastic MITC4 shells require INFINITESIMAL material kinematics'
      else if( fstrparam%solution_type /= kststatic .or. .not.fstrparam%nlgeom ) then
        message = 'Elastoplastic MITC4 shells require !SOLUTION, TYPE=STATIC, NONLINEAR'
      else if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
        message = 'Temperature loading is not supported in models with elastoplastic MITC4 shells'
      endif
      if( len_trim(message) == 0 .and. associated(material%shell_var) ) then
        ! The J2 return mapping uses the common elastic table, not layer-specific constants.
        call fetch_tabledata(mc_isoelastic, material%dict, elastic, missing)
        do ilayer = 1, material%totallyr
          if( .not.ieee_is_finite(material%shell_var(ilayer)%ee) .or. &
              .not.ieee_is_finite(material%shell_var(ilayer)%pp) ) then
            message = 'Invalid elastic properties in J2 shell layer'
            exit
          endif
          if( abs(material%shell_var(ilayer)%ee-elastic(1)) > 1.0d-12*elastic(1) .or. &
              abs(material%shell_var(ilayer)%pp-elastic(2)) > 1.0d-12 ) then
            message = 'J2 shell layers must share the !ELASTIC Young modulus and Poisson ratio'
            exit
          endif
        enddo
      endif
      if( len_trim(message) == 0 ) then
        isect = hecmesh%section_ID(i)
        index = hecmesh%section%sect_I_index(isect-1)
        if( hecmesh%section%sect_I_index(isect)-index /= 1 ) then
          message = 'Elastoplastic MITC4 shells require a SHELL section with 2 thickness integration points'
        else if( hecmesh%section%sect_I_item(index+1) /= numofshellthicknessquadpoints(741) ) then
          message = 'Elastoplastic MITC4 shells require 2 thickness integration points in !SECTION'
        endif
      endif
      if( len_trim(message) == 0 ) cycle
      write(*,*) '###ERROR### : ', trim(message)
      write(ilog,*) '###ERROR### : ', trim(message)
      call hecmw_abort(hecmw_comm_get_comm())
    enddo
  end subroutine fstr_check_shell_j2_inputs
 
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_setup_post_phys_alloc(phys, NDOF, n_node, n_elem)
    implicit none
    type(fstr_solid_physic_val), pointer :: phys
    integer(kind=kint) :: NDOF, n_node, n_elem, mdof, istat
    mdof = (ndof*ndof+ndof)/2;
    allocate ( phys%STRAIN  (mdof*n_node), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%STRAIN"
    allocate ( phys%STRESS  (mdof*n_node), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%STRESS"
    allocate ( phys%MISES   (     n_node), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%MISES"
    allocate ( phys%ESTRAIN (mdof*n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%ESTRAIN"
    allocate ( phys%ESTRESS (mdof*n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%ESTRESS"
    allocate ( phys%EMISES  (     n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%EMISES"
    allocate ( phys%EPLSTRAIN (   n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%EPLSTRAIN"
    allocate ( phys%ENQM    (12*n_elem), stat=istat)
    if( istat /= 0 ) stop "Allocation error: phys%ENQM"
  end subroutine fstr_setup_post_phys_alloc
 
  subroutine fstr_setup_post( P )
    implicit none
    integer(kind=kint) :: i
    type(fstr_param_pack) :: P
    type(fstr_solid_physic_val), pointer :: phys => null()
 
    if( p%PARAM%solution_type == kststatic &
        .or. p%PARAM%solution_type == ksteigen  &
        .or. p%PARAM%solution_type == kstdynamic &
        .or. p%PARAM%solution_type == kststaticeigen ) then
      ! Memory Allocation for Result Vectors ------------
      if( p%MESH%n_dof == 6 .or. p%SOLID%is_33shell == 1 ) then
        allocate ( p%SOLID%SHELL )
        call fstr_setup_post_phys_alloc(p%SOLID%SHELL,3, p%MESH%n_node,p%MESH%n_elem)
        allocate ( p%SOLID%SHELL%LAYER(p%SOLID%max_lyr) )
        do i=1,p%SOLID%max_lyr
          allocate ( p%SOLID%SHELL%LAYER(i)%PLUS )
          allocate ( p%SOLID%SHELL%LAYER(i)%MINUS )
          call fstr_setup_post_phys_alloc(p%SOLID%SHELL%LAYER(i)%PLUS , 3, p%MESH%n_node, p%MESH%n_elem)
          call fstr_setup_post_phys_alloc(p%SOLID%SHELL%LAYER(i)%MINUS, 3, p%MESH%n_node, p%MESH%n_elem)
        enddo
        phys => p%SOLID%SHELL
      else
        allocate ( p%SOLID%SOLID )
        phys => p%SOLID%SOLID
        call fstr_setup_post_phys_alloc(phys, p%MESH%n_dof, p%MESH%n_node,  p%MESH%n_elem)
      end if
      p%SOLID%STRAIN => phys%STRAIN
      p%SOLID%STRESS => phys%STRESS
      p%SOLID%MISES  => phys%MISES
      p%SOLID%ESTRAIN => phys%ESTRAIN
      p%SOLID%ESTRESS => phys%ESTRESS
      p%SOLID%EMISES  => phys%EMISES
      p%SOLID%EPLSTRAIN  => phys%EPLSTRAIN
      p%SOLID%ENQM    => phys%ENQM
      allocate( p%SOLID%REACTION( p%MESH%n_dof*p%MESH%n_node ), stat=i )
      if( i /= 0 ) stop "Allocation error: REACTION"
      p%SOLID%REACTION = 0.d0
    end if
 
    call hecmw_barrier( p%MESH ) ! JP-7
 
    if( p%HEAT%STEPtot == 0 ) then ! No !HEAT Input
      if( p%PARAM%analysis_n == 0 ) then  ! No !STATIC Input
        call reallocate_real( p%PARAM%dtime, 1)
        call reallocate_real( p%PARAM%etime, 1)
        call reallocate_real( p%PARAM%dtmin, 1)
        call reallocate_real( p%PARAM%delmax,1)
        call reallocate_integer( p%PARAM%itmax, 1)
        call reallocate_real( p%PARAM%eps,   1)
        p%PARAM%analysis_n = 1
        p%PARAM%dtime = 0
        p%PARAM%etime = 0
        p%PARAM%dtmin = 0
        p%PARAM%delmax = 0
        p%PARAM%itmax = 20
        p%PARAM%eps = 1.0e-6
      end if
      p%HEAT%STEPtot = 1
      call reallocate_real( p%HEAT%STEP_DLTIME, 1)
      call reallocate_real( p%HEAT%STEP_EETIME, 1)
      call reallocate_real( p%HEAT%STEP_DELMIN, 1)
      call reallocate_real( p%HEAT%STEP_DELMAX, 1)
      p%HEAT%STEP_DLTIME = 0
      p%HEAT%STEP_EETIME = 0
      p%HEAT%STEP_DELMIN = 0
      p%HEAT%STEP_DELMAX = 0
    end if
  end subroutine fstr_setup_post
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_convert_contact_type( hecMESH )
    implicit none
    type(hecmwst_local_mesh), pointer :: hecMESH
    integer(kind=kint) :: n, i, sgrp_id, ngrp_id
    ! convert SURF_SURF to NODE_SURF
    n = hecmesh%contact_pair%n_pair
    do i = 1,n
      if( hecmesh%contact_pair%type(i) /= hecmw_contact_type_surf_surf ) cycle
      sgrp_id = hecmesh%contact_pair%slave_grp_id(i)
      call append_node_grp_from_surf_grp( hecmesh, sgrp_id, ngrp_id )
      ! change type of contact and slave group ID
      hecmesh%contact_pair%type(i) = hecmw_contact_type_node_surf
      hecmesh%contact_pair%slave_grp_id(i) = ngrp_id
    enddo
  end subroutine fstr_convert_contact_type
 
  !-----------------------------------------------------------------------------!
  !-----------------------------------------------------------------------------!
 
  subroutine fstr_convert_rotation_dof( P )
    use hecmw_setup_util, only : append_new_groups, get_grp_id
    implicit none
    type(fstr_param_pack) :: P
    type(hecmwst_local_mesh), pointer :: hecMESH
    integer(kind=kint), allocatable :: rot_node(:), n_rot(:), n_rot_all(:), n_explicit(:), rot_grp(:), src(:)
    integer(kind=kint), allocatable :: new_index(:), new_list(:)
    character(len=HECMW_NAME_LEN), allocatable :: new_name(:)
    logical, allocatable :: is_rot(:), in_all(:), is_explicit(:)
    integer(kind=kint) :: icel, is, ie, nn, i, j, k, m, n, ig, n_new, first_id, ids, ide
 
    hecmesh => p%MESH
    if( hecmesh%n_dof /= 3 ) return
 
    ! explicit dummy nodes are in ALL, and the generated ones are not
    allocate( in_all(hecmesh%n_node) )
    in_all = .false.
    ig = get_grp_id( hecmesh, 'node_grp', 'ALL' )
    if( ig > 0 ) then
      is = hecmesh%node_group%grp_index(ig-1)
      in_all(hecmesh%node_group%grp_item(is+1:hecmesh%node_group%grp_index(ig))) = .true.
    endif
 
    allocate( rot_node(hecmesh%n_node), is_explicit(hecmesh%n_node) )
    rot_node = 0
    is_explicit = .false.
    do icel = 1, hecmesh%n_elem
      if( .not. hecmw_is_etype_33struct(hecmesh%elem_type(icel)) ) cycle
      nn = hecmw_get_max_node(hecmesh%elem_type(icel))/2
      is = hecmesh%elem_node_index(icel-1)
      do j = 1, nn
        rot_node(hecmesh%elem_node_item(is+j)) = hecmesh%elem_node_item(is+nn+j)
        if( in_all(hecmesh%elem_node_item(is+nn+j)) ) is_explicit(hecmesh%elem_node_item(is+j)) = .true.
      enddo
    enddo
    ! a rank may hold an external translational node without its 641/761/781 element,
    ! so the conditions are converted on every rank by the counts over all the ranks
    n = count( rot_node > 0 )
    call hecmw_allreduce_i1( hecmesh, n, hecmw_max )
    if( n == 0 ) return
 
    allocate( n_rot(hecmesh%node_group%n_grp), n_rot_all(hecmesh%node_group%n_grp) )
    allocate( n_explicit(hecmesh%node_group%n_grp) )
    do ig = 1, hecmesh%node_group%n_grp
      is = hecmesh%node_group%grp_index(ig-1)
      n_rot(ig) = count( rot_node(hecmesh%node_group%grp_item(is+1:hecmesh%node_group%grp_index(ig))) > 0 )
      n_explicit(ig) = count( is_explicit(hecmesh%node_group%grp_item(is+1:hecmesh%node_group%grp_index(ig))) )
    enddo
    n_rot_all = n_rot
    call hecmw_allreduce_i( hecmesh, n_rot_all, hecmesh%node_group%n_grp, hecmw_max )
 
    do i = 1, p%SOLID%SPRING_ngrp_tot
      if( p%SOLID%SPRING_ngrp_DOF(i) > 3 .and. n_rot_all(p%SOLID%SPRING_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!SPRING' )
    enddo
    do i = 1, p%SOLID%VELOCITY_ngrp_tot
      if( p%SOLID%VELOCITY_ngrp_rotID(i) > 0 ) cycle
      if( mod(p%SOLID%VELOCITY_ngrp_type(i), 10) > 3 .and. n_rot_all(p%SOLID%VELOCITY_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!VELOCITY' )
    enddo
    do i = 1, p%SOLID%ACCELERATION_ngrp_tot
      if( mod(p%SOLID%ACCELERATION_ngrp_type(i), 10) > 3 .and. n_rot_all(p%SOLID%ACCELERATION_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!ACCELERATION' )
    enddo
    do i = 1, p%FREQ%FLOAD_ngrp_tot
      if( p%FREQ%FLOAD_ngrp_TYPE(i) /= kfloadtype_node ) cycle
      if( p%FREQ%FLOAD_ngrp_DOF(i) > 3 .and. n_rot_all(p%FREQ%FLOAD_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!FLOAD' )
    enddo
    !the rotational DOFs of !FACCELERATION are also not allowed (the end DOF is mod(type,10))
    do i = 1, p%FREQ%FACC_ngrp_tot
      if( mod(p%FREQ%FACC_ngrp_type(i), 10) > 3 .and. n_rot_all(p%FREQ%FACC_ngrp_ID(i)) > 0 ) &
        call rotation_dof_err_stop( '!FACCELERATION' )
    enddo
    do i = 1, hecmesh%mpc%n_mpc
      do j = hecmesh%mpc%mpc_index(i-1)+1, hecmesh%mpc%mpc_index(i)
        if( hecmesh%mpc%mpc_dof(j) > 3 .and. rot_node(hecmesh%mpc%mpc_item(j)) > 0 ) &
          call rotation_dof_err_stop( '!EQUATION' )
      enddo
    enddo
 
    ! node groups whose rotation is moved: rot_grp(ig) = -1, then the id of the dummy node group
    allocate( rot_grp(hecmesh%node_group%n_grp) )
    rot_grp = 0
    do i = 1, p%SOLID%BOUNDARY_ngrp_tot
      if( p%SOLID%BOUNDARY_ngrp_rotID(i) > 0 ) cycle
      ig = p%SOLID%BOUNDARY_ngrp_ID(i)
      if( mod(p%SOLID%BOUNDARY_ngrp_type(i), 10) <= 3 .or. n_rot_all(ig) == 0 ) cycle
      if( n_explicit(ig) > 0 ) call explicit_dummy_node_err_stop( '!BOUNDARY' )
      rot_grp(ig) = -1
    enddo
    do i = 1, p%SOLID%CLOAD_ngrp_tot
      if( p%SOLID%CLOAD_ngrp_rotID(i) > 0 .or. p%SOLID%CLOAD_ngrp_DOF(i) <= 3 ) cycle
      ig = p%SOLID%CLOAD_ngrp_ID(i)
      is = hecmesh%node_group%grp_index(ig-1)
      ie = hecmesh%node_group%grp_index(ig)
      if( any( rot_node(hecmesh%node_group%grp_item(is+1:ie)) == 0 .and. &
        hecmesh%node_group%grp_item(is+1:ie) <= hecmesh%nn_internal ) ) call fstr_setup_util_err_stop( &
        '### Error: !CLOAD : dof 4-6 is given to nodes without rotational dof' )
      if( n_explicit(ig) > 0 ) call explicit_dummy_node_err_stop( '!CLOAD' )
      rot_grp(ig) = -1
    enddo
 
    n_new = count( rot_grp == -1 )
    if( n_new == 0 ) return
    allocate( new_name(n_new), new_index(0:n_new), new_list(sum(n_rot, mask=rot_grp == -1)) )
    new_index(0) = 0
    k = 0
    m = 0
    do ig = 1, hecmesh%node_group%n_grp
      if( rot_grp(ig) /= -1 ) cycle
      k = k + 1
      write( new_name(k), '(a,i0)' ) 'FSTR_ROT_', ig
      do j = hecmesh%node_group%grp_index(ig-1)+1, hecmesh%node_group%grp_index(ig)
        if( rot_node(hecmesh%node_group%grp_item(j)) == 0 ) cycle
        m = m + 1
        new_list(m) = rot_node(hecmesh%node_group%grp_item(j))
      enddo
      new_index(k) = m
      rot_grp(ig) = k
    enddo
    call append_new_groups( hecmesh, 'node_grp', n_new, new_name, new_index, new_list, first_id )
    where( rot_grp > 0 ) rot_grp = rot_grp + first_id - 1
 
    ! split each !BOUNDARY entry into dof 1-3 of its nodes and dof 1-3 of their dummy nodes
    n = p%SOLID%BOUNDARY_ngrp_tot
    if( n > 0 ) then
      allocate( src(2*n), is_rot(2*n) )
      m = 0
      do i = 1, n
        ide = mod(p%SOLID%BOUNDARY_ngrp_type(i), 10)
        if( p%SOLID%BOUNDARY_ngrp_rotID(i) <= 0 .and. ide > 3 .and. rot_grp(p%SOLID%BOUNDARY_ngrp_ID(i)) > 0 ) then
          if( p%SOLID%BOUNDARY_ngrp_type(i)/10 <= 3 ) then
            m = m + 1
            src(m) = i
            is_rot(m) = .false.
          endif
          m = m + 1
          src(m) = i
          is_rot(m) = .true.
        else
          m = m + 1
          src(m) = i
          is_rot(m) = .false.
        endif
      enddo
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_GRPID, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_ID, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_type, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_amp, src, m )
      call remap_real_array( p%SOLID%BOUNDARY_ngrp_val, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_istot, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_rotID, src, m )
      call remap_integer_array( p%SOLID%BOUNDARY_ngrp_centerID, src, m )
      p%SOLID%BOUNDARY_ngrp_tot = m
      do k = 1, m
        if( p%SOLID%BOUNDARY_ngrp_rotID(k) > 0 ) cycle
        ig = p%SOLID%BOUNDARY_ngrp_ID(k)
        ids = p%SOLID%BOUNDARY_ngrp_type(k)/10
        ide = mod(p%SOLID%BOUNDARY_ngrp_type(k), 10)
        if( ide <= 3 .or. rot_grp(ig) <= 0 ) cycle
        if( is_rot(k) ) then
          p%SOLID%BOUNDARY_ngrp_ID(k) = rot_grp(ig)
          p%SOLID%BOUNDARY_ngrp_type(k) = 10 * (max(ids, 4) - 3) + (ide - 3)
        else
          p%SOLID%BOUNDARY_ngrp_type(k) = 10 * ids + 3
        endif
      enddo
    endif
 
    do i = 1, p%SOLID%CLOAD_ngrp_tot
      if( p%SOLID%CLOAD_ngrp_rotID(i) > 0 .or. p%SOLID%CLOAD_ngrp_DOF(i) <= 3 ) cycle
      p%SOLID%CLOAD_ngrp_ID(i) = rot_grp(p%SOLID%CLOAD_ngrp_ID(i))
      p%SOLID%CLOAD_ngrp_DOF(i) = p%SOLID%CLOAD_ngrp_DOF(i) - 3
    enddo
  end subroutine fstr_convert_rotation_dof
 
  subroutine rotation_dof_err_stop( header_name )
    implicit none
    character(len=*) :: header_name
 
    call fstr_setup_util_err_stop( '### Error: '//header_name// &
      ' : dof 4-6 of beam/shell nodes in a mesh with solid elements is not supported' )
  end subroutine rotation_dof_err_stop
 
  subroutine explicit_dummy_node_err_stop( header_name )
    implicit none
    character(len=*) :: header_name
 
    call fstr_setup_util_err_stop( '### Error: '//header_name// &
      ' : dof 4-6 is given to beam/shell nodes with explicit dummy nodes; give it as dof 1-3 of the dummy nodes' )
  end subroutine explicit_dummy_node_err_stop
 
  subroutine remap_integer_array( array, src, n )
    implicit none
    integer(kind=kint), pointer :: array(:)
    integer(kind=kint), intent(in) :: src(:), n
    integer(kind=kint), pointer :: new_array(:)
 
    allocate( new_array(n) )
    new_array(1:n) = array(src(1:n))
    deallocate( array )
    array => new_array
  end subroutine remap_integer_array
 
  subroutine remap_real_array( array, src, n )
    implicit none
    real(kind=kreal), pointer :: array(:)
    integer(kind=kint), intent(in) :: src(:), n
    real(kind=kreal), pointer :: new_array(:)
 
    allocate( new_array(n) )
    new_array(1:n) = array(src(1:n))
    deallocate( array )
    array => new_array
  end subroutine remap_real_array
 
end module m_fstr_setup_finish