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
 
  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 :: candidate_intp = 3
 
  integer, parameter :: kcatfree = 1
  integer, parameter :: kcatnear = 2
  integer, parameter :: kcatcont = 3
 
  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 :: contactn2s = 1
  integer, parameter :: contacts2s = 2
 
  integer, parameter :: sparsity_none     = 0
  integer, parameter :: sparsity_neighbor = 1
 
  integer, parameter :: max_n_intp = 128
 
  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
    integer          :: id_lagrange
    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 tcontactsurf
    integer(kind=kint)              :: eid
    integer(kind=kint)              :: etype
    integer(kind=kint), pointer     :: nodes(:)=>null()     
    integer(kind=kint)              :: state = contactfree       
    integer(kind=kint)              :: state_prev = contactfree  
    integer(kind=kint)              :: n_intp
    type(tcontactstate), pointer    :: states(:)=>null()    
    integer(kind=kint), pointer     :: nslave_index(:)=>null() 
    ! --- lambda transaction buffers ---
    ! begin: previous substep's COMMIT product, immutable during a substep (warm-start source)
    ! working: current substep's accumulator, cleared at BEGIN, written by the augmentation update
    integer(kind=kint), pointer     :: lam_begin_id(:)=>null()  
    real(kind=kreal),   pointer     :: lam_begin_val(:,:)=>null() 
    integer(kind=kint)              :: lam_begin_n = 0
    integer(kind=kint), pointer     :: lam_work_id(:)=>null()   
    real(kind=kreal),   pointer     :: lam_work_val(:,:)=>null()  
    integer(kind=kint)              :: lam_work_n = 0
    ! --- tangential friction parallel arrays ---
    ! per-node covariant tangent multiplier (2 components, node a, rank r) and per-node
    ! stick/slip state. Rank r is keyed by lam_*_id (master), node a is the slave-surf node;
    ! both mirror lam_*_val.
    real(kind=kreal),   pointer     :: lam_begin_t(:,:,:)=>null()   
    real(kind=kreal),   pointer     :: lam_work_t(:,:,:)=>null()    
    integer(kind=kint), pointer     :: lam_begin_fstate(:,:)=>null()
    integer(kind=kint), pointer     :: lam_work_fstate(:,:)=>null() 
  end type tcontactsurf
 
  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                       :: n_master_owned = 0
    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               
 
    type(tcontactsurf), pointer   :: slave_surf(:)=>null()   
 
    ! !CONTACT_INTERFERENCE data; default-initialized because check_apply_Contact_IF only
    ! writes them when that card is present, while if_type is read unconditionally
    ! (if_flag = contact%if_type /= 0) by the contact force / search paths.
    real(kind=kreal)    :: ctime = 0.d0
    integer(kind=kint)  :: if_type = 0
    real(kind=kreal)    :: if_etime = 0.d0
    real(kind=kreal)    :: initial_pos = 0.d0
    real(kind=kreal)    :: end_pos = 0.d0
    ! 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
    integer                       :: method = contactn2s     
    integer                       :: sparsity_expansion = sparsity_none  
 
    logical                       :: mpced
    logical                       :: symmetric
    real(kind=kreal)              :: eps_fric_band = 0.d0    
 
    ! 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) :: n_statechange(3,3)
    integer(kind=kint) :: contact2neighbor
    integer(kind=kint) :: contact2beyond
    integer(kind=kint) :: contact2difflpos
    integer(kind=kint) :: free2contact_new
    integer(kind=kint) :: contactnode_previous
    integer(kind=kint) :: contactnode_current
  end type fstr_info_contactchange
 
  private :: is_active_contact
 
