c
c User subroutine VUCHARLENGTH for user-defined element characteristic length
c
      subroutine vucharlength(
c Read only -
     *     nblock, nfieldv, nprops, ncomp, ndim, nnode, nstatev,
     *     kSecPt, kLayer, kIntPt, jElType, jElem,
     *     totalTime, stepTime, dt,
     *     cmname, coordMp, coordNode, direct, T, props, 
     *     field, stateOld,
c Write only -
     *     charLength )
*
      include 'vaba_param.inc'
*     
      dimension jElType(3), jElem(nblock), coordMp(nblock, ndim), 
     *     coordNode(nblock,nnode,ndim),
     *     direct(nblock, 3, 3), T(nblock,3,3), props(nprops),
     *     stateOld(nblock, nstatev), charLength(nblock,ncomp), 
     *     field(nblock, nfieldv)
      character*80 cmname
*      
*T2D2, T3D2, SAX1
      if( jElType(1) .eq. 1 ) then
         do k = 1, nblock
            charLength(k,1) = props(1) * 
     *           abs (coordNode(k,2,2) - coordNode(k,1,2) )
         end do
      end if
* 
*S4R, S4RS, S4RSW, S4, CPS4R, CPE4R, CAX4R, M3D4R, M3D4
      if( jElType(1) .eq. 3 ) then
         do k = 1, nblock
            diagonal_1 = sqrt( (coordNode(k,1,1)-coordNode(k,3,1) )**2 +
     *           (coordNode(k,1,2)-coordNode(k,3,2))**2 )
            diagonal_2 = sqrt( (coordNode(k,2,1)-coordNode(k,4,1) )**2 +
     *           (coordNode(k,2,2)-coordNode(k,4,2))**2 )
            charLength(k,1)=max(diagonal_1,diagonal_2)
         end do
      end if   
*
*C3D8R, C3D8, SC8R
      if( jElType(1) .eq. 6 ) then
         do k = 1, nblock
            diagonal_1 = sqrt( (coordNode(k,1,1)-coordNode(k,7,1) )**2 +
     *           (coordNode(k,1,2)-coordNode(k,7,2))**2 +
     *           (coordNode(k,1,3)-coordNode(k,7,3))**2 )
            diagonal_2 = sqrt( (coordNode(k,2,1)-coordNode(k,8,1) )**2 +
     *           (coordNode(k,2,2)-coordNode(k,8,2))**2 +
     *           (coordNode(k,2,3)-coordNode(k,8,3))**2 )
            diagonal_3 = sqrt( (coordNode(k,3,1)-coordNode(k,5,1) )**2 +
     *           (coordNode(k,3,2)-coordNode(k,5,2))**2 +
     *           (coordNode(k,3,3)-coordNode(k,5,3))**2 )
            diagonal_4 = sqrt( (coordNode(k,4,1)-coordNode(k,6,1) )**2 +
     *           (coordNode(k,4,2)-coordNode(k,6,2))**2 +
     *           (coordNode(k,4,3)-coordNode(k,6,3))**2 )
            charLength(k,1)=max(diagonal_1,diagonal_2,diagonal_3,
     *           diagonal_4)
         end do
      end if   
*
      return
      end
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), 
     *          charLength(nblock),
     *          stateOld(nblock, nstatev), 
     *          stateNew(nblock, nstatev),
     *          field(nblock, nfieldv)
      character*80 cmname
c Properties array
c     props(1) -> eqpsCrit
c     props(2) -> failDisp
c     props(3) -> Dmax
c
      character*3 cData(maxblk)
      dimension jData(maxblk)
      dimension eqps(maxblk)
c
      parameter ( zero = 0.d0 )
c
c Multi-process information (if needed)
      call vgetnumcpus( numProcesses )
      call vgetrank( kProcessNum )
c
c Read properties
      eqpsCrit = props(1)
      failDisp = props(2)
      Dmax = props(3)
c Get stresses and strains from previous increment
      jStatus = 1
      call vgetvrm( 'PEEQ', eqps, jData, cData, jStatus )
c
      do k = 1, nblock
         eqpsOld = stateOld(k,1)
         damageOld = stateOld(k,2)
         damageNew = damageOld
         if ( eqps(k) .gt. eqpsCrit ) then
            deqps = eqps(k) - eqpsOld
            if ( damageOld .eq. zero ) ! First time criterion is met
     *           deqps = eqps(k) - eqpsCrit
            dUeq = charLength(k) * deqps 
            damageNew = damageOld + dUeq / failDisp
         end if
         if ( damageNew .ge. Dmax ) then
c     Element deletion
            stateNew(k,3) = zero
            damageNew = Dmax
         end if
         stateNew(k,1) = eqps(k)
         stateNew(k,2) = damageNew
         field(k,1) = damageNew
         stateNew(k,4) = charLength(k)
      end do
c
      return
      end