fstr_contact.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 mcontact
 
  use mcontactdef
  use hecmw
  use m_fstr
  use m_fstr_contact_assembly
  use m_fstr_contact_element
  use m_fstr_contact_geom
  use m_fstr_contact_search
  use m_fstr_contact_damping, only: is_damping_enabled
  use m_fstr_contact_mpc
  use m_fstr_contact_output
  implicit none
 
  logical, private :: active
 
  ! Re-export functions from sub-modules
  public :: initialize_contact_output_vectors
  public :: initialize_embed_vectors
  public :: setup_contact_elesurf_for_area
  public :: calc_contact_area
  public :: fstr_setup_parancon_contactvalue
  public :: update_contact_state_vectors
  public :: fstr_update_contact_state_vectors
  public :: fstr_calc_contact_output_force_exp
 
contains
 
  subroutine fstr_addcontactstiffness(cstep,ctAlgo, iter,hecMESH,conMAT,hecLagMAT,fstrSOLID)
 
    integer(kind=kint)                   :: cstep
    integer(kind=kint)                   :: ctAlgo
    integer(kind=kint)                   :: iter
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: conMAT
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    integer(kind=kint)                   :: i, grpid
 
    if( associated(heclagmat%AL_lagrange) ) heclagmat%AL_lagrange = 0.0d0
    if( associated(heclagmat%AU_lagrange) ) heclagmat%AU_lagrange = 0.0d0
 
    do i = 1, fstrsolid%n_contacts
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      if( fstrsolid%contacts(i)%method == contacts2s ) then
        call calcu_contact_stiffness_surfsurf( ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
          fstrsolid%dunode(:), conmat )
      else
        call calcu_contact_stiffness_nodesurf( ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
          fstrsolid%dunode(:), iter, heclagmat%Lagrange(:), conmat, heclagmat)
      endif
    enddo
 
    do i = 1, fstrsolid%n_embeds
 
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
 
      call calcu_contact_stiffness_nodesurf( ctalgo, fstrsolid%embeds(i), hecmesh%node(:), fstrsolid%unode(:), &
        fstrsolid%dunode(:), iter, heclagmat%Lagrange(:), conmat, heclagmat)
 
    enddo
 
  end subroutine fstr_addcontactstiffness
 
  subroutine fstr_update_ndforce_contact(cstep,ctAlgo,hecMESH,hecLagMAT,fstrSOLID,conMAT)
    integer(kind=kint), intent(in)       :: cstep
    integer(kind=kint), intent(in)       :: ctAlgo
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(hecmwst_matrix)                 :: conMAT
 
    call fstr_contact_ndforce_core(kctforresidual,cstep,ctalgo,hecmesh,heclagmat,fstrsolid,conmat)
 
  end subroutine fstr_update_ndforce_contact
 
  subroutine fstr_calc_contact_output_force(cstep,ctAlgo,hecMESH,hecLagMAT,fstrSOLID,conMAT)
    integer(kind=kint), intent(in)       :: cstep
    integer(kind=kint), intent(in)       :: ctAlgo
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(hecmwst_matrix)                 :: conMAT
 
    call fstr_contact_ndforce_core(kctforoutput,cstep,ctalgo,hecmesh,heclagmat,fstrsolid,conmat)
 
  end subroutine fstr_calc_contact_output_force
 
  subroutine fstr_calc_contact_output_force_exp( hecMESH, fstrSOLID )
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_solid), intent(inout)      :: fstrsolid
 
    integer(kind=kint) :: i
 
    if( .not. associated(fstrsolid%CONT_NFORCE) ) return
    fstrsolid%CONT_NFORCE(:) = 0.d0
    if( associated(fstrsolid%CONT_FRIC) ) fstrsolid%CONT_FRIC(:) = 0.d0
 
    do i = 1, fstrsolid%n_contacts
      call calcu_contact_ndforce_exp( fstrsolid%contacts(i), hecmesh%node(:), &
        fstrsolid%unode(:), fstrsolid%dunode(:), fstrsolid%CONT_NFORCE, fstrsolid%CONT_FRIC )
    enddo
 
  end subroutine fstr_calc_contact_output_force_exp
 
  subroutine fstr_contact_ndforce_core(purpose,cstep,ctAlgo,hecMESH,hecLagMAT,fstrSOLID,conMAT)
    integer(kind=kint), intent(in)       :: purpose
    integer(kind=kint), intent(in)       :: cstep
    integer(kind=kint), intent(in)       :: ctAlgo
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(hecmwst_matrix)                 :: conMAT
    integer(kind=kint) :: i, grpid
 
    ! Zero-clear output arrays only when computing for output
    if( purpose == kctforoutput ) then
      if( associated(fstrsolid%CONT_NFORCE) ) fstrsolid%CONT_NFORCE(:) = 0.d0
      if( associated(fstrsolid%CONT_FRIC) )   fstrsolid%CONT_FRIC(:) = 0.d0
      if( associated(fstrsolid%EMBED_NFORCE) ) fstrsolid%EMBED_NFORCE(:) = 0.d0
    endif
 
    do i = 1, fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
      if( fstrsolid%contacts(i)%method == contacts2s ) then
        call calcu_contact_ndforce_surfsurf( purpose, ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
          fstrsolid%dunode(:), conmat, fstrsolid%CONT_NFORCE, fstrsolid%CONT_FRIC )
      else
        call calcu_contact_ndforce_nodesurf( purpose, ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
          fstrsolid%dunode(:), heclagmat%Lagrange(:), conmat, &
          fstrsolid%CONT_NFORCE, fstrsolid%CONT_FRIC )
      endif
    enddo
 
    do i = 1, fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
      call calcu_contact_ndforce_nodesurf( purpose, ctalgo, fstrsolid%embeds(i), hecmesh%node(:), fstrsolid%unode(:), &
        fstrsolid%dunode(:), heclagmat%Lagrange(:), conmat, &
        fstrsolid%EMBED_NFORCE, fstrsolid%EMBED_NFORCE )
    enddo
 
  end subroutine fstr_contact_ndforce_core
 
  subroutine fstr_calc_contact_refstiff(cstep, hecMESH, hecMAT, fstrSOLID)
    use hecmw_matrix_misc
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: hecMAT
    type(fstr_solid)                     :: fstrSOLID
    
    integer(kind=kint) :: i, grpid, ndof
    real(kind=kreal), pointer :: diag(:)
    
    ndof = hecmat%NDOF
    
    ! Extract diagonal components once (efficient, called only once)
    diag => hecmw_mat_diag(hecmat)
    
    ! Loop over contact pairs
    do i = 1, fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if(.not. fstr_iscontactactive(fstrsolid, grpid, cstep)) cycle
      
      call calc_contact_pair_refstiff(fstrsolid%contacts(i), diag, ndof, hecmesh)
    enddo
    
    ! Same for embeds
    do i = 1, fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if(.not. fstr_isembedactive(fstrsolid, grpid, cstep)) cycle
      
      call calc_contact_pair_refstiff(fstrsolid%embeds(i), diag, ndof, hecmesh)
    enddo
    
    deallocate(diag)
  end subroutine fstr_calc_contact_refstiff
 
  subroutine fstr_scan_contact_state( cstep, sub_step, cont_step, ctAlgo, hecMESH, fstrSOLID, infoCTChange )
    integer(kind=kint), intent(in)         :: cstep
    integer(kind=kint), intent(in)         :: sub_step
    integer(kind=kint), intent(in)         :: cont_step
    integer(kind=kint), intent(in)         :: ctAlgo
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    type(fstr_info_contactchange), intent(inout):: infoCTChange
    character(len=9)                       :: flag_ctAlgo
    integer(kind=kint) :: i, j, k, grpid
    integer(kind=kint) :: s_chg(3,3), s_emov, s_ebey, s_islid, s_act
    integer(kind=kint), allocatable :: in_contact(:)
    logical :: iactive, is_init
 
    if( associated( fstrsolid%CONT_RELVEL ) ) fstrsolid%CONT_RELVEL(:) = 0.d0
    if( associated( fstrsolid%CONT_STATE ) ) fstrsolid%CONT_STATE(:) = 0.d0
 
    if( ctalgo == kcaslagrange ) then
      flag_ctalgo = 'SLagrange'
    elseif( ctalgo == kcaalagrange ) then
      flag_ctalgo = 'ALagrange'
    endif
 
    ! P.A. We redefine fstrSOLID%ddunode as current coordinate of every nodes
    !  fstrSOLID%ddunode(:) = fstrSOLID%unode(:) + fstrSOLID%dunode(:)
    do i = 1, size(fstrsolid%unode)
      fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
    enddo
    active = .false.
 
    infoctchange%n_statechange = 0
    infoctchange%contact2neighbor = 0
    infoctchange%contact2beyond = 0
    infoctchange%contact2diffLpos = 0
    infoctchange%free2contact_new = 0
    infoctchange%contactNode_current = 0
 
    is_init = ( cstep == 1 .and. sub_step == 1 .and. cont_step == 0 )
 
    do i=1,fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) then
        call clear_contact_state(fstrsolid%contacts(i));  cycle
      endif
      if( fstrsolid%contacts(i)%method == contacts2s ) then
        call scan_contact_state_ss( fstrsolid%contacts(i), fstrsolid%ddunode(:), infoctchange, is_init, iactive )
      else
        call scan_contact_state( fstrsolid%contacts(i), &
        & fstrsolid%ddunode(:), fstrsolid%dunode(:), infoctchange, &
        & hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), &
        & is_init, iactive, hecmesh=hecmesh, &
        & flag_ctalgo=flag_ctalgo, ndforce=fstrsolid%QFORCE(:) )
      endif
      if( .not. active ) active = iactive
    enddo
 
    do i=1,fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) then
        call clear_contact_state(fstrsolid%embeds(i));  cycle
      endif
      call scan_embed_state( fstrsolid%embeds(i), fstrsolid%ddunode(:), &
      & infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive )
      if( .not. active ) active = iactive
    enddo
    ! kept in infoCTChange as well, so that a cutback restores it together with the contact state
    infoctchange%active = active
 
    if( is_init .and. ctalgo == kcaslagrange .and. fstrsolid%n_contacts > 0 ) &
    &  call remove_duplication_tiedcontact( cstep, hecmesh, fstrsolid, infoctchange )
 
    infoctchange%contactNode_current = infoctchange%contactNode_previous &
      + sum(infoctchange%n_statechange(:,kcatcont)) - sum(infoctchange%n_statechange(kcatcont,:))
    infoctchange%contactNode_previous = infoctchange%contactNode_current
 
    ! Output summary of contact state changes (always on; per-node detail requires CONTACT_LOG_LEVEL>=1)
    ! allreduce must be called by all ranks regardless of local counts
    s_chg   = infoctchange%n_statechange
    s_chg(kcatfree,kcatcont) = s_chg(kcatfree,kcatcont) + infoctchange%free2contact_new
    s_emov  = infoctchange%contact2neighbor
    s_ebey  = infoctchange%contact2beyond
    s_islid = infoctchange%contact2diffLpos
    s_act   = infoctchange%contactNode_current
    ! SURF-SURF pairs add the slave nodes of their segments in contact, as update_contact_state_vectors marks them.
    ! A slave segment is kept by one rank only, so the marks are summed onto the owner of each node and counted there.
    allocate( in_contact(hecmesh%n_node) )
    in_contact(:) = 0
    do i = 1, fstrsolid%n_contacts
      if( fstrsolid%contacts(i)%method /= contacts2s ) cycle
      do j = 1, size(fstrsolid%contacts(i)%slave_surf)
        if( fstrsolid%contacts(i)%slave_surf(j)%state <= 0 ) cycle
        do k = 1, size(fstrsolid%contacts(i)%slave_surf(j)%nodes)
          in_contact(fstrsolid%contacts(i)%slave_surf(j)%nodes(k)) = 1
        enddo
      enddo
    enddo
    call hecmw_assemble_i(hecmesh, in_contact, hecmesh%n_node, 1)
    do i = 1, hecmesh%nn_internal
      if( in_contact(i) /= 0 ) s_act = s_act + 1
    enddo
    deallocate( in_contact )
    call hecmw_allreduce_i(hecmesh, s_chg, 9, hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_emov,  hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_ebey,  hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_islid, hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_act,   hecmw_sum)
    if (hecmw_comm_get_rank() == 0) then
      write(*,'(A,i7,A,i7,A,i7,A,i7,A,i7,A,i7)') ' Contact change: Free2Cont=', s_chg(kcatfree,kcatcont), &
        ' Cont2Free=', s_chg(kcatcont,kcatfree), &
        ' InnerSlid=', s_islid, ' MoveNeibr=', s_emov, ' MoveBynd=', s_ebey, ' ActiveNodes=', s_act
      if( sum(s_chg(kcatnear,:)) + sum(s_chg(:,kcatnear)) > 0 ) &
        write(*,'(A,i7,A,i7,A,i7,A,i7)') '   near state change: Free2Near=', s_chg(kcatfree,kcatnear), &
        ' Near2Free=', s_chg(kcatnear,kcatfree), &
        ' Near2Cont=', s_chg(kcatnear,kcatcont), ' Cont2Near=', s_chg(kcatcont,kcatnear)
    end if
 
    if( .not. active ) then
      if( associated( fstrsolid%CONT_NFORCE ) ) fstrsolid%CONT_NFORCE(:) = 0.d0
      if( associated( fstrsolid%CONT_FRIC ) ) fstrsolid%CONT_FRIC(:) = 0.d0
    end if
 
  end subroutine
 
  subroutine fstr_scan_contact_state_exp( cstep, is_init, hecMESH, fstrSOLID, infoCTChange )
    integer(kind=kint), intent(in)               :: cstep
    logical, intent(in)                          :: is_init
    type( hecmwst_local_mesh ), intent(in)       :: hecMESH
    type(fstr_solid), intent(inout)              :: fstrSOLID
    type(fstr_info_contactchange), intent(inout) :: infoCTChange
 
    integer(kind=kint) :: i, grpid
    integer(kind=kint) :: s_chg(3,3), s_emov, s_ebey, s_islid, s_act
    logical :: iactive
 
 
    ! P.A. We redefine fstrSOLID%ddunode as current coordinate of every nodes
    !  fstrSOLID%ddunode(:) = fstrSOLID%unode(:) + fstrSOLID%dunode(:)
    do i = 1, size(fstrsolid%unode)
      fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
    enddo
    infoctchange%active = .false.
 
    infoctchange%n_statechange = 0
    infoctchange%contact2neighbor = 0
    infoctchange%contact2beyond = 0
    infoctchange%contact2diffLpos = 0
    infoctchange%free2contact_new = 0
    infoctchange%contactNode_current = 0
 
    do i=1,fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) then
        call clear_contact_state(fstrsolid%contacts(i));  cycle
      endif
 
      call scan_contact_state( fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
      & infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive, hecmesh )
 
      infoctchange%active = infoctchange%active .or. iactive
    enddo
    call hecmw_allreduce_l1(hecmesh, infoctchange%active, hecmw_lor)
 
    infoctchange%contactNode_current = infoctchange%contactNode_previous &
      + sum(infoctchange%n_statechange(:,kcatcont)) - sum(infoctchange%n_statechange(kcatcont,:))
    infoctchange%contactNode_previous = infoctchange%contactNode_current
 
    ! Output summary of contact state changes (always on; per-node detail requires CONTACT_LOG_LEVEL>=1)
    ! allreduce must be called by all ranks regardless of local counts
    s_chg   = infoctchange%n_statechange
    s_chg(kcatfree,kcatcont) = s_chg(kcatfree,kcatcont) + infoctchange%free2contact_new
    s_emov  = infoctchange%contact2neighbor
    s_ebey  = infoctchange%contact2beyond
    s_islid = infoctchange%contact2diffLpos
    s_act   = infoctchange%contactNode_current
    call hecmw_allreduce_i(hecmesh, s_chg, 9, hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_emov,  hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_ebey,  hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_islid, hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_act,   hecmw_sum)
    if (hecmw_comm_get_rank() == 0) then
      write(*,'(A,i7,A,i7,A,i7,A,i7,A,i7,A,i7)') ' Contact change: Free2Cont=', s_chg(kcatfree,kcatcont), &
        ' Cont2Free=', s_chg(kcatcont,kcatfree), &
        ' InnerSlid=', s_islid, ' MoveNeibr=', s_emov, ' MoveBynd=', s_ebey, ' ActiveNodes=', s_act
      if( sum(s_chg(kcatnear,:)) + sum(s_chg(:,kcatnear)) > 0 ) &
        write(*,'(A,i7,A,i7,A,i7,A,i7)') '   near state change: Free2Near=', s_chg(kcatfree,kcatnear), &
        ' Near2Free=', s_chg(kcatnear,kcatfree), &
        ' Near2Cont=', s_chg(kcatnear,kcatcont), ' Cont2Near=', s_chg(kcatcont,kcatnear)
    end if
 
    fstrsolid%ddunode = 0.d0
  end subroutine
 
  logical function fstr_is_contact_active()
    fstr_is_contact_active = active
  end function
 
  subroutine fstr_set_contact_active( a )
    logical, intent(in) :: a
    active = a
  end subroutine
 
  subroutine fstr_set_contact_active_restart( cstep, fstrSOLID, infoCTChange )
    integer(kind=kint), intent(in)               :: cstep
    type(fstr_solid), intent(in)                 :: fstrSOLID
    type(fstr_info_contactchange), intent(inout) :: infoCTChange
    integer(kind=kint) :: i, grpid
 
    active = infoctchange%contactNode_previous > 0
    do i = 1, fstrsolid%n_contacts
      if( active ) exit
      if( fstrsolid%contacts(i)%method /= contactn2s ) cycle
      if( .not. is_damping_enabled(fstrsolid%contacts(i)) ) cycle
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
      active = any( fstrsolid%contacts(i)%states(:)%state == contactnear )
    enddo
    infoctchange%active = active
  end subroutine
 
  logical function fstr_is_contact_conv(infoCTChange,hecMESH)
    type (fstr_info_contactchange), intent(in) :: infoctchange
    type (hecmwst_local_mesh), intent(in) :: hecmesh
 
    fstr_is_contact_conv = .false.
    ! The two SURF-SURF counters are incremented on the mortar path only, so the sum is
    ! unchanged for a NODE-SURF analysis.
    if( sum(infoctchange%n_statechange) + infoctchange%contact2neighbor + infoctchange%contact2beyond + &
        infoctchange%free2contact_new + infoctchange%contact2diffLpos == 0 ) &
      fstr_is_contact_conv = .true.
 
    call hecmw_allreduce_l1(hecmesh, fstr_is_contact_conv, hecmw_land)
  end function
 
  logical function fstr_is_matrixstructure_changed(infoCTChange)
    type (fstr_info_contactchange)   :: infoctchange
    ! Matrix sparsity changes when a node connection appears or disappears. Each pair reports
    ! through the counters of its sparsity_expansion mode, so a plain sum covers any mix.
    ! contact2difflpos is excluded (no structural change).
    fstr_is_matrixstructure_changed = .false.
    if( sum(infoctchange%n_statechange) + infoctchange%contact2neighbor + infoctchange%contact2beyond + &
        infoctchange%free2contact_new > 0 ) &
      fstr_is_matrixstructure_changed = .true.
  end function
 
  subroutine fstr_update_contact_multiplier( cstep, hecMESH, hecLagMAT, fstrSOLID, ctchanged )
    integer(kind=kint), intent(in)        :: cstep
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(hecmwst_matrix_lagrange), intent(in) :: hecLagMAT
    type(fstr_solid), intent(inout)        :: fstrSOLID
    logical, intent(out)                   :: ctchanged
 
    integer(kind=kint) :: i, nc, algtype, grpid
 
    gnt = 0.d0;  ctchanged = .false.
    nc = fstrsolid%n_contacts+fstrsolid%n_embeds
    do i=1, fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      if( fstrsolid%contacts(i)%method == contacts2s ) then
        call update_contact_multiplier_surfsurf( fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
          fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff )
      else
        algtype = fstrsolid%contacts(i)%algtype
        if( algtype == contactsslid .or. algtype == contactfslid ) then
          call update_contact_multiplier( fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
            fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff, &
            hecmesh, heclagmat, gnt, ctchanged )
        else if( algtype == contacttied ) then
          call update_tied_multiplier( fstrsolid%contacts(i), fstrsolid%unode(:), fstrsolid%dunode(:), &
          &  ctchanged )
        endif
      endif
    enddo
 
    do i=1, fstrsolid%n_embeds
      call update_tied_multiplier( fstrsolid%embeds(i), fstrsolid%unode(:), fstrsolid%dunode(:), &
      &  ctchanged )
    enddo
 
    call hecmw_allreduce_l1(hecmesh, ctchanged, hecmw_lor)
 
    if( nc>0 ) gnt = gnt/nc
  end subroutine
 
  subroutine fstr_begin_lambda_txn( fstrSOLID, cstep )
    type(fstr_solid), intent(inout) :: fstrSOLID
    integer(kind=kint), intent(in)  :: cstep
    integer(kind=kint) :: i, s, grpid
 
    do i = 1, fstrsolid%n_contacts
      if( fstrsolid%contacts(i)%method /= contacts2s ) cycle
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
      if( .not. associated(fstrsolid%contacts(i)%slave_surf) ) cycle
      do s = 1, size(fstrsolid%contacts(i)%slave_surf)
        fstrsolid%contacts(i)%slave_surf(s)%lam_work_n = 0
        fstrsolid%contacts(i)%slave_surf(s)%lam_work_id(:)  = 0
        fstrsolid%contacts(i)%slave_surf(s)%lam_work_val(:,:) = 0.d0
        ! the tangential working accumulator is cleared alongside lam_work_val; the
        ! stick/slip warm-start is read from begin at resolve time, not restored here
        fstrsolid%contacts(i)%slave_surf(s)%lam_work_t(:,:,:)  = 0.d0
        fstrsolid%contacts(i)%slave_surf(s)%lam_work_fstate(:,:) = contactstick
      enddo
    enddo
  end subroutine
 
  subroutine fstr_commit_lambda_txn( fstrSOLID, cstep )
    type(fstr_solid), intent(inout) :: fstrSOLID
    integer(kind=kint), intent(in)  :: cstep
    integer(kind=kint) :: i, s, grpid, r, m
    type(tcontactsurf), pointer :: surf
 
    do i = 1, fstrsolid%n_contacts
      if( fstrsolid%contacts(i)%method /= contacts2s ) cycle
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
      if( .not. associated(fstrsolid%contacts(i)%slave_surf) ) cycle
      do s = 1, size(fstrsolid%contacts(i)%slave_surf)
        surf => fstrsolid%contacts(i)%slave_surf(s)
        m = 0
        ! The master is kept if at least one slave node is active against it.
        do r = 1, surf%lam_work_n
          if( maxval(surf%lam_work_val(1:size(surf%nodes),r)) > 0.d0 ) then
            m = m + 1
            surf%lam_begin_id(m)  = surf%lam_work_id(r)
            surf%lam_begin_val(1:size(surf%nodes),m) = surf%lam_work_val(1:size(surf%nodes),r)
            ! the tangential multiplier rides along the master that passed the lambda_n filter
            surf%lam_begin_t(:,1:size(surf%nodes),m)    = surf%lam_work_t(:,1:size(surf%nodes),r)
            surf%lam_begin_fstate(1:size(surf%nodes),m) = surf%lam_work_fstate(1:size(surf%nodes),r)
          endif
        enddo
        surf%lam_begin_n = m
      enddo
    enddo
  end subroutine
 
  subroutine fstr_update_contact_tangentforce( cstep, fstrSOLID )
    integer(kind=kint), intent(in)        :: cstep
    type(fstr_solid), intent(inout)        :: fstrSOLID
 
    integer(kind=kint) :: i, grpid
 
    do i=1, fstrsolid%n_contacts
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      call update_contact_tangentforce( fstrsolid%contacts(i) )
    enddo
  end subroutine
 
  subroutine fstr_update_contact_state_vectors( fstrSOLID, dt )
    type(fstr_solid), intent(inout)        :: fstrsolid
    real(kind=kreal), intent(in)           :: dt  
 
    integer(kind=kint) :: i
 
    if( associated( fstrsolid%CONT_RELVEL ) ) fstrsolid%CONT_RELVEL(:) = 0.d0
    if( associated( fstrsolid%CONT_STATE ) ) fstrsolid%CONT_STATE(:) = 0.d0
 
    do i=1, fstrsolid%n_contacts
      call update_contact_state_vectors( fstrsolid%contacts(i), dt, fstrsolid%CONT_RELVEL, fstrsolid%CONT_STATE )
    enddo
  end subroutine
 
end module mcontact