fstr_contact_search.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_contact_search
  use m_fstr
  use hecmw
  use elementinfo
  use mcontactdef
  use m_fstr_contact_geom
  use m_fstr_contact_element
  use m_fstr_contact_interference
  use m_fstr_contact_smoothing
  implicit none
 
  integer(kind=kint), parameter :: contact_log_level = 0
 
contains
 
  subroutine track_contact_position( nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID, flag_ctAlgo )
    integer, intent(in)                              :: nslave
    type( tcontact ), intent(inout)                  :: contact
    type( fstr_info_contactchange ), intent(inout)   :: infoCTChange
    real(kind=kreal), intent(in)                     :: currpos(:)   
    real(kind=kreal), intent(in)                     :: currdisp(:)  
    integer(kind=kint), intent(in)                   :: nodeID(:)
    integer(kind=kint), intent(in)                   :: elemID(:)
    character(len=9), intent(in), optional           :: flag_ctAlgo
 
    integer(kind=kint) :: slave, sid0, sid, etype
    integer(kind=kint) :: nn, i, j, iSS
    real(kind=kreal)    :: coord(3), elem(3, l_max_elem_node), elem0(3, l_max_elem_node)
    logical            :: isin
    real(kind=kreal)    :: opos(2)
    integer(kind=kint) :: bktID, nCand, idm
    integer(kind=kint), allocatable :: indexCand(:)
    logical            :: is_implicit
 
    is_implicit = present(flag_ctalgo)
 
    sid = 0
 
    slave = contact%slave(nslave)
    coord(:) = currpos(3*slave-2:3*slave)
    sid0 = contact%states(nslave)%surface
    opos = contact%states(nslave)%lpos(1:2)
    etype = contact%master(sid0)%etype
    nn = getnumberofnodes( etype )
    do j=1,nn
      iss = contact%master(sid0)%nodes(j)
      elem0(1:3,j)=currpos(3*iss-2:3*iss)-currdisp(3*iss-2:3*iss)
    enddo
    call project_point2surfelement( coord, contact%master(sid0), currpos, &
      contact%states(nslave), isin, contact%cparam%DISTCLR_NOCHECK, &
      contact%states(nslave)%lpos(1:2), contact%cparam%CLR_SAME_ELEM, smoothing=contact%smoothing )
    if( .not. isin ) then
      do i=1, contact%master(sid0)%n_neighbor
        sid = contact%master(sid0)%neighbor(i)
        call project_point2surfelement( coord, contact%master(sid), currpos, &
          contact%states(nslave), isin, contact%cparam%DISTCLR_NOCHECK, &
          localclr=contact%cparam%CLEARANCE, smoothing=contact%smoothing )
        if( isin ) then
          contact%states(nslave)%surface = sid
          exit
        endif
      enddo
    endif
 
    if( .not. isin ) then   ! such case is considered to rarely or never occur
      write(*,*) 'Warning: contact moved beyond neighbor elements'
      ! get master candidates from bucketDB
      bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
      ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
      if (ncand > 0) then
        allocate(indexcand(ncand))
        call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
        do idm= 1, ncand
          sid = indexcand(idm)
          if( sid==sid0 ) cycle
          if( associated(contact%master(sid0)%neighbor) ) then
            if( any(sid==contact%master(sid0)%neighbor(:)) ) cycle
          endif
          if (.not. is_in_surface_box( contact%master(sid), coord(1:3), contact%cparam%BOX_EXP_RATE )) cycle
          call project_point2surfelement( coord, contact%master(sid), currpos, &
            contact%states(nslave), isin, contact%cparam%DISTCLR_NOCHECK, &
            localclr=contact%cparam%CLEARANCE, smoothing=contact%smoothing )
          if( isin ) then
            contact%states(nslave)%surface = sid
            exit
          endif
        enddo
        deallocate(indexcand)
      endif
    endif
 
    if( isin ) then
      if( contact%states(nslave)%surface==sid0 ) then
        if(any(dabs(contact%states(nslave)%lpos(1:2)-opos(:)) >= contact%cparam%CLR_DIFFLPOS))  then
          !$omp atomic
          infoctchange%contact2difflpos = infoctchange%contact2difflpos + 1
        endif
      else
        if (contact_log_level >= 1) write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3)') "Node",nodeid(slave)," move to contact with", &
          elemid(contact%master(sid)%eid), " with distance ",      &
          contact%states(nslave)%distance," at ",contact%states(nslave)%lpos(1:2)
        !$omp atomic
        infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
      endif
      if( is_implicit .and. flag_ctalgo=='SLagrange' ) then
        ! Setup elem array for update_TangentForce
        etype = contact%master(contact%states(nslave)%surface)%etype
        nn = size(contact%master(contact%states(nslave)%surface)%nodes)
        do j=1,nn
          iss = contact%master(contact%states(nslave)%surface)%nodes(j)
          elem(1:3,j)=currpos(3*iss-2:3*iss)
        enddo
        call update_tangentforce(etype,nn,elem0,elem,contact%states(nslave))
      endif
      iss = isinsideelement( etype, contact%states(nslave)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
      if( iss>0 ) &
        call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
        contact%states(nslave)%lpos(1:2), contact%states(nslave)%direction(:) )
    else if( .not. isin ) then
      if (contact_log_level >= 1) write(*,'(A,i10,A)') "Node",nodeid(slave)," move out of contact"
      contact%states(nslave)%state = contactfree
      contact%states(nslave)%multiplier(:) = 0.d0
    endif
 
  end subroutine track_contact_position
 
  subroutine scan_contact_state( contact, currpos, currdisp, infoCTChange, &
    nodeID, elemID, is_init, active, hecMESH, flag_ctAlgo, ndforce )
    type( tcontact ), intent(inout)                  :: contact
    type( fstr_info_contactchange ), intent(inout)   :: infoCTChange
    real(kind=kreal), intent(in)                     :: currpos(:)   
    real(kind=kreal), intent(in)                     :: currdisp(:)  
    integer(kind=kint), intent(in)                   :: nodeID(:)
    integer(kind=kint), intent(in)                   :: elemID(:)
    logical, intent(in)                              :: is_init
    logical, intent(out)                             :: active
    type( hecmwst_local_mesh ), intent(in), optional :: hecMESH
    character(len=9), intent(in), optional           :: flag_ctAlgo
    real(kind=kreal), intent(in), optional           :: ndforce(:)   
 
    real(kind=kreal)    :: distclr
    integer(kind=kint)  :: slave, id, etype
    integer(kind=kint)  :: i, iSS, nactive
    real(kind=kreal)    :: coord(3)
    real(kind=kreal)    :: nlforce
    logical             :: isin
    integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
    real(kind=kreal)    :: effective_near_dist, distclr_use
    !
    integer, pointer :: indexMaster(:),indexCand(:)
    integer   ::  nMaster,idm,nMasterMax,bktID,nCand
    logical :: is_implicit
 
    is_implicit = present(flag_ctalgo)
 
    if( is_init ) then
      distclr = contact%cparam%DISTCLR_INIT
    else
      distclr = contact%cparam%DISTCLR_FREE
      if( contact%algtype == contacttied ) then
        active = .false.
        do i= 1, size(contact%slave)
          if( contact%states(i)%state==contactstick ) then
            active = .true.
            exit
          endif
        enddo
      endif
    endif
 
    allocate(contact_surf(size(nodeid)))
    allocate(states_prev(size(contact%slave)))
    contact_surf(:) = huge(0)
    do i = 1, size(contact%slave)
      states_prev(i) = contact%states(i)%state
    enddo
 
    if( contact%smoothing == kcsnagata ) then
      call update_surface_normal( contact%master, currpos, hecmesh )
      call create_intermediate_points( contact%master, currpos, contact%pair_name )
    endif
 
    call update_surface_box_info( contact%master, currpos )
    call update_surface_bucket_info( contact%master, contact%master_bktDB )
 
    ! Compute effective near distance for CONTACTNEAR detection
    ! If damping is configured (damp_gact > 0), use it as the near distance
    effective_near_dist = contact%cparam%NEAR_DIST
    if (contact%damp_gact > 0.0d0) then
      effective_near_dist = max(effective_near_dist, contact%damp_gact)
    endif
 
    !$omp parallel do &
    !$omp& default(none) &
    !$omp& private(i,slave,id,nlforce,coord,indexMaster,nMaster,iSS,idm,etype,isin, &
    !$omp&         bktID,nCand,indexCand,distclr_use) &
    !$omp& firstprivate(nMasterMax,is_implicit,effective_near_dist) &
    !$omp& shared(contact,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,nodeID,elemID,distclr,contact_surf,is_init) &
    !$omp& reduction(.or.:active) &
    !$omp& schedule(dynamic,1)
    do i= 1, size(contact%slave)
      if(contact%if_type /= 0) call set_shrink_factor(contact%ctime, contact%states(i), contact%if_etime, contact%if_type)
      slave = contact%slave(i)
      if( contact%states(i)%state==contactstick .or. contact%states(i)%state==contactslip ) then
        id = contact%states(i)%surface
        nlforce = contact%states(i)%multiplier(1)
 
        ! update direction of TIED contact
        if( contact%algtype == contacttied ) then
          call update_direction( i, contact, currpos )
          if (.not.is_init) cycle
        endif
 
        if( nlforce < contact%cparam%TENSILE_FORCE ) then
          contact%states(i)%state = contactfree
          contact%states(i)%multiplier(:) = 0.d0
          if (contact_log_level >= 1) write(*,'(A,i10,A,i10,A,e12.3)') "Node",nodeid(slave)," free from contact with element", &
            elemid(contact%master(id)%eid), " with tensile force ", nlforce
          cycle
        endif
        if( contact%algtype /= contactfslid ) then   ! small slide problem
          contact_surf(contact%slave(i)) = -elemid(contact%master(id)%eid)
        else
          if( is_implicit ) then
            call track_contact_position( i, contact, currpos, currdisp, infoctchange, nodeid, elemid, &
              flag_ctalgo=flag_ctalgo )
          else
            call track_contact_position( i, contact, currpos, currdisp, infoctchange, nodeid, elemid )
          endif
          if( contact%states(i)%state /= contactfree ) then
            id = contact%states(i)%surface
            contact_surf(contact%slave(i)) = -elemid(contact%master(id)%eid)
          endif
        endif
 
      else if( contact%states(i)%state==contactnear ) then
        ! NEAR tracking: re-project on current master surface
        coord(:) = currpos(3*slave-2:3*slave)
        id = contact%states(i)%surface
        if (effective_near_dist > 0.0d0) then
          distclr_use = max(distclr, effective_near_dist / contact%master(id)%reflen)
        else
          distclr_use = distclr
        end if
        call project_point2surfelement( coord, contact%master(id), currpos, &
          contact%states(i), isin, distclr_use, &
          contact%states(i)%lpos(1:2), contact%cparam%CLR_SAME_ELEM, smoothing=contact%smoothing )
        if (isin) then
          if (contact%states(i)%distance <= distclr * contact%master(id)%reflen) then
            ! Within contact threshold -> upgrade to STICK
            contact%states(i)%state = contactstick
            contact%states(i)%multiplier(:) = 0.d0
            if (contact_log_level >= 1) write(*,'(A,i10,A,i10,A,i6)') "Node",nodeid(slave)," upgraded NEAR->STICK on element", &
              elemid(contact%master(id)%eid)," rank=",hecmw_comm_get_rank()
          else if (contact%states(i)%distance > effective_near_dist) then
            ! Beyond NEAR range -> free
            contact%states(i)%state = contactfree
            contact%states(i)%multiplier(:) = 0.d0
          end if
          ! else: stay NEAR
          if (.not. is_contact_free(contact%states(i)%state)) then
            etype = contact%master(id)%etype
            iss = isinsideelement( etype, contact%states(i)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
            if( iss>0 ) &
              call cal_node_normal( id, iss, contact%master, currpos, &
              contact%states(i)%lpos(1:2), contact%states(i)%direction(:) )
            contact_surf(contact%slave(i)) = elemid(contact%master(id)%eid)
          end if
        else
          ! Lost projection -> free
          contact%states(i)%state = contactfree
          contact%states(i)%multiplier(:) = 0.d0
        end if
 
      else if( contact%states(i)%state==contactfree ) then
        if( contact%algtype == contacttied .and. .not. is_init ) cycle
        coord(:) = currpos(3*slave-2:3*slave)
 
        ! get master candidates from bucketDB
        bktid = bucketdb_getbucketid(contact%master_bktDB, coord)
        ncand = bucketdb_getnumcand(contact%master_bktDB, bktid)
        if (ncand == 0) cycle
        allocate(indexcand(ncand))
        call bucketdb_getcand(contact%master_bktDB, bktid, ncand, indexcand)
 
        nmastermax = ncand
        allocate(indexmaster(nmastermax))
        nmaster = 0
 
        ! narrow down candidates
        do idm= 1, ncand
          id = indexcand(idm)
          if (.not. is_in_surface_box( contact%master(id), coord(1:3), contact%cparam%BOX_EXP_RATE )) cycle
          nmaster = nmaster + 1
          indexmaster(nmaster) = id
        enddo
        deallocate(indexcand)
 
        if(nmaster == 0) then
          deallocate(indexmaster)
          cycle
        endif
 
        do idm = 1,nmaster
          id = indexmaster(idm)
          ! Expand distclr for NEAR detection
          if (effective_near_dist > 0.0d0) then
            distclr_use = max(distclr, effective_near_dist / contact%master(id)%reflen)
          else
            distclr_use = distclr
          end if
          call project_point2surfelement( coord, contact%master(id), currpos, &
            contact%states(i), isin, distclr_use, localclr=contact%cparam%CLEARANCE, smoothing=contact%smoothing )
          if( .not. isin ) cycle
          ! Classify: STICK or NEAR
          if (effective_near_dist > 0.0d0 .and. &
              contact%states(i)%distance > distclr * contact%master(id)%reflen) then
            contact%states(i)%state = contactnear
          end if
          contact%states(i)%surface = id
          contact%states(i)%multiplier(:) = 0.d0
          etype = contact%master(id)%etype
          iss = isinsideelement( etype, contact%states(i)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
          if( iss>0 ) &
            call cal_node_normal( id, iss, contact%master, currpos, contact%states(i)%lpos(1:2), &
            contact%states(i)%direction(:) )
          contact_surf(contact%slave(i)) = elemid(contact%master(id)%eid)
          if (contact_log_level >= 1) then
            if (contact%states(i)%state == contactnear) then
              write(*,'(A,i10,A,i10,A,f7.3,A,i6)') "Node",nodeid(slave)," near element", &
                elemid(contact%master(id)%eid), &
                " with distance ", contact%states(i)%distance," rank=",hecmw_comm_get_rank()
            else
              write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3,A,3f7.3,A,i6)') "Node",nodeid(slave)," contact with element", &
                elemid(contact%master(id)%eid),       &
                " with distance ", contact%states(i)%distance," at ",contact%states(i)%lpos(1:2), &
                " along direction ", contact%states(i)%direction," rank=",hecmw_comm_get_rank()
            end if
          end if
          exit
        enddo
        deallocate(indexmaster)
      endif
    enddo
    !$omp end parallel do
 
    if( contact%algtype == contacttied .and. .not. is_init ) then
      deallocate(contact_surf)
      deallocate(states_prev)
      return
    endif
 
    call hecmw_contact_comm_allreduce_i(contact%comm, contact_surf, hecmw_min)
    nactive = 0
    do i = 1, size(contact%slave)
      if (.not. is_contact_free(contact%states(i)%state)) then             ! any slave in contact or near
        id = contact%states(i)%surface
        if (abs(contact_surf(contact%slave(i))) /= elemid(contact%master(id)%eid)) then ! that is in contact with other surface
          contact%states(i)%state = contactfree                           ! should be freed
          if (contact_log_level >= 1) &
          &  write(*,'(A,i10,A,i10,A,i6,A,i6,A)') "Node",nodeid(contact%slave(i)), &
          &  " contact with element",elemid(contact%master(id)%eid), &
          &  " in rank",hecmw_comm_get_rank()," freed due to duplication"
        else if (is_contact_active(contact%states(i)%state)) then
          nactive = nactive + 1
        endif
      endif
      if (is_contact_free(states_prev(i)) .and. .not. is_contact_free(contact%states(i)%state)) then
        infoctchange%free2contact = infoctchange%free2contact + 1
      elseif (.not. is_contact_free(states_prev(i)) .and. is_contact_free(contact%states(i)%state)) then
        infoctchange%contact2free = infoctchange%contact2free + 1
      endif
    enddo
    active = (nactive > 0)
    deallocate(contact_surf)
    deallocate(states_prev)
  end subroutine scan_contact_state
 
  subroutine scan_embed_state( flag_ctAlgo, embed, currpos, currdisp, ndforce, infoCTChange, &
    nodeID, elemID, is_init, active, B )
    character(len=9), intent(in)                     :: flag_ctAlgo
    type( tcontact ), intent(inout)                  :: embed
    type( fstr_info_contactchange ), intent(inout)   :: infoCTChange
    real(kind=kreal), intent(in)                     :: currpos(:)   
    real(kind=kreal), intent(in)                     :: currdisp(:)  
    real(kind=kreal), intent(in)                     :: ndforce(:)   
    integer(kind=kint), intent(in)                   :: nodeID(:)
    integer(kind=kint), intent(in)                   :: elemID(:)
    logical, intent(in)                              :: is_init
    logical, intent(out)                             :: active
    real(kind=kreal), optional, target               :: b(:)         
 
    real(kind=kreal)    :: distclr
    integer(kind=kint)  :: slave, id, etype
    integer(kind=kint)  :: nn, i, j, iSS, nactive
    real(kind=kreal)    :: coord(3), elem(3, l_max_elem_node )
    logical             :: isin
    integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
    !
    integer, pointer :: indexMaster(:),indexCand(:)
    integer   ::  nMaster,idm,nMasterMax,bktID,nCand
    logical :: is_present_B
    real(kind=kreal), pointer :: bp(:)
 
    if( is_init ) then
      distclr = embed%cparam%DISTCLR_INIT
    else
      distclr = embed%cparam%DISTCLR_FREE
      active = .false.
      do i= 1, size(embed%slave)
        if( embed%states(i)%state==contactstick ) then
          active = .true.
          exit
        endif
      enddo
      return
    endif
 
    allocate(contact_surf(size(nodeid)))
    allocate(states_prev(size(embed%slave)))
    contact_surf(:) = huge(0)
    do i = 1, size(embed%slave)
      states_prev(i) = embed%states(i)%state
    enddo
 
    call update_surface_box_info( embed%master, currpos )
    call update_surface_bucket_info( embed%master, embed%master_bktDB )
 
    ! for gfortran-10: optional parameter seems not allowed within omp parallel
    is_present_b = present(b)
    if( is_present_b ) bp => b
 
    !$omp parallel do &
    !$omp& default(none) &
    !$omp& private(i,slave,id,coord,indexMaster,nMaster,nn,j,iSS,elem,idm,etype,isin, &
    !$omp&         bktID,nCand,indexCand) &
    !$omp& firstprivate(nMasterMax,is_present_B) &
    !$omp& shared(embed,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,nodeID,elemID,Bp,distclr,contact_surf) &
    !$omp& reduction(.or.:active) &
    !$omp& schedule(dynamic,1)
    do i= 1, size(embed%slave)
      slave = embed%slave(i)
      if( embed%states(i)%state==contactfree ) then
        coord(:) = currpos(3*slave-2:3*slave)
 
        ! get master candidates from bucketDB
        bktid = bucketdb_getbucketid(embed%master_bktDB, coord)
        ncand = bucketdb_getnumcand(embed%master_bktDB, bktid)
        if (ncand == 0) cycle
        allocate(indexcand(ncand))
        call bucketdb_getcand(embed%master_bktDB, bktid, ncand, indexcand)
 
        nmastermax = ncand
        allocate(indexmaster(nmastermax))
        nmaster = 0
 
        ! narrow down candidates
        do idm= 1, ncand
          id = indexcand(idm)
          if (.not. is_in_surface_box( embed%master(id), coord(1:3), embed%cparam%BOX_EXP_RATE )) cycle
          nmaster = nmaster + 1
          indexmaster(nmaster) = id
        enddo
        deallocate(indexcand)
 
        if(nmaster == 0) then
          deallocate(indexmaster)
          cycle
        endif
 
        do idm = 1,nmaster
          id = indexmaster(idm)
          etype = embed%master(id)%etype
          if( mod(etype,10) == 2 ) etype = etype - 1 !search by 1st-order shape function
          nn = getnumberofnodes(etype)
          do j=1,nn
            iss = embed%master(id)%nodes(j)
            elem(1:3,j)=currpos(3*iss-2:3*iss)
          enddo
          call project_point2solidelement( coord,etype,nn,elem,embed%master(id)%reflen,embed%states(i), &
            isin,distclr,localclr=embed%cparam%CLEARANCE )
          if( .not. isin ) cycle
          embed%states(i)%surface = id
          embed%states(i)%multiplier(:) = 0.d0
          contact_surf(embed%slave(i)) = elemid(embed%master(id)%eid)
          if (contact_log_level >= 1) write(*,'(A,i10,A,i10,A,3f7.3,A,i6)') "Node",nodeid(slave)," embeded to element", &
            elemid(embed%master(id)%eid), " at ",embed%states(i)%lpos(:)," rank=",hecmw_comm_get_rank()
          exit
        enddo
        deallocate(indexmaster)
      endif
    enddo
    !$omp end parallel do
 
    call hecmw_contact_comm_allreduce_i(embed%comm, contact_surf, hecmw_min)
    nactive = 0
    do i = 1, size(embed%slave)
      if (embed%states(i)%state /= contactfree) then                    ! any slave in contact
        id = embed%states(i)%surface
        if (abs(contact_surf(embed%slave(i))) /= elemid(embed%master(id)%eid)) then ! that is in contact with other surface
          embed%states(i)%state = contactfree                           ! should be freed
          if (contact_log_level >= 1) write(*,'(A,i10,A,i10,A,i6,A,i6,A)') "Node",nodeid(embed%slave(i))," contact with element", &
          &  elemid(embed%master(id)%eid), " in rank",hecmw_comm_get_rank()," freed due to duplication"
        else
          nactive = nactive + 1
        endif
      endif
      if (states_prev(i) == contactfree .and. embed%states(i)%state /= contactfree) then
        infoctchange%free2contact = infoctchange%free2contact + 1
      elseif (states_prev(i) /= contactfree .and. embed%states(i)%state == contactfree) then
        infoctchange%contact2free = infoctchange%contact2free + 1
      endif
    enddo
    active = (nactive > 0)
    deallocate(contact_surf)
    deallocate(states_prev)
  end subroutine scan_embed_state
 
  subroutine remove_duplication_tiedcontact( cstep, hecMESH, fstrSOLID, infoCTChange )
    integer(kind=kint), intent(in)         :: cstep
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    type(fstr_info_contactchange), intent(inout):: infoCTChange
 
    integer(kind=kint) :: i, j, grpid, slave
    integer(kind=kint) :: k, id, iSS
    integer(kind=kint) :: ig0, ig, iS0, iE0
    integer(kind=kint), allocatable :: states(:)
 
    allocate(states(hecmesh%n_node))
    states(:) = contactfree
 
    ! if a boundary condition is given, the slave
    do ig0= 1, fstrsolid%BOUNDARY_ngrp_tot
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
      ig= fstrsolid%BOUNDARY_ngrp_ID(ig0)
      is0= hecmesh%node_group%grp_index(ig-1) + 1
      ie0= hecmesh%node_group%grp_index(ig  )
      do k= is0, ie0
        iss = hecmesh%node_group%grp_item(k)
        !states(iSS) = CONTACTSTICK
      enddo
    enddo
 
    do i=1,fstrsolid%n_contacts
      if( fstrsolid%contacts(i)%algtype /= contacttied ) cycle
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      do j=1, size(fstrsolid%contacts(i)%slave)
        if( is_contact_free(fstrsolid%contacts(i)%states(j)%state) ) cycle   ! free
        slave = fstrsolid%contacts(i)%slave(j)
        if( is_contact_free(states(slave)) ) then
          states(slave) = fstrsolid%contacts(i)%states(j)%state
          id = fstrsolid%contacts(i)%states(j)%surface
          do k=1,size( fstrsolid%contacts(i)%master(id)%nodes )
            iss = fstrsolid%contacts(i)%master(id)%nodes(k)
            states(iss) = fstrsolid%contacts(i)%states(j)%state
          enddo
        else !found duplicate tied contact slave node
          fstrsolid%contacts(i)%states(j)%state = contactfree
          infoctchange%free2contact = infoctchange%free2contact - 1
          if (contact_log_level >= 1) write(*,'(A,i10,A,i6,A,i6,A)') "Node",hecmesh%global_node_ID(slave), &
            " in rank",hecmw_comm_get_rank()," freed due to duplication"
        endif
      enddo
    enddo
 
  end subroutine
 
end module m_fstr_contact_search