*
* 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, fluxSlv, fluxMst, 
     *   state, sed, 
* Write only -
     *   sfd, scd, spd, svd,     
* Read only - 
     *   nBlock, nBlockAnal, nBlockEdge, 
     *   nNodState, nNodSlv, nNodMst, nDir, 
     *   nStates, nProps, nTemp, nFields, 
     *   jFlags, rData, 
     *   surfInt, surfSlv, surfMst, 
     *   jSlvUid, jMstUid, props,
     *   penetration, drDisp, dRot, dircos, 
     *   stiffDef, conductDef, 
     *   coordSlv, coordMst, areaProx, shapeSlv, shapeMst, 
     *   tempSlv, tempMst, dTempSlv, dTempMst, 
     *   fieldSlv, fieldMst, dFieldSlv, dFieldMst )
c
c
c Surface dependent variables
c
c      nBlockAnal = 1         (analytical rigid master surface)
c      nBlockAnal = nBlock    (non-analytical-rigid master surface)
c      nBlockEdge = 1         (non-edge-type slave surface)
c      nBlockEdge = nBlock    (edge-type slave surface)
c      nNodState = 1          (node-type slave surface)
c      nNodState = 4          (edge-type slave surface)
c      nNodSlv = 1            (node-type slave surface)
c      nNodSlv = 2            (edge-type slave surface)
c      nNodMst = 1            (analytical rigid master surface)
c      nNodMst = 2            (edge-type master surface)
c      nNodMst = 4            (facet-type master surface)
c
c Surface names surfSlv and surfMst are not available for
c general contact (set to blank).
c
      include 'vaba_param.inc'
c
      dimension stress(nDir,nBlock),
     *   fluxSlv(nBlock),
     *   fluxMst(nBlock),
     *   state(nStates,nNodState,nBlock),
     *   sed(nBlock), 
     *   sfd(nBlock),
     *   scd(nBlock), 
     *   spd(nBlock), 
     *   svd(nBlock),
     *   jSlvUid(nNodSlv,nBlock), 
     *   jMstUid(nNodMst,nBlockAnal), 
     *   props(nProps), 
     *   penetration(nBlock), 
     *   drDisp(nDir,nBlock), 
     *   dRot(2,2,nBlock), 
     *   dircos(nDir,nDir,nBlock), 
     *   stiffDef(nBlock), 
     *   conductDef(nBlock),
     *   coordSlv(nDir,nNodSlv,nBlock), 
     *   coordMst(nDir,nNodMst,nBlockAnal), 
     *   areaProx(nBlock),
     *   shapeSlv(nNodSlv,nBlockEdge), 
     *   shapeMst(nNodMst,nBlockAnal), 
     *   tempSlv(nBlock), 
     *   tempMst(nBlockAnal),
     *   dTempSlv(nBlock), 
     *   dTempMst(nBlockAnal),
     *   fieldSlv(nFields,nBlock),
     *   fieldMst(nFields,nBlockAnal), 
     *   dFieldSlv(nFields,nBlock),
     *   dFieldMst(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, surfSlv, surfMst
      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
            fluxSlv(k) = Cond*(tempMst(k) - tempSlv(k))
            fluxMst(k) = -fluxSlv(k)
         end do
      else
         do k = 1, nBlock
            fluxSlv(k) = zero
            fluxMst(k) = zero
         end do
      end if
c
      return
      end