fstr_dynamic_nlimplicit.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_dynamic_nlimplicit
  use m_fstr
  use m_dynamic_output
  use m_fstr_eig_setmass
  use m_dynamic_mat_ass_bc_ac
  use m_dynamic_mat_ass_bc
  use m_dynamic_mat_ass_bc_vl
  use m_dynamic_mat_ass_load
  use m_fstr_creatematrix_and_dampingforce
  use m_fstr_update
  use m_fstr_restart
  use fstr_matrix_con_contact
  use m_fstr_residual
  use m_fstr_iterationcontrol
  use mcontact
  use m_solve_lineq_contact
  use m_dynamic_init_variables
  use m_fstr_timeinc
  use m_fstr_cutback
  use m_fstr_spring
 
  !-------- for couple -------
  use m_dynamic_mat_ass_couple
  use m_fstr_rcap_io
 
  real(kind=kreal), parameter :: pi = 3.14159265358979323846d0
 
contains
 
  subroutine fstr_solve_nlgeom_dynamic_implicit_contactslag(hecMESH,hecMAT,fstrSOLID,fstrEIG   &
      ,fstrDYNAMIC,fstrRESULT,fstrPARAM &
      ,fstrCPL,hecLagMAT,restart_step_num,restart_substep_num,infoCTChange  &
      ,conMAT,restart_step_count )
    implicit none
    !C-- global variable
    type(hecmwst_local_mesh)             :: hecmesh
    type(hecmwst_matrix)                 :: hecMAT
    type(hecmwst_matrix), pointer        :: hecMAT0
    type(fstr_eigen)                     :: fstrEIG
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_result_data)            :: fstrRESULT ! unused - kept for interface compatibility
    type(fstr_param)                     :: fstrPARAM
    type(fstr_dynamic)                   :: fstrDYNAMIC
    type(fstr_couple)                    :: fstrCPL !for COUPLE
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(fstr_info_contactchange)        :: infoCTChange
    type(hecmwst_matrix)                 :: conMAT
 
    !C-- local variable
    integer(kind=kint) :: nnod, ndof, nn
    integer(kind=kint) :: i, tot_step_print, CBbound
    real(kind=kreal) :: time_1, time_2, factor
    integer(kind=kint) :: sub_step
 
    integer(kind=kint) :: restart_step_num, restart_substep_num, restart_step_count, tot_step, step_count
    integer(kind=kint) :: ctAlgo
    integer(kind=kint) :: max_iter_contact
    real(kind=kreal) :: converg_dlag
    type(fstr_info_contactchange)        :: infoctchange_bak
 
    logical :: is_outpoint
    integer(kind=kint) :: n_node_global
    logical :: is_mat_symmetric, is_interaction_active
 
    if(hecmesh%my_rank==0) call fstr_timeinc_printstatus_init
 
    is_interaction_active = ( associated( fstrsolid%contacts ) .or. associated( fstrsolid%embeds ) )
 
    nullify(hecmat0)
 
    ! sum of n_node among all subdomains (to be used to calc res)
    n_node_global = hecmesh%nn_internal
    call hecmw_allreduce_i1(hecmesh,n_node_global,hecmw_sum)
 
    ctalgo = fstrparam%contact_algo
 
    if( hecmat%Iarray(99)==4 .and. .not.fstr_is_matrixstruct_symmetric(fstrsolid,hecmesh) ) then
      write(*,*) ' This type of direct solver is not yet available in such case ! '
      write(*,*) ' Please use intel MKL direct solver !'
      call  hecmw_abort(hecmw_comm_get_comm())
    endif
 
    hecmat%NDOF=hecmesh%n_dof
 
    nnod=hecmesh%n_node
    ndof=hecmat%NDOF
    nn=ndof*ndof
 
    if( associated( fstrsolid%contacts ) ) then
      call initialize_contact_output_vectors(fstrsolid,hecmat)
      call setup_contact_elesurf_for_area( 1, hecmesh, fstrsolid )
    endif
 
    !!-- initial value
    time_1 = hecmw_wtime()
 
    !C-- check parameters
    if(dabs(fstrdynamic%beta) < 1.0e-20) then
      if( hecmesh%my_rank == 0 ) then
        write(imsg,*) 'stop due to Newmark-beta = 0'
      endif
      call hecmw_abort( hecmw_comm_get_comm())
    endif
 
    !C-- matrix [M] lumped mass matrix
    if(fstrdynamic%idx_mas == 1) then
      call setmass(fstrsolid,hecmesh,hecmat,fstreig)
 
    !C-- consistent mass matrix
    else if(fstrdynamic%idx_mas == 2) then
      if( hecmesh%my_rank .eq. 0 ) then
        write(imsg,*) 'stop: consistent mass matrix is not yet available !'
      endif
      call hecmw_abort( hecmw_comm_get_comm())
    endif
 
    hecmat%Iarray(98) = 1   !Assembly complete
    hecmat%Iarray(97) = 1   !Need numerical factorization
 
    !C-- initialize variables
    if( restart_step_num == 1 .and. fstrdynamic%VarInitialize .and. abs(fstrdynamic%ray_m) > 1.0d-15 ) &
      call dynamic_init_varibles( hecmesh, hecmat, fstrsolid, fstreig, fstrdynamic, fstrparam )
 
    !C-- output of initial state
    if( restart_step_num == 1 ) then
      call fstr_dynamic_output(1, 0, 0.d0, hecmesh, fstrsolid, fstrdynamic, fstrparam, .true.)
      call dynamic_output_monit(1, 0, 0.d0, hecmesh, fstrparam, fstrdynamic, fstreig, fstrsolid)
    endif
 
    fstrdynamic%VEC3(:) =0.d0
    hecmat%X(:) =0.d0
 
    call fstr_save_originalmatrixstructure(hecmat)
    call fstr_scan_contact_state(restart_step_num, restart_step_num, 0, fstrdynamic%t_delta, &
          ctalgo, hecmesh, fstrsolid, infoctchange)
 
    call hecmw_mat_copy_profile( hecmat, conmat )
 
    if ( fstr_is_contact_active() ) then
      call fstr_mat_con_contact( restart_step_num, ctalgo, hecmat, fstrsolid, heclagmat, &
          infoctchange, conmat, fstr_is_contact_active())
    elseif( hecmat%Iarray(99)==4 ) then
      write(*,*) ' This type of direct solver is not yet available in such case ! '
      write(*,*) ' Please change solver type to intel MKL direct solver !'
      call  hecmw_abort(hecmw_comm_get_comm())
    endif
    is_mat_symmetric = fstr_is_matrixstruct_symmetric(fstrsolid,hecmesh)
    call solve_lineq_contact_init(hecmesh,hecmat,heclagmat,is_mat_symmetric)
 
    fstrsolid%FACTOR = 0.0d0
    call fstr_cutback_init( hecmesh, fstrsolid, fstrparam )
    call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak )
    
    step_count = restart_step_count
    do tot_step=1, fstrsolid%nstep_tot
      tot_step_print = tot_step+restart_step_num-1
      if(hecmesh%my_rank==0) write(*,*) ''
      if(hecmesh%my_rank==0) write(*,'(a,i5)') ' loading step=',tot_step_print
 
      sub_step = restart_substep_num
      do while(.true.)
        if (ndof == 4 .and. hecmesh%my_rank==0) write(*,'(a,i5)')"iter: ",sub_step
        ! ----- time history of factor
        call fstr_timeinc_settimeincrement( fstrsolid%step_ctrl(tot_step), fstrparam, sub_step, &
          &  fstrsolid%NRstat_i, fstrsolid%NRstat_r, fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
 
        fstrdynamic%t_curr = fstr_get_time()
        fstrdynamic%t_delta = fstr_get_timeinc()
 
        call fstr_newton_dynamic_contactslag(tot_step, hecmesh, hecmat, fstrsolid, fstreig, &
            fstrdynamic, fstrparam, fstrcpl, heclagmat, infoctchange, conmat, &
            restart_step_num, hecmat0, sub_step, fstrdynamic%t_curr, fstrdynamic%t_delta)
 
        ! Time Increment
        if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus( fstrsolid%step_ctrl(tot_step), fstrparam, &
          &  tot_step_print, sub_step, fstrsolid%NRstat_i, fstrsolid%NRstat_r,   &
          &  fstrsolid%AutoINC_stat, fstrsolid%CutBack_stat )
        if( fstr_cutback_active() ) then
 
          if( fstrsolid%CutBack_stat == 0 ) then ! converged
            call fstr_cutback_save( fstrsolid, infoctchange, infoctchange_bak )  ! save analysis state
            call fstr_proceed_time()             ! current time += time increment
            fstrdynamic%t_curr = fstr_get_time()
 
          else                                   ! not converged
            cbbound = fstrparam%ainc(fstrsolid%step_ctrl(tot_step)%AincParam_id)%CBbound
            if( fstrsolid%CutBack_stat == cbbound ) then
              if( hecmesh%my_rank == 0 ) then
                write(*,*) 'Number of successive cutback reached max number: ',cbbound
                call fstr_timeinc_printstatus_final(.false.)
              endif
              call hecmw_abort( hecmw_comm_get_comm() )
            endif
            call fstr_cutback_load( fstrsolid, infoctchange, infoctchange_bak )  ! load analysis state
            call fstr_set_contact_active( infoctchange%contactNode_current > 0 )
 
            ! restore matrix structure for slagrange contact analysis
            if( is_interaction_active .and. fstrparam%contact_algo == kcaslagrange ) then
              call fstr_mat_con_contact( tot_step, fstrparam%contact_algo, hecmat, fstrsolid, heclagmat, &
                &  infoctchange, conmat, fstr_is_contact_active())
              conmat%B(:) = 0.0d0
              call solve_lineq_contact_init(hecmesh, hecmat, heclagmat, fstr_is_matrixstruct_symmetric(fstrsolid, hecmesh))
            endif
            if( hecmesh%my_rank == 0 ) write(*,*) '### State has been restored at time =',fstr_get_time()
 
            !stop if # of substeps reached upper bound.
            if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
              if( hecmesh%my_rank == 0 ) then
                write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
                  & tot_step_print, '  time=', fstr_get_time()
              endif
              call hecmw_abort( hecmw_comm_get_comm())
            endif
 
            ! output time
            time_2 = hecmw_wtime()
            if( hecmesh%my_rank==0) write(imsg,'(a,",",2(I8,","),f10.2)') &
              &  'step, substep, solve (sec) :', tot_step_print, sub_step, time_2 - time_1
            cycle
          endif
        else
          if( fstrsolid%CutBack_stat > 0 ) stop
          call fstr_proceed_time() ! current time += time increment
          fstrdynamic%t_curr = fstr_get_time()
        endif
 
        step_count = step_count + 1
 
        ! ----- Result output (include visualize output)
        is_outpoint = fstr_timeinc_istimepoint( fstrsolid%step_ctrl(tot_step), fstrparam ) &
          & .or. fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) )
        !C-- output new displacement, velocity and acceleration
        call fstr_dynamic_output(tot_step, step_count, fstrdynamic%t_curr, hecmesh, fstrsolid, fstrdynamic, fstrparam, is_outpoint)
 
        !C-- output result of monitoring node
        call dynamic_output_monit(tot_step, i, fstrdynamic%t_curr, hecmesh, fstrparam, fstrdynamic, fstreig, fstrsolid)
 
        !---  Restart info
        if( fstrdynamic%restart_nout > 0 ) then
          if( mod(step_count,fstrdynamic%restart_nout).eq.0 ) then
            call fstr_write_restart_dyna_nl(tot_step,sub_step,hecmesh,fstrsolid,fstrdynamic,fstrparam,&
              .false.,infoctchange%contactNode_current,step_count)
          endif
        endif
 
        !if time reached the end time of step, exit loop.
        if( fstr_timeinc_isstepfinished( fstrsolid%step_ctrl(tot_step) ) ) exit
 
        if( sub_step == fstrsolid%step_ctrl(tot_step)%num_substep ) then
          if( hecmesh%my_rank == 0 ) then
            write(*,'(a,i5,a,f6.3)') '### Number of substeps reached max number: at total_step=', &
              & tot_step_print, '  time=', fstr_get_time()
          endif
          if( hecmesh%my_rank == 0 ) call fstr_timeinc_printstatus_final(.false.)
          stop !stop if # of substeps reached upper bound.
        endif
 
        sub_step = sub_step + 1
      enddo
 
      !--- Restart at the end of step
      if( fstrdynamic%restart_nout > 0 ) then
        call fstr_write_restart_dyna_nl(tot_step,sub_step,hecmesh,fstrsolid,fstrdynamic,fstrparam,&
          .true.,infoctchange%contactNode_current,step_count)
      endif
 
      restart_substep_num = 1
      !C-- end of time step loop
    enddo
 
    if (associated(hecmat0)) then
      call hecmw_mat_finalize(hecmat0)
      deallocate(hecmat0)
    endif
 
    !  message
    if( hecmesh%my_rank == 0 ) then
      call fstr_timeinc_printstatus_final(.true.)
      write(imsg,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
      write(*,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
    endif
 
  end subroutine fstr_solve_nlgeom_dynamic_implicit_contactslag
 
  subroutine fstr_newton_dynamic_contactslag(cstep, hecMESH, hecMAT, fstrSOLID, fstrEIG, &
      fstrDYNAMIC, fstrPARAM, fstrCPL, hecLagMAT, infoCTChange, conMAT, &
      restart_step_num, hecMAT0, istep, t_curr, t_delta)
    implicit none
    !C-- arguments
    integer(kind=kint), intent(in)       :: cstep, restart_step_num, istep
    real(kind=kreal), intent(in)         :: t_curr, t_delta
    type(hecmwst_local_mesh)             :: hecmesh
    type(hecmwst_matrix)                 :: hecMAT
    type(hecmwst_matrix), pointer        :: hecMAT0
    type(fstr_solid)                     :: fstrSOLID
    type(fstr_eigen)                     :: fstrEIG
    type(fstr_dynamic)                   :: fstrDYNAMIC
    type(fstr_param)                     :: fstrPARAM
    type(fstr_couple)                    :: fstrCPL
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(fstr_info_contactchange)        :: infoCTChange
    type(hecmwst_matrix)                 :: conMAT
 
    !C-- local variables
    integer(kind=kint) :: j, kk, idm, imm
    integer(kind=kint) :: iter
    real(kind=kreal) :: a1, a2, a3, b1, b2, b3, c1, c2
    real(kind=kreal) :: coef(6)
    real(kind=kreal) :: res, res1, res0, relres
    integer(kind=kint) :: count_step, stepcnt
    real(kind=kreal) :: maxdlag
    integer(kind=kint) :: contact_changed_global
    logical :: is_mat_symmetric
    integer :: istat
    logical :: is_cycle
    integer(kind=kint) :: ctAlgo, max_iter_contact
    integer(kind=kint) :: nnod, ndof, nn
    real(kind=kreal) :: converg_dlag
    real(kind=kreal), allocatable :: coord(:)
    integer(kind=kint) :: iterStatus, nresid, n_node_global
    real(kind=kreal), allocatable :: resid_work(:)
 
    fstrsolid%NRstat_i(:) = 0 ! logging newton iteration(init)
 
    !C-- initialize local variables
    n_node_global = hecmesh%nn_internal
    call hecmw_allreduce_i1(hecmesh, n_node_global, hecmw_sum)
    ctalgo = fstrparam%contact_algo
    max_iter_contact = fstrsolid%step_ctrl(cstep)%max_contiter
    converg_dlag = fstrsolid%step_ctrl(cstep)%converg_lag
    nnod = hecmesh%n_node
    ndof = hecmat%NDOF
    nn = ndof*ndof
    allocate(coord(hecmesh%n_node*ndof))
    allocate(resid_work(hecmesh%n_node*ndof + conmat%NP*ndof))
 
    a1 = .5d0/fstrdynamic%beta - 1.d0
    a2 = 1.d0/(fstrdynamic%beta*t_delta)
    a3 = 1.d0/(fstrdynamic%beta*t_delta*t_delta)
    b1 = ( .5d0*fstrdynamic%gamma/fstrdynamic%beta - 1.d0 )*t_delta
    b2 = fstrdynamic%gamma/fstrdynamic%beta - 1.d0
    b3 = fstrdynamic%gamma/(fstrdynamic%beta*t_delta)
    c1 = 1.d0 + fstrdynamic%ray_k*b3
    c2 = a3 + fstrdynamic%ray_m*b3
 
    coef(1) = a1; coef(2) = a2; coef(3) = a3
    coef(4) = b1; coef(5) = b2; coef(6) = b3
 
    if(hecmesh%my_rank==0) then
      write(*,'(A)')'-------------------------------------------------'
      write(*,'('' time step='',i10,'' time='',1pe13.4e3)') istep,t_curr
    endif
 
    do j = 1 ,ndof*nnod
      fstrdynamic%VEC1(j) = a1*fstrdynamic%ACC(j,1) + a2*fstrdynamic%VEL(j,1)
      fstrdynamic%VEC2(j) = b1*fstrdynamic%ACC(j,1) + b2*fstrdynamic%VEL(j,1)
    enddo
 
    count_step = 0
    stepcnt = 0
 
    !C for couple analysis
    do
      fstrsolid%dunode(:) =0.d0
      ! call fstr_UpdateEPState( hecMESH, fstrSOLID )
      call fstr_solve_dynamic_nlimplicit_couple_init(fstrparam, fstrcpl)
 
    loopforcontactanalysis: do while( .true. )
    count_step = count_step + 1
 
      ! ----- Inner Iteration
      res0   = 0.d0
      res1   = 0.d0
      relres = 1.d0
 
      do iter = 1, fstrsolid%step_ctrl(cstep)%max_iter
        stepcnt=stepcnt+1
        call fstr_creatematrix_and_dampingforce( hecmesh, hecmat, fstrsolid, t_curr, t_delta, fstrdynamic, coef )
 
        call fstr_addspring(cstep, hecmesh, hecmat, fstrsolid, fstrparam)
 
        !C-- mechanical boundary condition
        call dynamic_mat_ass_load (cstep, t_curr+t_delta, hecmesh, hecmat, fstrsolid, fstrdynamic, fstrparam, iter )
        do j=1, hecmesh%n_node*  hecmesh%n_dof
          hecmat%B(j) = hecmat%B(j) - fstrsolid%QFORCE(j) + fstrsolid%DFORCE(j)
        enddo
 
        call fstr_update_ndforce_spring( cstep, hecmesh, fstrsolid, hecmat%B )
 
        !C for couple analysis
        call fstr_solve_dynamic_nlimplicit_couple_pre(hecmesh, hecmat, fstrsolid, &
          & fstrparam, fstrdynamic, fstrcpl, restart_step_num, istep)
 
        call hecmw_mat_clear( conmat )
        call hecmw_mat_clear_b( conmat )
        conmat%X = 0.0d0
 
        if( fstr_is_contact_active() ) then
          call fstr_update_ndforce_contact(cstep,ctalgo,hecmesh,heclagmat,fstrsolid,conmat)
          call fstr_addcontactstiffness(cstep,ctalgo,iter,hecmesh,conmat,heclagmat,fstrsolid)
        endif
 
        !C-- geometrical boundary condition
        call dynamic_mat_ass_bc   (cstep, hecmesh, hecmat, fstrsolid, fstrdynamic, &
          &  fstrparam, heclagmat, t_curr+t_delta, stepcnt, conmat=conmat)
        call dynamic_mat_ass_bc_vl(cstep, hecmesh, hecmat, fstrsolid, fstrdynamic, &
          &  fstrparam, heclagmat, t_curr+t_delta, stepcnt, conmat=conmat)
        call dynamic_mat_ass_bc_ac(cstep, hecmesh, hecmat, fstrsolid, fstrdynamic, &
          &  fstrparam, heclagmat, t_curr+t_delta, stepcnt, conmat=conmat)
 
        ! ----- check convergence
        call fstr_assemble_residual_contact(hecmat, heclagmat, conmat, hecmesh, resid_work, nresid)
 
        if( .not.fstr_is_contact_active() ) then
          maxdlag = 0.0d0
        elseif( abs(maxdlag) < 1.0d-15) then
          maxdlag = 1.0d0
        endif
        call hecmw_allreduce_r1(hecmesh, maxdlag, hecmw_max)
 
        call fstr_check_convergence(hecmesh, hecmat, fstrsolid, fstrpr, &
            ndof, iter, istep, cstep, &
            resid_work, nresid, &
            res0, res, &
            n_node_global, &
            iterstatus, &
            maxdlag, converg_dlag)
        if (iterstatus == kitrconverged) exit
        if (iterstatus == kitrdiverged .or. iterstatus == kitrfloatingerror) then
          fstrsolid%NRstat_i(knstciter) = count_step
          return
        endif
 
        !   ----  For Parallel Contact with Multi-Partition Domains
        hecmat%X = 0.0d0
        call fstr_set_current_config_to_mesh(hecmesh,fstrsolid,coord)
        call solve_lineq_contact(hecmesh,hecmat,heclagmat,conmat,istat,1.0d0,fstr_is_contact_active())
        call fstr_recover_initial_config_to_mesh(hecmesh,fstrsolid,coord)
        ! ----- check matrix solver error
        if( istat /= 0 ) then
          if( hecmesh%my_rank == 0) then
            write(   *,'(a,i5,a,i5)') '     ### Fail to Converge  : at total_step=', cstep, '  sub_step=', istep
          end if
          fstrsolid%NRstat_i(knstdresn) = 4
          fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
          return
        end if
 
        ! ----- update external nodal displacement increments
        call hecmw_update_r (hecmesh, hecmat%X, hecmat%NP, hecmat%NDOF)
 
        ! ----- update the strain, stress, and internal force
        do j=1,hecmesh%n_node*ndof
          fstrsolid%dunode(j)  = fstrsolid%dunode(j)+hecmat%X(j)
        enddo
        call fstr_updatenewton( hecmesh, hecmat, fstrsolid, t_curr, &
          &   t_delta,iter, fstrdynamic%strainEnergy )
 
        ! ----- update reaction force at constrained DOFs using converged QFORCE
        call fstr_update_reaction_spc( cstep, hecmesh, fstrsolid )
 
        ! ----- update the Lagrange multipliers
        ! Done before the nlgeom early-exit so the multipliers are accumulated
        ! even for infinitesimal analysis. Otherwise hecLagMAT%lagrange stays
        ! zero when NONLINEAR is not specified, and the contact output forces
        ! CONT_NFORCE/CONT_FRIC (computed from the multipliers) are always zero
        ! even though contact is enforced. Matches the unconditional update in
        ! the static path (fstr_solve_NonLinear).
        if( fstr_is_contact_active() ) then
          maxdlag = 0.0d0
          do j=1,heclagmat%num_lagrange
            heclagmat%lagrange(j) = heclagmat%lagrange(j) + hecmat%X(hecmesh%n_node*ndof+j)
            if(dabs(hecmat%X(hecmesh%n_node*ndof+j))>maxdlag) maxdlag=dabs(hecmat%X(hecmesh%n_node*ndof+j))
          enddo
        endif
 
        if(.not. fstrparam%nlgeom) exit
      enddo
 
      fstrsolid%NRstat_i(knstmaxit) = max(fstrsolid%NRstat_i(knstmaxit),iter) ! logging newton iteration(maxtier)
      fstrsolid%NRstat_i(knstsumit) = fstrsolid%NRstat_i(knstsumit) + iter    ! logging newton iteration(sum of iter)
 
      ! ----- compute CONT_NFORCE/CONT_FRIC for output
      if( fstr_is_contact_active() ) &
        call fstr_calc_contact_output_force(cstep,ctalgo,hecmesh,heclagmat,fstrsolid,conmat)
 
      call fstr_scan_contact_state(cstep, istep, count_step, t_delta, ctalgo, &
                                   hecmesh, fstrsolid, infoctchange)
 
      if( hecmat%Iarray(99)==4 .and. .not. fstr_is_contact_active() ) then
        write(*,*) ' This type of direct solver is not yet available in such case ! '
        write(*,*) ' Please use intel MKL direct solver !'
        call  hecmw_abort(hecmw_comm_get_comm())
      endif
 
      is_mat_symmetric = fstr_is_matrixstruct_symmetric(fstrsolid,hecmesh)
      contact_changed_global=0
      if( fstr_is_contact_conv(ctalgo,infoctchange,hecmesh) ) then
        exit loopforcontactanalysis
      elseif( fstr_is_matrixstructure_changed(infoctchange) ) then
        call fstr_mat_con_contact( cstep, ctalgo, hecmat, fstrsolid, heclagmat, infoctchange, conmat, fstr_is_contact_active())
        contact_changed_global=1
      endif
      call hecmw_allreduce_i1(hecmesh,contact_changed_global,hecmw_max)
      if (contact_changed_global > 0) then
        call hecmw_mat_clear_b( hecmat )
        call hecmw_mat_clear_b( conmat )
        call solve_lineq_contact_init(hecmesh,hecmat,heclagmat,is_mat_symmetric)
      endif
 
      if( count_step > max_iter_contact ) exit loopforcontactanalysis
 
      ! ----- check divergence
      if( count_step >= fstrsolid%step_ctrl(cstep)%max_contiter ) then
        if( hecmesh%my_rank == 0) then
          write(   *,'(a,i5,a,i5)') '     ### Contact failed to Converge  : at total_step=', cstep, '  sub_step=', istep
        end if
        fstrsolid%NRstat_i(knstciter) = count_step                              ! logging contact iteration
        fstrsolid%CutBack_stat = fstrsolid%CutBack_stat + 1
        fstrsolid%NRstat_i(knstdresn) = 3
        return
      end if
 
    enddo loopforcontactanalysis
 
    fstrsolid%NRstat_i(knstciter) = count_step ! logging contact iteration
 
    !C for couple analysis
      call fstr_solve_dynamic_nlimplicit_couple_post(hecmesh, hecmat, fstrsolid, &
        & fstrparam, fstrdynamic, fstrcpl, a1, a2, a3, b1, b2, b3, istep, is_cycle)
      if(is_cycle) cycle
      exit
    enddo
 
    !C-- new displacement, velocity and acceleration
    fstrdynamic%kineticEnergy = 0.0d0
    do j = 1 ,ndof*nnod
      fstrdynamic%ACC (j,2) = -a1*fstrdynamic%ACC(j,1) - a2*fstrdynamic%VEL(j,1) + &
        a3*fstrsolid%dunode(j)
      fstrdynamic%VEL (j,2) = -b1*fstrdynamic%ACC(j,1) - b2*fstrdynamic%VEL(j,1) + &
        b3*fstrsolid%dunode(j)
      fstrdynamic%ACC (j,1) = fstrdynamic%ACC (j,2)
      fstrdynamic%VEL (j,1) = fstrdynamic%VEL (j,2)
 
      fstrsolid%unode(j)  = fstrsolid%unode(j)+fstrsolid%dunode(j)
      fstrdynamic%DISP(j,2) = fstrsolid%unode(j)
 
      fstrdynamic%kineticEnergy = fstrdynamic%kineticEnergy + &
        0.5d0*fstreig%mass(j)*fstrdynamic%VEL(j,2)*fstrdynamic%VEL(j,2)
    enddo
 
    call fstr_updatestate( hecmesh, fstrsolid, t_delta )
    call fstr_update_contact_tangentforce( cstep, fstrsolid )
 
    deallocate(coord)
    deallocate(resid_work)
    fstrsolid%CutBack_stat = 0
  end subroutine fstr_newton_dynamic_contactslag
 
  subroutine fstr_solve_dynamic_nlimplicit_couple_init(fstrPARAM, fstrCPL)
    implicit none
    type(fstr_param)         :: fstrparam
    type(fstr_couple)        :: fstrCPL
    if( fstrparam%fg_couple == 1) then
      if( fstrparam%fg_couple_type==1 .or. &
          fstrparam%fg_couple_type==3 .or. &
          fstrparam%fg_couple_type==5 ) call fstr_rcap_get( fstrcpl )
    endif
  end subroutine fstr_solve_dynamic_nlimplicit_couple_init
 
  subroutine fstr_solve_dynamic_nlimplicit_couple_pre(hecMESH, hecMAT, fstrSOLID, &
    & fstrPARAM, fstrDYNAMIC, fstrCPL, restart_step_num, i)
    implicit none
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_solid)         :: fstrSOLID
    type(fstr_param)         :: fstrPARAM
    type(fstr_dynamic)       :: fstrDYNAMIC
    type(fstr_couple)        :: fstrCPL
    integer(kint) :: kkk0, kkk1, j, kk, i, restart_step_num
    real(kreal) :: bsize
 
    if( fstrparam%fg_couple == 1) then
      if( fstrparam%fg_couple_first /= 0 ) then
        bsize = dfloat( i ) / dfloat( fstrparam%fg_couple_first )
        if( bsize > 1.0 ) bsize = 1.0
        do kkk0 = 1, fstrcpl%coupled_node_n
          kkk1 = 3 * kkk0
          fstrcpl%trac(kkk1-2) = bsize * fstrcpl%trac(kkk1-2)
          fstrcpl%trac(kkk1-1) = bsize * fstrcpl%trac(kkk1-1)
          fstrcpl%trac(kkk1  ) = bsize * fstrcpl%trac(kkk1  )
        enddo
      endif
      if( fstrparam%fg_couple_window > 0 ) then
        j = i - restart_step_num + 1
        kk = fstrdynamic%n_step - restart_step_num + 1
        bsize = 0.5*(1.0-cos(2.0*pi*dfloat(j)/dfloat(kk)))
        do kkk0 = 1, fstrcpl%coupled_node_n
          kkk1 = 3 * kkk0
          fstrcpl%trac(kkk1-2) = bsize * fstrcpl%trac(kkk1-2)
          fstrcpl%trac(kkk1-1) = bsize * fstrcpl%trac(kkk1-1)
          fstrcpl%trac(kkk1  ) = bsize * fstrcpl%trac(kkk1  )
        enddo
      endif
      call dynamic_mat_ass_couple( hecmesh, hecmat, fstrsolid, fstrcpl )
    endif
  end subroutine fstr_solve_dynamic_nlimplicit_couple_pre
 
  subroutine fstr_solve_dynamic_nlimplicit_couple_post(hecMESH, hecMAT, fstrSOLID, &
    & fstrPARAM, fstrDYNAMIC, fstrCPL, a1, a2, a3, b1, b2, b3, i, is_cycle)
    implicit none
    type(hecmwst_local_mesh) :: hecMESH ! unused - kept for interface compatibility
    type(hecmwst_matrix)     :: hecMAT ! unused - kept for interface compatibility
    type(fstr_solid)         :: fstrSOLID
    type(fstr_param)         :: fstrPARAM
    type(fstr_dynamic)       :: fstrDYNAMIC
    type(fstr_couple)        :: fstrCPL
    integer(kint) :: kkk0, kkk1, j, i, revocap_flag
    real(kreal) :: a1, a2, a3, b1, b2, b3
    logical :: is_cycle
 
    is_cycle = .false.
 
    if( fstrparam%fg_couple == 1 ) then
      if( fstrparam%fg_couple_type>1 ) then
        do j=1, fstrcpl%coupled_node_n
          if( fstrcpl%dof == 3 ) then
            kkk0 = j*3
            kkk1 = fstrcpl%coupled_node(j)*3
 
            fstrcpl%disp (kkk0-2) = fstrsolid%unode(kkk1-2) + fstrsolid%dunode(kkk1-2)
            fstrcpl%disp (kkk0-1) = fstrsolid%unode(kkk1-1) + fstrsolid%dunode(kkk1-1)
            fstrcpl%disp (kkk0  ) = fstrsolid%unode(kkk1  ) + fstrsolid%dunode(kkk1  )
 
            fstrcpl%velo (kkk0-2) = -b1*fstrdynamic%ACC(kkk1-2,1) - b2*fstrdynamic%VEL(kkk1-2,1) + &
              b3*fstrsolid%dunode(kkk1-2)
            fstrcpl%velo (kkk0-1) = -b1*fstrdynamic%ACC(kkk1-1,1) - b2*fstrdynamic%VEL(kkk1-1,1) + &
              b3*fstrsolid%dunode(kkk1-1)
            fstrcpl%velo (kkk0  ) = -b1*fstrdynamic%ACC(kkk1,1) - b2*fstrdynamic%VEL(kkk1,1) + &
              b3*fstrsolid%dunode(kkk1)
            fstrcpl%accel(kkk0-2) = -a1*fstrdynamic%ACC(kkk1-2,1) - a2*fstrdynamic%VEL(kkk1-2,1) + &
              a3*fstrsolid%dunode(kkk1-2)
            fstrcpl%accel(kkk0-1) = -a1*fstrdynamic%ACC(kkk1-1,1) - a2*fstrdynamic%VEL(kkk1-1,1) + &
              a3*fstrsolid%dunode(kkk1-1)
            fstrcpl%accel(kkk0  ) = -a1*fstrdynamic%ACC(kkk1,1) - a2*fstrdynamic%VEL(kkk1,1) + &
              a3*fstrsolid%dunode(kkk1)
          else
            kkk0 = j*2
            kkk1 = fstrcpl%coupled_node(j)*2
 
            fstrcpl%disp (kkk0-1) = fstrsolid%unode(kkk1-1) + fstrsolid%dunode(kkk1-1)
            fstrcpl%disp (kkk0  ) = fstrsolid%unode(kkk1  ) + fstrsolid%dunode(kkk1  )
 
            fstrcpl%velo (kkk0-1) = -b1*fstrdynamic%ACC(kkk1-1,1) - b2*fstrdynamic%VEL(kkk1-1,1) + &
              b3*fstrsolid%dunode(kkk1-1)
            fstrcpl%velo (kkk0  ) = -b1*fstrdynamic%ACC(kkk1,1) - b2*fstrdynamic%VEL(kkk1,1) + &
              b3*fstrsolid%dunode(kkk1)
            fstrcpl%accel(kkk0-1) = -a1*fstrdynamic%ACC(kkk1-1,1) - a2*fstrdynamic%VEL(kkk1-1,1) + &
              a3*fstrsolid%dunode(kkk1-1)
            fstrcpl%accel(kkk0  ) = -a1*fstrdynamic%ACC(kkk1,1) - a2*fstrdynamic%VEL(kkk1,1) + &
              a3*fstrsolid%dunode(kkk1)
          endif
        enddo
        call fstr_rcap_send( fstrcpl )
      endif
 
      select case ( fstrparam%fg_couple_type )
        case (4)
          call fstr_rcap_get( fstrcpl )
        case (5)
          call fstr_get_convergence( revocap_flag )
          if( revocap_flag==0 ) is_cycle = .true.
        case (6)
          call fstr_get_convergence( revocap_flag )
          if( revocap_flag==0 ) then
            call fstr_rcap_get( fstrcpl )
            is_cycle = .true.
          else
            if( i /= fstrdynamic%n_step ) call fstr_rcap_get( fstrcpl )
          endif
      end select
    endif
  end subroutine fstr_solve_dynamic_nlimplicit_couple_post
 
end module fstr_dynamic_nlimplicit