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
  use m_fstr_contact_elem_alag, only: get_unique_map
  use m_fstr_contact_damping, only: is_damping_enabled
  implicit none
 
  integer(kind=kint), parameter :: contact_log_level = 0
 
contains
 
  subroutine is_known_master_in_segment(contact, known_masters, n_known, id, is_known)
    type(tcontact), intent(in)         :: contact
    integer(kind=kint), intent(inout)  :: known_masters(:)
    integer(kind=kint), intent(inout)  :: n_known
    integer(kind=kint), intent(in)     :: id
    logical, intent(out)               :: is_known
 
    integer(kind=kint) :: k
 
    is_known = any(known_masters(1:n_known) == id)
    if( .not. is_known .and. contact%sparsity_expansion == sparsity_neighbor ) then
      ! neighbors of known masters are pre-registered in sparsity
      do k = 1, n_known
        if( associated(contact%master(known_masters(k))%neighbor) ) then
          if( any(contact%master(known_masters(k))%neighbor(1:contact%master(known_masters(k))%n_neighbor) == id) ) then
            is_known = .true.; exit
          endif
        endif
      enddo
    endif
    if( .not. is_known ) then
      n_known = n_known + 1
      known_masters(n_known) = id
    endif
  end subroutine is_known_master_in_segment
 
  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, id_best
    integer(kind=kint), allocatable :: indexCand(:)
    logical            :: is_implicit, update_tangent
    type(tcontactstate) :: cstate_free, cstate_try, cstate_best
 
    is_implicit = present(flag_ctalgo)
    if( is_implicit ) then
      update_tangent = (flag_ctalgo == 'SLagrange')
    else
      update_tangent = .true.
    endif
 
    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%cparam%PENCLR_NOCHECK, &
      contact%states(nslave)%lpos(1:2), contact%cparam%CLR_SAME_ELEM, smoothing=contact%smoothing )
    if( .not. isin ) then
      ! The neighbors include the faces across a convex edge, which the loose checks of
      ! DISTCLR_NOCHECK and PENCLR_NOCHECK accept at a large depth, so the nearest surface
      ! is taken here as well as in the search of free nodes.
      cstate_free = contact%states(nslave)
      id_best = 0
      do i=1, contact%master(sid0)%n_neighbor
        sid = contact%master(sid0)%neighbor(i)
        cstate_try = cstate_free
        call project_point2surfelement( coord, contact%master(sid), currpos, &
          cstate_try, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_NOCHECK, &
          localclr=contact%cparam%CLEARANCE, smoothing=contact%smoothing )
        if( .not. isin ) cycle
        if( id_best /= 0 ) then
          if( dabs(cstate_try%distance) > dabs(cstate_best%distance) ) cycle
          if( dabs(cstate_try%distance) == dabs(cstate_best%distance) .and. sid > id_best ) cycle
        endif
        id_best = sid
        cstate_best = cstate_try
      enddo
      isin = ( id_best /= 0 )
      if( isin ) then
        sid = id_best
        contact%states(nslave) = cstate_best
        contact%states(nslave)%surface = sid
      endif
    endif
 
    if( .not. isin ) then   ! such case is considered to rarely or never occur
      !$omp critical
      write(*,*) 'Warning: contact moved beyond neighbor elements'
      !$omp end critical
      cstate_free = contact%states(nslave)
      ! 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)
        id_best = 0
        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(1:contact%master(sid0)%n_neighbor)) ) cycle
          endif
          cstate_try = cstate_free
          call project_point2surfelement( coord, contact%master(sid), currpos, &
            cstate_try, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_FREE, &
            localclr=contact%cparam%CLEARANCE, smoothing=contact%smoothing )
          if( .not. isin ) cycle
          if( id_best /= 0 ) then
            if( dabs(cstate_try%distance) > dabs(cstate_best%distance) ) cycle
            if( dabs(cstate_try%distance) == dabs(cstate_best%distance) .and. sid > id_best ) cycle
          endif
          id_best = sid
          cstate_best = cstate_try
        enddo
        deallocate(indexcand)
        isin = ( id_best /= 0 )
        if( isin ) then
          sid = id_best
          contact%states(nslave) = cstate_best
          contact%states(nslave)%surface = sid
        endif
      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) then
          !$omp critical
          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 end critical
        endif
        !$omp atomic
        infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
      endif
      if( update_tangent .and. contact%fcoeff /= 0.d0 ) 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 .and. contact%smoothing /= kcsnagata ) &
        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) then
        !$omp critical
        write(*,'(A,i10,A)') "Node",nodeid(slave)," move out of contact"
        !$omp end critical
      endif
      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, icat_prev, icat_curr, surf_key
    real(kind=kreal)    :: coord(3)
    real(kind=kreal)    :: nlforce
    logical             :: isin
    integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
    integer(kind=kint), allocatable :: contact_prio(:), contact_prio_min(:)
    real(kind=kreal)    :: effective_near_dist, distclr_use
    !
    integer, pointer :: indexCand(:)
    integer   ::  idm,bktID,nCand,id_best
    type(tcontactstate) :: cstate_free, cstate_try, cstate_best
    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(contact_prio(size(nodeid)))
    allocate(states_prev(size(contact%slave)))
    contact_surf(:) = huge(0)
    contact_prio(:) = 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,iSS,idm,etype,isin, &
    !$omp&         bktID,nCand,indexCand,distclr_use,id_best,cstate_free,cstate_try,cstate_best) &
    !$omp& firstprivate(is_implicit,effective_near_dist) &
    !$omp& shared(contact,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,nodeID,elemID,distclr,contact_surf,contact_prio, &
    !$omp&        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)%multiplier(:) = 0.d0
          if( effective_near_dist > 0.0d0 ) then
            ! keep the projection so that the NEAR branch re-projects on the next scan and
            ! decides between NEAR and FREE with an up-to-date distance
            contact%states(i)%state = contactnear
            contact_surf(contact%slave(i)) = elemid(contact%master(id)%eid)
            contact_prio(contact%slave(i)) = 1
            if (contact_log_level >= 1) then
              !$omp critical
              write(*,'(A,i10,A,i10,A,e12.3)') "Node",nodeid(slave), &
                " released to near contact with element", &
                elemid(contact%master(id)%eid), " with tensile force ", nlforce
              !$omp end critical
            endif
          else
            contact%states(i)%state = contactfree
            if (contact_log_level >= 1) then
              !$omp critical
              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
              !$omp end critical
            endif
          endif
          cycle
        endif
        if( contact%algtype /= contactfslid ) then   ! small slide problem
          contact_surf(contact%slave(i)) = -elemid(contact%master(id)%eid)
          contact_prio(contact%slave(i)) = 0
        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)
            contact_prio(contact%slave(i)) = 0
          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%cparam%PENCLR_FREE, &
          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) then
              !$omp critical
              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()
              !$omp end critical
            endif
          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 .and. contact%smoothing /= kcsnagata ) &
              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)
            contact_prio(contact%slave(i)) = 1
          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)
 
        ! The bucket collects candidates in scan order, not in distance order, so every
        ! candidate is projected on a trial state and the closest one is adopted. Equal
        ! distances are settled by the master index so that the choice never depends on
        ! the order the bucket happened to return.
        id_best = 0
        cstate_free = contact%states(i)
        do idm = 1,ncand
          id = indexcand(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
          cstate_try = cstate_free
          call project_point2surfelement( coord, contact%master(id), currpos, &
            cstate_try, isin, distclr_use, contact%cparam%PENCLR_FREE, &
            localclr=contact%cparam%CLEARANCE, smoothing=contact%smoothing )
          if( .not. isin ) cycle
          if( id_best /= 0 ) then
            if( dabs(cstate_try%distance) > dabs(cstate_best%distance) ) cycle
            if( dabs(cstate_try%distance) == dabs(cstate_best%distance) .and. id > id_best ) cycle
          endif
          id_best = id
          cstate_best = cstate_try
        enddo
        deallocate(indexcand)
        if( id_best == 0 ) cycle
 
        id = id_best
        contact%states(i) = cstate_best
        ! 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 .and. contact%smoothing /= kcsnagata ) &
          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)
        contact_prio(contact%slave(i)) = 1
        if (contact_log_level >= 1) then
          !$omp critical
          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
          !$omp end critical
        end if
      endif
    enddo
    !$omp end parallel do
 
    if( contact%algtype == contacttied .and. .not. is_init ) then
      deallocate(contact_surf)
      deallocate(contact_prio)
      deallocate(states_prev)
      return
    endif
 
    ! continued contacts take precedence; the sign of contact_surf cannot tell them apart
    ! because a refined mesh has negative element IDs
    allocate(contact_prio_min(size(nodeid)))
    contact_prio_min(:) = contact_prio(:)
    call hecmw_contact_comm_allreduce_i(contact%comm, contact_prio_min, hecmw_min)
    where( contact_prio(:) > contact_prio_min(:) ) contact_surf(:) = huge(0)
    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
        surf_key = elemid(contact%master(id)%eid)
        if (contact_prio(contact%slave(i)) == 0) surf_key = -surf_key
        if (contact_prio(contact%slave(i)) /= contact_prio_min(contact%slave(i)) .or. &
          & contact_surf(contact%slave(i)) /= surf_key) 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) .or. is_damping_enabled(contact)) then
          nactive = nactive + 1   ! a NEAR node with damping also contributes stiffness and residual
        endif
      endif
      icat_prev = contact_state_category(states_prev(i))
      icat_curr = contact_state_category(contact%states(i)%state)
      if (icat_prev /= icat_curr) infoctchange%n_statechange(icat_prev,icat_curr) = &
        infoctchange%n_statechange(icat_prev,icat_curr) + 1
    enddo
    active = (nactive > 0)
    deallocate(contact_surf)
    deallocate(contact_prio, contact_prio_min)
    deallocate(states_prev)
  end subroutine scan_contact_state
 
  subroutine scan_embed_state( embed, currpos, infoCTChange, &
    nodeID, elemID, is_init, active, B )
    type( tcontact ), intent(inout)                  :: embed
    type( fstr_info_contactchange ), intent(inout)   :: infoCTChange
    real(kind=kreal), intent(in)                     :: currpos(:)   
    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, icat_prev, icat_curr
    real(kind=kreal)    :: coord(3), elem(3, l_max_elem_node )
    logical             :: isin
    integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
    !
    integer, pointer :: indexCand(:)
    integer   ::  idm,bktID,nCand,id_best
    type(tcontactstate) :: cstate_free, cstate_try, cstate_best
    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,nn,j,iSS,elem,idm,etype,isin, &
    !$omp&         bktID,nCand,indexCand,id_best,cstate_free,cstate_try,cstate_best) &
    !$omp& firstprivate(is_present_B) &
    !$omp& shared(embed,infoCTChange,currpos,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)
 
        ! As in scan_contact_state, the bucket order is not a distance order, so every
        ! candidate is evaluated on a trial state and the closest element is adopted.
        id_best = 0
        cstate_free = embed%states(i)
        do idm = 1,ncand
          id = indexcand(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
          cstate_try = cstate_free
          call project_point2solidelement( coord,etype,nn,elem,cstate_try, &
            isin,distclr,localclr=embed%cparam%CLEARANCE )
          if( .not. isin ) cycle
          if( id_best /= 0 ) then
            if( cstate_try%distance > cstate_best%distance ) cycle
            if( cstate_try%distance == cstate_best%distance .and. id > id_best ) cycle
          endif
          id_best = id
          cstate_best = cstate_try
        enddo
        deallocate(indexcand)
        if( id_best == 0 ) cycle
 
        id = id_best
        embed%states(i) = cstate_best
        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) then
          !$omp critical
          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()
          !$omp end critical
        endif
      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 (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
      icat_prev = contact_state_category(states_prev(i))
      icat_curr = contact_state_category(embed%states(i)%state)
      if (icat_prev /= icat_curr) infoctchange%n_statechange(icat_prev,icat_curr) = &
        infoctchange%n_statechange(icat_prev,icat_curr) + 1
    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, icat
    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
          icat = contact_state_category(fstrsolid%contacts(i)%states(j)%state)
          fstrsolid%contacts(i)%states(j)%state = contactfree
          infoctchange%n_statechange(kcatfree,icat) = infoctchange%n_statechange(kcatfree,icat) - 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
 
  subroutine scan_contact_state_ss( contact, currpos, infoCTChange, is_init, active )
    type( tcontact ), intent(inout)                  :: contact
    type( fstr_info_contactchange ), intent(inout)   :: infoCTChange
    real(kind=kreal), intent(in)                     :: currpos(:)   
    logical, intent(in)                              :: is_init
    logical, intent(out)                             :: active
 
    real(kind=kreal)    :: distclr
    integer(kind=kint)  :: slave, id
    integer(kind=kint)  :: nnode_s, i, j
    real(kind=kreal)    :: coord(3), ncoord(2)
    real(kind=kreal)    :: surf_node_pos(3,4), sfunc(4)
    logical             :: isin
    !
    integer, pointer :: indexCand(:)
    integer   ::  idm,bktID,nCand,id_best
    type(tcontactstate) :: state_free, state_try, state_best
    ! per-segment known-master tracking for free2contact_new
    integer(kind=kint) :: maplist(MAX_N_INTP), master_idxs(MAX_N_INTP)
    integer(kind=kint) :: unique_count
    integer(kind=kint) :: known_masters(MAX_N_INTP), n_known
    logical :: is_known
    ! each kind of change is counted once per slave segment; the convergence and rebuild checks only test for zero
    logical :: counted_f2c, counted_nbr, counted_byd
    if( is_init ) then
      distclr = contact%cparam%DISTCLR_INIT
    else
      distclr = contact%cparam%DISTCLR_FREE
    endif
    call update_surface_box_info( contact%master, currpos )
    call update_surface_bucket_info( contact%master, contact%master_bktDB )
    ! the reduction below starts from the value active has on entry, so it must be defined
    active = .false.
    !$omp parallel do &
    !$omp& default(none) &
    !$omp& private(i,slave,id,coord,ncoord,surf_node_pos,sfunc,nnode_s, &
    !$omp&         j,idm,isin,bktID,nCand,indexCand,id_best,state_free,state_try,state_best, &
    !$omp&         maplist,master_idxs,unique_count,known_masters,n_known,is_known, &
    !$omp&         counted_f2c,counted_nbr,counted_byd) &
    !$omp& shared(contact,infoCTChange,currpos,distclr,is_init) &
    !$omp& reduction(.or.:active) &
    !$omp& schedule(dynamic,1)
 
    ! loop element; if any node(nslave_index) in contact, call intp point projection
    do i = 1, size(contact%slave_surf)
      surf_node_pos(:,:) = 0.d0
      nnode_s = size(contact%slave_surf(i)%nodes)
      ! Snapshot the previous-step active master set for this slave segment.
      ! Used by is_known_master_in_segment to detect genuinely new masters
      ! (= those that change matrix sparsity).
      call get_unique_map(contact%slave_surf(i), maplist, master_idxs, unique_count)
      n_known = unique_count
      if( n_known > 0 ) known_masters(1:n_known) = master_idxs(1:n_known)
      counted_f2c = .false.
      counted_nbr = .false.
      counted_byd = .false.
 
      ! ALagrange SS only (SLAG+MORTAR rejected at fstr_setup): set all IPs as
      ! candidates unconditionally (do not depend on NTS node-level contact state)
      do j = 1, contact%slave_surf(i)%n_intp
        if (contact%slave_surf(i)%states(j)%state == contactfree) then
          contact%slave_surf(i)%states(j)%state = candidate_intp
        endif
      enddo
 
      do j = 1, nnode_s
        slave = contact%slave_surf(i)%nodes(j)
        surf_node_pos(:,j) = currpos(3*slave-2:3*slave)
      enddo
 
      ! Snapshot this segment's previous state before reset (for seg on/off detection)
      contact%slave_surf(i)%state_prev = contact%slave_surf(i)%state
      contact%slave_surf(i)%state = contactfree
      do j = 1, contact%slave_surf(i)%n_intp
        ! All FREE IPs were promoted to CANDIDATE_INTP above; no IP can be CONTACTFREE here.
        ! CONTACT_INTERFERENCE (if_type) not supported with MORTAR=YES (rejected at check_apply_Contact_IF)
 
        ! IP natural coord (ncoord) and physical position (coord): used by both the
        ! already-active (track) and candidate (bucket search + projection) paths.
        call getintpoint4ss( contact%slave_surf(i)%etype, j, ncoord, contact%slave_surf(i)%n_intp, sfunc )
        coord(:) = matmul( surf_node_pos(:,1:nnode_s), sfunc(1:nnode_s) )
 
        !! already active contact
        if( contact%slave_surf(i)%states(j)%state==contactstick .or. contact%slave_surf(i)%states(j)%state==contactslip ) then
          active = .true.
          id = contact%slave_surf(i)%states(j)%surface
 
          ! INTERACTION=TIED not supported with MORTAR=YES (rejected at fstr_setup)
          if( contact%algtype /= contactfslid ) cycle  ! small slide problem
 
          call track_contact_position_ss( coord, contact%slave_surf(i)%states(j),&
            contact%slave_surf(i), ncoord, contact, currpos, infoctchange, counted_nbr, counted_byd )
 
          if( contact%slave_surf(i)%states(j)%state /= contactfree ) contact%slave_surf(i)%state = contactstick
 
        !! contact candidate integration point
        else if( contact%slave_surf(i)%states(j)%state==candidate_intp ) then
          contact%slave_surf(i)%states(j)%state=contactfree
 
          ! 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)
 
          ! The bucket returns the candidates in scan order, so every candidate is
          ! projected on a trial state and the nearest one is adopted, as in the
          ! NODE-SURF scan. A surface across a convex edge is accepted as well and
          ! would otherwise bind the integration point at a large depth.
          state_free = contact%slave_surf(i)%states(j)
          id_best = 0
          do idm = 1,ncand
            id = indexcand(idm)
            ! OFF->ON uses is_init-dependent distclr (DISTCLR_INIT on initial scan,
            ! DISTCLR_FREE on re-scan) to form a hysteresis band with DISTCLR_C2F.
            ! The _ss wrapper sets direction to the slave inward normal (mortar normal).
            state_try = state_free
            call project_point2surfelement_ss( coord, contact%master(id), contact%slave_surf(i), &
              ncoord, currpos, state_try, isin, distclr=distclr, &
              penclr=contact%cparam%PENCLR_FREE, localclr=contact%cparam%CLEARANCE )
            if( .not. isin ) cycle
            if( id_best /= 0 ) then
              if( dabs(state_try%distance) > dabs(state_best%distance) ) cycle
              if( dabs(state_try%distance) == dabs(state_best%distance) .and. id > id_best ) cycle
            endif
            id_best = id
            state_best = state_try
          enddo
          if( id_best /= 0 ) then
            id = id_best
            contact%slave_surf(i)%states(j) = state_best
            contact%slave_surf(i)%states(j)%surface = id
            contact%slave_surf(i)%states(j)%multiplier(:) = 0.d0
 
            contact%slave_surf(i)%states(j)%state = contactstick
            active = .true.
            contact%slave_surf(i)%state = contactstick
            ! Always maintain the known-master set. SPARSITY_NEIGHBOR pre-registers the neighbors, so only a
            ! genuinely-new master rebuilds / breaks convergence; current-master-only: every new contact rebuilds.
            call is_known_master_in_segment(contact, known_masters, n_known, id, is_known)
            if( (contact%sparsity_expansion == sparsity_none .or. .not. is_known) .and. .not. counted_f2c ) then
              !$omp atomic
              infoctchange%free2contact_new = infoctchange%free2contact_new + 1
              counted_f2c = .true.
            endif
          endif
          deallocate(indexcand)
        endif
      enddo
 
    enddo
 
  end subroutine scan_contact_state_ss
 
  subroutine track_contact_position_ss( coord, state, sSurf, ncoord, contact, currpos, infoCTChange, counted_nbr, counted_byd )
    type( tcontact ), intent(inout)                  :: contact
    type( tcontactstate ), intent(inout)             :: state
    type( tcontactsurf ), intent(in)                 :: sSurf
    real(kind=kreal), intent(in)                     :: ncoord(2)  
    type( fstr_info_contactchange ), intent(inout)   :: infoctchange
    logical, intent(inout)                           :: counted_nbr
    logical, intent(inout)                           :: counted_byd
    real(kind=kreal), intent(in)                     :: currpos(:)   
    real(kind=kreal), intent(inout)                  :: coord(:)         
 
    integer(kind=kint) :: sid0, sid
    integer(kind=kint) :: i
    logical            :: isin, found_in_neighbor
    real(kind=kreal)    :: opos(2)
    integer(kind=kint) :: bktid, ncand, idm, id_best
    integer(kind=kint), allocatable :: indexcand(:)
    type(tcontactstate) :: state_free, state_try, state_best
 
    sid = 0
    found_in_neighbor = .false.
 
    sid0 = state%surface
    opos = state%lpos(1:2)
 
    call project_point2surfelement_ss( coord, contact%master(sid0), ssurf, ncoord, currpos, &
      state, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_NOCHECK, &
      state%lpos(1:2), contact%cparam%CLR_SAME_ELEM )
    if( .not. isin ) then ! not contact previous master, search neighbor master surf
      ! A neighbor across a convex edge is accepted at a large depth, so the nearest
      ! surface is taken instead of the first accepted one.
      state_free = state
      id_best = 0
      do i=1, contact%master(sid0)%n_neighbor
        sid = contact%master(sid0)%neighbor(i)
        state_try = state_free
        call project_point2surfelement_ss( coord, contact%master(sid), ssurf, ncoord, currpos, &
          state_try, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_NOCHECK, &
          localclr=contact%cparam%CLEARANCE )
        if( .not. isin ) cycle
        if( id_best /= 0 ) then
          if( dabs(state_try%distance) > dabs(state_best%distance) ) cycle
          if( dabs(state_try%distance) == dabs(state_best%distance) .and. sid > id_best ) cycle
        endif
        id_best = sid
        state_best = state_try
      enddo
      isin = ( id_best /= 0 )
      if( isin ) then
        sid = id_best
        state = state_best
        state%surface = sid
        found_in_neighbor = .true.
      endif
    endif
 
    if( .not. isin ) then   ! such case is considered to rarely or never occur
      state_free = state
      ! 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)
        id_best = 0
        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(1:contact%master(sid0)%n_neighbor)) ) cycle
          endif
          state_try = state_free
          call project_point2surfelement_ss( coord, contact%master(sid), ssurf, ncoord, currpos, &
            state_try, isin, contact%cparam%DISTCLR_NOCHECK, contact%cparam%PENCLR_FREE, &
            localclr=contact%cparam%CLEARANCE )
          if( .not. isin ) cycle
          if( id_best /= 0 ) then
            if( dabs(state_try%distance) > dabs(state_best%distance) ) cycle
            if( dabs(state_try%distance) == dabs(state_best%distance) .and. sid > id_best ) cycle
          endif
          id_best = sid
          state_best = state_try
        enddo
        deallocate(indexcand)
        isin = ( id_best /= 0 )
        if( isin ) then
          sid = id_best
          state = state_best
          state%surface = sid
        endif
      endif
    endif
 
    if( isin ) then
      if( state%distance > contact%cparam%DISTCLR_C2F * contact%master(state%surface)%reflen ) then
        state%state = contactfree
        state%multiplier(:) = 0.d0
        return
      endif
      ! An IP that stays on the same master is not counted: unlike NODE-SURF, its local position
      ! never settles under sliding, so the jitter is not a state change.
      if( state%surface /= sid0 ) then
        if( found_in_neighbor ) then
          ! SPARSITY_NEIGHBOR: neighbor is pre-registered, tolerated (no rebuild / no relance)
          if( contact%sparsity_expansion == sparsity_none .and. .not. counted_nbr ) then
            !$omp atomic
            infoctchange%contact2neighbor = infoctchange%contact2neighbor + 1
            counted_nbr = .true.
          endif
        else if( .not. counted_byd ) then
          !$omp atomic
          infoctchange%contact2beyond = infoctchange%contact2beyond + 1
          counted_byd = .true.
        endif
      endif
      ! direction already set to the slave inward normal by project_Point2SurfElement_ss
    else if( .not. isin ) then
      state%state = contactfree
      state%multiplier(:) = 0.d0
    endif
 
  end subroutine track_contact_position_ss
 
end module m_fstr_contact_search