fstr_contact_elem_common.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_elem_common
  use hecmw
  use elementinfo
  use mcontactdef
  use m_fstr_contact_geom
  use m_fstr_contact_smoothing
  implicit none
 
  public :: computetm_tt
  public :: computecontactmaps_alag
  public :: computefrictionforce_alag
  public :: gettrialfricforceandcheckfricstate
  public :: update_tangentforce
 
contains
 
  subroutine computetm_tt(ctState, tSurf, fcoeff, Tm, Tt, smoothing_type, Bn)
    implicit none
 
    ! Input arguments
    type(tcontactstate), intent(in) :: ctstate
    type(tsurfelement), intent(in)  :: tsurf
    real(kind=kreal), intent(in)    :: fcoeff          
    integer(kind=kint), optional, intent(in) :: smoothing_type
 
    ! Output arguments
    real(kind=kreal), intent(out)   :: tm(3, 3*(l_max_surface_node+1)) 
    real(kind=kreal), intent(out)   :: tt(3, 3*(l_max_surface_node+1)) 
    real(kind=kreal), optional, intent(out) :: bn(:)   
 
    ! Local variables
    integer(kind=kint) :: i, j
    integer(kind=kint) :: nnode
    integer(kind=kint) :: smoothing
    real(kind=kreal)   :: shapefunc(l_max_surface_node) 
    real(kind=kreal)   :: p_matrix(3, 3*l_max_surface_node)  
    real(kind=kreal)   :: normal(3)       
    real(kind=kreal)   :: pt(3,3)         
 
    tm = 0.0d0
    tt = 0.0d0
 
    nnode = size(tsurf%nodes)
    normal(1:3) = ctstate%direction(1:3)
 
    ! Determine smoothing type
    smoothing = kcsnone
    if (present(smoothing_type)) smoothing = smoothing_type
 
    ! Construct Tm
    if ( smoothing == kcsnone ) then
      ! Standard linear shape functions
      call getshapefunc(tsurf%etype, ctstate%lpos, shapefunc)
 
      ! First block (slave node): Identity matrix I_3
      tm(1,1) = 1.0d0
      tm(2,2) = 1.0d0
      tm(3,3) = 1.0d0
 
      ! Remaining blocks (master nodes): -N_i * I_3
      do i = 1, nnode
        tm(1, i*3+1) = -shapefunc(i)
        tm(2, i*3+2) = -shapefunc(i)
        tm(3, i*3+3) = -shapefunc(i)
      enddo
 
    else if ( smoothing == kcsnagata ) then
      ! Use Nagata patch interpolation
      call compute_interpolation_matrix_p(tsurf%etype, nnode, ctstate%lpos, &
                                          tsurf%vertex_normals, p_matrix)
 
      ! Tm = [I_3; -P_matrix]
      tm(1,1) = 1.0d0
      tm(2,2) = 1.0d0
      tm(3,3) = 1.0d0
      tm(1:3, 4:3*(nnode+1)) = -p_matrix(1:3, 1:3*nnode)
 
    endif
 
    ! Optionally compute normal distribution vector: Bn = Tm^T * n
    if (present(bn)) then
      bn(:) = 0.0d0
      bn(1:3*(nnode+1)) = matmul(transpose(tm(1:3, 1:3*(nnode+1))), normal(1:3))
    endif
 
    ! Compute Tt only if friction coefficient is non-zero
    if (fcoeff /= 0.0d0) then
      ! Construct tangential projection operator Pt = I_3 - n x n
      do i = 1, 3
        do j = 1, 3
          if (i == j) then
            pt(i,j) = 1.0d0 - normal(i)*normal(j)
          else
            pt(i,j) = -normal(i)*normal(j)
          endif
        enddo
      enddo
 
      ! Compute Tt = Pt * Tm using matrix multiplication
      tt(1:3, 1:3*(l_max_surface_node+1)) = matmul(pt(1:3,1:3), tm(1:3, 1:3*(l_max_surface_node+1)))
    endif
 
  end subroutine computetm_tt
 
  subroutine computecontactmaps_alag(ctState, tSurf, ele, &
                                      Bn, metric, Ht, Gt, smoothing_type)
    implicit none
 
    ! Input arguments
    type(tcontactstate), intent(in) :: ctstate
    type(tsurfelement), intent(in)  :: tsurf
    real(kind=kreal), intent(in)    :: ele(:,:)        
    integer(kind=kint), optional, intent(in) :: smoothing_type
 
    ! Output arguments
    real(kind=kreal), intent(out)   :: bn(:)           
    real(kind=kreal), intent(out)   :: metric(2,2)     
    real(kind=kreal), intent(out)   :: ht(:,:)         
    real(kind=kreal), intent(out)   :: gt(:,:)         
 
    ! Local variables
    integer(kind=kint) :: nnode, ndof
    real(kind=kreal)   :: tangent(3,2)    
    real(kind=kreal)   :: tm(3, 3*(l_max_surface_node+1)) 
    real(kind=kreal)   :: tt(3, 3*(l_max_surface_node+1)) 
 
    ! Initialize arrays to zero for safety
    bn(:) = 0.0d0
    metric(:,:) = 0.0d0
    ht(:,:) = 0.0d0
    gt(:,:) = 0.0d0
 
    nnode = size(tsurf%nodes)
    ndof = nnode*3 + 3
 
    ! Call computeTm_Tt as the core routine (fcoeff=0 since we don't need Tt for ALag)
    ! Bn is computed inside computeTm_Tt as Tm^T * normal
    call computetm_tt(ctstate, tsurf, 0.0d0, tm, tt, smoothing_type, bn=bn)
 
    ! Call DispIncreMatrix to get tangent displacement map and metric
    ! DispIncreMatrix returns: tangent basis, metric tensor, and displacement increment matrix
    call dispincrematrix(ctstate%lpos(1:2), tsurf%etype, nnode, ele(1:3,1:nnode), &
                         tangent, metric, ht(1:2,1:ndof))
 
    ! Compute covariant displacement map: Gt = metric * Ht
    ! This gives: dxy = Gt * edisp = metric * (Ht * edisp)
    gt(1:2,1:ndof) = matmul(metric(1:2,1:2), ht(1:2,1:ndof))
 
  end subroutine computecontactmaps_alag
 
  subroutine computefrictionforce_alag(ctState, fcoeff, lambda_n, metric, Ht, Gt, edisp, edof, ctTForce, &
                                        mut, update_multiplier, slave_id, ctchanged, &
                                        norm_trial, alpha, that, jump_ratio)
    implicit none
 
    type(tcontactstate), intent(inout) :: ctstate
    real(kind=kreal), intent(in)       :: fcoeff        
    real(kind=kreal), intent(in)       :: lambda_n      
    real(kind=kreal), intent(in)       :: metric(2,2)   
    real(kind=kreal), intent(in)       :: ht(:,:)       
    real(kind=kreal), intent(in)       :: gt(:,:)       
    real(kind=kreal), intent(in)       :: edisp(:)      
    integer(kind=kint), intent(in)     :: edof
    real(kind=kreal), intent(out)      :: cttforce(:)   
    real(kind=kreal), intent(in)       :: mut           
    logical, optional, intent(in)      :: update_multiplier
    integer(kind=kint), optional, intent(in) :: slave_id
    logical, optional, intent(inout)   :: ctchanged
    real(kind=kreal), optional, intent(out) :: norm_trial 
    real(kind=kreal), optional, intent(out) :: alpha      
    real(kind=kreal), optional, intent(out) :: that(2)    
    real(kind=kreal), optional, intent(out) :: jump_ratio 
 
    real(kind=kreal) :: dxy(2)       
    real(kind=kreal) :: fric(2)      
    real(kind=kreal) :: f3(edof)     
    real(kind=kreal) :: invmetric(2,2) 
    real(kind=kreal) :: det          
    real(kind=kreal) :: norm_fric    
    real(kind=kreal) :: tmp(2)       
    real(kind=kreal) :: alpha_local  
    logical          :: do_update
 
    cttforce = 0.0d0
    do_update = .false.
    if(present(jump_ratio)) jump_ratio = 0.0d0
    if(present(update_multiplier)) do_update = update_multiplier
 
    ! Compute inverse of metric tensor (2x2)
    det = metric(1,1)*metric(2,2) - metric(1,2)*metric(2,1)
    if( abs(det) < 1.0d-20 ) then
      ! Degenerate metric: set friction to zero
      cttforce = 0.0d0
      return
    endif
    invmetric(1,1) = metric(2,2) / det
    invmetric(2,2) = metric(1,1) / det
    invmetric(1,2) = -metric(1,2) / det
    invmetric(2,1) = -metric(2,1) / det
 
    ! Compute tangent displacement using Gt map: dxy = Gt * edisp
    dxy(1:2) = matmul( gt(1:2,1:edof), edisp(1:edof) )
 
    ! Trial friction force: fric = λ_t + ρ_t * dxy
    fric(1:2) = ctstate%multiplier(2:3) + mut*dxy(1:2)
 
    ! Compute 2D local norm: ||fric||_g = sqrt(fric^T * inv(metric) * fric)
    tmp(1:2) = matmul(invmetric(1:2,1:2), fric(1:2))
    norm_fric = dsqrt( dot_product(fric(1:2), tmp(1:2)) )
 
    ! Return projection information for stiffness (optional outputs)
    if(present(norm_trial)) norm_trial = norm_fric
    if(present(that) .and. norm_fric > 1.0d-20) then
      that(1:2) = fric(1:2) / norm_fric
    else if(present(that)) then
      that(1:2) = 0.0d0
    endif
 
    ! Compute projection ratio: alpha = min(1, r/norm_trial)
    if( lambda_n > 0.0d0 .and. norm_fric > 1.0d-20 ) then
      alpha_local = min(1.0d0, fcoeff*lambda_n / norm_fric)
    else
      alpha_local = 0.0d0
    endif
    if(present(alpha)) alpha = alpha_local
 
    ! Check stick/slip state using projection ratio
    if( do_update .and. lambda_n > 0.0d0 ) then
      ! Update mode: check state and update multiplier
      if( alpha_local < 1.0d0 ) then
        ! Slip state: projected to friction cone
        if( ctstate%state == contactstick ) then
          if(present(ctchanged)) ctchanged = .true.
          if(present(slave_id)) print *, "Node", slave_id, "to slip state", norm_fric, fcoeff*lambda_n
          if(present(jump_ratio) .and. lambda_n > 1.0e-20 .and. norm_fric > 1.0e-20) &
            &  jump_ratio = norm_fric / (fcoeff*lambda_n)
        endif
        ctstate%state = contactslip
        fric(1:2) = fric(1:2) * alpha_local
      else
        ! Stick state
        if( ctstate%state == contactslip ) then
          if(present(ctchanged)) ctchanged = .true.
          if(present(slave_id)) print *, "Node", slave_id, "to stick state", norm_fric, fcoeff*lambda_n
        endif
        ctstate%state = contactstick
      endif
      ! Update multiplier
      ctstate%multiplier(2:3) = fric(1:2)
    else if( lambda_n <= 0.0d0 ) then
      ! No normal force: no friction
      fric(1:2) = 0.0d0
    else
      ! Read-only mode: just scale if alpha < 1
      if( alpha_local < 1.0d0 ) then
        fric(1:2) = fric(1:2) * alpha_local
      endif
    endif
 
    ! Compute friction force vector using Ht map: f3 = Ht^T * fric
    f3(1:edof) = matmul( transpose(ht(1:2,1:edof)), fric(1:2) )
 
    ! Distribute friction force (negative for equilibrium)
    cttforce(1:edof) = -f3(1:edof)
 
  end subroutine computefrictionforce_alag
 
  subroutine gettrialfricforceandcheckfricstate(ctstate,fcoeff,tPenalty,lagrange)
 
    type(tcontactstate) :: ctstate
    real(kind=kreal)   :: fcoeff, tpenalty 
    real(kind=kreal)   :: lagrange 
 
    integer(kind=kint) :: i
    real(kind=kreal)   :: tf_yield
 
    do i = 1, 3
      ctstate%tangentForce_trial(i) = ctstate%tangentForce1(i) + tpenalty*ctstate%reldisp(i)
    enddo
 
    tf_yield = fcoeff*dabs(lagrange)
    if(ctstate%state == contactslip) tf_yield =0.99d0*tf_yield
    if( dsqrt(dot_product(ctstate%tangentForce_trial,ctstate%tangentForce_trial)) <= tf_yield ) then
      ctstate%state = contactstick
    else
      ctstate%state = contactslip
    endif
 
  end subroutine gettrialfricforceandcheckfricstate
 
  subroutine update_tangentforce(etype,nn,elemt0,elemt,cstate)
    integer, intent(in)                 :: etype
    integer, intent(in)                 :: nn
    real(kind=kreal),intent(in)         :: elemt0(3,nn)  
    real(kind=kreal),intent(in)         :: elemt(3,nn)   
    type(tcontactstate), intent(inout)  :: cstate
 
    integer           ::  i
    real(kind=kreal)  ::  tangent0(3,2), tangent(3,2)    ! base vectors in tangent space
    real(kind=kreal)  ::  coeff(2), norm, norm_tmp
    real(kind=kreal)  ::  tangentforce_tmp(3)
 
    call tangentbase( etype, nn, cstate%lpos(1:2), elemt0, tangent0 )
    call tangentbase( etype, nn, cstate%lpos(1:2), elemt, tangent )
 
    !project tangentforce to base vector tangent0
    do i=1,2
      coeff(i) = dot_product(cstate%tangentForce(1:3),tangent0(1:3,i))
      coeff(i) = coeff(i)/dot_product(tangent0(1:3,i),tangent0(1:3,i))
    enddo
    tangentforce_tmp(1:3) = coeff(1)*tangent0(1:3,1) + coeff(2)*tangent0(1:3,2)
    norm_tmp = dsqrt(dot_product(tangentforce_tmp,tangentforce_tmp))
    !adjust tangent force of slave point which moved over element boundary
    if( norm_tmp > 1.d-6 ) then
      norm = dsqrt(dot_product(cstate%tangentForce,cstate%tangentForce))
      coeff(1:2) = (norm/norm_tmp)*coeff(1:2)
    end if
 
    !set rotated tangentforce to tangentforce1
    cstate%tangentForce1(1:3) = coeff(1)*tangent(1:3,1) + coeff(2)*tangent(1:3,2)
 
  end subroutine update_tangentforce
 
end module m_fstr_contact_elem_common