fstr_contact_elem_slag.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_slag
  use hecmw
  use elementinfo
  use mcontactdef
  use msurfelement
  use m_fstr_contact_geom
  use m_fstr_contact_elem_common
  implicit none
 
  public :: getcontactstiffness_slag
  public :: getcontactnodalforce_slag
  public :: gettiedstiffness_slag
  public :: gettiednodalforce_slag
 
contains
 
  subroutine getcontactstiffness_slag(ctState,tSurf,iter,tPenalty,fcoeff,lagrange,stiffness)
 
    use msurfelement
    type(tcontactstate) :: ctstate
    type(tsurfelement)  :: tsurf
    integer(kind=kint)  :: iter
    integer(kind=kint)  :: nnode
    integer(kind=kint)  :: i, j
    real(kind=kreal)    :: normal(3)
    real(kind=kreal)    :: ntm((l_max_surface_node + 1)*3) 
    real(kind=kreal)    :: fcoeff, tpenalty 
    real(kind=kreal)    :: lagrange 
    real(kind=kreal)    :: tf_trial(3), length_tft
    real(kind=kreal)    :: tangent(3)
    real(kind=kreal)    :: stiffness(:,:) 
    real(kind=kreal)    :: tm(3, 3*(l_max_surface_node+1)), tt(3, 3*(l_max_surface_node+1)) 
    real(kind=kreal)    :: q((l_max_surface_node+1)*3)  
    real(kind=kreal)    :: dot_tn  
 
    nnode = size(tsurf%nodes)
 
    stiffness = 0.0d0
 
    ! Compute Tm and Tt matrices using standard shape functions
    call computetm_tt(ctstate, tsurf, fcoeff, tm, tt)
 
    normal(1:3) = ctstate%direction(1:3)
 
    ! Compute nTm = Tm^T * n for KKT constraint (Lagrange multiplier row/column)
    ntm(1:(nnode+1)*3) = matmul(transpose(tm(1:3, 1:3*(nnode+1))), normal)
 
    i = (nnode+1)*3 + 1
    do j = 1, (nnode+1)*3
      stiffness(i,j) = ntm(j);  stiffness(j,i) = ntm(j)
    enddo
 
    if( fcoeff /= 0.0d0 ) then
      if( lagrange>0.0d0 .or. iter==1 ) then
 
        ! Stick friction: K_uu^stick = ρ * Tt^T * Tt
        do i = 1, (nnode+1)*3
          do j = 1, i
            stiffness(i,j) = stiffness(i,j) + tpenalty * dot_product(tt(1:3,i), tt(1:3,j))
            ! Fill symmetric part
            if (i /= j) then
              stiffness(j,i) = stiffness(i,j)
            endif
          enddo
        enddo
 
        if( ctstate%state == contactslip ) then
 
          tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
          length_tft = dsqrt(dot_product(tf_trial,tf_trial))
          tangent(1:3) = tf_trial(1:3)/length_tft
 
          ! Compute q = Tm^T * tangent
          q(1:(nnode+1)*3) = matmul(transpose(tm(1:3, 1:3*(nnode+1))), tangent)
 
          ! Slip correction:
          dot_tn = dot_product(tangent, normal)
 
          do i = 1, (nnode+1)*3
            do j = 1, i
              stiffness(i,j) = stiffness(i,j) + tpenalty * (-q(i)*q(j) + q(i)*ntm(j)*dot_tn)
              if (i /= j) then
                stiffness(j,i) = stiffness(i,j)
              endif
            enddo
          enddo
 
          ! Scale (u,u) block by (μλ/|t_tr|)
          stiffness(1:(nnode+1)*3,1:(nnode+1)*3) = (fcoeff*lagrange/length_tft) * stiffness(1:(nnode+1)*3,1:(nnode+1)*3)
 
          ! (u,λ) cross term:)
          stiffness(1:(nnode+1)*3, (nnode+1)*3+1) = stiffness(1:(nnode+1)*3, (nnode+1)*3+1) + fcoeff * q(1:(nnode+1)*3)
 
        endif
      endif
    endif
 
  end subroutine getcontactstiffness_slag
 
  subroutine getcontactnodalforce_slag(ctState,tSurf,ndCoord,ndDu,tPenalty,fcoeff,lagrange,ctNForce,ctTForce,cflag)
 
    use msurfelement
    type(tcontactstate) :: ctstate
    type(tsurfelement)  :: tsurf
    integer(kind=kint) :: nnode
    integer(kind=kint) :: j
    real(kind=kreal)   :: fcoeff, tpenalty 
    real(kind=kreal)   :: lagrange 
    real(kind=kreal)   :: ndcoord(:), nddu(:) 
    real(kind=kreal)   :: normal(3) 
    real(kind=kreal)   :: ntm((l_max_surface_node + 1)*3) 
    real(kind=kreal)   :: tf_trial(3), length_tft, tangent(3), tf_final(3)
    real(kind=kreal)       :: ctnforce(:)                     
    real(kind=kreal)       :: cttforce(:)                     
    logical            :: cflag
 
    real(kind=kreal)   :: tm(3, 3*(l_max_surface_node+1)), tt(3, 3*(l_max_surface_node+1)) 
 
    ctstate%multiplier(1) = lagrange
 
    nnode = size(tsurf%nodes)
 
    ctnforce = 0.0d0
    cttforce = 0.0d0
 
    ! Compute Tm matrix (Tt not needed, pass fcoeff=0 to skip)
    call computetm_tt(ctstate, tsurf, fcoeff, tm, tt)
 
    normal(1:3) = ctstate%direction(1:3)
 
    ! Compute nTm = Tm^T * normal
    ntm(1:(nnode+1)*3) = -matmul(transpose(tm(1:3, 1:3*(nnode+1))), normal)
 
    do j = 1, (nnode+1)*3
      ctnforce(j) = lagrange*ntm(j)
    enddo
    j = (nnode+1)*3 + 1
    ctnforce(j) = dot_product(ntm(1:(nnode+1)*3),ndcoord(1:(nnode+1)*3))
 
    if(fcoeff /= 0.0d0 .and. lagrange > 0.0d0)then
 
      if( cflag ) then !calc tf_final and set it to tangentForce_final
        ctstate%reldisp(1:3) = matmul(tt(1:3,1:3*(nnode+1)), nddu(1:3*(nnode+1)))
 
        call gettrialfricforceandcheckfricstate(ctstate,fcoeff,tpenalty,lagrange)
 
        if( ctstate%state == contactstick ) then
          tf_final(1:3) = ctstate%tangentForce_trial(1:3)
        elseif( ctstate%state == contactslip ) then
          tf_trial(1:3) = ctstate%tangentForce_trial(1:3)
          length_tft = dsqrt(dot_product(tf_trial,tf_trial))
          tangent(1:3) = tf_trial(1:3)/length_tft
          tf_final(1:3) = fcoeff*dabs(lagrange)*tangent(1:3)
        endif
        ctstate%tangentForce_final(1:3) = tf_final(1:3)
      else ! just set tangentForce_final to tf_final (used for fstr_ass_load_contact)
        tf_final(1:3) = ctstate%tangentForce_final(1:3)
      endif
 
      ! Distribute friction force using Tm: ctTForce = -Tm^T * tf_final
      cttforce(1:(nnode+1)*3) = -matmul(transpose(tm(1:3, 1:3*(nnode+1))), tf_final)
 
    endif
 
  end subroutine getcontactnodalforce_slag
 
  subroutine gettiedstiffness_slag(ctState,tSurf,idof,stiffness)
    use msurfelement
    type(tcontactstate) :: ctstate
    type(tsurfelement)  :: tsurf
    integer(kind=kint)  :: nnode
    integer(kind=kint)  :: idof
    real(kind=kreal)    :: stiffness(:,:) 
 
    integer(kind=kint)  :: i, j
    real(kind=kreal)    :: ntm((l_max_surface_node+1)*3)
    real(kind=kreal)    :: tm(3, 3*(l_max_surface_node+1)), tt(3, 3*(l_max_surface_node+1)) 
    real(kind=kreal)    :: e_idof(3)  
 
    nnode = size(tsurf%nodes)
    stiffness = 0.0d0
 
    ! Compute Tm matrix
    call computetm_tt(ctstate, tsurf, 0.0d0, tm, tt)
 
    ! Create unit vector in idof direction
    e_idof = 0.0d0
    e_idof(idof) = 1.0d0
 
    ! Compute nTm = Tm^T * e_idof
    ! This gives constraint direction for tied contact
    ntm(1:(nnode+1)*3) = matmul(transpose(tm(1:3, 1:3*(nnode+1))), e_idof)
 
    i = (nnode+1)*3 + 1
    do j = 1, (nnode+1)*3
      stiffness(i,j) = ntm(j);  stiffness(j,i) = ntm(j)
    enddo
 
  end subroutine gettiedstiffness_slag
 
  subroutine gettiednodalforce_slag(ctState,tSurf,idof,ndu,lagrange,ctNForce,ctTForce)
    use msurfelement
    type(tcontactstate) :: ctstate
    type(tsurfelement)  :: tsurf
    integer(kind=kint)  :: nnode
    integer(kind=kint)  :: idof
    real(kind=kreal)   :: ndu(:) 
    real(kind=kreal)   :: lagrange 
    real(kind=kreal)       :: ctnforce(:)  
    real(kind=kreal)       :: cttforce(:)  
 
    integer(kind=kint) :: j
    real(kind=kreal)   :: normal(3) 
    real(kind=kreal)   :: ntm((l_max_surface_node + 1)*3) 
    real(kind=kreal)   :: ttm((l_max_surface_node + 1)*3) 
    real(kind=kreal)   :: tm(3, 3*(l_max_surface_node+1)), tt(3, 3*(l_max_surface_node+1)) 
    real(kind=kreal)   :: e_idof(3), e_normal(3), e_tangent(3)  
 
    nnode = size(tsurf%nodes)
 
    ctnforce = 0.0d0
    cttforce = 0.0d0
 
    ! Compute Tm matrix
    call computetm_tt(ctstate, tsurf, 0.0d0, tm, tt)
 
    normal(1:3) = ctstate%direction(1:3)
 
    ! Create unit vector in idof direction
    e_idof = 0.0d0
    e_idof(idof) = 1.0d0
 
    ! Normal component: -normal(idof) * normal (OLD: nTm(1:3) = -normal(idof)*normal(1:3))
    e_normal = -normal(idof) * normal
 
    ! Tangential component: -e_idof - e_normal (OLD: tTm(idof)=-1.0; tTm(1:3)-=nTm(1:3))
    e_tangent = -e_idof - e_normal
 
    ! Compute nTm and tTm using Tm^T
    ntm(1:(nnode+1)*3) = matmul(transpose(tm(1:3, 1:3*(nnode+1))), e_normal)
    ttm(1:(nnode+1)*3) = matmul(transpose(tm(1:3, 1:3*(nnode+1))), e_tangent)
 
    do j = 1, (nnode+1)*3
      ctnforce(j) = lagrange*ntm(j)
      cttforce(j) = lagrange*ttm(j)
    enddo
    j = (nnode+1)*3 + 1
    ctnforce(j) = dot_product(ntm(1:(nnode+1)*3),ndu(1:(nnode+1)*3))
    cttforce(j) = dot_product(ttm(1:(nnode+1)*3),ndu(1:(nnode+1)*3))
  end subroutine gettiednodalforce_slag
 
end module m_fstr_contact_elem_slag