fstr_contact_def.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 mcontactdef
 
  use hecmw
  use elementinfo
  use msurfelement
  use m_hecmw_contact_comm
  use bucket_search
  use mcontactparam
  use m_fstr_contact_smoothing
 
  implicit none
 
  real(kind=kreal), save :: cgn=1.d-5 
  real(kind=kreal), save :: cgt=1.d-3 
 
  real(kind=kreal), save :: gnt(2)    
  real(kind=kreal), save :: bakgnt(2) 
 
  integer, parameter :: contactunknown = -1
  integer, parameter :: contactfree = -1
  integer, parameter :: contactnear = 0
  integer, parameter :: contactstick = 1
  integer, parameter :: contactslip = 2
 
  integer, parameter :: contacttied = 1
  integer, parameter :: contactglued = 2
  integer, parameter :: contactsslid = 3
  integer, parameter :: contactfslid = 4
 
  integer, parameter :: kcsnone   = 0
  integer, parameter :: kcsnagata = 1
 
  integer, parameter :: c_if_slave = 1
  integer, parameter :: c_if_master = 2
 
  integer, parameter :: kctforresidual = 1
  integer, parameter :: kctforoutput   = 2
 
  type tcontactstate
    integer          :: state
    integer          :: surface
    real(kind=kreal) :: distance 
    real(kind=kreal) :: wkdist 
    real(kind=kreal) :: lpos(3) 
    real(kind=kreal) :: gpos(3) 
    real(kind=kreal) :: direction(3) 
    real(kind=kreal) :: multiplier(3) 
    real(kind=kreal) :: tangentforce(3) 
    real(kind=kreal) :: tangentforce1(3) 
    real(kind=kreal) :: tangentforce_trial(3) 
    real(kind=kreal) :: tangentforce_final(3) 
    real(kind=kreal)    :: reldisp(3)
    !
    real(kind=kreal)    :: shrink_factor
    real(kind=kreal)    :: time_factor
    real(kind=kreal)    :: init_pos
    real(kind=kreal)    :: end_pos
    integer             :: interference_flag
  end type
 
  type tcontact
    ! following contact definition
    character(len=HECMW_NAME_LEN) :: name
    integer                       :: ctype
    integer                       :: group
    character(len=HECMW_NAME_LEN) :: pair_name
    integer                       :: surf_id1, surf_id2
    integer                       :: surf_id1_sgrp
    type(tsurfelement), pointer   :: master(:)=>null()       
    integer, pointer              :: slave(:)=>null()        
    real(kind=kreal)              :: fcoeff                  
    real(kind=kreal)              :: npenalty                
    real(kind=kreal)              :: tpenalty                
    real(kind=kreal)              :: refstiff                
    real(kind=kreal)              :: damp_alpha              
    real(kind=kreal)              :: damp_gact               
    
    real(kind=kreal)    :: ctime
    integer(kind=kint)  :: if_type
    real(kind=kreal)    :: if_etime
    real(kind=kreal)    :: initial_pos
    real(kind=kreal)    :: end_pos
    ! following algorithm
    ! -1: not initialized
    ! 1: TIED-Just rigidly fixed the two surfaces
    ! 2: GLUED-Distance between the two surfaces to zero and glue them
    ! 3: SSLID-Small sliding contact( no position but with contact state change)
    ! 4: FSLID-Finite sliding contact (both changes in contact state and position possible)
    integer                       :: algtype
    integer                       :: smoothing
 
    logical                       :: mpced
    logical                       :: symmetric
 
    ! following contact state
    type(tcontactstate), pointer  :: states(:)=>null()       
 
    type(hecmwst_contact_comm)    :: comm
    type(bucketdb)                :: master_bktdb
 
    type(tcontactparam), pointer  :: cparam=>null()          
  end type tcontact
 
  type fstr_info_contactchange
    logical            :: active
    integer(kind=kint) :: contact2free
    integer(kind=kint) :: contact2neighbor
    integer(kind=kint) :: contact2difflpos
    integer(kind=kint) :: free2contact
    integer(kind=kint) :: contactnode_previous
    integer(kind=kint) :: contactnode_current
  end type fstr_info_contactchange
 
  private :: is_mpc_available
  private :: is_active_contact
 
