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
  implicit none
 
contains
 
  subroutine reset_contact_force( contact, currpos, lslave, omaster, opos, odirec, B )
    type( tcontact ), intent(inout)   :: contact
    real(kind=kreal), intent(in)      :: currpos(:)     
    integer, intent(in)               :: lslave
    integer, intent(in)               :: omaster
    real(kind=kreal), intent(in)      :: opos(2)        
    real(kind=kreal), intent(in)      :: odirec(3)      
    real(kind=kreal), intent(inout)   :: b(:)           
 
    integer(kind=kint)  :: slave,  etype, master
    integer(kind=kint)  :: nn, j, iSS
    real(kind=kreal)    :: nrlforce, fcoeff, tangent(3,2)
    real(kind=kreal)    :: elemcrd(3, l_max_elem_node )
    real(kind=kreal)    :: shapefunc(l_max_surface_node)
    real(kind=kreal)    :: metric(2,2), dispmat(2,l_max_elem_node*3+3)
    real(kind=kreal)    :: fric(2), f3(l_max_elem_node*3+3)
    integer(kind=kint)  :: i, idx0
 
    slave = contact%slave(lslave)
    fcoeff = contact%fcoeff
    ! clear contact force in former contacted element
    nrlforce = contact%states(lslave)%multiplier(1)
    b(3*slave-2:3*slave) = b(3*slave-2:3*slave)+nrlforce*odirec
    nn = size( contact%master(omaster)%nodes )
    etype = contact%master(omaster)%etype
    call getshapefunc( etype, opos(:), shapefunc )
    do j=1,nn
      iss = contact%master(omaster)%nodes(j)
      ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-nrlforce*shapefunc(j)*odirec
      idx0 = 3*(iss-1)
      do i=1,3
        !$omp atomic
        b(idx0+i) = b(idx0+i)-nrlforce*shapefunc(j)*odirec(i)
      enddo
    enddo
    if( fcoeff/=0.d0 ) then
      do j=1,nn
        iss = contact%master(omaster)%nodes(j)
        elemcrd(:,j) = currpos(3*iss-2:3*iss)
      enddo
      call dispincrematrix( opos(:), etype, nn, elemcrd, tangent,   &
        metric, dispmat )
      fric(1:2) = contact%states(lslave)%multiplier(2:3)
      f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
      b(3*slave-2:3*slave) = b(3*slave-2:3*slave)+f3(1:3)
      do j=1,nn
        iss = contact%master(omaster)%nodes(j)
        ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)+f3(3*j+1:3*j+3)
        idx0 = 3*(iss-1)
        do i=1,3
          !$omp atomic
          b(idx0+i) = b(idx0+i)+f3(3*j+i)
        enddo
      enddo
    endif
 
    ! reset contact force in new contacting element
    master = contact%states(lslave)%surface
    nn = size( contact%master(master)%nodes )
    etype = contact%master(master)%etype
    call getshapefunc( etype, contact%states(lslave)%lpos(1:2), shapefunc )
    b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-nrlforce*contact%states(lslave)%direction
    do j=1,nn
      iss = contact%master(master)%nodes(j)
      ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)+nrlforce*        &
      !            shapefunc(j)*contact%states(lslave)%direction
      idx0 = 3*(iss-1)
      do i=1,3
        !$omp atomic
        b(idx0+i) = b(idx0+i)+nrlforce*        &
          shapefunc(j)*contact%states(lslave)%direction(i)
      enddo
    enddo
    if( fcoeff/=0.d0 ) then
      do j=1,nn
        iss = contact%master(master)%nodes(j)
        elemcrd(:,j) = currpos(3*iss-2:3*iss)
      enddo
      call dispincrematrix( contact%states(lslave)%lpos(1:2), etype, nn, elemcrd, tangent,   &
        metric, dispmat )
      fric(1:2) = contact%states(lslave)%multiplier(2:3)
      f3(:) = fric(1)*dispmat(1,:)+fric(2)*dispmat(2,:)
      b(3*slave-2:3*slave) = b(3*slave-2:3*slave)-f3(1:3)
      do j=1,nn
        iss = contact%master(master)%nodes(j)
        ! B(3*iSS-2:3*iSS) = B(3*iSS-2:3*iSS)-f3(3*j+1:3*j+3)
        idx0 = 3*(iss-1)
        do i=1,3
          !$omp atomic
          b(idx0+i) = b(idx0+i)-f3(3*j+i)
        enddo
      enddo
    endif
 
  end subroutine reset_contact_force
 
  subroutine track_contact_position( flag_ctAlgo, nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID, B )
    character(len=9), intent(in)                    :: 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(:)
    real(kind=kreal), intent(inout)                  :: b(:)         
 
    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), odirec(3)
    integer(kind=kint) :: bktID, nCand, idm
    integer(kind=kint), allocatable :: indexCand(:)
 
    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)
    odirec = contact%states(nslave)%direction
    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 )
    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 )
        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 )
          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
        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
        if( flag_ctalgo=='ALagrange' )  &
          call reset_contact_force( contact, currpos, nslave, sid0, opos, odirec, b )
      endif
      if( 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
      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 track_contact_position_exp( nslave, contact, currpos, currdisp, infoCTChange, nodeID, elemID )
    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(:)
 
    integer(kind=kint) :: slave, sid0, sid, etype
    integer(kind=kint) :: nn, i, j, iSS
    real(kind=kreal)    :: coord(3), elem0(3, l_max_elem_node )
    logical            :: isin
    real(kind=kreal)    :: opos(2), odirec(3)
    integer(kind=kint) :: bktID, nCand, idm
    integer(kind=kint), allocatable :: indexCand(:)
 
    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)
    odirec = contact%states(nslave)%direction
    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 )
    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 )
        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( any(sid==contact%master(sid0)%neighbor(:)) ) cycle
          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 )
          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
        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
      iss = isinsideelement( etype, contact%states(nslave)%lpos(1:2), contact%cparam%CLR_CAL_NORM )
      if( iss>0 ) then
        call cal_node_normal( contact%states(nslave)%surface, iss, contact%master, currpos, &
          contact%states(nslave)%lpos(1:2), contact%states(nslave)%direction(:) )
      endif
    else if( .not. isin ) then
      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_exp
 
  subroutine scan_contact_state( flag_ctAlgo, contact, currpos, currdisp, ndforce, infoCTChange, &
    nodeID, elemID, is_init, active, B )
    character(len=9), intent(in)                     :: flag_ctAlgo
    type( tcontact ), intent(inout)                  :: contact
    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)  :: i, iSS, nactive
    real(kind=kreal)    :: coord(3)
    real(kind=kreal)    :: nlforce, slforce(3)
    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 = 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
 
    call update_surface_box_info( contact%master, currpos )
    call update_surface_bucket_info( contact%master, contact%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,slforce,id,nlforce,coord,indexMaster,nMaster,iSS,idm,etype,isin, &
    !$omp&         bktID,nCand,indexCand) &
    !$omp& firstprivate(nMasterMax,is_present_B) &
    !$omp& shared(contact,ndforce,flag_ctAlgo,infoCTChange,currpos,currdisp,nodeID,elemID,Bp,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
        slforce(1:3)=ndforce(3*slave-2:3*slave)
        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
          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 .or. (.not. is_present_b) ) then   ! small slide problem
          contact_surf(contact%slave(i)) = -elemid(contact%master(id)%eid)
        else
          call track_contact_position( flag_ctalgo, i, contact, currpos, currdisp, infoctchange, nodeid, elemid, bp )
          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==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)
          call project_point2surfelement( coord, contact%master(id), currpos, &
            contact%states(i), isin, distclr, localclr=contact%cparam%CLEARANCE )
          if( .not. isin ) cycle
          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)
          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()
          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 (contact%states(i)%state /= contactfree) then                    ! any slave in contact
        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
          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
          nactive = nactive + 1
        endif
      endif
      if (states_prev(i) == contactfree .and. contact%states(i)%state /= contactfree) then
        infoctchange%free2contact = infoctchange%free2contact + 1
      elseif (states_prev(i) /= contactfree .and. contact%states(i)%state == contactfree) then
        infoctchange%contact2free = infoctchange%contact2free + 1
      endif
    enddo
    active = (nactive > 0)
    deallocate(contact_surf)
    deallocate(states_prev)
  end subroutine scan_contact_state
 
  subroutine scan_contact_state_exp( contact, currpos, currdisp, infoCTChange, &
    nodeID, elemID, is_init, active )
    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
 
    real(kind=kreal)    :: distclr
    integer(kind=kint)  :: slave, id, etype
    integer(kind=kint)  :: i, iSS, nactive
    real(kind=kreal)    :: coord(3)
    logical             :: isin
    integer(kind=kint), allocatable :: contact_surf(:), states_prev(:)
    !
    integer, pointer :: indexMaster(:),indexCand(:)
    integer   ::  nMaster,idm,nMasterMax,bktID,nCand
 
    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
        return
      endif
    endif
 
    allocate(contact_surf(size(nodeid)))
    allocate(states_prev(size(contact%slave)))
    contact_surf(:) = size(elemid)+1
    do i = 1, size(contact%slave)
      states_prev(i) = contact%states(i)%state
    enddo
 
    call update_surface_box_info( contact%master, currpos )
    call update_surface_bucket_info( contact%master, contact%master_bktDB )
 
    !$omp parallel do &
    !$omp& default(none) &
    !$omp& private(i,slave,id,coord,indexMaster,nMaster,iSS,idm,etype,isin, &
    !$omp&         bktID,nCand,indexCand) &
    !$omp& firstprivate(nMasterMax) &
    !$omp& shared(contact,infoCTChange,currpos,currdisp,nodeID,elemID,distclr,contact_surf) &
    !$omp& reduction(.or.:active) &
    !$omp& schedule(dynamic,1)
    do i= 1, size(contact%slave)
      slave = contact%slave(i)
      if( contact%states(i)%state==contactstick .or. contact%states(i)%state==contactslip ) then
        call track_contact_position_exp( i, contact, currpos, currdisp, infoctchange, nodeid, elemid )
        if( contact%states(i)%state /= contactfree ) then
          contact_surf(contact%slave(i)) = -contact%states(i)%surface
        endif
      else if( contact%states(i)%state==contactfree ) then
        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)
          call project_point2surfelement( coord, contact%master(id), currpos, &
            contact%states(i), isin, distclr, localclr=contact%cparam%CLEARANCE )
          if( .not. isin ) cycle
          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)) = id
          write(*,'(A,i10,A,i10,A,f7.3,A,2f7.3,A,3f7.3)') "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
          exit
        enddo
        deallocate(indexmaster)
      endif
    enddo
    !$omp end parallel do
 
    call hecmw_contact_comm_allreduce_i(contact%comm, contact_surf, hecmw_min)
    nactive = 0
    do i = 1, size(contact%slave)
      if (contact%states(i)%state /= contactfree) then                    ! any slave in contact
        if (abs(contact_surf(contact%slave(i))) /= contact%states(i)%surface) then ! that is in contact with other surface
          contact%states(i)%state = contactfree                           ! should be freed
          write(*,'(A,i10,A,i6,A,i6,A)') "Node",nodeid(contact%slave(i)), &
            " in rank",hecmw_comm_get_rank()," freed due to duplication"
        else
          nactive = nactive + 1
        endif
      endif
      if (states_prev(i) == contactfree .and. contact%states(i)%state /= contactfree) then
        infoctchange%free2contact = infoctchange%free2contact + 1
      elseif (states_prev(i) /= contactfree .and. contact%states(i)%state == contactfree) then
        infoctchange%contact2free = infoctchange%contact2free + 1
      endif
    enddo
    active = (nactive > 0)
    deallocate(contact_surf)
    deallocate(states_prev)
  end subroutine scan_contact_state_exp
 
  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)
          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
          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( fstrsolid%contacts(i)%states(j)%state==contactfree ) cycle   ! free
        slave = fstrsolid%contacts(i)%slave(j)
        if( states(slave) == contactfree ) 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
          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