

!=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~=~
!                                          numgeo
!                              Copyright (C) 2017 Jan Machacek
!=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~=~
!
! SUBROUTINE: user_contact_properties
!
!> author: Patrick Staubach, patrick.staubach@yahoo.de
!> date: 28.04.2020
!
!> Interface to user subroutine to define contact model properties 
!
!>### History
!> 28.04.2020, P. Staubach - Initial version
!> 29.06.2020, P. Staubach - Added arguments to subroutine
!
!>@todo
!> After program is completed, this subroutine should be excluded again and made accessible to the
!> main program numgeo.exe as a dynamic linked library (*.so). This should be done to prevent excessive 
!> recompilation before every calculation with changed initial conditions.
!>@endtodo
!=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~~=~=~
subroutine user_contact_properties(istep,node,slave,nprops,interaction_type,step_time,coords,coords_connected,disp,disp_connected,props) 
								   ! bind(c,name='user_contact_properties')
  use, intrinsic :: iso_c_binding
  implicit none
  !DEC$ ATTRIBUTES DLLEXPORT, STDCALL, REFERENCE :: user_contact_properties
  integer(c_int)                    , intent(in)    :: istep               !! step number
  integer(c_int)                    , intent(in)    :: node                !! node label
  logical                           , intent(in)    :: slave               !! true if current node belongs to slave surface
  integer(c_int)                    , intent(in)    :: nprops              !! number of contact properties
  character(c_char)                 , intent(in)    :: interaction_type(*) !! constitutive contact name
  real(c_double)                    , intent(in)    :: step_time           !! current step time cdb%step(istep)%time%is
  real(c_double), dimension(3)      , intent(in)    :: coords              !! reference coordinates of the contact point
  real(c_double), dimension(3)      , intent(in)    :: coords_connected    !! reference coordinates of the connected contact point
  real(c_double), dimension(3)      , intent(in)    :: disp                !! displacement of the contact point
  real(c_double), dimension(3)      , intent(in)    :: disp_connected      !! displacement of the connected contact point
  real(c_double), dimension(nprops) , intent(inout) :: props               !! properties of this contact
  
  real(8) :: yref1, yref2, y, y_connected, yref1_connected

  ! user coding to define the contact properties props
  props(:) = 0.0d0
  props(1) = 30d4
  
  yref2 = 0.0
  
  ! The "real" tip ends at approx 0.027 made
  ! we therfore neglect friction below 0.03 m (accounting for some rounding of the tip)
  ! before it was 0.05 which is the end of the tip
  yref1 = -0.030 + disp(2)
  
  y = coords(2) + disp(2)
  y_connected = coords_connected(2) + disp_connected(2)
  yref1_connected = -0.030 + disp_connected(2)
  
  if(slave .eqv. .false.) then
    if( y >= yref1 .and. y <= yref2 ) then
	  ! friction coefficient corresponds to tan(20.88°)
      props(2) = 0.381
    endif
  else
    if( y_connected >= yref1_connected ) then
      props(2) = 0.381
    endif
  endif
  
end subroutine user_contact_properties


