*
* User subroutine VUINTERACTION example to model a uniform thickness interface
* bonded to both surfaces. 
* The interface material has "uniaxial" plasticity in the normal direction 
* with linear hardening; the shear behavior remains elastic. 
* Membrane straining of the interface does not affect the stress transmitted 
* across the interface. 
* Simple heat transfer across the interface is modeled using a conductance 
* that is independent of gap distance or pressure. 
* Heat generation at the interface is not considered.
*
      subroutine vuinteraction( 
* Read/Write - 
     *   stress, fluxSec, fluxMain, 
     *   state, sed, 
* Write only -
     *   sfd, scd, spd, svd,     
* Read only - 
     *   nBlock, nBlockAnal, nBlockEdge, 
     *   nNodState, nNodSec, nNodMain, nDir, 
     *   nStates, nProps, nTemp, nFields, 
     *   jFlags, rData, 
     *   surfInt, surfSec, surfMain, 
     *   jSecUid, jMainUid, props,
     *   penetration, drDisp, dRot, dircos, 
     *   stiffDef, conductDef, 
     *   coordSec, coordMain, areaProx, shapeSec, shapeMain, 
     *   tempSec, tempMain, dTempSec, dTempMain, 
     *   fieldSec, fieldMain, dFieldSec, dFieldMain )
c
c
c Surface dependent variables
c
c      nBlockAnal = 1         (analytical rigid main surface)
c      nBlockAnal = nBlock    (non-analytical-rigid main surface)
c      nBlockEdge = 1         (non-edge-type secondary surface)
c      nBlockEdge = nBlock    (edge-type secondary surface)
c      nNodState = 1          (node-type secondary surface)
c      nNodState = 4          (edge-type secondary surface)
c      nNodSec = 1            (node-type secondary surface)
c      nNodSec = 2            (edge-type secondary surface)
c      nNodMain = 1            (analytical rigid main surface)
c      nNodMain = 2            (edge-type main surface)
c      nNodMain = 4            (facet-type main surface)
c
c Surface names surfSec and surfMain are not available for
c general contact (set to blank).
c
      include 'vaba_param.inc'
c
      dimension stress(nDir,nBlock),
     *   fluxSec(nBlock),
     *   fluxMain(nBlock),
     *   state(nStates,nNodState,nBlock),
     *   sed(nBlock), 
     *   sfd(nBlock),
     *   scd(nBlock), 
     *   spd(nBlock), 
     *   svd(nBlock),
     *   jSecUid(nNodSec,nBlock), 
     *   jMainUid(nNodMain,nBlockAnal), 
     *   props(nProps), 
     *   penetration(nBlock), 
     *   drDisp(nDir,nBlock), 
     *   dRot(2,2,nBlock), 
     *   dircos(nDir,nDir,nBlock), 
     *   stiffDef(nBlock), 
     *   conductDef(nBlock),
     *   coordSec(nDir,nNodSec,nBlock), 
     *   coordMain(nDir,nNodMain,nBlockAnal), 
     *   areaProx(nBlock),
     *   shapeSec(nNodSec,nBlockEdge), 
     *   shapeMain(nNodMain,nBlockAnal), 
     *   tempSec(nBlock), 
     *   tempMain(nBlockAnal),
     *   dTempSec(nBlock), 
     *   dTempMain(nBlockAnal),
     *   fieldSec(nFields,nBlock),
     *   fieldMain(nFields,nBlockAnal), 
     *   dFieldSec(nFields,nBlock),
     *   dFieldMain(nFields,nBlockAnal)
c
      parameter( iKStep    = 1,
     *           iKInc     = 2,
     *           iLConType = 3,
     *           nFlags    = 3 )
c
      parameter( iTimStep    = 1,
     *           iTimGlb     = 2,
     *           iDTimCur    = 3,
     *           iTrackThick = 4, 
     *           nData       = 4 )
c
      dimension jFlags(nFlags), rData(nData)
      character*80 surfInt, surfSec, surfMain
      parameter( zero=0.d0, half=0.5d0, one=1.d0 )
c
c     Parameters for the properties array (nProps=7) are taken from vuinter3d_n.
c     For the current test i_prp_GapCutOff is not used.
c
      parameter ( i_prp_GapCutOff = 1,
     *            i_prp_YoungsModulus = 2,
     *            i_prp_PoissonsRatio = 3,
     *            i_prp_InitYield = 4,
     *            i_prp_HardenMod = 5,
     *            i_prp_ThickInter = 6,
     *            i_prp_IfcCond = 7 )
c
c     Use state(1,1,k) to store the current yield stress for odd increments
c     Use state(2,1,k) to store the current yield stress for even increments
c
c     At increment k: retieve old state with "state(iGet,1,k)"
c                     store   new state with "state(iPut,1,k)"
c
      kInc = jFlags(iKInc)
      iGet = mod( kInc  , 2 ) + 1
      iPut = mod( kInc+1, 2 ) + 1
c
c     Initialize states to initial yield stress
c
      if (kInc .eq. 0) then
         do i = 1, nBlock
            state(1,1,i) = props(i_prp_InitYield)
            state(2,1,i) = props(i_prp_InitYield)
         end do
      end if
c
c     Compute interface properties
c
      E    = props(i_prp_YoungsModulus)
      ET   = props(i_prp_HardenMod)
      EH   = E * ET / (E - ET)
      G    = half * E / (one + props(i_prp_PoissonsRatio))
      hInv = one / props(i_prp_ThickInter)
c
      do k = 1, nBlock
c
         dE11 =  drDisp(1,k) * hInv
         dE12 = -drDisp(2,k) * hInv
         dE13 = -drDisp(3,k) * hInv
c
c        Compute elastic trial stress
c
         stress(1,k) = stress(1,k) + E*dE11
         stress(2,k) = stress(2,k) + G*dE12
         stress(3,k) = stress(3,k) + G*dE13
c
c        Determine how much to scale S11 due to yielding
c
         sTrial = abs( stress(1,k) )
         sYield = state(iGet,1,k)
         if( sTrial .gt. sYield ) then
c
c           Compute plastic strain increment
            dEP11 = ( sTrial - sYield ) * ( one/ET - one/E )
c
c           Update yield stress
            sYield  = sYield + dEP11 * EH
c
c           Compute normal stress and store the new state
            stress(1,k) = sign( sYield, stress(1,k) )
            state(iPut,1,k) = sYield
         end if
c
      end do
c
c     Heat flux convention: positive for flux going into a surface
c
      if( nTemp .gt. 0 ) then
c
c        Compute interface conductance
         Cond = props(i_prp_IfcCond) * hInv
c
         do k = 1, nBlock
            fluxSec(k) = Cond*(tempMain(k) - tempSec(k))
            fluxMain(k) = -fluxSec(k)
         end do
      else
         do k = 1, nBlock
            fluxSec(k) = zero
            fluxMain(k) = zero
         end do
      end if
c
      return
      end