fstr_mat_con_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
!-------------------------------------------------------------------------------
!         stiffness matrix structure for the contact analysis
!         employing standard Lagrange multiplier algorithm
 
module fstr_matrix_con_contact
 
  use m_fstr
  use elementinfo
  use m_fstr_contact_damping, only: is_damping_enabled
  use m_fstr_contact_elem_alag, only: get_unique_map
  use hecmw_matrix_ass, only: hecmw_mat_profile_has_node
  use m_fstr_finiterotationkinematics, only: fstr_uses_finite_rotation_kinematics
 
  implicit none
  private
  public :: fstr_get_num_lagrange_pernode
  public :: hecmwst_matrix_lagrange
  public :: fstr_save_originalmatrixstructure
  public :: fstr_mat_con_contact
  public :: fstr_s2s_profile_needs_refresh
  public :: fstr_is_matrixstruct_symmetric
  public :: fstr_is_contactalag_symmetric
  public :: fstr_is_tangent_symmetric
  public :: fstr_set_lagrange_diagonal
  public :: fstr_get_lagrange_diagonal
 
  integer(kind=kint), save         :: NPL_org, NPU_org
  type(nodeRelated), pointer, save :: list_nodeRelated_org(:) => null() 
 
  type(nodeRelated), pointer       :: list_nodeRelated(:) => null() 
 
  logical                          :: permission = .false.
 
