fstr_mat_con_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
!-------------------------------------------------------------------------------
!         stiffness matrix structure for the contact analysis
!         employing standard Lagrange multiplier algorithm
 
module fstr_matrix_con_contact
 
  use m_fstr
  use elementinfo
  use m_fstr_contact_damping, only: is_damping_enabled
 
  implicit none
  private
  public :: fstr_get_num_lagrange_pernode
  public :: hecmwst_matrix_lagrange
  public :: fstr_save_originalmatrixstructure
  public :: fstr_mat_con_contact
  public :: fstr_is_matrixstruct_symmetric
  public :: fstr_set_lagrange_diagonal
  public :: fstr_get_lagrange_diagonal
 
  integer(kind=kint), save         :: NPL_org, NPU_org
  type(nodeRelated), pointer, save :: list_nodeRelated_org(:) => null() 
 
  type(nodeRelated), pointer       :: list_nodeRelated(:) => null() 
 
  logical                          :: permission = .false.
 
contains
 
  integer(kind=kint) function fstr_get_num_lagrange_pernode(algtype)
    integer(kind=kint) :: algtype
    if( algtype == contactsslid .or. algtype == contactfslid ) then
      fstr_get_num_lagrange_pernode = 1
    else if( algtype == contacttied ) then
      fstr_get_num_lagrange_pernode = 3
    endif
  end function
 
  subroutine fstr_save_originalmatrixstructure(hecMAT)
 
    type(hecmwst_matrix) :: hecmat
 
    if( associated(list_noderelated_org) ) return
    call hecmw_construct_noderelated_from_hecmat(hecmat, npl_org, npu_org, list_noderelated_org)
 
  end subroutine fstr_save_originalmatrixstructure
 
  subroutine fstr_mat_con_contact(cstep,contact_algo,hecMAT,fstrSOLID,hecLagMAT,infoCTChange,conMAT,is_contact_active_flag)
 
    integer(kind=kint)                   :: cstep
    integer(kind=kint)                   :: contact_algo
    type(hecmwst_matrix)                 :: hecmat
    type(fstr_solid)                     :: fstrsolid
    type(hecmwst_matrix_lagrange) :: heclagmat
    type(fstr_info_contactchange)        :: infoctchange
 
    integer(kind=kint)                   :: num_lagrange
    integer(kind=kint)                   :: countnon0lu_node, countnon0lu_lagrange
    integer(kind=kint)                   :: numnon0_node, numnon0_lagrange
    type (hecmwst_matrix)                :: conmat
    logical, intent(in)                  :: is_contact_active_flag
 
    integer(kind=kint)                   :: i, j, grpid
    integer(kind=kint)                   :: nlag
 
    num_lagrange = 0
    if( contact_algo == kcaslagrange ) then
      do i = 1, fstrsolid%n_contacts
        grpid = fstrsolid%contacts(i)%group
        if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
        nlag = fstr_get_num_lagrange_pernode(fstrsolid%contacts(i)%algtype)
        do j = 1, size(fstrsolid%contacts(i)%slave)
          if( .not. is_contact_active(fstrsolid%contacts(i)%states(j)%state) ) cycle
          num_lagrange = num_lagrange + nlag
        enddo
      enddo
 
      do i = 1, fstrsolid%n_embeds
        grpid = fstrsolid%embeds(i)%group
        if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
        nlag = 3
        do j = 1, size(fstrsolid%embeds(i)%slave)
          if( .not. is_contact_active(fstrsolid%embeds(i)%states(j)%state) ) cycle
          num_lagrange = num_lagrange + nlag
        enddo
      enddo
    endif
 
    ! Get original list of related nodes
    call hecmw_init_noderelated_from_org(hecmat%NP,num_lagrange,is_contact_active_flag,list_noderelated_org,list_noderelated)
 
    ! Construct new list of related nodes and Lagrange multipliers
    countnon0lu_node = npl_org + npu_org
    countnon0lu_lagrange = 0
    if( is_contact_active_flag ) call getnewlistofrelatednodesandlagrangemultipliers(cstep,contact_algo, &
    &  hecmat%NP,fstrsolid,countnon0lu_node,countnon0lu_lagrange,list_noderelated)
 
    ! Construct new matrix structure(hecMAT&hecLagMAT)
    numnon0_node = countnon0lu_node/2
    numnon0_lagrange = countnon0lu_lagrange/2
    call hecmw_construct_hecmat_from_noderelated(hecmat%N, hecmat%NP, hecmat%NDOF, &
    & numnon0_node, num_lagrange, list_noderelated, hecmat)
    call hecmw_construct_hecmat_from_noderelated(hecmat%N, hecmat%NP, hecmat%NDOF, &
    & numnon0_node, num_lagrange, list_noderelated, conmat)
    if( contact_algo == kcaslagrange ) call hecmw_construct_heclagmat_from_noderelated(hecmat%NP, &
    & hecmat%NDOF, num_lagrange, numnon0_lagrange, is_contact_active_flag, list_noderelated, heclagmat)
    call hecmw_finalize_noderelated(list_noderelated)
 
    ! Copy Lagrange multipliers
    if( is_contact_active_flag .and. contact_algo == kcaslagrange ) &
      call fstr_copy_lagrange_contact(fstrsolid,heclagmat)
 
  end subroutine fstr_mat_con_contact
 
  subroutine getnewlistofrelatednodesandlagrangemultipliers( &
  & cstep, contact_algo, np, fstrSOLID, countNon0LU_node, countNon0LU_lagrange, list_nodeRelated )
    integer(kind=kint),intent(in)             :: cstep
    integer(kind=kint),intent(in)             :: contact_algo
    integer(kind=kint),intent(in)             :: np
    type(fstr_solid),intent(in)               :: fstrsolid
    integer(kind=kint), intent(inout)         :: countnon0lu_node, countnon0lu_lagrange
    type(noderelated), pointer, intent(inout) :: list_noderelated(:)
 
    integer(kind=kint)            :: grpid
    integer(kind=kint)            :: count_lagrange
    integer(kind=kint)            :: ctsurf, etype, nnode, ndlocal(l_max_surface_node + 1)
    integer(kind=kint)            :: i, j, k, nlag, algtype
    real(kind=kreal)              :: fcoeff 
    logical                       :: necessary_to_insert_node, necessary_to_insert_node_pair
    logical                       :: is_contact_active_flag, is_damping_active_flag
 
    count_lagrange = 0
    do i = 1, fstrsolid%n_contacts
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      fcoeff = fstrsolid%contacts(i)%fcoeff
      necessary_to_insert_node = ( fcoeff /= 0.0d0 .or. contact_algo == kcaalagrange )
 
      algtype = fstrsolid%contacts(i)%algtype
      nlag = fstr_get_num_lagrange_pernode(algtype)
      if( contact_algo == kcaalagrange ) nlag = 1
      if( algtype == contacttied ) permission = .true.
 
      do j = 1, size(fstrsolid%contacts(i)%slave)
        ! stick or sliding contact is active
        is_contact_active_flag = is_contact_active(fstrsolid%contacts(i)%states(j)%state)
        ! damping is active
        is_damping_active_flag = fstrsolid%contacts(i)%states(j)%state == contactnear .and. &
          &  is_damping_enabled(fstrsolid%contacts(i))
        
        if( is_contact_active_flag .or. is_damping_active_flag ) then
 
          ctsurf = fstrsolid%contacts(i)%states(j)%surface
          etype = fstrsolid%contacts(i)%master(ctsurf)%etype
          if( etype/=fe_tri3n .and. etype/=fe_quad4n ) &
            stop " ##Error: This element type is not supported in contact analysis !!! "
          nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
          ndlocal(1) = fstrsolid%contacts(i)%slave(j)
          ndlocal(2:nnode+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(1:nnode)
 
          ! For CONTACTNEAR damping (especially S-Lagrange + frictionless),
          ! we still need slave-master connectivity to assemble damping terms.
          necessary_to_insert_node_pair = necessary_to_insert_node .or. is_damping_active_flag
 
          if( is_contact_active_flag ) then
            do k=1,nlag
              if( contact_algo == kcaslagrange ) count_lagrange = count_lagrange + 1
              call hecmw_ass_noderelated_from_contact_pair(np, nnode, ndlocal, count_lagrange, permission, &
              & necessary_to_insert_node_pair, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
            enddo
          else
            ! NEAR damping only: no Lagrange multiplier, insert connectivity once
            call hecmw_ass_noderelated_from_contact_pair(np, nnode, ndlocal, 0, permission, &
            & necessary_to_insert_node_pair, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
          endif
              
        end if
 
      enddo
 
    enddo
 
    do i = 1, fstrsolid%n_embeds
 
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_isembedactive( fstrsolid, grpid, cstep ) ) cycle
 
      necessary_to_insert_node = ( contact_algo == kcaalagrange )
 
      nlag = 3
      if( contact_algo == kcaalagrange ) nlag = 1
      permission = .true.
 
      do j = 1, size(fstrsolid%embeds(i)%slave)
 
        if( .not. is_contact_active(fstrsolid%embeds(i)%states(j)%state) ) cycle
        ctsurf = fstrsolid%embeds(i)%states(j)%surface
        etype = fstrsolid%embeds(i)%master(ctsurf)%etype
        nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
        ndlocal(1) = fstrsolid%embeds(i)%slave(j)
        ndlocal(2:nnode+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(1:nnode)
 
        do k=1,nlag
          if( contact_algo == kcaslagrange ) count_lagrange = count_lagrange + 1
          call hecmw_ass_noderelated_from_contact_pair(np, nnode, ndlocal, count_lagrange, permission, &
          & necessary_to_insert_node, list_noderelated_org, list_noderelated, countnon0lu_node, countnon0lu_lagrange )
        enddo
      enddo
 
    enddo
 
  end subroutine getnewlistofrelatednodesandlagrangemultipliers
 
  subroutine fstr_copy_lagrange_contact(fstrSOLID,hecLagMAT)
 
    type(fstr_solid)                        :: fstrsolid
    type(hecmwst_matrix_lagrange)          :: heclagmat
    integer (kind=kint)                    :: id_lagrange, algtype, i, j, k, nlag, slave_node
 
    id_lagrange = 0
 
    do i = 1, fstrsolid%n_contacts
 
      algtype = fstrsolid%contacts(i)%algtype
      nlag = fstr_get_num_lagrange_pernode(algtype)
 
      do j = 1, size(fstrsolid%contacts(i)%slave)
        if( .not. is_contact_active(fstrsolid%contacts(i)%states(j)%state) ) cycle
        slave_node = fstrsolid%contacts(i)%slave(j)
        heclagmat%lag_node_table(slave_node) = id_lagrange + 1
        do k=1,nlag
          id_lagrange = id_lagrange + 1
          heclagmat%Lagrange(id_lagrange)=fstrsolid%contacts(i)%states(j)%multiplier(k)
        enddo
      enddo
    enddo
 
    do i = 1, fstrsolid%n_embeds
      nlag = 3
      do j = 1, size(fstrsolid%embeds(i)%slave)
        if( .not. is_contact_active(fstrsolid%embeds(i)%states(j)%state) ) cycle
        slave_node = fstrsolid%embeds(i)%slave(j)
        heclagmat%lag_node_table(slave_node) = id_lagrange + 1
        do k=1,nlag
          id_lagrange = id_lagrange + 1
          heclagmat%Lagrange(id_lagrange)=fstrsolid%embeds(i)%states(j)%multiplier(k)
        enddo
      enddo
    enddo
 
  end subroutine fstr_copy_lagrange_contact
 
  logical function fstr_is_matrixstruct_symmetric(fstrSOLID,hecMESH)
 
    type(fstr_solid )        :: fstrsolid
    type(hecmwst_local_mesh) :: hecmesh
    integer (kind=kint)      :: is_in_contact
 
    is_in_contact = 0
    if( fstrsolid%n_contacts>0 ) then
      if( any(fstrsolid%contacts(:)%fcoeff /= 0.0d0) )  is_in_contact = 1
    endif
    call hecmw_allreduce_i1(hecmesh, is_in_contact, hecmw_max)
    if( is_in_contact == 0 .and. hecmesh%n_dof /= 4 ) then
      fstr_is_matrixstruct_symmetric = .true.
    else
      fstr_is_matrixstruct_symmetric = .false.
    endif
 
  end function fstr_is_matrixstruct_symmetric
 
  subroutine fstr_set_lagrange_diagonal(hecLagMAT, ilag, value)
    type(hecmwst_matrix_lagrange), intent(inout) :: heclagmat
    integer(kind=kint), intent(in) :: ilag
    real(kind=kreal), intent(in) :: value 
 
    if (ilag < 1 .or. ilag > heclagmat%num_lagrange) then
      write(*,*) 'Error in fstr_set_lagrange_diagonal: invalid Lagrange multiplier index', ilag
      stop
    endif
 
    if (.not. associated(heclagmat%D_lagrange)) then
      write(*,*) 'Error in fstr_set_lagrange_diagonal: D_lagrange not allocated'
      stop
    endif
 
    heclagmat%D_lagrange(ilag) = value
 
  end subroutine fstr_set_lagrange_diagonal
 
  real(kind=kreal) function fstr_get_lagrange_diagonal(hecLagMAT, ilag)
    type(hecmwst_matrix_lagrange), intent(in) :: heclagmat
    integer(kind=kint), intent(in) :: ilag
 
    if (ilag < 1 .or. ilag > heclagmat%num_lagrange) then
      write(*,*) 'Error in fstr_get_lagrange_diagonal: invalid Lagrange multiplier index', ilag
      fstr_get_lagrange_diagonal = 0.0d0
      return
    endif
 
    if (.not. associated(heclagmat%D_lagrange)) then
      write(*,*) 'Error in fstr_get_lagrange_diagonal: D_lagrange not allocated'
      fstr_get_lagrange_diagonal = 0.0d0
      return
    endif
 
    fstr_get_lagrange_diagonal = heclagmat%D_lagrange(ilag)
 
  end function fstr_get_lagrange_diagonal
 
 
end module fstr_matrix_con_contact