contains
 
  subroutine contact_state_init(cstate)
    type(tcontactstate), intent(inout) :: cstate
    cstate%state = -1
    cstate%surface = -1
    cstate%distance = 0.0d0
    cstate%wkdist = 0.0d0
    cstate%lpos(:) = 0.0d0
    cstate%gpos(:) = 0.0d0
    cstate%direction(:) = 0.0d0
    cstate%multiplier(:) = 0.0d0
    cstate%tangentForce(:) = 0.0d0
    cstate%tangentForce1(:) = 0.0d0
    cstate%tangentForce_trial(:) = 0.0d0
    cstate%tangentForce_final(:) = 0.0d0
    cstate%reldisp(:) = 0.0d0
    cstate%shrink_factor = 0.0d0
    cstate%time_factor = 0.0d0
    cstate%init_pos = 0.0d0
    cstate%end_pos = 0.0d0
    cstate%interference_flag = 0
  end subroutine
 
  subroutine contact_state_copy(cstate1, cstate2)
    type(tcontactstate), intent(in)    :: cstate1
    type(tcontactstate), intent(inout) :: cstate2
    cstate2 = cstate1
  end subroutine
 
  pure logical function is_contact_active(state)
    integer, intent(in) :: state
    is_contact_active = (state >= contactstick)
  end function
 
  pure logical function has_projection(state)
    integer, intent(in) :: state
    has_projection = (state >= contactnear)
  end function
 
  pure logical function is_contact_free(state)
    integer, intent(in) :: state
    is_contact_free = (state == contactfree)
  end function
 
  subroutine print_contact_state(fnum, cstate)
    integer, intent(in)             :: fnum
    type(tcontactstate), intent(in) :: cstate
    write(fnum, *) "--Contact state=",cstate%state
    write(fnum, *) cstate%surface, cstate%distance
    write(fnum, *) cstate%lpos
    write(fnum, *) cstate%direction
    write(fnum, *) cstate%multiplier
  end subroutine
 
  logical function is_mpc_available( contact )
    type(tcontact), intent(in)        :: contact
    is_mpc_available = .true.
    if( contact%fcoeff/=0.d0 ) is_mpc_available = .false.
  end function
 
  subroutine fstr_write_contact( file, contact )
    integer(kind=kint), intent(in)    :: file
    type(tcontact), intent(in)        :: contact
    integer :: i
    write(file,*) "CONTACT:", contact%ctype,contact%group,trim(contact%pair_name),contact%fcoeff
    write(file,*) "---Slave----"
    write(file,*) "num.slave",size(contact%slave)
    if( associated(contact%slave) ) then
      do i=1,size(contact%slave)
        write(file, *) contact%slave(i)
      enddo
    endif
    write(file,*) "----master---"
    write(file,*) "num.master",size(contact%master)
    if( associated(contact%master) ) then
      do i=1,size(contact%master)
        call write_surf( file, contact%master(i) )
      enddo
    endif
  end subroutine
 
  subroutine fstr_contact_finalize( contact )
    type(tcontact), intent(inout)     :: contact
    integer  :: i
    if( associated( contact%slave ) ) deallocate(contact%slave)
    if( associated( contact%master ) ) then
      do i=1,size( contact%master )
        call  finalize_surf( contact%master(i) )
      enddo
      deallocate(contact%master)
    endif
    if( associated(contact%states) ) deallocate(contact%states)
    call hecmw_contact_comm_finalize(contact%comm)
    call bucketdb_finalize( contact%master_bktDB )
  end subroutine
 
  logical function fstr_contact_check( contact, hecMESH )
    type(tcontact), intent(inout)     :: contact
    type(hecmwst_local_mesh), pointer :: hecmesh
 
    integer  :: i
    logical  :: isfind
 
    fstr_contact_check = .false.
 
    ! if contact pair exist?
    isfind = .false.
    do i=1,hecmesh%contact_pair%n_pair
      if( hecmesh%contact_pair%name(i) == contact%pair_name ) then
        contact%ctype = hecmesh%contact_pair%type(i)
        contact%surf_id1 = hecmesh%contact_pair%slave_grp_id(i)
        contact%surf_id2 = hecmesh%contact_pair%master_grp_id(i)
        contact%surf_id1_sgrp = hecmesh%contact_pair%slave_orisgrp_id(i)
        isfind = .true.
      endif
    enddo
    if( .not. isfind ) return;
    if( contact%fcoeff<=0.d0 ) contact%fcoeff=0.d0
    if( contact%ctype < 1 .and. contact%ctype > 3 ) return
    if( contact%group<=0 ) return
 
    fstr_contact_check = .true.
  end function
 
  logical function fstr_contact_init( contact, hecMESH, cparam, myrank )
    type(tcontact), intent(inout)     :: contact
    type(hecmwst_local_mesh), pointer :: hecmesh
    type(tcontactparam), target       :: cparam
    integer(kint),intent(in),optional :: myrank
 
    integer  :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
    integer  :: count,nodeid
 
    fstr_contact_init = .false.
 
    contact%cparam => cparam
 
    !  master surface
    cgrp = contact%surf_id2
    if( cgrp<=0 ) return
    is= hecmesh%surf_group%grp_index(cgrp-1) + 1
    ie= hecmesh%surf_group%grp_index(cgrp  )
 
    if(present(myrank)) then
      ! PARA_CONTACT
      count = 0
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        if(hecmesh%elem_ID(ic*2) /= myrank) cycle
        count = count + 1
      enddo
      allocate( contact%master(count) )
      count = 0
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        if(hecmesh%elem_ID(ic*2) /= myrank) cycle
        count = count + 1
        nsurf = hecmesh%surf_group%grp_item(2*i)
        ic_type = hecmesh%elem_type(ic)
        call initialize_surf( ic, ic_type, nsurf, contact%master(count) )
        iss = hecmesh%elem_node_index(ic-1)
        do j=1, size( contact%master(count)%nodes )
          nn = contact%master(count)%nodes(j)
          contact%master(count)%nodes(j) = hecmesh%elem_node_item( iss+nn )
        enddo
      enddo
 
    else
      ! not PARA_CONTACT
      allocate( contact%master(ie-is+1) )
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        nsurf = hecmesh%surf_group%grp_item(2*i)
        ic_type = hecmesh%elem_type(ic)
        call initialize_surf( ic, ic_type, nsurf, contact%master(i-is+1) )
        iss = hecmesh%elem_node_index(ic-1)
        do j=1, size( contact%master(i-is+1)%nodes )
          nn = contact%master(i-is+1)%nodes(j)
          contact%master(i-is+1)%nodes(j) = hecmesh%elem_node_item( iss+nn )
        enddo
      enddo
 
    endif
 
    call update_surface_reflen( contact%master, hecmesh%node )
 
    cgrp = contact%surf_id1
    if( cgrp<=0 ) return
    is= hecmesh%node_group%grp_index(cgrp-1) + 1
    ie= hecmesh%node_group%grp_index(cgrp  )
    nslave = 0
    do i=is,ie
      nodeid = hecmesh%global_node_ID(hecmesh%node_group%grp_item(i))
      if(present(myrank)) then
        ! PARA_CONTACT
        nslave = nslave + 1
      else
        ! not PARA_CONTACT
        if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal) then
          nslave = nslave + 1
        endif
      endif
    enddo
    allocate( contact%slave(nslave) )
    ii = 0
    do i=is,ie
      if(.not.present(myrank)) then
        ! not PARA_CONTACT
        if( hecmesh%node_group%grp_item(i) > hecmesh%nn_internal) cycle
      endif
      ii = ii + 1
      contact%slave(ii) = hecmesh%node_group%grp_item(i)
    enddo
 
    ! contact state
    allocate( contact%states(nslave) )
    do i=1,nslave
      call contact_state_init( contact%states(i) )
    enddo
 
    ! neighborhood of surface group
    call update_surface_box_info( contact%master, hecmesh%node )
    call bucketdb_init( contact%master_bktDB )
    call update_surface_bucket_info( contact%master, contact%master_bktDB )
    call find_surface_neighbor( contact%master, contact%master_bktDB )
 
    ! initialize contact communication table
    call hecmw_contact_comm_init( contact%comm, hecmesh, 1, nslave, contact%slave )
 
    contact%symmetric = .true.
    fstr_contact_init = .true.
  end function
 
  logical function fstr_embed_init( embed, hecMESH, cparam, myrank )
    type(tcontact), intent(inout)     :: embed
    type(hecmwst_local_mesh), pointer :: hecmesh
    type(tcontactparam), target       :: cparam
    integer(kint),intent(in),optional :: myrank
 
    integer  :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
    integer  :: count,nodeid
 
    fstr_embed_init = .false.
 
    embed%cparam => cparam
 
    !  master surface
    cgrp = embed%surf_id2
    if( cgrp<=0 ) return
    is= hecmesh%elem_group%grp_index(cgrp-1) + 1
    ie= hecmesh%elem_group%grp_index(cgrp  )
 
    if(present(myrank)) then
      ! PARA_CONTACT
      count = 0
      do i=is,ie
        ic   = hecmesh%elem_group%grp_item(i)
        if(hecmesh%elem_ID(ic*2) /= myrank) cycle
        count = count + 1
      enddo
      allocate( embed%master(count) )
      count = 0
      do i=is,ie
        ic   = hecmesh%elem_group%grp_item(i)
        if(hecmesh%elem_ID(ic*2) /= myrank) cycle
        count = count + 1
        ic_type = hecmesh%elem_type(ic)
        call initialize_surf( ic, ic_type, 0, embed%master(count) )
        iss = hecmesh%elem_node_index(ic-1)
        do j=1, size( embed%master(count)%nodes )
          nn = embed%master(count)%nodes(j)
          embed%master(count)%nodes(j) = hecmesh%elem_node_item( iss+nn )
        enddo
      enddo
 
    else
      ! not PARA_CONTACT
      allocate( embed%master(ie-is+1) )
      do i=is,ie
        ic   = hecmesh%elem_group%grp_item(i)
        ic_type = hecmesh%elem_type(ic)
        call initialize_surf( ic, ic_type, 0, embed%master(i-is+1) )
        iss = hecmesh%elem_node_index(ic-1)
        do j=1, size( embed%master(i-is+1)%nodes )
          nn = embed%master(i-is+1)%nodes(j)
          embed%master(i-is+1)%nodes(j) = hecmesh%elem_node_item( iss+nn )
        enddo
      enddo
 
    endif
 
    ! slave surface
    cgrp = embed%surf_id1
    if( cgrp<=0 ) return
    is= hecmesh%node_group%grp_index(cgrp-1) + 1
    ie= hecmesh%node_group%grp_index(cgrp  )
    nslave = 0
    do i=is,ie
      nodeid = hecmesh%global_node_ID(hecmesh%node_group%grp_item(i))
      if(present(myrank)) then
        ! PARA_CONTACT
        nslave = nslave + 1
      else
        ! not PARA_CONTACT
        if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal) then
          nslave = nslave + 1
        endif
      endif
    enddo
    allocate( embed%slave(nslave) )
    ii = 0
    do i=is,ie
      if(.not.present(myrank)) then
        ! not PARA_CONTACT
        if( hecmesh%node_group%grp_item(i) > hecmesh%nn_internal) cycle
      endif
      ii = ii + 1
      embed%slave(ii) = hecmesh%node_group%grp_item(i)
    enddo
 
    ! embed state
    allocate( embed%states(nslave) )
    do i=1,nslave
      call contact_state_init( embed%states(i) )
    enddo
 
    ! neighborhood of surface group
    call update_surface_box_info( embed%master, hecmesh%node )
    call bucketdb_init( embed%master_bktDB )
    call update_surface_bucket_info( embed%master, embed%master_bktDB )
    call find_surface_neighbor( embed%master, embed%master_bktDB )
 
    ! initialize contact communication table
    call hecmw_contact_comm_init( embed%comm, hecmesh, 1, nslave, embed%slave )
 
    ! initialize penalty coefficients
    embed%nPenalty = 1.0d0     ! default normal penalty coefficient
    embed%tPenalty = 0.1d0     ! default tangential penalty coefficient
    embed%refStiff = 0.0d0     ! will be calculated after first stiffness assembly
    embed%damp_alpha = 0.0d0
    embed%damp_gact = 0.0d0
 
    embed%symmetric = .true.
    fstr_embed_init = .true.
  end function
 
  function check_apply_contact_if( contact_if, contacts )
    type(tcontactinterference), intent(inout)     :: contact_if
    type(tcontact)     :: contacts(:)
    
    integer  :: i, j
    logical  :: isfind
    integer(kind=kint)            :: check_apply_contact_if
 
    check_apply_contact_if = -1
    ! if contact pair exist?
    isfind = .false.
    do i = 1, size(contacts)
      if( contacts(i)%pair_name == contact_if%cp_name ) then
        contacts(i)%if_type     = contact_if%if_type
        contacts(i)%if_etime    = contact_if%etime
        contacts(i)%initial_pos = contact_if%initial_pos
        contacts(i)%end_pos     = contact_if%end_pos
        do j = 1, size(contacts(i)%states)
          contacts(i)%states(j)%interference_flag = contact_if%if_type
          contacts(i)%states(j)%init_pos = contact_if%initial_pos
          contacts(i)%states(j)%end_pos  = contact_if%end_pos
          if( contact_if%if_type /= c_if_slave )then
            contacts(i)%states(j)%time_factor = (contact_if%end_pos - contact_if%initial_pos) / contact_if%etime
          else
            contacts(i)%states(j)%time_factor = contact_if%etime
          end if
        end do
        isfind = .true.
        check_apply_contact_if = 0; return
      endif
    enddo
    if( .not. isfind ) return;
    check_apply_contact_if = 0
 
  end function
 
  subroutine clear_contact_state( contact )
    type(tcontact), intent(inout) :: contact
    integer :: i
    if( .not. associated(contact%states) ) return
    do i=1,size( contact%states )
      contact%states(i)%state = -1
    enddo
  end subroutine
 
  logical function is_active_contact( acgrp, contact )
    integer, intent(in)        :: acgrp(:)
    type(tcontact), intent(in) :: contact
    if( any( acgrp==contact%group ) ) then
      is_active_contact = .true.
    else
      is_active_contact = .false.
    endif
  end function
 
  subroutine print_contatct_pair( file, pair )
    integer(kind=kint), intent(in)           :: file
    type( hecmwst_contact_pair ), intent(in) :: pair
 
    integer(kind=kint) :: i
    write(file,*) "Number of contact pair", pair%n_pair
    do i=1,pair%n_pair
      write(file,*) trim(pair%name(i)), pair%type(i), pair%slave_grp_id(i)  &
        ,pair%master_grp_id(i), pair%slave_orisgrp_id(i)
    enddo
  end subroutine
 
 
end module mcontactdef