contains
 
  integer(kind=kint) function fstr_get_num_lagrange_pernode(algtype)
    integer(kind=kint) :: algtype
    if( algtype == contactsslid .or. algtype == contactfslid ) then
      fstr_get_num_lagrange_pernode = 1
    else if( algtype == contacttied ) then
      fstr_get_num_lagrange_pernode = 3
    endif
  end function
 
  subroutine fstr_save_originalmatrixstructure(hecMAT)
 
    type(hecmwst_matrix) :: hecmat
 
    if( associated(list_noderelated_org) ) return
    call hecmw_construct_noderelated_from_hecmat(hecmat, npl_org, npu_org, list_noderelated_org)
 
  end subroutine fstr_save_originalmatrixstructure
 
  subroutine fstr_mat_con_contact(cstep,contact_algo,hecMAT,fstrSOLID,hecLagMAT,conMAT,is_contact_active_flag)
 
    integer(kind=kint)                   :: cstep
    integer(kind=kint)                   :: contact_algo
    type(hecmwst_matrix)                 :: hecmat
    type(fstr_solid)                     :: fstrsolid
    type(hecmwst_matrix_lagrange) :: heclagmat
 
    integer(kind=kint)                   :: num_lagrange
    integer(kind=kint)                   :: countnon0lu_node, countnon0lu_lagrange
    integer(kind=kint)                   :: numnon0_node, numnon0_lagrange
    type (hecmwst_matrix)                :: conmat
    logical, intent(in)                  :: is_contact_active_flag
 
    integer(kind=kint)                   :: i, j, grpid
    integer(kind=kint)                   :: count_n2s, count_s2s
    integer(kind=kint)                   :: nlag
 
    count_n2s = 0
    count_s2s = 0
    do i = 1, fstrsolid%n_contacts
      if( fstrsolid%contacts(i)%method == contactn2s ) count_n2s = count_n2s + 1
      if( fstrsolid%contacts(i)%method == contacts2s ) count_s2s = count_s2s + 1
    enddo
    count_n2s = count_n2s + fstrsolid%n_embeds
 
    num_lagrange = 0
    if( contact_algo == kcaslagrange ) then
      do i = 1, fstrsolid%n_contacts
        do j = 1, size(fstrsolid%contacts(i)%slave)
          fstrsolid%contacts(i)%states(j)%id_lagrange = 0
        enddo
      enddo
      do i = 1, fstrsolid%n_embeds
        do j = 1, size(fstrsolid%embeds(i)%slave)
          fstrsolid%embeds(i)%states(j)%id_lagrange = 0
        enddo
      enddo
 
      do i = 1, fstrsolid%n_contacts
        grpid = fstrsolid%contacts(i)%group
        if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
        nlag = fstr_get_num_lagrange_pernode(fstrsolid%contacts(i)%algtype)
        do j = 1, size(fstrsolid%contacts(i)%slave)
          if( .not. is_contact_active(fstrsolid%contacts(i)%states(j)%state) ) cycle
          fstrsolid%contacts(i)%states(j)%id_lagrange = num_lagrange + 1
          num_lagrange = num_lagrange + nlag
        enddo
      enddo
 
      do i = 1, fstrsolid%n_embeds
        grpid = fstrsolid%embeds(i)%group
        if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
        nlag = 3
        do j = 1, size(fstrsolid%embeds(i)%slave)
          if( .not. is_contact_active(fstrsolid%embeds(i)%states(j)%state) ) cycle
          fstrsolid%embeds(i)%states(j)%id_lagrange = num_lagrange + 1
          num_lagrange = num_lagrange + nlag
        enddo
      enddo
    endif
 
    ! Get original list of related nodes
    call hecmw_init_noderelated_from_org(hecmat%NP,num_lagrange,is_contact_active_flag,list_noderelated_org,list_noderelated)
 
    ! Construct new list of related nodes and Lagrange multipliers
    countnon0lu_node = npl_org + npu_org
    countnon0lu_lagrange = 0
    if( is_contact_active_flag )then
      if( count_n2s > 0 ) &
        call getnewlistofrelatednodesandlagrangemultipliers(cstep,contact_algo, &
        &  hecmat%NP,fstrsolid,countnon0lu_node,countnon0lu_lagrange,list_noderelated)
      if( count_s2s > 0 ) &
        call getnewlistofrelatednodesandlagrangemultipliers_ss(cstep,contact_algo, &
        &  hecmat%NP,fstrsolid,countnon0lu_node,countnon0lu_lagrange,list_noderelated)
    endif
 
    ! Construct new matrix structure(hecMAT&hecLagMAT)
    numnon0_node = countnon0lu_node/2
    numnon0_lagrange = countnon0lu_lagrange/2
    call hecmw_construct_hecmat_from_noderelated(hecmat%N, hecmat%NP, hecmat%NDOF, &
    & numnon0_node, num_lagrange, list_noderelated, hecmat)
    call hecmw_construct_hecmat_from_noderelated(hecmat%N, hecmat%NP, hecmat%NDOF, &
    & numnon0_node, num_lagrange, list_noderelated, conmat)
    if( contact_algo == kcaslagrange ) call hecmw_construct_heclagmat_from_noderelated(hecmat%NP, &
    & hecmat%NDOF, num_lagrange, numnon0_lagrange, is_contact_active_flag, list_noderelated, heclagmat)
    call hecmw_finalize_noderelated(list_noderelated)
 
    ! Copy Lagrange multipliers
    if( is_contact_active_flag .and. contact_algo == kcaslagrange ) &
      call fstr_copy_lagrange_contact(fstrsolid,heclagmat)
 
  end subroutine fstr_mat_con_contact
 
  subroutine getnewlistofrelatednodesandlagrangemultipliers( &
  & cstep, contact_algo, np, fstrSOLID, countNon0LU_node, countNon0LU_lagrange, list_nodeRelated )
    integer(kind=kint),intent(in)             :: cstep
    integer(kind=kint),intent(in)             :: contact_algo
    integer(kind=kint),intent(in)             :: np
    type(fstr_solid),intent(in)               :: fstrsolid
    integer(kind=kint), intent(inout)         :: countnon0lu_node, countnon0lu_lagrange
    type(noderelated), pointer, intent(inout) :: list_noderelated(:)
 
    integer(kind=kint)            :: grpid
    integer(kind=kint)            :: count_lagrange
    integer(kind=kint)            :: ctsurf, etype, nnode, ndlocal(l_max_surface_node + 1)
    integer(kind=kint)            :: i, j, k, nlag, algtype
    real(kind=kreal)              :: fcoeff 
    logical                       :: necessary_to_insert_node, necessary_to_insert_node_pair
    logical                       :: is_contact_active_flag, is_damping_active_flag
 
    count_lagrange = 0
    do i = 1, fstrsolid%n_contacts
      ! Process only NODE-SURF pairs here (mortar pairs are handled by getNewListOF..._ss).
      ! For a mortar pair states(j)%surface is -1, so master(ctsurf)%etype below would be an
      ! out-of-bounds read on a mixed deck.
      if( fstrsolid%contacts(i)%method /= contactn2s ) cycle
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      fcoeff = fstrsolid%contacts(i)%fcoeff
      necessary_to_insert_node = ( fcoeff /= 0.0d0 .or. contact_algo == kcaalagrange )
 
      algtype = fstrsolid%contacts(i)%algtype
      nlag = fstr_get_num_lagrange_pernode(algtype)
      if( contact_algo == kcaalagrange ) nlag = 1
      if( algtype == contacttied ) permission = .true.
 
      do j = 1, size(fstrsolid%contacts(i)%slave)
        ! stick or sliding contact is active
        is_contact_active_flag = is_contact_active(fstrsolid%contacts(i)%states(j)%state)
        ! damping is active
        is_damping_active_flag = fstrsolid%contacts(i)%states(j)%state == contactnear .and. &
          &  is_damping_enabled(fstrsolid%contacts(i))
        
        if( is_contact_active_flag .or. is_damping_active_flag ) then
 
          ctsurf = fstrsolid%contacts(i)%states(j)%surface
          etype = fstrsolid%contacts(i)%master(ctsurf)%etype
          if( etype/=fe_tri3n .and. etype/=fe_quad4n ) &
            stop " ##Error: This element type is not supported in contact analysis !!! "
          nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
          ndlocal(1) = fstrsolid%contacts(i)%slave(j)
          ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
 
          ! For CONTACTNEAR damping (especially S-Lagrange + frictionless),
          ! we still need slave-master connectivity to assemble damping terms.
          necessary_to_insert_node_pair = necessary_to_insert_node .or. is_damping_active_flag
 
          if( is_contact_active_flag ) then
            do k=1,nlag
              if( contact_algo == kcaslagrange ) &
                count_lagrange = fstrsolid%contacts(i)%states(j)%id_lagrange + k - 1
              call hecmw_ass_noderelated_from_contact_pair(np, nnode, ndlocal, count_lagrange, permission, &
              & necessary_to_insert_node_pair, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
            enddo
          else
            ! NEAR damping only: no Lagrange multiplier, insert connectivity once
            call hecmw_ass_noderelated_from_contact_pair(np, nnode, ndlocal, 0, permission, &
            & necessary_to_insert_node_pair, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
          endif
              
        end if
 
      enddo
 
    enddo
 
    do i = 1, fstrsolid%n_embeds
 
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
 
      necessary_to_insert_node = ( contact_algo == kcaalagrange )
 
      nlag = 3
      if( contact_algo == kcaalagrange ) nlag = 1
      permission = .true.
 
      do j = 1, size(fstrsolid%embeds(i)%slave)
 
        if( .not. is_contact_active(fstrsolid%embeds(i)%states(j)%state) ) cycle
        ctsurf = fstrsolid%embeds(i)%states(j)%surface
        etype = fstrsolid%embeds(i)%master(ctsurf)%etype
        nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
        ndlocal(1) = fstrsolid%embeds(i)%slave(j)
        ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
 
        do k=1,nlag
          if( contact_algo == kcaslagrange ) &
            count_lagrange = fstrsolid%embeds(i)%states(j)%id_lagrange + k - 1
          call hecmw_ass_noderelated_from_contact_pair(np, nnode, ndlocal, count_lagrange, permission, &
          & necessary_to_insert_node, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
        enddo
      enddo
 
    enddo
 
  end subroutine getnewlistofrelatednodesandlagrangemultipliers
 
  subroutine fstr_copy_lagrange_contact(fstrSOLID,hecLagMAT)
 
    type(fstr_solid)                        :: fstrsolid
    type(hecmwst_matrix_lagrange)          :: heclagmat
    integer (kind=kint)                    :: id_lagrange, algtype, i, j, k, nlag
 
    do i = 1, fstrsolid%n_contacts
 
      algtype = fstrsolid%contacts(i)%algtype
      nlag = fstr_get_num_lagrange_pernode(algtype)
 
      do j = 1, size(fstrsolid%contacts(i)%slave)
        id_lagrange = fstrsolid%contacts(i)%states(j)%id_lagrange
        if( id_lagrange == 0 ) cycle
        do k=1,nlag
          heclagmat%Lagrange(id_lagrange+k-1)=fstrsolid%contacts(i)%states(j)%multiplier(k)
        enddo
      enddo
    enddo
 
    do i = 1, fstrsolid%n_embeds
      nlag = 3
      do j = 1, size(fstrsolid%embeds(i)%slave)
        id_lagrange = fstrsolid%embeds(i)%states(j)%id_lagrange
        if( id_lagrange == 0 ) cycle
        do k=1,nlag
          heclagmat%Lagrange(id_lagrange+k-1)=fstrsolid%embeds(i)%states(j)%multiplier(k)
        enddo
      enddo
    enddo
 
  end subroutine fstr_copy_lagrange_contact
 
  logical function fstr_is_matrixstruct_symmetric(fstrSOLID,hecMESH)
 
    type(fstr_solid )        :: fstrsolid
    type(hecmwst_local_mesh) :: hecmesh
    integer (kind=kint)      :: is_in_contact
 
    is_in_contact = 0
    if( fstrsolid%n_contacts>0 ) then
      if( any(fstrsolid%contacts(:)%fcoeff /= 0.0d0) )  is_in_contact = 1
    endif
    call hecmw_allreduce_i1(hecmesh, is_in_contact, hecmw_max)
    if( is_in_contact == 0 .and. hecmesh%n_dof /= 4 .and. fstr_is_tangent_symmetric(fstrsolid,hecmesh) ) then
      fstr_is_matrixstruct_symmetric = .true.
    else
      fstr_is_matrixstruct_symmetric = .false.
    endif
 
  end function fstr_is_matrixstruct_symmetric
 
  logical function fstr_is_contactalag_symmetric(fstrSOLID,hecMESH)
 
    type(fstr_solid )        :: fstrsolid
    type(hecmwst_local_mesh) :: hecmesh
    integer (kind=kint)      :: is_unsymmetric
 
    ! the ALagrange contact terms are symmetric while the friction cone radius stays frozen
    ! at the multiplier; !CONTACT_ALGO, FRICTION_CONE=FOLLOW clears contact%symmetric and
    ! the tangent then carries the coupling block of getContactStiffness_Alag
    is_unsymmetric = 0
    if( fstrsolid%n_contacts>0 ) then
      if( any( fstrsolid%contacts(:)%fcoeff /= 0.0d0 .and. .not.fstrsolid%contacts(:)%symmetric ) ) is_unsymmetric = 1
    endif
    call hecmw_allreduce_i1(hecmesh, is_unsymmetric, hecmw_max)
    fstr_is_contactalag_symmetric = ( is_unsymmetric == 0 ) .and. fstr_is_tangent_symmetric(fstrsolid,hecmesh)
 
  end function fstr_is_contactalag_symmetric
 
  logical function fstr_is_tangent_symmetric(fstrSOLID,hecMESH)
 
    type(fstr_solid )        :: fstrsolid
    type(hecmwst_local_mesh) :: hecmesh
    integer (kind=kint)      :: is_unsymmetric, i, ytype, nn
 
    ! non-associated flow (dilatancy angle psi /= friction angle phi) makes the consistent
    ! tangent unsymmetric; when psi is omitted the parser stores psi=phi, so the exact
    ! comparison below classifies that case as associated (symmetric)
    is_unsymmetric = 0
    if( associated(fstrsolid%materials) ) then
      do i = 1, size(fstrsolid%materials)
        ytype = getyieldfunction( fstrsolid%materials(i)%mtype )
        if( ytype == 1 ) then      ! Mohr-Coulomb: PLCONST3=phi, PLCONST4=psi
          if( fstrsolid%materials(i)%variables(m_plconst3) /= fstrsolid%materials(i)%variables(m_plconst4) ) &
            is_unsymmetric = 1
        elseif( ytype == 2 ) then  ! Drucker-Prager: PLCONST3=eta(phi), PLCONST5=etabar(psi)
          if( fstrsolid%materials(i)%variables(m_plconst3) /= fstrsolid%materials(i)%variables(m_plconst5) ) &
            is_unsymmetric = 1
        endif
      enddo
    endif
    ! Finite-rotation UL shells have an unsymmetric tangent even for elastic materials.
    if( associated(fstrsolid%elements) ) then
      do i = 1, hecmesh%n_elem
        if( .not. associated(fstrsolid%elements(i)%gausses) ) cycle
        if( fstrsolid%elements(i)%gausses(1)%pMaterial%nlgeom_flag /= updatelag ) cycle
        nn = hecmesh%elem_node_index(i) - hecmesh%elem_node_index(i-1)
        if( fstr_uses_finite_rotation_kinematics(hecmesh%elem_type(i), nn, &
            fstrsolid%elements(i)%gausses(1)%pMaterial) ) then
          is_unsymmetric = 1
          exit
        endif
      enddo
    endif
    call hecmw_allreduce_i1(hecmesh, is_unsymmetric, hecmw_max)
    fstr_is_tangent_symmetric = (is_unsymmetric == 0)
 
  end function fstr_is_tangent_symmetric
 
  subroutine fstr_set_lagrange_diagonal(hecLagMAT, ilag, value)
    type(hecmwst_matrix_lagrange), intent(inout) :: heclagmat
    integer(kind=kint), intent(in) :: ilag
    real(kind=kreal), intent(in) :: value 
 
    if (ilag < 1 .or. ilag > heclagmat%num_lagrange) then
      write(*,*) 'Error in fstr_set_lagrange_diagonal: invalid Lagrange multiplier index', ilag
      stop
    endif
 
    if (.not. associated(heclagmat%D_lagrange)) then
      write(*,*) 'Error in fstr_set_lagrange_diagonal: D_lagrange not allocated'
      stop
    endif
 
    heclagmat%D_lagrange(ilag) = value
 
  end subroutine fstr_set_lagrange_diagonal
 
  real(kind=kreal) function fstr_get_lagrange_diagonal(hecLagMAT, ilag)
    type(hecmwst_matrix_lagrange), intent(in) :: heclagmat
    integer(kind=kint), intent(in) :: ilag
 
    if (ilag < 1 .or. ilag > heclagmat%num_lagrange) then
      write(*,*) 'Error in fstr_get_lagrange_diagonal: invalid Lagrange multiplier index', ilag
      fstr_get_lagrange_diagonal = 0.0d0
      return
    endif
 
    if (.not. associated(heclagmat%D_lagrange)) then
      write(*,*) 'Error in fstr_get_lagrange_diagonal: D_lagrange not allocated'
      fstr_get_lagrange_diagonal = 0.0d0
      return
    endif
 
    fstr_get_lagrange_diagonal = heclagmat%D_lagrange(ilag)
 
  end function fstr_get_lagrange_diagonal
 
 
  subroutine getnewlistofrelatednodesandlagrangemultipliers_ss( &
      & cstep, contact_algo, np, fstrSOLID, countNon0LU_node, countNon0LU_lagrange, list_nodeRelated )
    integer(kind=kint),intent(in)             :: cstep
    integer(kind=kint),intent(in)             :: contact_algo
    integer(kind=kint),intent(in)             :: np
    type(fstr_solid),intent(in)               :: fstrsolid
    integer(kind=kint), intent(inout)         :: countnon0lu_node, countnon0lu_lagrange
    type(noderelated), pointer, intent(inout) :: list_noderelated(:)
    integer(kind=kint)            :: grpid
    integer(kind=kint)            :: ctsurf, nsurf
    integer(kind=kint)            :: i, j, m
    integer(kind=kint)            :: g, unique_count
    integer(kind=kint)            :: maplist(max_n_intp), master_idxs(max_n_intp)
    real(kind=kreal)              :: fcoeff 
    logical                       :: necessary_to_insert_node
 
    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
 
      fcoeff = fstrsolid%contacts(i)%fcoeff
      necessary_to_insert_node = ( fcoeff /= 0.0d0 .or. contact_algo == kcaalagrange )
 
      do j = 1, size(fstrsolid%contacts(i)%slave_surf)
        if( fstrsolid%contacts(i)%slave_surf(j)%state == contactfree ) cycle
 
        if( contact_algo == kcaalagrange ) then
          ! ALagrange SS: register sparsity once per unique (slave_surf, master) pair.
          ! sparsity_expansion == SPARSITY_NONE:     current master element only (matrix rebuilt on contact2neighbor)
          ! sparsity_expansion == SPARSITY_NEIGHBOR: current + neighbor master elements (matrix rebuilt only on contact2beyond)
          call get_unique_map( fstrsolid%contacts(i)%slave_surf(j), maplist, master_idxs, unique_count )
          do g = 1, unique_count
            ctsurf = master_idxs(g)
            ! Register current master element
            call register_pair_to_sparsity( np, fstrsolid%contacts(i)%slave_surf(j)%nodes, &
              fstrsolid%contacts(i)%master(ctsurf), necessary_to_insert_node, &
              countnon0lu_node, countnon0lu_lagrange, list_noderelated )
 
            ! SPARSITY_NEIGHBOR: also register neighbor master elements
            if( fstrsolid%contacts(i)%sparsity_expansion == sparsity_neighbor ) then
              do m = 1, fstrsolid%contacts(i)%master(ctsurf)%n_neighbor
                nsurf = fstrsolid%contacts(i)%master(ctsurf)%neighbor(m)
                call register_pair_to_sparsity( np, fstrsolid%contacts(i)%slave_surf(j)%nodes, &
                  fstrsolid%contacts(i)%master(nsurf), necessary_to_insert_node, &
                  countnon0lu_node, countnon0lu_lagrange, list_noderelated )
              enddo
            endif
          enddo
        endif
      enddo
    enddo
 
  end subroutine getnewlistofrelatednodesandlagrangemultipliers_ss
 
  subroutine register_pair_to_sparsity( np, slave_nodes, master_surf, &
      & necessary_to_insert_node, countNon0LU_node, countNon0LU_lagrange, list_nodeRelated )
    integer(kind=kint), intent(in)            :: np
    integer(kind=kint), intent(in)            :: slave_nodes(:)
    type(tsurfelement), intent(in)            :: master_surf
    logical, intent(in)                       :: necessary_to_insert_node
    integer(kind=kint), intent(inout)         :: countnon0lu_node, countnon0lu_lagrange
    type(noderelated), pointer, intent(inout) :: list_noderelated(:)
 
    integer(kind=kint) :: nnode_s, nnode_m, l, m, idx, nnode_pair, etype
    integer(kind=kint) :: ndlocal(2*l_max_surface_node + 1)
 
    etype = master_surf%etype
    if( etype/=fe_tri3n .and. etype/=fe_quad4n ) stop " ##Error: This element type is not supported in contact analysis !!! "
    nnode_s = size(slave_nodes)
    nnode_m = size(master_surf%nodes)
    do l = 1, nnode_s
      ndlocal(1) = slave_nodes(l)
      ! ndLocal: slave_l + other_slave_nodes + master_nodes
      idx = 1
      do m = 1, nnode_s
        if( m == l ) cycle
        idx = idx + 1
        ndlocal(idx) = slave_nodes(m)
      enddo
      ndlocal(idx+1:idx+nnode_m) = master_surf%nodes(:)
      nnode_pair = nnode_s - 1 + nnode_m
      ! mortar pairs have no Lagrange row: 0 keeps insert_lagrange out of the reservation
      call hecmw_ass_noderelated_from_contact_pair( np, nnode_pair, ndlocal, 0, permission, &
        & necessary_to_insert_node, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
    enddo
  end subroutine register_pair_to_sparsity
 
  logical function fstr_s2s_profile_needs_refresh( cstep, contact_algo, fstrSOLID, conMAT )
    integer(kind=kint), intent(in)   :: cstep
    integer(kind=kint), intent(in)   :: contact_algo
    type(fstr_solid), intent(in)     :: fstrsolid
    type(hecmwst_matrix), intent(in) :: conmat
 
    integer(kind=kint) :: i, j, g, ctsurf, grpid, unique_count
    integer(kind=kint) :: maplist(max_n_intp), master_idxs(max_n_intp)
 
    fstr_s2s_profile_needs_refresh = .false.
    if( contact_algo /= kcaalagrange ) return
 
    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
 
      do j = 1, size(fstrsolid%contacts(i)%slave_surf)
        if( fstrsolid%contacts(i)%slave_surf(j)%state == contactfree ) cycle
 
        call get_unique_map( fstrsolid%contacts(i)%slave_surf(j), maplist, master_idxs, unique_count )
        do g = 1, unique_count
          ctsurf = master_idxs(g)
          if( .not. pair_cross_profile_complete( fstrsolid%contacts(i)%slave_surf(j)%nodes, &
              &  fstrsolid%contacts(i)%master(ctsurf), conmat ) ) then
            fstr_s2s_profile_needs_refresh = .true.
            return
          endif
        enddo
      enddo
    enddo
  end function fstr_s2s_profile_needs_refresh
 
  logical function pair_cross_profile_complete( slave_nodes, master_surf, conMAT )
    integer(kind=kint), intent(in)   :: slave_nodes(:)
    type(tsurfelement), intent(in)   :: master_surf
    type(hecmwst_matrix), intent(in) :: conmat
 
    integer(kind=kint) :: s, m
 
    pair_cross_profile_complete = .true.
    do s = 1, size(slave_nodes)
      do m = 1, size(master_surf%nodes)
        if( .not. hecmw_mat_profile_has_node( conmat, slave_nodes(s), master_surf%nodes(m) ) ) then
          pair_cross_profile_complete = .false.
          return
        endif
      enddo
    enddo
  end function pair_cross_profile_complete
 
end module fstr_matrix_con_contact