c
c User subroutine VUINTERACTION to define surface interaction
c
      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, areaSlv, 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), 
     *   areaSlv(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     Indices to user-defined properties (nprops=7):
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     Elastic stresses
c
      G = half * props(i_prp_YoungsModulus) /
     *    (one + props(i_prp_PoissonsRatio))
c
      do k = 1, nBlock
c
         if( penetration(k) .gt. zero ) then
c
            E = stiffDef(k) * props(i_prp_ThickInter)
            dE1 =  drDisp(1,k) / props(i_prp_ThickInter)
            dE2 = -drDisp(2,k) / props(i_prp_ThickInter)
            dE3 = -drDisp(3,k) / props(i_prp_ThickInter)
            stress(1,k) = stress(1,k) + E*dE1
            stress(2,k) = stress(2,k) + G*dE2
            stress(3,k) = stress(3,k) + G*dE3
c
         else
c
            stress(1,k) = zero
            stress(2,k) = zero
            stress(3,k) = zero
c
         end if
c
      end do
c
c     Heat flux convention: positive for flux goibng into a surface
c
      if( nTemp .gt. 0 ) then
         do k = 1, nBlock
            fluxSlv(k) = props(i_prp_IfcCond)*(tempMst(k) - tempSlv(k))
     *                 / props(i_prp_ThickInter)
            fluxMst(k) = -fluxSlv(k)
         end do
      else
         do k = 1, nBlock
            fluxSlv(k) = zero
            fluxMst(k) = zero
         end do
      end if
c
      return
      end