! This is a VUMAT subroutine implementing a continuum damage mechanics
! framework for composite materials in Abaqus.
! Copyright (C) 2021 Rutger Wouter Kok
! Modified by Tianyu Zhang


      subroutine vumat(
        ! Read only - input variables from Abaqus
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
        ! Write only - outputs Abaqus needs from subroutine
     7  stressNew, stateNew, enerInternNew, enerInelasNew)

        include 'vaba_param.inc'
        
        ! Dimension the Abaqus input variables
        dimension props(nprops), density(nblock), coordMp(nblock),
     1    charLength(nblock), strainInc(nblock, ndir+nshr),
     2    relSpinInc(nblock, nshr), tempOld(nblock),
     3    stretchOld(nblock, ndir+nshr), fieldOld(nblock, nfieldv),
     4    defgradOld(nblock,ndir+nshr+nshr), stateOld(nblock, nstatev),
     5    stressOld(nblock, ndir+nshr), enerInternOld(nblock),
     6    enerInelasOld(nblock), tempNew(nblock),
     7    stretchNew(nblock, ndir+nshr), fieldNew(nblock, nfieldv),
     8    defgradNew(nblock,ndir+nshr+nshr), enerInelasNew(nblock),
     9    stressNew(nblock,ndir+nshr), stateNew(nblock, nstatev),
     1    enerInternNew(nblock)     
        character*80 cmname
  
        
      ! Declare variables used in main part of script
      real*8, dimension(6) :: stress, strain
      real*8, dimension(6,6) :: C
      real*8 :: e11,e22,e33,nu12,nu13,nu23,g12,g13,g23
      real*8 :: nu21,nu31,nu32,psi
      integer :: k,i                                                  

      ! Elastic constants orthotropic ply
      e11 = props(1) ! stiffness fiber 
      e22 = props(2) ! stiffness transverse  (in-plane)
      e33 = props(3) ! stiffness transverse  (out-of-plane)
      nu12 = props(4) ! Poisson's ratio 12 
      nu13 = props(5) ! Poisson's ratio 13 
      nu23 = props(6) ! Poisson's ratio 23 
      ! shear moduli multiplied by two (tensorial-engineering strain)
      g12 = 2.0d0*props(7) ! shear modulus 12  
      g13 = 2.0d0*props(8) ! shear modulus 13 
      g23 = 2.0d0*props(9) ! shear modulus 23 

      ! Calculate remaining Poisson's ratios
      nu21 = nu12*(e22/e11) ! Poisson's ratio 21
      nu31 = nu13*(e33/e11) ! Poisson's ratio 31
      nu32 = nu23*(e33/e22) ! Poisson's ratio 32

      ! Calculate undamaged stiffness matrix
      psi = (1.0d0 - nu12*nu21 - nu23*nu32 - nu31*nu13 -
     1       2.0d0*nu12*nu23*nu31)/(e11*e22*e33)

      C = 0.0d0 ! initialize stiffness matrix, set all elements = 0
      C(1,1) = (1.0d0 - nu23*nu32)/(e22*e33*psi)
      C(1,2) = (nu21 + nu31*nu23)/(e22*e33*psi)
      C(1,3) = (nu31 + nu21*nu32)/(e22*e33*psi)
      C(2,1) = (nu12 + nu13*nu32)/(e11*e33*psi)
      C(2,2) = (1.0d0 - nu31*nu13)/(e11*e33*psi)
      C(2,3) = (nu32 + nu31*nu12)/(e11*e33*psi)
      C(3,1) = (nu13 + nu12*nu23)/(e11*e22*psi)
      C(3,2) = (nu23 + nu13*nu21)/(e11*e22*psi)
      C(3,3) = (1.0d0 - nu12*nu21)/(e11*e22*psi)
      C(4,4) = g12
      C(5,5) = g23
      C(6,6) = g13

      ! Initial elastic step, for Abaqus tests 
      if (stepTime == 0) then
          do k = 1, nblock
              ! Initialize state variables
              stateNew(k,:) = 0.d0
              ! Initial strain
              strain = strainInc(k,:)
              ! Calculate initial stress (multiply stiffness and strain)
              stress = matmul(C, strain)
              ! Assign calculated stresses to stressNew array
              stressNew(k,:) = stress
          end do
      else
          ! If not initial step, calc. stresses according to CDM model
          do k = 1,nblock
              ! Calc. new strain and assign to state variable array (1-6)
              stateNew(k,1:6) = stateOld(k,1:6) + strainInc(k,:)
          
              ! Assign strain values to strain array
              strain = stateNew(k,1:6)

              ! Update stressNew and stateNew
              stressNew(k,:) = matmul(C,strain)

              ! Calculate internal energy
              stressPower = 0.5d0*((stressNew(k,1)+stressOld(k,1))*
     1                        strainInc(k,1)+(stressNew(k,2)+
     2                        stressOld(k,2))*strainInc(k,2)+
     3                        (stressNew(k,3)+stressOld(k,3))*
     4                        strainInc(k,3)+2.0d0*(stressNew(k,4)+
     5                        stressOld(k,4))*strainInc(k,4)+
     6                        2.0d0*(stressNew(k,5)+stressOld(k,5))*
     7                        strainInc(k,5)+2.0d0*(stressNew(k,6)+
     8                        stressOld(k,6))*strainInc(k,6))
              ! Assign energies to enerInternNew array
              enerInternNew(k) = enerInternOld(k)+stressPower/density(k)
          end do
      end if
          
      close(95)
      end subroutine vumat
*======================================================================================================================