c
c User subroutine VUSDFLD for user-defined fields
c
      subroutine vusdfld(
c Read only -
     *   nblock, nstatev, nfieldv, nprops, ndir, nshr, 
     *   jElemUid, kIntPt, kLayer, kSecPt, 
     *   stepTime, totalTime, dt, cmname, 
     *   coordMp, direct, T, charLength, props, 
     *   stateOld, 
c Write only -
     *   stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension props(nprops),
     *          jElemUid(nblock), coordMp(nblock, *), 
     *          direct(nblock, 3, 3), T(nblock,3,3), 
     *          stateOld(nblock, nstatev), 
     *          stateNew(nblock, nstatev),
     *          field(nblock, nfieldv)
      character*80 cmname
c Properties array
c     props(1) -> Transverse tensile strength, Yt
c     props(2) -> Matreix compressive strength, Yc
c     props(3) -> Ply shear strength, Sc
c     props(4) -> Fiber buckling strength, Xc
c     props(5) -> Initial shear modulus, G12
c     props(6) -> Nonlinear shear factor, alpha
c
      character*3 cData(maxblk*6)
      dimension jData(maxblk*6)
      dimension stress(maxblk*6),strain(maxblk*6)
c Read properties
      yt    = props(1) 
      yc    = props(2)
      sc    = props(3)
      xc    = props(4)
      g12   = props(5) 
      alpha = props(6)
c Get stresses and strains from previous increment
      jStatus = 1
      call vgetvrm( 'S', stress, jData, cData, jStatus )
      jStatus = 1
      call vgetvrm( 'LE', strain, jData, cData, jStatus )
c
      call evaluateDamage( nblock, nstatev,
     *     nfieldv, ndir, nshr, 
     *     yt, yc, sc, xc, g12, alpha,
     *     stress, strain,
     *     stateOld, 
     *     stateNew, field )
c
      return 
      end
c
      subroutine evaluateDamage ( nblock, nstatev, 
     *     nfieldv, ndir, nshr, 
     *     yt, yc, sc, xc, g12, alpha,
     *     stress, strain,
     *     stateOld, 
     *     stateNew, field )
c
      include 'vaba_param.inc'
c
      dimension stress(nblock,ndir+nshr),
     *     strain(nblock,ndir+nshr),
     *     stateOld(nblock,nstatev), 
     *     stateNew(nblock,nstatev),
     *     field(nblock,nfieldv)
c
c initialize failure flags from statev. 
      do k = 1, nblock
         stateNew(k,1) = stateOld(k,1)
         stateNew(k,2) = stateOld(k,2)
         stateNew(k,3) = stateOld(k,3)
c
         em     = stateOld(k,1)
         efs    = stateOld(k,2)
         damage = stateOld(k,3)
c     
         s11 = stress(k,1)
         s22 = stress(k,2)
         s12 = stress(k,4)
*
         e12 = 2.0*strain(k,4) !e12 is engineering strain
c
c damage index: = 0 if no strain to prevent divide by zero
c
         damage = 0.d0
         if (e12.ne.0) 
     *        damage = (3.d0*alpha*g12*s12**2-2.d0*alpha*(s12**3)/e12)/ 
     *                 (1.d0+3.d0*alpha*g12*s12**2)
c
         f1 = s12**2/(2.d0*g12) + 0.75d0*alpha*s12**4
         f2 = sc**2 /(2.d0*g12) + 0.75d0*alpha*sc**4
c
c matrix tensile/compressive failure
         if (em .lt. 1.d0) then
            if (s22 .lt. 0.d0) then
               em = sqrt((s22/yc)**2 + f1/f2)
            else 
               em = sqrt((s22/yt)**2 + f1/f2)
            endif
            stateNew(k,1) = em
         endif
c
c fiber-matrix shear failure
         if (efs .lt. 1.d0) then
            if (s11 .lt. 0.d0) then
               efs = sqrt((s11/xc)**2 + f1/f2)
            else
               efs = 0.d0
            endif
            stateNew(k,2) = efs
         endif
c
c state transition diagram
c
c fv1: matrix compr/tens failure
c fv2: fiber/matrix shear failure
c fv3: material damage (shear nonlinearity)
c                            fv1  fv2  fv3      e1  e2  nu12  g12
c
c (0) no failure              0    0    0  -->  e1  e2  nu12  g12
c (1) matrix (compr/tens)     1    0    0  -->  e1   0   0    g12
c (2) fib/mtx shear           0    1    0  -->  e1  e2   0     0
c (3) matrix & f/m shear      1    1    0  -->  e1   0   0     0
c (4) pure damage             0    0    1  -->  e1  e2  nu12   0
c (5) mtrx & damage           1    0    1  -->  e1   0   0     0
c (6) f/m shear & damage      0    1    1  -->  e1  e2   0     0
c (7) mtrx, f/m shr & damage  1    1    1  -->  e1   0   0     0
c
c     update field variables 
c          
         field(k,1) = 0.d0
         field(k,2) = 0.d0
         if (em .ge. 1.d0) then 
            field(k,1) = 1.d0
         end if
         if (efs .ge. 1.d0) then
            field(k,2) = 1.d0
         end if
         field(k,3) = damage
         stateNew(k,3) = field(k,3)
      end do
c
      return
      end