fstr_contact_elem_alag.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_alag
  use hecmw
  use elementinfo
  use mcontactdef
  use msurfelement
  use m_fstr_contact_geom
  use m_fstr_contact_elem_common
  use m_fstr_contact_smoothing
  implicit none
 
  public :: getcontactstiffness_alag
  public :: getcontactnodalforce_alag
  public :: gettiedstiffness_alag
  public :: gettiednodalforce_alag
  public :: updatecontactmultiplier_alag
 
contains
 
  subroutine getcontactstiffness_alag(cstate, tSurf, ele, mu, mut, fcoeff, symm, stiff, force, smoothing_type, edisp, iter)
 
    type(tcontactstate), intent(inout) :: cstate
    type(tsurfelement), intent(in)  :: tsurf
    real(kind=kreal), intent(in)    :: ele(:,:)        
    real(kind=kreal), intent(in)    :: mu, mut         
    real(kind=kreal), intent(in)    :: fcoeff          
    logical, intent(in)             :: symm
    real(kind=kreal), intent(out)   :: stiff(:,:)      
    real(kind=kreal), intent(out)   :: force(:)        
    integer(kind=kint), optional, intent(in) :: smoothing_type
    real(kind=kreal), optional, intent(in) :: edisp(:)  
    integer(kind=kint), intent(in) :: iter
 
    integer          :: i, j, nnode
    real(kind=kreal) :: bn(size(tsurf%nodes)*3+3), ht(2,size(tsurf%nodes)*3+3), gt(2,size(tsurf%nodes)*3+3)
    real(kind=kreal) :: metric(2,2)
    real(kind=kreal) :: a(2,2)       
    real(kind=kreal) :: alpha_proj, that_dir(2)
    real(kind=kreal) :: k_fric(size(tsurf%nodes)*3+3,size(tsurf%nodes)*3+3)  
    real(kind=kreal) :: tmp_vec(2)
    real(kind=kreal) :: dummy_force(size(tsurf%nodes)*3+3)  
    real(kind=kreal) :: eval_disp(size(tsurf%nodes)*3+3)   
 
    nnode = size(tsurf%nodes)
 
    ! Use common mapping routine to compute Bn, metric, Ht, Gt
    call computecontactmaps_alag(cstate, tsurf, ele, bn, metric, ht, gt, smoothing_type)
 
    ! Normal stiffness: stiff = mu * Bn * Bn^T
    do j = 1, nnode*3+3
      do i = 1, nnode*3+3
        stiff(i,j) = mu * bn(i) * bn(j)
      enddo
    enddo
    force(1:nnode*3+3) = bn(:)
 
    ! frictional component
    if( fcoeff /= 0.d0 ) then
      ! Evaluate trial friction at current displacement for consistent tangent
      if( present(edisp) ) then
        eval_disp(1:nnode*3+3) = edisp(1:nnode*3+3)
      else
        eval_disp = 0.0d0
      endif
      call computefrictionforce_alag(cstate, fcoeff, cstate%multiplier(1), metric, &
                                      ht, gt, eval_disp, nnode*3+3, dummy_force, &
                                      mut, alpha=alpha_proj, that=that_dir)
 
      ! Friction tangent operator A in 2D metric space:  K_fric = Ht^T * A * Ht
      if( cstate%multiplier(1) <= 0.0d0 .or. alpha_proj <= 1.0d-20 ) then
        ! No normal contact force: no friction contribution
        a = 0.0d0
      else if( alpha_proj >= 0.999d0 ) then
        ! Stick: A = mu_t * M (exact)
        a(1,1) = mut * metric(1,1)
        a(1,2) = mut * metric(1,2)
        a(2,1) = mut * metric(2,1)
        a(2,2) = mut * metric(2,2)
      else
        ! Slip: switch tangent by NR iteration count for stability.
        !   iter <= 2: A = alpha*mu_t*M  (stable, no directional correction)
        !   iter >= 3: A = alpha*mu_t*(M - t_hat x t_hat)  (consistent tangent)
        ! The first two NR steps use M to let t_hat stabilize; after that the consistent tangent is used to regain quadratic convergence.
        if( iter <= 2 ) then
          a(1,1) = alpha_proj * mut * metric(1,1)
          a(1,2) = alpha_proj * mut * metric(1,2)
          a(2,1) = alpha_proj * mut * metric(2,1)
          a(2,2) = alpha_proj * mut * metric(2,2)
        else
          a(1,1) = alpha_proj * mut * (metric(1,1) - that_dir(1)*that_dir(1))
          a(1,2) = alpha_proj * mut * (metric(1,2) - that_dir(1)*that_dir(2))
          a(2,1) = alpha_proj * mut * (metric(2,1) - that_dir(2)*that_dir(1))
          a(2,2) = alpha_proj * mut * (metric(2,2) - that_dir(2)*that_dir(2))
        endif
      endif
 
      ! Compute friction stiffness: K_fric = Ht^T * A * Ht (consistent tangent)
      do j = 1, nnode*3+3
        tmp_vec = matmul(a, ht(1:2,j))
        do i = 1, nnode*3+3
          k_fric(i,j) = dot_product(ht(1:2,i), tmp_vec)
        enddo
      enddo
 
      stiff(1:nnode*3+3,1:nnode*3+3) = stiff(1:nnode*3+3,1:nnode*3+3) + k_fric(1:nnode*3+3,1:nnode*3+3)
    endif
 
  end subroutine getcontactstiffness_alag
 
  subroutine getcontactnodalforce_alag(ctState,tSurf,ndCoord,ndDu,mu,mut,fcoeff,lagrange,ctNForce,ctTForce,cflag,smoothing_type)
 
    use msurfelement
    type(tcontactstate) :: ctstate
    type(tsurfelement)  :: tsurf
    integer(kind=kint) :: nnode
    integer(kind=kint) :: j
    real(kind=kreal), intent(in) :: mu, mut 
    real(kind=kreal)   :: fcoeff 
    real(kind=kreal)   :: lagrange 
    real(kind=kreal)   :: ndcoord(:), nddu(:) 
    real(kind=kreal)   :: ctnforce(:) 
    real(kind=kreal)   :: cttforce(:) 
    logical            :: cflag
    integer(kind=kint), optional, intent(in) :: smoothing_type
 
    real(kind=kreal)   :: normal(3) 
    real(kind=kreal)   :: bn(3*l_max_elem_node+3) 
    real(kind=kreal)   :: ht(2,3*l_max_elem_node+3), gt(2,3*l_max_elem_node+3) 
    real(kind=kreal)   :: elemcrd(3, l_max_elem_node) 
    real(kind=kreal)   :: edisp(3*l_max_elem_node+3) 
    real(kind=kreal)   :: dgn, nrlforce 
    real(kind=kreal)   :: metric(2,2)
    integer(kind=kint) :: edof
 
    nnode = size(tsurf%nodes)
    edof = nnode*3+3
 
    ctnforce = 0.0d0
    cttforce = 0.0d0
 
    normal(1:3) = ctstate%direction(1:3)
 
    ! Prepare elemcrd = ndCoord - ndDu (i.e., coord + disp) for computeContactMaps_ALag
    do j = 1, nnode
      elemcrd(1:3, j) = ndcoord(j*3+1:j*3+3) - nddu(j*3+1:j*3+3)
    enddo
 
    ! Use common mapping routine to compute Bn, metric, Ht, Gt
    call computecontactmaps_alag(ctstate, tsurf, elemcrd(:,1:nnode), &
                                  bn, metric, ht, gt, smoothing_type)
 
    ! Normal gap: dgn = Bn^T * ndCoord (using normal distribution vector)
    dgn = dot_product( bn(1:edof), ndcoord(1:edof) )
 
    ! Normal force: multiplier + penalty * gap
    nrlforce = ctstate%multiplier(1) + mu*dgn
 
    ! Distribute normal force using Bn: ctNForce = -nrlforce * Bn
    ctnforce(1:edof) = -nrlforce * bn(1:edof)
 
    ! Lagrange row (not used in ALagrange, set to 0)
    ctnforce((nnode+1)*3+1) = 0.d0
 
    if( fcoeff == 0.d0 ) return
 
    ! --- Tangent component ---
 
    ! Prepare edisp from ndDu
    edisp(1:3) = nddu(1:3)  ! slave
    do j = 1, nnode
      edisp(j*3+1:j*3+3) = nddu(j*3+1:j*3+3)  ! master nodes
    enddo
 
    ! Compute friction force using common routine
    call computefrictionforce_alag(ctstate, fcoeff, ctstate%multiplier(1), metric, &
                                    ht, gt, edisp, edof, cttforce, &
                                    mut)
 
    ! Lagrange row (not used in ALagrange, set to 0)
    cttforce((nnode+1)*3+1) = 0.d0
 
  end subroutine getcontactnodalforce_alag
 
  subroutine updatecontactmultiplier_alag(ctState,ndLocal,coord,disp,ddisp,&
     &  mu,mut,fcoeff,tSurf,lgnt,ctchanged,ctNForce,ctTForce,jump_ratio,smoothing_type)
 
    type(tcontactstate), intent(inout)   :: ctstate
    integer(kind=kint), intent(in)       :: ndlocal(:)
    real(kind=kreal), intent(in)         :: coord(:)            
    real(kind=kreal), intent(in)         :: disp(:)             
    real(kind=kreal), intent(in)         :: ddisp(:)            
    real(kind=kreal), intent(in)         :: mu, mut             
    real(kind=kreal), intent(in)         :: fcoeff              
    type(tsurfelement), intent(in)       :: tsurf
    real(kind=kreal), intent(inout)      :: lgnt(2)             
    logical, intent(inout)               :: ctchanged
    real(kind=kreal), intent(out)        :: ctnforce(:)         
    real(kind=kreal), intent(out)        :: cttforce(:)         
    real(kind=kreal), intent(out)        :: jump_ratio          
    integer(kind=kint), optional, intent(in) :: smoothing_type
 
    integer(kind=kint)  :: nnode
    integer(kind=kint)  :: slave, j
    real(kind=kreal)    :: bn(3*l_max_elem_node+3) 
    real(kind=kreal)    :: ht(2,3*l_max_elem_node+3), gt(2,3*l_max_elem_node+3) 
    real(kind=kreal)    :: elemcrd(3,l_max_elem_node) 
    real(kind=kreal)    :: curpos(3*l_max_elem_node+3) 
    real(kind=kreal)    :: edisp(3*l_max_elem_node+3) 
    real(kind=kreal)    :: dgn, nrlforce 
    real(kind=kreal)    :: metric(2,2)
    real(kind=kreal)    :: dxy(2)        
    integer(kind=kint)  :: edof
 
    nnode = size(ndlocal) - 1
    slave = ndlocal(1)
    edof = nnode*3+3
 
    ctnforce = 0.0d0
    cttforce = 0.0d0
 
    ! Prepare elemcrd (coord+disp) and current positions (coord+disp+ddisp)
    curpos(1:3) = coord(3*slave-2:3*slave) + disp(3*slave-2:3*slave) + ddisp(3*slave-2:3*slave)
    edisp(1:3) = ddisp(3*slave-2:3*slave)
    do j = 1, nnode
      elemcrd(1:3,j) = coord(3*ndlocal(j+1)-2:3*ndlocal(j+1)) + disp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
      curpos(j*3+1:j*3+3) = elemcrd(1:3,j) + ddisp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
      edisp(j*3+1:j*3+3) = ddisp(3*ndlocal(j+1)-2:3*ndlocal(j+1))
    enddo
 
    ! Use common mapping routine to compute Bn, metric, Ht, Gt
    call computecontactmaps_alag(ctstate, tsurf, elemcrd(:,1:nnode), &
                                  bn, metric, ht, gt, smoothing_type)
 
    ! Normal gap: dgn = Bn^T * curpos (using normal distribution vector)
    dgn = dot_product( bn(1:edof), curpos(1:edof) )
 
    ! Update multiplier and working distance
    ctstate%wkdist = -dgn
    ctstate%multiplier(1) = ctstate%multiplier(1) - mu*ctstate%wkdist
    ctstate%distance = ctstate%wkdist
    lgnt(1) = lgnt(1) - ctstate%wkdist
 
    ! Normal force: use updated multiplier
    nrlforce = ctstate%multiplier(1)
 
    ! Distribute normal force using Bn: ctNForce = -nrlforce * Bn
    ctnforce(1:edof) = -nrlforce * bn(1:edof)
 
    if( fcoeff == 0.d0 ) return
 
    ! --- Tangent component ---
 
    ! Compute friction force with multiplier update and state check
    call computefrictionforce_alag(ctstate, fcoeff, ctstate%multiplier(1), metric, &
                                    ht, gt, edisp, edof, cttforce, &
                                    mut, &
                                    update_multiplier=.true., slave_id=slave, ctchanged=ctchanged, &
                                    jump_ratio=jump_ratio)
 
    ! Tangent displacement for convergence check: use Gt to project curpos directly
    dxy = matmul( gt(:,1:edof), curpos(1:edof) )
    lgnt(2) = lgnt(2) + dsqrt( dxy(1)*dxy(1) + dxy(2)*dxy(2) )
 
  end subroutine updatecontactmultiplier_alag
 
  subroutine gettiedstiffness_alag(cstate, tSurf, mu, stiff, force)
 
    type(tcontactstate), intent(in) :: cstate
    type(tsurfelement), intent(in)  :: tsurf
    real(kind=kreal), intent(in)    :: mu              
    real(kind=kreal), intent(out)   :: stiff(:,:)      
    real(kind=kreal), intent(out)   :: force(:)        
 
    integer          :: i, j, nnode, edof
    real(kind=kreal) :: tm(3, 3*(l_max_surface_node+1))
    real(kind=kreal) :: tt(3, 3*(l_max_surface_node+1))  
 
    nnode = size(tsurf%nodes)
    edof = nnode*3+3
 
    stiff = 0.d0
 
    ! Use common mapping routine to compute Tm
    call computetm_tt(cstate, tsurf, 0.0d0, tm, tt)
 
    ! Tied stiffness: stiff = mu * Tm^T * Tm
    do j = 1, edof
      do i = 1, edof
        stiff(i,j) = mu * dot_product(tm(1:3,i), tm(1:3,j))
      enddo
    enddo
    force(1:edof) = 0.d0  ! not used for tied (3-direction constraint)
 
  end subroutine gettiedstiffness_alag
 
  subroutine gettiednodalforce_alag(ctState,tSurf,ndu,mu,ctNForce,ctTForce)
 
    use msurfelement
    type(tcontactstate) :: ctstate
    type(tsurfelement)  :: tsurf
    integer(kind=kint)  :: nnode
    real(kind=kreal)   :: ndu(:) 
    real(kind=kreal), intent(in) :: mu 
    real(kind=kreal)   :: ctnforce(:)  
    real(kind=kreal)   :: cttforce(:)  
 
    integer(kind=kint) :: edof
    real(kind=kreal)   :: tm(3, 3*(l_max_surface_node+1))  
    real(kind=kreal)   :: tt(3, 3*(l_max_surface_node+1))  
    real(kind=kreal)   :: dg(3) 
    real(kind=kreal)   :: nrlforce(3) 
 
    nnode = size(tsurf%nodes)
    edof = nnode*3+3
 
    ctnforce = 0.0d0
    cttforce = 0.0d0
 
    ! Use common mapping routine to compute Tm
    call computetm_tt(ctstate, tsurf, 0.0d0, tm, tt)
 
    ! Gap vector: dg = Tm * ndu (3-component relative displacement)
    dg(1:3) = matmul(tm(1:3, 1:edof), ndu(1:edof))
 
    ! Force: multiplier + penalty * gap (3 components)
    nrlforce(1:3) = ctstate%multiplier(1:3) + mu*dg(1:3)
 
    ! Distribute force: ctNForce = -Tm^T * nrlforce
    ctnforce(1:edof) = -matmul(transpose(tm(1:3, 1:edof)), nrlforce(1:3))
 
    ! Lagrange row (not used in ALagrange, set to 0)
    ctnforce((nnode+1)*3+1) = 0.d0
    cttforce((nnode+1)*3+1) = 0.d0
 
  end subroutine gettiednodalforce_alag
 
end module m_fstr_contact_elem_alag