contains
 
  subroutine contact_state_init(cstate)
    type(tcontactstate), intent(inout) :: cstate
    cstate%state = -1
    cstate%surface = -1
    cstate%id_lagrange = 0
    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 is_contact_free(state)
    integer, intent(in) :: state
    is_contact_free = (state == contactfree)
  end function
 
  pure integer function contact_state_category(state)
    integer, intent(in) :: state
    if( is_contact_free(state) ) then
      contact_state_category = kcatfree
    else if( is_contact_active(state) ) then
      contact_state_category = kcatcont
    else
      contact_state_category = kcatnear
    endif
  end function
 
  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
 
  integer(kind=kint) function fstr_count_internal_slaves( contact, hecMESH )
    type(tcontact), intent(in)        :: contact
    type(hecmwst_local_mesh), pointer :: hecmesh
 
    integer  :: i, is, ie, cgrp
 
    fstr_count_internal_slaves = 0
 
    cgrp = contact%surf_id1
    if( cgrp<=0 ) return
    is= hecmesh%node_group%grp_index(cgrp-1) + 1
    ie= hecmesh%node_group%grp_index(cgrp  )
    do i=is,ie
      if( hecmesh%node_group%grp_item(i) <= hecmesh%nn_internal ) then
        fstr_count_internal_slaves = fstr_count_internal_slaves + 1
      endif
    enddo
  end function
 
  logical function fstr_contact_init( contact, hecMESH, cparam )
    type(tcontact), intent(inout)     :: contact
    type(hecmwst_local_mesh), pointer :: hecmesh
    type(tcontactparam), target       :: cparam
 
    integer  :: i, j, is, ie, cgrp, nsurf, nslave, ic, ic_type, iss, nn, ii
    integer  :: count, id_area
    logical  :: slave_owner, take_master
    integer, allocatable  :: slave_index(:)
 
    fstr_contact_init = .false.
 
    contact%cparam => cparam
 
    slave_owner = hecmw_partcontact_get_owner( hecmesh%hecmw_flag_partcontact ) == hecmw_flag_partcontact_owner_slave
 
    ! update_surface_normal normalises vertex normals only after the cross-rank assembly, so a uniform factor cancels
    ! but a partial one does not: a rank owning no slave must take no master surface at all, even though the
    ! partitioner does leave master elements here to keep the communication tables symmetric
    take_master = .true.
    if( slave_owner ) take_master = fstr_count_internal_slaves( contact, hecmesh ) > 0
 
    !  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  )
 
    ! Owned-only master face count, independent of take_master: the SURF-SURF visibility
    ! guard in fstr_setup sums it over ranks to get the global unique master face count.
    contact%n_master_owned = 0
    do i=is,ie
      ic   = hecmesh%surf_group%grp_item(2*i-1)
      if( hecmesh%elem_ID(ic*2) == hecmesh%my_rank ) contact%n_master_owned = contact%n_master_owned + 1
    enddo
 
    count = 0
    if( take_master ) then
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        id_area = hecmesh%elem_ID(ic*2)
        if( .not. slave_owner .and. id_area /= hecmesh%my_rank ) cycle
        count = count + 1
      enddo
    endif
    allocate( contact%master(count) )
    count = 0
    if( take_master ) then
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        id_area = hecmesh%elem_ID(ic*2)
        if( .not. slave_owner .and. id_area /= hecmesh%my_rank ) 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
    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
      if( slave_owner .and. hecmesh%node_group%grp_item(i) > hecmesh%nn_internal ) cycle
      nslave = nslave + 1
    enddo
    allocate( contact%slave(nslave) )
    allocate( slave_index(hecmesh%n_node) )
    slave_index(:) = 0
    ii = 0
    do i=is,ie
      if( slave_owner .and. hecmesh%node_group%grp_item(i) > hecmesh%nn_internal ) cycle
      ii = ii + 1
      contact%slave(ii) = hecmesh%node_group%grp_item(i)
      slave_index(contact%slave(ii)) = ii
    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 )
 
    if(contact%method == contacts2s) then
      ! The mortar integral passes the master shape functions through a length-4 array and
      ! sizes its element vectors for a first-order pair, so a second-order master face
      ! overruns them. initialize_csurf rejects the slave side for the same reason.
      do i=1, size( contact%master )
        if( contact%master(i)%etype /= fe_quad4n .and. contact%master(i)%etype /= fe_tri3n ) then
          write(*,*) '### Error: MORTAR=YES supports first-order surfaces only (quad4/tri3) : etype=', &
            contact%master(i)%etype
          stop hecmw_exit_model
        endif
      enddo
 
      !  slave surface
      cgrp = contact%surf_id1_sgrp
      if( cgrp<=0 ) return
      is= hecmesh%surf_group%grp_index(cgrp-1) + 1
      ie= hecmesh%surf_group%grp_index(cgrp  )
 
      ! Slave segments are taken owned-only (a serial mesh owns every element, so this is
      ! the full slave surface there). The master surface must be visible in full on every
      ! slave-owning rank; that is what !PARTITION, CONTACT_OWNER=SLAVE gives and what the
      ! visibility guard in fstr_setup checks.
      count = 0
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        if( hecmesh%elem_ID(ic*2) /= hecmesh%my_rank ) cycle
        count = count + 1
      enddo
      allocate( contact%slave_surf(count) )
      count = 0
      do i=is,ie
        ic   = hecmesh%surf_group%grp_item(2*i-1)
        if( hecmesh%elem_ID(ic*2) /= hecmesh%my_rank ) cycle
        count = count + 1
        nsurf = hecmesh%surf_group%grp_item(2*i)
        ic_type = hecmesh%elem_type(ic)
        call initialize_csurf( ic, ic_type, nsurf, contact%slave_surf(count) )
        iss = hecmesh%elem_node_index(ic-1)
        do j=1, size( contact%slave_surf(count)%nodes )
          nn = contact%slave_surf(count)%nodes(j)
          contact%slave_surf(count)%nodes(j) = hecmesh%elem_node_item( iss+nn )
          contact%slave_surf(count)%nslave_index(j) = slave_index(hecmesh%elem_node_item( iss+nn ))
        enddo
      enddo
 
      ! state for each integration points
      do i=1, size( contact%slave_surf )
        nn = contact%slave_surf(i)%n_intp
        allocate( contact%slave_surf(i)%states(nn) )
        do j = 1, contact%slave_surf(i)%n_intp
          contact%slave_surf(i)%states(j)%state = -1
          contact%slave_surf(i)%states(j)%multiplier(:) = 0.d0
          contact%slave_surf(i)%states(j)%tangentForce(:) = 0.d0
          contact%slave_surf(i)%states(j)%tangentForce1(:) = 0.d0
          contact%slave_surf(i)%states(j)%tangentForce_trial(:) = 0.d0
          contact%slave_surf(i)%states(j)%tangentForce_final(:) = 0.d0
          contact%slave_surf(i)%states(j)%reldisp(:) = 0.d0
          contact%slave_surf(i)%states(j)%time_factor = 0.d0
          contact%slave_surf(i)%states(j)%interference_flag = 0
        enddo
      enddo
    endif
 
    ! 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 )
    type(tcontact), intent(inout)     :: embed
    type(hecmwst_local_mesh), pointer :: hecmesh
    type(tcontactparam), target       :: cparam
 
    integer  :: i, j, is, ie, cgrp, nslave, ic, ic_type, iss, nn, ii
    integer  :: count, id_area
    logical  :: slave_owner, take_master
 
    fstr_embed_init = .false.
 
    embed%cparam => cparam
 
    slave_owner = hecmw_partcontact_get_owner( hecmesh%hecmw_flag_partcontact ) == hecmw_flag_partcontact_owner_slave
 
    take_master = .true.
    if( slave_owner ) take_master = fstr_count_internal_slaves( embed, hecmesh ) > 0
 
    !  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  )
 
    count = 0
    if( take_master ) then
      do i=is,ie
        ic   = hecmesh%elem_group%grp_item(i)
        id_area = hecmesh%elem_ID(ic*2)
        if( .not. slave_owner .and. id_area /= hecmesh%my_rank ) cycle
        count = count + 1
      enddo
    endif
    allocate( embed%master(count) )
    count = 0
    if( take_master ) then
      do i=is,ie
        ic   = hecmesh%elem_group%grp_item(i)
        id_area = hecmesh%elem_ID(ic*2)
        if( .not. slave_owner .and. id_area /= hecmesh%my_rank ) 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
    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
      if( slave_owner .and. hecmesh%node_group%grp_item(i) > hecmesh%nn_internal ) cycle
      nslave = nslave + 1
    enddo
    allocate( embed%slave(nslave) )
    ii = 0
    do i=is,ie
      if( slave_owner .and. hecmesh%node_group%grp_item(i) > hecmesh%nn_internal ) cycle
      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
        ! !CONTACT_INTERFERENCE is not implemented for the SURF-SURF mortar path
        if( contacts(i)%method == contacts2s ) then
          write(*,*) '### Error: CONTACT_INTERFERENCE is not supported with MORTAR=YES'
          stop hecmw_exit_model
        endif
        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 initialize_csurf( eid, etype, nsurf, surf )
    use elementinfo
    integer(kind=kint), intent(in)    :: eid
    integer(kind=kint), intent(in)    :: etype
    integer(kind=kint), intent(in)    :: nsurf
    type(tcontactsurf), intent(inout) :: surf
    integer(kind=kint) :: n, outtype, nodes(100)
    surf%eid = eid
 
    call getsubface( etype, nsurf, outtype, nodes )
    surf%etype = outtype
    n=getnumberofnodes( outtype )
    if(surf%etype == fe_quad4n )then
      surf%n_intp = 16
    else if(surf%etype == fe_tri3n )then
      surf%n_intp = 27
    else
      write(*,*) '### Error: MORTAR=YES supports first-order surfaces only (quad4/tri3) : etype=', surf%etype
      stop hecmw_exit_model
    end if
    allocate( surf%nodes(n) )
    allocate( surf%nslave_index(n) )
    surf%nodes(1:n)=nodes(1:n)
    surf%nslave_index(:)= 0
    ! lam_*_val: dim1 = slave-surf node a (size n), dim2 = rank r (master, size n_intp)
    allocate( surf%lam_begin_id(surf%n_intp), surf%lam_begin_val(n,surf%n_intp) )
    allocate( surf%lam_work_id (surf%n_intp), surf%lam_work_val (n,surf%n_intp) )
    surf%lam_begin_id(:) = 0; surf%lam_begin_val(:,:) = 0.0d0; surf%lam_begin_n = 0
    surf%lam_work_id (:) = 0; surf%lam_work_val (:,:) = 0.0d0; surf%lam_work_n  = 0
    ! tangential friction parallel arrays: (2, node a, rank r) / (node a, rank r)
    allocate( surf%lam_begin_t(2,n,surf%n_intp), surf%lam_work_t(2,n,surf%n_intp) )
    allocate( surf%lam_begin_fstate(n,surf%n_intp), surf%lam_work_fstate(n,surf%n_intp) )
    surf%lam_begin_t(:,:,:) = 0.0d0; surf%lam_work_t(:,:,:) = 0.0d0
    surf%lam_begin_fstate(:,:) = contactstick; surf%lam_work_fstate(:,:) = contactstick
  end subroutine
 
 
end module mcontactdef