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
  implicit none
 
  private :: l_contact2mpc, l_tied2mpc
  integer(kind=kint), save :: n_contact_mpc
  logical, private :: active
 
  real(kind=kreal), save :: mu=1.d10  
  real(kind=kreal), save :: mut=1.d6  
  real(kind=kreal), save :: cdotp=1.d3  
 
  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) 
 
contains
 
  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
 
  subroutine fstr_set_contact_penalty( maxv )
    real(kind=kreal), intent(in) :: maxv
    mu = cdotp * maxv
    if( gnt(1)<1.d-3 ) mu=cdotp*10.d0* maxv
    bakgnt = 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
    integer(kind=kint) :: is_conv
 
    fstr_is_contact_conv = .false.
    if( infoctchange%contact2free+infoctchange%contact2neighbor+      &
    infoctchange%contact2difflpos+infoctchange%free2contact == 0 ) &
    fstr_is_contact_conv = .true.
 
    is_conv = 0
    if (fstr_is_contact_conv) is_conv = 1
    call hecmw_allreduce_i1(hecmesh, is_conv, hecmw_min)
    if (is_conv == 0) then
      fstr_is_contact_conv = .false.
    else
      fstr_is_contact_conv = .true.
    endif
  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 l_contact2mpc( contact, mpcs, nmpc )
    use fstr_ctrl_modifier
    type( tcontact ), intent(in)         :: contact
    type( hecmwst_mpc ), intent(inout)   :: mpcs
    integer(kind=kint), intent(out)      :: nmpc
    integer(kind=kint), parameter :: ndof = 3        ! 3D problem only, currently
    real(kind=kreal), parameter   :: tol =1.d-10
    integer(kind=kint) :: i, j, k, nn, csurf, nenode, etype, tdof
    integer(kind=kint) :: nodes(l_max_surface_node*ndof+ndof), dofs(l_max_surface_node*ndof+ndof)
    real(kind=kreal) :: values(l_max_surface_node*ndof+ndof+1),val(l_max_surface_node*ndof+ndof+1)
    nmpc=0
    do i=1,size(contact%states)
      if( contact%states(i)%state == -1 ) cycle    ! in free
      csurf = contact%states(i)%surface
      if( csurf<=0 ) stop "error in contact state"
      etype =  contact%master(csurf)%etype
      nenode = size(contact%master(csurf)%nodes)
      tdof = nenode*ndof+ndof
      call contact2mpcval( contact%states(i), etype, nenode, values(1:tdof+1) )
      tdof = 0
      do j=1,ndof
        if( dabs(values(j))<tol ) cycle
        tdof = tdof+1
        nodes(tdof) = contact%slave(i)
        dofs(tdof) = j
        val(tdof) = values(j)
      enddo
      do j=1,nenode
        nn =  contact%master(csurf)%nodes(j)
        nodes( j*ndof+1:j*ndof+ndof ) =  nn
        do k=1,ndof
          if( dabs(values(j*ndof+k)) < tol ) cycle
          tdof=tdof+1
          nodes(tdof)=nn
          dofs(tdof ) = k
          val(tdof)=values(j*ndof+k)
        enddo
      enddo
      val(tdof+1) = values(nenode*ndof+ndof+1)
 
      call fstr_append_mpc( tdof, nodes(1:tdof), dofs(1:tdof), val(1:tdof+1), mpcs )
      nmpc=nmpc+1
    enddo
  end subroutine l_contact2mpc
 
  subroutine l_tied2mpc( contact, mpcs, nmpc )
    use fstr_ctrl_modifier
    type( tcontact ), intent(in)         :: contact
    type( hecmwst_mpc ), intent(inout)   :: mpcs
    integer(kind=kint), intent(out)      :: nmpc
    integer(kind=kint) :: i, j, csurf, nenode, etype, tdof
    integer(kind=kint) :: nodes(l_max_surface_node+1), dofs(l_max_surface_node+1)
    real(kind=kreal) :: values(l_max_surface_node+2)
    nmpc=0
    do i=1,size(contact%slave)
      csurf = contact%states(i)%surface
      if( csurf<=0 ) cycle                         ! contactor not exists
      nenode = size(contact%master(csurf)%nodes)
      tdof = nenode+1
      nodes(1) = contact%slave(i)
      nodes( 2:tdof ) =  contact%master(csurf)%nodes(:)
      values(1) = -1.d0
      values(2:tdof) = 1.d0
      values(tdof+1) = 0.d0
      etype =  contact%master(csurf)%etype
      do j=1,3
        dofs(1:tdof) = j
        call fstr_append_mpc( tdof, nodes(1:tdof), dofs(1:tdof), values(1:tdof+1), mpcs )
        nmpc=nmpc+1
      enddo
    enddo
  end subroutine l_tied2mpc
 
  subroutine fstr_contact2mpc( contacts, mpcs )
    type( tcontact ), intent(in)       :: contacts(:)
    type( hecmwst_mpc ), intent(inout) :: mpcs
    integer(kind=kint) :: i, nmpc
    n_contact_mpc = 0
    do i=1,size(contacts)
      if( contacts(i)%algtype == contactunknown ) cycle     ! not initialized
      if( contacts(i)%algtype == contactfslid ) then
        print *, "Cannot deal with finit slip problems by MPC!"
        cycle
      endif
      if( contacts(i)%algtype == contactsslid ) then
        call l_contact2mpc( contacts(i), mpcs, nmpc )
        n_contact_mpc = n_contact_mpc + nmpc
      elseif( contacts(i)%algtype == contacttied ) then
        call l_tied2mpc( contacts(i), mpcs, nmpc )
        n_contact_mpc = n_contact_mpc + nmpc
      endif
    enddo
  end subroutine
 
  subroutine fstr_del_contactmpc( mpcs )
    use fstr_ctrl_modifier
    type( hecmwst_mpc ), intent(inout) :: mpcs
    call fstr_delete_mpc( n_contact_mpc, mpcs )
  end subroutine
 
  subroutine fstr_write_mpc( file, mpcs )
    integer(kind=kint), intent(in)  :: file
    type( hecmwst_mpc ), intent(in) :: mpcs
 
    integer(kind=kint) :: i,j,n0,n1
    write(file, *) "Number of equation", mpcs%n_mpc
    do i=1,mpcs%n_mpc
      write(file,*) "--Equation",i
      n0=mpcs%mpc_index(i-1)+1
      n1=mpcs%mpc_index(i)
      write(file, *) n0,n1
      write(file,'(30i5)') (mpcs%mpc_item(j),j=n0,n1)
      write(file,'(30i5)') (mpcs%mpc_dof(j),j=n0,n1)
      write(file,'(30f7.2)') (mpcs%mpc_val(j),j=n0,n1),mpcs%mpc_const(i)
    enddo
  end subroutine
 
  subroutine fstr_scan_contact_state( cstep, sub_step, cont_step, dt, ctAlgo, hecMESH, fstrSOLID, infoCTChange, B )
    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
    !      logical, intent(inout)                 :: changed     !< if contact state changed
    real(kind=kreal), optional             :: b(:)        
    character(len=9)                       :: flag_ctAlgo
    integer(kind=kint) :: i, grpid
    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
      if( present(b) ) then
        call scan_contact_state( flag_ctalgo, fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
           & fstrsolid%QFORCE(:), infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive, mu, b )
      else
        call scan_contact_state( flag_ctalgo, fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
           & fstrsolid%QFORCE(:), infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive, mu )
      endif
      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, mu )
      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 )
 
    !for output contact state
    do i=1,fstrsolid%n_contacts
      call set_contact_state_vector( fstrsolid%contacts(i), dt, fstrsolid%CONT_RELVEL, fstrsolid%CONT_STATE )
    enddo
 
    infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
    infoctchange%contactNode_previous = infoctchange%contactNode_current
 
    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 remove_duplication_tiedcontact( 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, j, grpid, slave
    integer(kind=kint) :: k, id, iSS
    integer(kind=kint) :: ig0, ig, iS0, iE0
    integer(kind=kint), allocatable :: states(:)
 
    allocate(states(hecmesh%n_node))
    states(:) = contactfree
 
    ! if a boundary condition is given, the slave
    do ig0= 1, fstrsolid%BOUNDARY_ngrp_tot
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
      ig= fstrsolid%BOUNDARY_ngrp_ID(ig0)
      is0= hecmesh%node_group%grp_index(ig-1) + 1
      ie0= hecmesh%node_group%grp_index(ig  )
      do k= is0, ie0
        iss = hecmesh%node_group%grp_item(k)
        !states(iSS) = CONTACTSTICK
      enddo
    enddo
 
    do i=1,fstrsolid%n_contacts
      if( fstrsolid%contacts(i)%algtype /= contacttied ) cycle
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      do j=1, size(fstrsolid%contacts(i)%slave)
        if( fstrsolid%contacts(i)%states(j)%state==contactfree ) cycle   ! free
        slave = fstrsolid%contacts(i)%slave(j)
        if( states(slave) == contactfree ) then
          states(slave) = fstrsolid%contacts(i)%states(j)%state
          id = fstrsolid%contacts(i)%states(j)%surface
          do k=1,size( fstrsolid%contacts(i)%master(id)%nodes )
            iss = fstrsolid%contacts(i)%master(id)%nodes(k)
            states(iss) = fstrsolid%contacts(i)%states(j)%state
          enddo
        else !found duplicate tied contact slave node
          fstrsolid%contacts(i)%states(j)%state = contactfree
          infoctchange%free2contact = infoctchange%free2contact - 1
          write(*,'(A,i10,A,i6,A,i6,A)') "Node",hecmesh%global_node_ID(slave), &
            " in rank",hecmw_comm_get_rank()," freed due to duplication"
        endif
      enddo
    enddo
 
  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
    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_exp( fstrsolid%contacts(i), fstrsolid%ddunode(:), fstrsolid%dunode(:), &
           & infoctchange, hecmesh%global_node_ID(:), hecmesh%global_elem_ID(:), is_init, iactive )
 
      infoctchange%active = infoctchange%active .or. iactive
    enddo
 
    infoctchange%contactNode_current = infoctchange%contactNode_previous+infoctchange%free2contact-infoctchange%contact2free
    infoctchange%contactNode_previous = infoctchange%contactNode_current
    fstrsolid%ddunode = 0.d0
  end subroutine
 
  subroutine fstr_update_contact0( hecMESH, fstrSOLID, B )
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    real(kind=kreal), intent(inout)        :: b(:)        
 
    integer(kind=kint) :: i, algtype
 
    do i=1, fstrsolid%n_contacts
      !   if( contacts(i)%mpced ) cycle
      algtype = fstrsolid%contacts(i)%algtype
      if( algtype == contactsslid .or. algtype == contactfslid ) then
        call calcu_contact_force0( fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:)  &
        , fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff, mu, mut, b )
      else if( algtype == contacttied ) then
        call calcu_tied_force0( fstrsolid%contacts(i), fstrsolid%unode(:), fstrsolid%dunode(:), mu, b )
      endif
    enddo
 
    do i=1, fstrsolid%n_embeds
      call calcu_tied_force0( fstrsolid%embeds(i), fstrsolid%unode(:), fstrsolid%dunode(:), mu, b )
    enddo
 
  end subroutine
 
  subroutine fstr_update_contact_multiplier( hecMESH, fstrSOLID, ctchanged )
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    logical, intent(out)                   :: ctchanged
 
    integer(kind=kint) :: i, nc, algtype
 
    gnt = 0.d0;  ctchanged = .false.
    nc = fstrsolid%n_contacts+fstrsolid%n_embeds
    do i=1, fstrsolid%n_contacts
      algtype = fstrsolid%contacts(i)%algtype
      if( algtype == contactsslid .or. algtype == contactfslid ) then
        call update_contact_multiplier( fstrsolid%contacts(i), hecmesh%node(:), fstrsolid%unode(:)  &
        , fstrsolid%dunode(:), fstrsolid%contacts(i)%fcoeff, mu, mut, gnt, ctchanged )
      else if( algtype == contacttied ) then
        call update_tied_multiplier( fstrsolid%contacts(i), fstrsolid%unode(:), fstrsolid%dunode(:), &
                                      &  mu, ctchanged )
      endif
    enddo
 
    do i=1, fstrsolid%n_embeds
      call update_tied_multiplier( fstrsolid%embeds(i), fstrsolid%unode(:), fstrsolid%dunode(:), &
      &  mu, ctchanged )
    enddo
 
    if( nc>0 ) gnt = gnt/nc
  end subroutine
 
  subroutine fstr_ass_load_contactalag( hecMESH, fstrSOLID, B )
    type( hecmwst_local_mesh ), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    real(kind=kreal), intent(inout)        :: b(:)        
 
    integer(kind=kint) :: i, algtype
 
    do i = 1, fstrsolid%n_contacts
      algtype = fstrsolid%contacts(i)%algtype
      if( algtype == contactsslid .or. algtype == contactfslid ) then
        call ass_contact_force( fstrsolid%contacts(i), hecmesh%node, fstrsolid%unode, b )
      else if( algtype == contacttied ) then
        call calcu_tied_force0( fstrsolid%contacts(i), fstrsolid%unode(:), fstrsolid%dunode(:), mu, b )
      endif
    enddo
 
    do i = 1, fstrsolid%n_embeds
      call calcu_tied_force0( fstrsolid%embeds(i), fstrsolid%unode(:), fstrsolid%dunode(:), mu, b )
    enddo
  end subroutine
 
  subroutine fstr_update_contact_tangentforce( fstrSOLID )
    type(fstr_solid), intent(inout)        :: fstrSOLID
 
    integer(kind=kint) :: i
 
    do i=1, fstrsolid%n_contacts
      call update_contact_tangentforce( fstrsolid%contacts(i) )
    enddo
  end subroutine
 
  subroutine fstr_contactbc( cstep, iter, hecMESH, hecMAT, fstrSOLID )
    use fstr_ctrl_modifier
    integer(kind=kint)                       :: cstep
    integer(kind=kint), intent(in)           :: iter
    type (hecmwST_local_mesh), intent(inout) :: hecMESH
    type (hecmwST_matrix), intent(inout)     :: hecMAT
    type(fstr_solid), intent(inout)          :: fstrSOLID
 
    integer(kind=kint), parameter :: NDOF=3
 
    integer(kind=kint) :: i, j, k, m, nnode, nd, etype, grpid
    integer(kind=kint) :: ctsurf, ndLocal(l_max_surface_node+1)
    integer(kind=kint) :: algtype    
    real(kind=kreal) :: factor, elecoord( 3, l_max_surface_node)
    real(kind=kreal) :: stiff(l_max_surface_node*3+3, l_max_surface_node*3+3)
    real(kind=kreal) :: nrlforce, force(l_max_surface_node*3+3)
 
    factor = fstrsolid%FACTOR(2)
 
    do i=1,fstrsolid%n_contacts
 
      grpid = fstrsolid%contacts(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
 
      algtype = fstrsolid%contacts(i)%algtype
 
      do j=1, size(fstrsolid%contacts(i)%slave)
        if( fstrsolid%contacts(i)%states(j)%state==contactfree ) cycle   ! free
        ctsurf = fstrsolid%contacts(i)%states(j)%surface          ! contacting surface
        etype = fstrsolid%contacts(i)%master(ctsurf)%etype
        nnode = size(fstrsolid%contacts(i)%master(ctsurf)%nodes)
        ndlocal(1) = fstrsolid%contacts(i)%slave(j)
        do k=1,nnode
          ndlocal(k+1) = fstrsolid%contacts(i)%master(ctsurf)%nodes(k)
          elecoord(1:3,k)=hecmesh%node(3*ndlocal(k+1)-2:3*ndlocal(k+1))
        enddo
        if( algtype == contactsslid .or. algtype == contactfslid ) then
          call contact2stiff( algtype, fstrsolid%contacts(i)%states(j),    &
          etype, nnode, elecoord(:,:), mu, mut, fstrsolid%contacts(i)%fcoeff,    &
          fstrsolid%contacts(i)%symmetric, stiff(:,:), force(:) )
        else if( algtype == contacttied ) then
          call tied2stiff( algtype, fstrsolid%contacts(i)%states(j),    &
          etype, nnode, mu, mut, stiff(:,:), force(:) )
        endif
        ! ----- CONSTRUCT the GLOBAL MATRIX STARTED
        call hecmw_mat_ass_elem(hecmat, nnode+1, ndlocal, stiff)
 
        if( iter>1 ) cycle
        !  if( fstrSOLID%contacts(i)%states(j)%multiplier(1)/=0.d0 ) cycle
        ! In case of new contact nodes, add enforced disp constraint
        fstrsolid%contacts(i)%states(j)%wkdist = fstrsolid%contacts(i)%states(j)%distance
        nrlforce = -mu*fstrsolid%contacts(i)%states(j)%distance
        force(1:nnode*ndof+ndof) = force(1:nnode*ndof+ndof)*nrlforce
        do m=1,nnode+1
          nd = ndlocal(m)
          do k=1,ndof
            hecmat%B(ndof*(nd-1)+k)=hecmat%B(ndof*(nd-1)+k)-force((m-1)*ndof+k)
          enddo
        enddo
 
      enddo
    enddo
 
    do i=1,fstrsolid%n_embeds
      grpid = fstrsolid%embeds(i)%group
      if( .not. fstr_iscontactactive( fstrsolid, grpid, cstep ) ) cycle
      do j=1, size(fstrsolid%embeds(i)%slave)
        if( fstrsolid%embeds(i)%states(j)%state==contactfree ) cycle   ! free
        ctsurf = fstrsolid%embeds(i)%states(j)%surface          ! contacting surface
        etype = fstrsolid%embeds(i)%master(ctsurf)%etype
        nnode = size(fstrsolid%embeds(i)%master(ctsurf)%nodes)
        ndlocal(1) = fstrsolid%embeds(i)%slave(j)
        do k=1,nnode
          ndlocal(k+1) = fstrsolid%embeds(i)%master(ctsurf)%nodes(k)
          elecoord(1:3,k)=hecmesh%node(3*ndlocal(k+1)-2:3*ndlocal(k+1))
        enddo
        call tied2stiff( algtype, fstrsolid%embeds(i)%states(j),    &
        etype, nnode, mu, mut, stiff(:,:), force(:) )
        ! ----- CONSTRUCT the GLOBAL MATRIX STARTED
        call hecmw_mat_ass_elem(hecmat, nnode+1, ndlocal, stiff)
      enddo
    enddo
 
  end subroutine
 
  subroutine initialize_contact_output_vectors(fstrSOLID,hecMAT)
    type(fstr_solid)       :: fstrsolid
    type(hecmwst_matrix)   :: hecMAT
 
    if( .not. associated(fstrsolid%CONT_NFORCE) ) then
      allocate( fstrsolid%CONT_NFORCE(hecmat%NP*3) )
      fstrsolid%CONT_NFORCE(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_FRIC) ) then
      allocate( fstrsolid%CONT_FRIC(hecmat%NP*3) )
      fstrsolid%CONT_FRIC(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_RELVEL) ) then
      allocate( fstrsolid%CONT_RELVEL(hecmat%NP*3) )
      fstrsolid%CONT_RELVEL(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_STATE) ) then
      allocate( fstrsolid%CONT_STATE(hecmat%NP*1) )
      fstrsolid%CONT_STATE(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_AREA) ) then
      allocate( fstrsolid%CONT_AREA(hecmat%NP) )
      fstrsolid%CONT_AREA(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_NTRAC) ) then
      allocate( fstrsolid%CONT_NTRAC(hecmat%NP*3) )
      fstrsolid%CONT_NTRAC(:) = 0.d0
    end if
 
    if( .not. associated(fstrsolid%CONT_FTRAC) ) then
      allocate( fstrsolid%CONT_FTRAC(hecmat%NP*3) )
      fstrsolid%CONT_FTRAC(:) = 0.d0
    end if
 
  end subroutine
 
  subroutine initialize_embed_vectors(fstrSOLID,hecMAT)
    type(fstr_solid)       :: fstrSOLID
    type(hecmwst_matrix)   :: hecMAT
 
    if( .not. associated(fstrsolid%EMBED_NFORCE) ) then
      allocate( fstrsolid%EMBED_NFORCE(hecmat%NP*3) )
      fstrsolid%EMBED_NFORCE(:) = 0.d0
    end if
  end subroutine
 
  subroutine setup_contact_elesurf_for_area( cstep, hecMESH, fstrSOLID )
    integer(kind=kint), intent(in)               :: cstep
    type( hecmwst_local_mesh ), intent(in)       :: hecMESH
    type(fstr_solid), intent(inout)              :: fstrSOLID
 
    integer(kind=kint) :: i, j, k, sgrp_id, iS, iE, eid, sid, n_cdsurfs
    logical, pointer :: cdef_surf(:,:)
    real(kind=kreal), pointer   :: coord(:)
 
    if( associated(fstrsolid%CONT_SGRP_ID) ) deallocate(fstrsolid%CONT_SGRP_ID)
 
    allocate(cdef_surf(l_max_elem_surf,hecmesh%n_elem))
    cdef_surf(:,:) = .false.
 
    ! label contact defined surfaces
    n_cdsurfs = 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
 
      do k=1,2 !slave,master
        if( k==1 ) then !slave
          sgrp_id = fstrsolid%contacts(i)%surf_id1_sgrp
        else if( k==2 ) then !master
          sgrp_id = fstrsolid%contacts(i)%surf_id2
        end if
 
        if( sgrp_id <= 0 ) cycle
 
        is = hecmesh%surf_group%grp_index(sgrp_id-1) + 1
        ie = hecmesh%surf_group%grp_index(sgrp_id  )
        do j=is,ie
          eid = hecmesh%surf_group%grp_item(2*j-1)
          sid = hecmesh%surf_group%grp_item(2*j)
          ! only internal and boundary element should be added
          if( .not. cdef_surf(sid,eid) ) n_cdsurfs = n_cdsurfs + 1
          cdef_surf(sid,eid) = .true.
        enddo
      end do
    enddo
 
    !gather info of contact defined surfaces
    allocate(fstrsolid%CONT_SGRP_ID(2*n_cdsurfs))
    n_cdsurfs = 0
    do i=1,hecmesh%n_elem
      do j=1,l_max_elem_surf
        if( cdef_surf(j,i) ) then
          n_cdsurfs = n_cdsurfs + 1
          fstrsolid%CONT_SGRP_ID(2*n_cdsurfs-1) = i
          fstrsolid%CONT_SGRP_ID(2*n_cdsurfs  ) = j
        endif
      end do
    end do
    deallocate(cdef_surf)
 
  end subroutine
 
  subroutine calc_contact_area( hecMESH, fstrSOLID, flag )
    type( hecmwst_local_mesh ), intent(in)       :: hecmesh
    type(fstr_solid), intent(inout)              :: fstrSOLID
    integer(kind=kint), intent(in)               :: flag
 
    integer(kind=kint), parameter :: NDOF=3
    integer(kind=kint) :: i, isuf, icel, sid, etype, nn, iS, stype, idx
    integer(kind=kint) :: ndlocal(l_max_elem_node)
    real(kind=kreal), pointer   :: coord(:)
    real(kind=kreal)   :: ecoord(ndof,l_max_elem_node), vect(l_max_elem_node)
 
    fstrsolid%CONT_AREA(:) = 0.d0
 
    if( .not. associated(fstrsolid%CONT_SGRP_ID) ) return
 
    allocate(coord(ndof*hecmesh%n_node))
    do i=1,ndof*hecmesh%n_node
      coord(i) = hecmesh%node(i)+fstrsolid%unode(i)
    end do
    if( flag == 1 ) then
      do i=1,ndof*hecmesh%n_node
        coord(i) = coord(i)+fstrsolid%dunode(i)
      end do
    end if
 
    do isuf=1,size(fstrsolid%CONT_SGRP_ID)/2
      icel = fstrsolid%CONT_SGRP_ID(2*isuf-1)
      sid  = fstrsolid%CONT_SGRP_ID(2*isuf  )
 
      etype = hecmesh%elem_type(icel)
      nn = hecmw_get_max_node(etype)
      is = hecmesh%elem_node_index(icel-1)
      ndlocal(1:nn) = hecmesh%elem_node_item (is+1:is+nn)
 
      do i=1,nn
        idx = ndof*(ndlocal(i)-1)
        ecoord(1:ndof,i) = coord(idx+1:idx+ndof)
      end do
 
      call calc_nodalarea_surfelement( etype, nn, ecoord, sid, vect )
 
      do i=1,nn
        idx = ndlocal(i)
        fstrsolid%CONT_AREA(idx) = fstrsolid%CONT_AREA(idx) + vect(i)
      end do
 
    end do
 
    deallocate(coord)
  end subroutine
 
  subroutine calc_nodalarea_surfelement( etype, nn, ecoord, sid, vect )
    integer(kind=kint), intent(in)   :: etype
    integer(kind=kint), intent(in)   :: nn
    real(kind=kreal), intent(in)     :: ecoord(:,:)
    integer(kind=kint), intent(in)   :: sid
    real(kind=kreal), intent(out)    :: vect(:)
 
    integer(kind=kint), parameter :: NDOF=3
    integer(kind=kint) :: nod(l_max_surface_node)
    integer(kind=kint) :: nsur, stype, ig0, i
    real(kind=kreal)   :: localcoord(2), normal(3), area, wg
    real(kind=kreal)   :: scoord(ndof,l_max_surface_node), h(l_max_surface_node)
 
    vect(:) = 0.d0
 
    call getsubface( etype, sid, stype, nod )
    nsur = getnumberofnodes( stype )
    do i=1,nsur
      scoord(1:ndof,i) = ecoord(1:ndof,nod(i))
    end do
 
    area = 0.d0
    do ig0=1,numofquadpoints( stype )
      call getquadpoint( stype, ig0, localcoord(1:2) )
      call getshapefunc( stype, localcoord(1:2), h(1:nsur) )
 
      wg=getweight( stype, ig0 )
      ! normal = dx/dr_1 \times dx/dr_2
      normal(1:3) = surfacenormal( stype, nsur, localcoord(1:2), scoord(1:ndof,1:nsur) )
      area = area + dsqrt(dot_product(normal,normal))*wg
    enddo
    area = area/dble(nsur)
    do i=1,nsur
      vect(nod(i)) = area
    end do
 
  end subroutine
 
  subroutine fstr_setup_parancon_contactvalue(hecMESH,ndof,vec,vtype)
    use m_fstr
    implicit none
    type(hecmwst_local_mesh), intent(in)      :: hecmesh
    integer(kind=kint), intent(in)            :: ndof
    real(kind=kreal), pointer, intent(inout)  :: vec(:)
    integer(kind=kint), intent(in)            :: vtype !1:value, 2:state
    !
    real(kind=kreal) ::  rhsb
    integer(kind=kint) ::  i,j,n,i0,n_loc,nndof
    integer(kind=kint) :: offset, pid, lid
    integer(kind=kint), allocatable :: displs(:)
    real(kind=kreal), allocatable   :: vec_all(:)
  
    !
    n_loc = hecmesh%nn_internal
    allocate(displs(0:nprocs))
    displs(:) = 0
    displs(myrank+1) = n_loc
    call hecmw_allreduce_i(hecmesh, displs, nprocs+1, hecmw_sum)
    do i=1,nprocs
      displs(i) = displs(i-1) + displs(i)
    end do
    offset = displs(myrank)
    n = displs(nprocs)
  
    allocate(vec_all(ndof*n))
  
    if( vtype == 1 ) then
      vec_all(:) = 0.d0
      do i= hecmesh%nn_internal+1,hecmesh%n_node
        pid = hecmesh%node_ID(i*2)
        lid = hecmesh%node_ID(i*2-1)
        i0 = (displs(pid) + (lid-1))*ndof
        vec_all(i0+1:i0+ndof) = vec((i-1)*ndof+1:i*ndof)
        vec((i-1)*ndof+1:i*ndof) = 0.d0
      enddo
  
      call hecmw_allreduce_r(hecmesh, vec_all, n*ndof, hecmw_sum)
  
      do i=1,ndof*n_loc
        vec(i) = vec(i) + vec_all(offset*ndof+i)
      end do
    else if( vtype == 2 ) then
      vec_all(:) = -1000.d0
      do i= hecmesh%nn_internal+1,hecmesh%n_node
        if( vec(i) == 0.d0 ) cycle
        pid = hecmesh%node_ID(i*2)
        lid = hecmesh%node_ID(i*2-1)
        i0 = displs(pid) + lid
        vec_all(i0) = vec(i)
      enddo
  
      call hecmw_allreduce_r(hecmesh, vec_all, n, hecmw_max)
  
      do i=1,n_loc
        if( vec_all(offset+i) == -1000.d0 ) cycle
        if( vec(i) < vec_all(offset+i) ) vec(i) = vec_all(offset+i)
      end do
    end if
  
    deallocate(displs,vec_all)
    end subroutine
 
end module mcontact