C
C User subroutine vugens
      subroutine vugens ( 
C Read only -
     *     nblock, ndir, nshr, nstatev, nfieldv,
     *     jElem, kIntPt, kIncr, 
     *     stepTime, totalTime, dt,
     *     cElsetName, coordMp,
     *     strainInc, curvInc, 
     *     tempOld, fieldOld, 
     *     defgradOld, direct,
     *     forceOld, xmomOld, stateOld, 
     *     enerInternOld,  enerInelasOld,
     *     props, nprops, jprops, njprops,
     *     density, thick,
     *     tempNew, fieldNew, 
     *     defgradNew, curvNew, gamma,
C Write only -
     *     ddsdde, 
     *     forceNew, xmomNew, stateNew, 
     *     enerInternNew, enerInelasNew )
*
      include 'vaba_param.inc'
*
      character*80 cElsetName
*
      dimension forceNew(nblock,ndir+nshr),
     *     forceOld(nblock,ndir+nshr),
     *     xmomNew(nblock,ndir+nshr),
     *     xmomOld(nblock,ndir+nshr),
     *     stateOld(nblock,nstatev),
     *     stateNew(nblock,nstatev),
     *     coordMp(nblock,3),
     *     strainInc(nblock,ndir+1+nshr),
     *     curvInc(nblock,ndir+nshr),
     *     ddsdde(nblock,*),
     *     defgradNew(nblock,ndir+1+nshr+nshr),
     *     defgradOld(nblock,ndir+1+nshr+nshr),
     *     direct(nblock,3,3),
     *     gamma(nblock,2),
     *     tempNew(nblock), tempOld(nblock),
     *     fieldNew(nblock,nfieldv),
     *     fieldOld(nblock,nfieldv),
     *     props(nprops),
     *     density(nblock),
     *     enerInternNew(nblock),  
     *     enerInternOld(nblock),
     *     enerInelasOld(nblock),
     *     enerInelasNew(nblock)
 
*
      parameter ( zero = 0.d0, half =0.5d0, two = 2.d0 )
      parameter ( tfail = 1.0d12 )
*
      D11 = props(1)
      D12 = props(2) 
      D22 = props(3)
      D13 = props(4)
      D23 = props(5)
      D33 = props(6)
      D14 = props(7)
      D24 = props(8)
      D34 = props(9)
      D44 = props(10)
      D15 = props(11)
      D25 = props(12)
      D35 = props(13)
      D45 = props(14)
      D55 = props(15)
      D16 = props(16)
      D26 = props(17)
      D36 = props(18)
      D46 = props(19)
      D56 = props(20)
      D66 = props(21)
*     
      ldebug = 0
      if ( ldebug.eq.1 .and. dt.eq. totalTime) then
         write(*,*) '*** VUGENS: Elset Name ',cElsetName
         write(*,*) ' direct = ', ((direct(1,i,j), i=1,3), j=1,3)
      end if
*
      do k= 1, nblock
*
         strainInc(k,3) = zero 
*
         de11 = strainInc(k,1)
         de22 = strainInc(k,2)
         de12 = two*strainInc(k,4)
*
         dc11 = curvInc(k,1)
         dc22 = curvInc(k,2)
         dc12 = curvInc(k,3)
*
         forceNew(k,1) = forceOld(k,1) 
     *        + D11 * de11
     *        + D12 * de22
     *        + D13 * de12
     *        + D14 * dc11
     *        + D15 * dc22
     *        + D16 * dc12
*
         forceNew(k,2) = forceOld(k,2) 
     *        + D12 * de11
     *        + D22 * de22
     *        + D23 * de12
     *        + D24 * dc11
     *        + D25 * dc22
     *        + D26 * dc12
*
         forceNew(k,3) = forceOld(k,3) 
     *        + D13 * de11
     *        + D23 * de22
     *        + D33 * de12
     *        + D34 * dc11
     *        + D35 * dc22
     *        + D36 * dc12
*
         xmomNew(k,1) = xmomOld(k,1)
     *        + D14 * de11
     *        + D24 * de22
     *        + D34 * de12
     *        + D44 * dc11
     *        + D45 * dc22
     *        + D46 * dc12
*
         xmomNew(k,2) = xmomOld(k,2)
     *        + D15 * de11
     *        + D25 * de22
     *        + D35 * de12
     *        + D45 * dc11
     *        + D55 * dc22
     *        + D56 * dc12
*
         xmomNew(k,3) = xmomOld(k,3)
     *        + D16 * de11
     *        + D26 * de22
     *        + D36 * de12
     *        + D46 * dc11
     *        + D56 * dc22
     *        + D66 * dc12
*
         ddsdde(k, 1) = D11 
         ddsdde(k, 2) = D12 
         ddsdde(k, 3) = D22 
         ddsdde(k, 4) = D13 
         ddsdde(k, 5) = D23 
         ddsdde(k, 6) = D33 
         ddsdde(k, 7) = D14 
         ddsdde(k, 8) = D24 
         ddsdde(k, 9) = D34 
         ddsdde(k,10) = D44 
         ddsdde(k,11) = D15 
         ddsdde(k,12) = D25 
         ddsdde(k,13) = D35 
         ddsdde(k,14) = D45 
         ddsdde(k,15) = D55 
         ddsdde(k,16) = D16 
         ddsdde(k,17) = D26 
         ddsdde(k,18) = D36 
         ddsdde(k,19) = D46 
         ddsdde(k,20) = D56 
         ddsdde(k,21) = D66 
*
         enerInternNew(k) = enerInternOld(k) + half * (
     *        +  (forceNew(k,1)+forceOld(k,1))*de11
     *        +  (forceNew(k,2)+forceOld(k,2))*de22
     *        +  (forceNew(k,3)+forceOld(k,3))*de12
     *        +  (xmomNew(k,1)+xmomOld(k,1))*dc11
     *        +  (xmomNew(k,2)+xmomOld(k,2))*dc22
     *        +  (xmomNew(k,3)+xmomOld(k,3))*dc12 ) / density(k)
*
      end do
      if ( nstatev .ge. 7 ) then
         do k = 1, nblock
            stateNew(k,1)=forceNew(k,1)
            stateNew(k,2)=forceNew(k,2)
            stateNew(k,3)=forceNew(k,3)
            stateNew(k,4)=xmomNew(k,1)
            stateNew(k,5)=xmomNew(k,2)
            stateNew(k,6)=xmomNew(k,3)
         end do
         if ( totalTime .ge. tfail ) then
            do k = 1, nblock
               stateNew(k,7)= zero
            end do
         end if
      end if
*     
      return
      end