fstr_contact.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 mcontact
 
  use mcontactdef
  use hecmw
  use m_fstr
  use m_fstr_contact_assembly
  use m_fstr_contact_element
  use m_fstr_contact_geom
  use m_fstr_contact_search
  use m_fstr_contact_mpc
  use m_fstr_contact_output
  implicit none
 
  logical, private :: active
 
  ! Re-export functions from sub-modules
  public :: initialize_contact_output_vectors
  public :: initialize_embed_vectors
  public :: setup_contact_elesurf_for_area
  public :: calc_contact_area
  public :: fstr_setup_parancon_contactvalue
  public :: update_contact_state_vectors
  public :: fstr_update_contact_state_vectors
 
contains
 
  subroutine fstr_addcontactstiffness(cstep,ctAlgo, iter,hecMESH,conMAT,hecLagMAT,fstrSOLID)
 
    integer(kind=kint)                   :: cstep
    integer(kind=kint)                   :: ctAlgo
    integer(kind=kint)                   :: iter
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: conMAT
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    integer(kind=kint)                   :: i, grpid
 
    if( associated(heclagmat%AL_lagrange) ) heclagmat%AL_lagrange = 0.0d0
    if( associated(heclagmat%AU_lagrange) ) heclagmat%AU_lagrange = 0.0d0
 
    do i = 1, fstrsolid%n_contacts
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      call calcu_contact_stiffness_nodesurf( ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
        fstrsolid%dunode(:), iter, heclagmat%Lagrange(:), conmat, heclagmat)
 
    enddo
 
    do i = 1, fstrsolid%n_embeds
 
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
 
      call calcu_contact_stiffness_nodesurf( ctalgo, fstrsolid%embeds(i), hecmesh%node(:), fstrsolid%unode(:), &
        fstrsolid%dunode(:), iter, heclagmat%Lagrange(:), conmat, heclagmat)
 
    enddo
 
  end subroutine fstr_addcontactstiffness
 
  subroutine fstr_update_ndforce_contact(cstep,ctAlgo,hecMESH,hecLagMAT,fstrSOLID,conMAT)
    integer(kind=kint), intent(in)       :: cstep
    integer(kind=kint), intent(in)       :: ctAlgo
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(hecmwst_matrix)                 :: conMAT
 
    call fstr_contact_ndforce_core(kctforresidual,cstep,ctalgo,hecmesh,heclagmat,fstrsolid,conmat)
 
  end subroutine fstr_update_ndforce_contact
 
  subroutine fstr_calc_contact_output_force(cstep,ctAlgo,hecMESH,hecLagMAT,fstrSOLID,conMAT)
    integer(kind=kint), intent(in)       :: cstep
    integer(kind=kint), intent(in)       :: ctAlgo
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(hecmwst_matrix)                 :: conMAT
 
    call fstr_contact_ndforce_core(kctforoutput,cstep,ctalgo,hecmesh,heclagmat,fstrsolid,conmat)
 
  end subroutine fstr_calc_contact_output_force
 
  subroutine fstr_contact_ndforce_core(purpose,cstep,ctAlgo,hecMESH,hecLagMAT,fstrSOLID,conMAT)
    integer(kind=kint), intent(in)       :: purpose
    integer(kind=kint), intent(in)       :: cstep
    integer(kind=kint), intent(in)       :: ctAlgo
    type(hecmwst_local_mesh)             :: hecMESH
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(hecmwst_matrix)                 :: conMAT
    integer(kind=kint) :: i, grpid
 
    ! Zero-clear output arrays only when computing for output
    if( purpose == kctforoutput ) then
      if( associated(fstrsolid%CONT_NFORCE) ) fstrsolid%CONT_NFORCE(:) = 0.d0
      if( associated(fstrsolid%CONT_FRIC) )   fstrsolid%CONT_FRIC(:) = 0.d0
      if( associated(fstrsolid%EMBED_NFORCE) ) fstrsolid%EMBED_NFORCE(:) = 0.d0
    endif
 
    do i = 1, fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
      call calcu_contact_ndforce_nodesurf( purpose, ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
        fstrsolid%dunode(:), heclagmat%Lagrange(:), conmat, &
        fstrsolid%CONT_NFORCE, fstrsolid%CONT_FRIC, heclagmat )
    enddo
 
    do i = 1, fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
      call calcu_contact_ndforce_nodesurf( purpose, ctalgo, fstrsolid%embeds(i), hecmesh%node(:), fstrsolid%unode(:), &
        fstrsolid%dunode(:), heclagmat%Lagrange(:), conmat, &
        fstrsolid%EMBED_NFORCE, fstrsolid%EMBED_NFORCE, heclagmat )
    enddo
 
  end subroutine fstr_contact_ndforce_core
 
  subroutine fstr_calc_contact_refstiff(cstep, hecMESH, hecMAT, fstrSOLID)
    use hecmw_matrix_misc
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: hecMAT
    type(fstr_solid)                     :: fstrSOLID
    
    integer(kind=kint) :: i, grpid, ndof
    real(kind=kreal), pointer :: diag(:)
    
    ndof = hecmat%NDOF
    
    ! Extract diagonal components once (efficient, called only once)
    diag => hecmw_mat_diag(hecmat)
    
    ! Loop over contact pairs
    do i = 1, fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if(.not. fstr_iscontactactive(fstrsolid, grpid, cstep)) cycle
      
      call calc_contact_pair_refstiff(fstrsolid%contacts(i), diag, ndof, hecmesh)
    enddo
    
    ! Same for embeds
    do i = 1, fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if(.not. fstr_isembedactive(fstrsolid, grpid, cstep)) cycle
      
      call calc_contact_pair_refstiff(fstrsolid%embeds(i), diag, ndof, hecmesh)
    enddo
    
    deallocate(diag)
  end subroutine fstr_calc_contact_refstiff
 
  subroutine fstr_scan_contact_state( cstep, sub_step, cont_step, dt, ctAlgo, hecMESH, fstrSOLID, infoCTChange )
    integer(kind=kint), intent(in)         :: cstep
    integer(kind=kint), intent(in)         :: sub_step
    integer(kind=kint), intent(in)         :: cont_step
    real(kind=kreal), intent(in)           :: dt
    integer(kind=kint), intent(in)         :: ctAlgo
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    type(fstr_info_contactchange), intent(inout):: infoCTChange
    character(len=9)                       :: flag_ctAlgo
    integer(kind=kint) :: i, grpid
    integer(kind=kint) :: s_f2c, s_c2f, s_emov, s_islid, s_act
    logical :: iactive, is_init
 
    if( associated( fstrsolid%CONT_RELVEL ) ) fstrsolid%CONT_RELVEL(:) = 0.d0
    if( associated( fstrsolid%CONT_STATE ) ) fstrsolid%CONT_STATE(:) = 0.d0
 
    if( ctalgo == kcaslagrange ) then
      flag_ctalgo = 'SLagrange'
    elseif( ctalgo == kcaalagrange ) then
      flag_ctalgo = 'ALagrange'
    endif
 
    ! P.A. We redefine fstrSOLID%ddunode as current coordinate of every nodes
    !  fstrSOLID%ddunode(:) = fstrSOLID%unode(:) + fstrSOLID%dunode(:)
    do i = 1, size(fstrsolid%unode)
      fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
    enddo
    active = .false.
 
    infoctchange%contact2free = 0
    infoctchange%contact2neighbor = 0
    infoctchange%contact2diffLpos = 0
    infoctchange%free2contact = 0
    infoctchange%contactNode_current = 0
 
    is_init = ( cstep == 1 .and. sub_step == 1 .and. cont_step == 0 )
 
    do i=1,fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) then
        call clear_contact_state(fstrsolid%contacts(i));  cycle
      endif
      call scan_contact_state( fstrsolid%contacts(i), &
      & fstrsolid%ddunode(:), fstrsolid%dunode(:), infoctchange, &
      & hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), &
      & is_init, iactive, hecmesh=hecmesh, &
      & flag_ctalgo=flag_ctalgo, ndforce=fstrsolid%QFORCE(:) )
      if( .not. active ) active = iactive
    enddo
 
    do i=1,fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) then
        call clear_contact_state(fstrsolid%embeds(i));  cycle
      endif
      call scan_embed_state( flag_ctalgo, fstrsolid%embeds(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
      & fstrsolid%QFORCE(:), infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive )
      if( .not. active ) active = iactive
    enddo
 
    if( is_init .and. ctalgo == kcaslagrange .and. fstrsolid%n_contacts > 0 ) &
    &  call remove_duplication_tiedcontact( cstep, hecmesh, fstrsolid, infoctchange )
 
    infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
    infoctchange%contactNode_previous = infoctchange%contactNode_current
 
    ! Output summary of contact state changes (always on; per-node detail requires CONTACT_LOG_LEVEL>=1)
    ! allreduce must be called by all ranks regardless of local counts
    s_f2c   = infoctchange%free2contact
    s_c2f   = infoctchange%contact2free
    s_emov  = infoctchange%contact2neighbor
    s_islid = infoctchange%contact2diffLpos
    s_act   = infoctchange%contactNode_current
    call hecmw_allreduce_i1(hecmesh, s_f2c,   hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_c2f,   hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_emov,  hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_islid, hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_act,   hecmw_sum)
    if (hecmw_comm_get_rank() == 0) then
      write(*,'(A,i0,A,i0,A,i0,A,i0,A,i0)') ' Contact change: Free2Cont=', s_f2c, ', Cont2Free=', s_c2f, &
        ', ElemMoved=', s_emov, ', InElemSlid=', s_islid, ', ActiveNodes=', s_act
    end if
 
    if( .not. active ) then
      if( associated( fstrsolid%CONT_NFORCE ) ) fstrsolid%CONT_NFORCE(:) = 0.d0
      if( associated( fstrsolid%CONT_FRIC ) ) fstrsolid%CONT_FRIC(:) = 0.d0
    end if
 
  end subroutine
 
  subroutine fstr_scan_contact_state_exp( 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
    integer(kind=kint) :: s_f2c, s_c2f, s_emov, s_islid, s_act
    logical :: iactive, is_init
 
 
    ! P.A. We redefine fstrSOLID%ddunode as current coordinate of every nodes
    !  fstrSOLID%ddunode(:) = fstrSOLID%unode(:) + fstrSOLID%dunode(:)
    do i = 1, size(fstrsolid%unode)
      fstrsolid%ddunode(i) = hecmesh%node(i) + fstrsolid%unode(i) + fstrsolid%dunode(i)
    enddo
    infoctchange%active = .false.
 
    infoctchange%contact2free = 0
    infoctchange%contact2neighbor = 0
    infoctchange%contact2diffLpos = 0
    infoctchange%free2contact = 0
    infoctchange%contactNode_current = 0
 
    is_init = ( cstep == 1 )
 
    do i=1,fstrsolid%n_contacts
      !   grpid = fstrSOLID%contacts(i)%group
      !   if( .not. fstr_isContactActive( fstrSOLID, grpid, cstep ) ) then
      !     call clear_contact_state(fstrSOLID%contacts(i));  cycle
      !   endif
 
      call scan_contact_state( fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
      & infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive, hecmesh )
 
      infoctchange%active = infoctchange%active .or. iactive
    enddo
 
    infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
    infoctchange%contactNode_previous = infoctchange%contactNode_current
 
    ! Output summary of contact state changes (always on; per-node detail requires CONTACT_LOG_LEVEL>=1)
    ! allreduce must be called by all ranks regardless of local counts
    s_f2c   = infoctchange%free2contact
    s_c2f   = infoctchange%contact2free
    s_emov  = infoctchange%contact2neighbor
    s_islid = infoctchange%contact2diffLpos
    s_act   = infoctchange%contactNode_current
    call hecmw_allreduce_i1(hecmesh, s_f2c,   hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_c2f,   hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_emov,  hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_islid, hecmw_sum)
    call hecmw_allreduce_i1(hecmesh, s_act,   hecmw_sum)
    if (hecmw_comm_get_rank() == 0) then
      write(*,'(A,i0,A,i0,A,i0,A,i0,A,i0)') ' Contact change: Free2Cont=', s_f2c, ', Cont2Free=', s_c2f, &
        ', ElemMoved=', s_emov, ', InElemSlid=', s_islid, ', ActiveNodes=', s_act
    end if
 
    fstrsolid%ddunode = 0.d0
  end subroutine
 
  logical function fstr_is_contact_active()
    fstr_is_contact_active = active
  end function
 
  subroutine fstr_set_contact_active( a )
    logical, intent(in) :: a
    active = a
  end subroutine
 
  logical function fstr_is_contact_conv(ctAlgo,infoCTChange,hecMESH)
    integer(kind=kint), intent(in)  :: ctalgo
    type (fstr_info_contactchange), intent(in) :: infoctchange
    type (hecmwst_local_mesh), intent(in) :: hecmesh
 
    fstr_is_contact_conv = .false.
    if( infoctchange%contact2free+infoctchange%contact2neighbor+      &
      infoctchange%contact2difflpos+infoctchange%free2contact == 0 ) &
      fstr_is_contact_conv = .true.
 
    call hecmw_allreduce_l1(hecmesh, fstr_is_contact_conv, hecmw_land)
  end function
 
  logical function fstr_is_matrixstructure_changed(infoCTChange)
    type (fstr_info_contactchange)   :: infoctchange
    fstr_is_matrixstructure_changed = .false.
    if( infoctchange%contact2free+infoctchange%contact2neighbor+infoctchange%free2contact > 0 ) &
      fstr_is_matrixstructure_changed = .true.
  end function
 
  subroutine fstr_update_contact_multiplier( cstep, ctAlgo, hecMESH, hecLagMAT, fstrSOLID, ctchanged )
    integer(kind=kint), intent(in)        :: cstep
    integer(kind=kint), intent(in)        :: ctAlgo
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(hecmwst_matrix_lagrange), intent(in) :: hecLagMAT
    type(fstr_solid), intent(inout)        :: fstrSOLID
    logical, intent(out)                   :: ctchanged
 
    integer(kind=kint) :: i, nc, algtype, grpid
 
    gnt = 0.d0;  ctchanged = .false.
    nc = fstrsolid%n_contacts+fstrsolid%n_embeds
    do i=1, fstrsolid%n_contacts
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      algtype = fstrsolid%contacts(i)%algtype
      if( algtype == contactsslid .or. algtype == contactfslid ) then
        call update_contact_multiplier( ctalgo, fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:), &
          fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff, &
          hecmesh, heclagmat, gnt, ctchanged )
      else if( algtype == contacttied ) then
        call update_tied_multiplier( fstrsolid%contacts(i), fstrsolid%unode(:), fstrsolid%dunode(:), &
        &  ctchanged )
      endif
    enddo
 
    do i=1, fstrsolid%n_embeds
      call update_tied_multiplier( fstrsolid%embeds(i), fstrsolid%unode(:), fstrsolid%dunode(:), &
      &  ctchanged )
    enddo
 
    call hecmw_allreduce_l1(hecmesh, ctchanged, hecmw_lor)
 
    if( nc>0 ) gnt = gnt/nc
  end subroutine
 
  subroutine fstr_update_contact_tangentforce( cstep, fstrSOLID )
    integer(kind=kint), intent(in)        :: cstep
    type(fstr_solid), intent(inout)        :: fstrSOLID
 
    integer(kind=kint) :: i, grpid
 
    do i=1, fstrsolid%n_contacts
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      call update_contact_tangentforce( fstrsolid%contacts(i) )
    enddo
  end subroutine
 
  subroutine fstr_update_contact_state_vectors( fstrSOLID, dt )
    type(fstr_solid), intent(inout)        :: fstrsolid
    real(kind=kreal), intent(in)           :: dt  
 
    integer(kind=kint) :: i
 
    if( associated( fstrsolid%CONT_RELVEL ) ) fstrsolid%CONT_RELVEL(:) = 0.d0
    if( associated( fstrsolid%CONT_STATE ) ) fstrsolid%CONT_STATE(:) = 0.d0
 
    do i=1, fstrsolid%n_contacts
      call update_contact_state_vectors( fstrsolid%contacts(i), dt, fstrsolid%CONT_RELVEL, fstrsolid%CONT_STATE )
    enddo
  end subroutine
 
end module mcontact