!DEC$ FREEFORM
!===============================================================================
! Material Routine for Abaqus for Mises plasticity with isotropic hardening
! in incremental formulation for usage under small and large deformations
! - required: hardening law as a UHARD routine
! - not supported yet: rate-dependent hardening laws
!
! distributed under CC BY-NC-SA 4.0 license
! (https://creativecommons.org/licenses/by-nc-sa/4.0/)
! Reference: 
!   N. Lange, G. Huetter, B. Kiefer: "An efficient monolithic solution scheme for FE2 problems",
!   https://arxiv.org/abs/2101.01802
!
! Geralf Hütter, 2020-1110: first version
! Nils Lange, 25.11.2020
! - handover of STATEVS to UHARD corrected
! - assigne PNEWDT with a value >1.0 if step was linear or only the half amount
!   of allowed Newton Steps were necessary
! - PROPPS(NTENS) -> PROPS(NPROPS)
! - changed convergence parameters to those of umatmst3
! Nils Lange, 10.12.2020
! - corrected Newton algorithmn
! - corrected handover of STATEVS to UHARD
!===============================================================================

include 'ABQinterface.f90'
include 'UHARD.f'

subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,DRPLDT,&
                STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,&
                NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,&
                CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
    use ABQINTERFACE
    implicit none
    !---------------------------------------------------------------------------------
    ! Declaration of interface of UMAT to call it afterwards (See paragraph 1.1.40 of 
    ! Abaqus User Subroutines Reference Manual).
    !---------------------------------------------------------------------------------
    ! Note that some arguments are already declared in UEL declaration section.
    integer(kind=AbqIK):: NSTATV,NPROPS, & ! Number of solution-dependent state variables
           ! that are associated with UMAT.
           NDI,NSHR,NTENS, & !Dimensions and size of tensors
           LAYER, & ! Layer number (for composite shells and layered solids).
           NOEL,NPT,KSPT,&  ! Element and GP no.,  Section point number within the current layer.
           KINC
    integer(kind=AbqIK),dimension(4):: JSTEP
    real(kind=AbqRK):: SSE, SPD, SCD, & ! Specific elastic strain energy, plastic 
           ! dissipation, and “creep” dissipation, respectively.
           RPL, & ! Volumetric heat generation per unit time.
           DRPLDT, & ! Variation of RPL with respect to the temperature.
           TEMP, & ! Temperature at the start of the increment.
           DTEMP, & ! Increment of temperature.
           CELENT,& ! Characteristic element length.
           DTIME, PNEWDT
    real(kind=AbqRK),dimension(2):: TIME
    real(kind=AbqRK):: STRESS(NTENS),& ! Classical Stress
           ! higher order strain 'K'. 'NMICROMORPHICGRAD' is a part of microstrain Bmat module.
           STRAN(NTENS), & ! Array containing the macroscopic strains at the beginning
           ! of the increment.
           STATEV(NSTATV), & ! Internal state variables
           DSTRAN(NTENS), & ! Array of macrostrain increments.
           DDSDDE(NTENS,NTENS), & ! Partial derivative of stress 'sigma' w.r.t. STRAN.
           DDSDDT(NTENS), & ! Variation of the stress increments with respect to the
           ! temperature.
           DRPLDE(NTENS), & ! Variation of RPL with respect to the strain increments.
           PROPS(NPROPS),& ! Constitutive parameters
           COORDS(3),&
           PREDEF(1), & ! Array of interpolated values of predefined field variables at
           ! this point. _UMAT added to avoid collision with PREDEF array of UEL subroutine.
           DPRED(1), & ! Array of increments of predefined field variables.
           DROT(3,3), & ! Rotation increment matrix.
           DFGRD0(3,3), & ! Array containing the deformation gradient at the beginning of 
           ! the increment.
           DFGRD1(3,3) ! Array containing the deformation gradient at the end of the
           ! increment.
    character*80 CMNAME ! User-defined material name, left justified.
    !---------------------------------------------------------------------------------
    ! Local variables
    !---------------------------------------------------------------------------------
    real(kind=AbqRK):: E, Nu,shearmod,shearmod2,shearmod3,compmod, & !elastic parameters
                       PEEQ,DilatStrain,MeanStrain,StressHydro,StressDevTrial(NTENS),StressMisesTrial, &
                           SYIELD, HARD(3),PEEQRATE,dPEEQ,allowedRes,SdevScaleFactor,coef1,coef2,R
    integer(kind=AbqIK)::n_iter,i,j
        integer(kind=AbqIK),parameter:: NSTATVUMAT=1, NPROPSUMAT=2, n_iter_max=10
        real(kind=AbqRK),parameter:: allowedResRel=0.000001
        real(kind=AbqRK),dimension(NSTATV-NSTATVUMAT):: STATEVHARD
        real(kind=AbqRK),dimension(6):: STRESS_out
        pointer(pt_STRESS_out,STRESS_out)
        pt_STRESS_out=SMAFloatArrayAccess(1)
    !---------------------------------------------------------------------------------
    ! Actual Subroutine
    !---------------------------------------------------------------------------------
        PNEWDT=1.0_AbqRK !initialize PNEWDT, lower it, when necessary (too many iterations when calculation plastic multiplier)
        E = PROPS(1) ! Young's modulus is the 1st property in input file.
        Nu = PROPS(2) ! Poisson's ratio is the 2nd property.
        ! elastic coefficients
        shearmod2=E/(1.0_AbqRK+Nu) ! 2*shear modulus = first elastic coefficient
        shearmod=shearmod2/2.0_AbqRK ! shear modulus
        compmod=E/(1.0_AbqRK-2.0_AbqRK*Nu)/3.0_AbqRK ! compression modulus
        PEEQ=STATEV(1)
        PEEQRATE=0.0_AbqRK
        ! Elastic predictor
        DilatStrain = sum( DSTRAN(1:NDI) )
        MeanStrain = DilatStrain/3
        StressHydro=sum(STRESS(1:NDI))/NDI ! just store old hydrostatic stress temporarily in order to avoid to compute it twice
        StressDevTrial(1:NDI) = STRESS(1:NDI)-StressHydro+ shearmod2*(DSTRAN(1:NDI)- MeanStrain)
        StressDevTrial(NDI+1:NTENS) = STRESS(NDI+1:NTENS)+ shearmod*DSTRAN(NDI+1:NTENS)
        StressHydro=StressHydro+compmod*DilatStrain !hydrostatic stress
        StressMisesTrial=sqrt(1.5_AbqRK*dot_product(StressDevTrial(1:NDI),StressDevTrial(1:NDI))+3*dot_product(StressDevTrial(NDI+1:NTENS),StressDevTrial(NDI+1:NTENS)))
        if (NPROPS>NPROPSUMAT) then
            if (size(STATEVHARD)>0) STATEVHARD(:)=STATEV(NSTATVUMAT+1:NSTATV)
           call UHARD(SYIELD,HARD,PEEQ,PEEQRATE,TIME,DTIME,TEMP, DTEMP,NOEL,NPT,LAYER,KSPT,JSTEP(1),KINC,CMNAME,size(STATEVHARD),STATEVHARD,size(PREDEF),PREDEF,DPRED,NPROPS-NPROPSUMAT,PROPS(NPROPSUMAT+1:NPROPS))
        end if
    !Elastic step
        if (StressMisesTrial<=SYIELD .OR. NPROPS==NPROPSUMAT) then
          STRESS(1:NDI)=StressHydro+StressDevTrial(1:NDI)
          STRESS(NDI+1:NTENS)=StressDevTrial(NDI+1:NTENS)
          coef2=compmod-shearmod2/3
          coef1=coef2+shearmod2
          ! Tangent stiffness
          DDSDDE = 0.0_AbqRK
          forall(i=1:NDI)         
            DDSDDE(i,i) = coef1
          end forall
          forall(i=1:NDI, j=1:NDI, i.ne.j) DDSDDE(i,j) = coef2
          forall(i=NDI+1:NTENS) DDSDDE(i,i) = shearmod
          PNEWDT=1.5_AbqRK
    !Plastic Corrector
    else
      !Solve for plastic multiplier
          shearmod3=shearmod*3.0_AbqRK
          allowedRes=allowedResRel*StressMisesTrial !allowed residual
          dPEEQ=(StressMisesTrial-SYIELD)/(shearmod3+HARD(1))
          n_iter=1
          do while (.TRUE.)
           PEEQ=STATEV(1)+dPEEQ
           if (size(STATEVHARD)>0) STATEVHARD(:)=STATEV(NSTATVUMAT+1:NSTATV)
           call UHARD(SYIELD,HARD,PEEQ,PEEQRATE,TIME,DTIME,TEMP, DTEMP,NOEL,NPT,LAYER,KSPT,JSTEP(1),KINC,CMNAME,SIZE(STATEVHARD),STATEVHARD,SIZE(PREDEF),PREDEF,DPRED,NPROPS-NPROPSUMAT,PROPS(NPROPSUMAT+1))
           R=StressMisesTrial-shearmod3*dPEEQ-SYIELD !Residual
           if (abs(R)<allowedRes) then
              if (n_iter<=n_iter_max*0.8_AbqRK) PNEWDT=1.5_AbqRK
              exit
           end if
           dPEEQ=dPEEQ+R/(shearmod3+HARD(1))
           n_iter=n_iter+1
           !abort if too many iterations required and request smaller time increment
           if (n_iter>n_iter_max) then
              PNEWDT=0.5_AbqRK
              exit
           end if
          end do
          !Stress Update
          SdevScaleFactor=SYIELD/StressMisesTrial
          STRESS(1:NDI)=StressHydro+SdevScaleFactor*StressDevTrial(1:NDI)
          STRESS(NDI+1:NTENS)=SdevScaleFactor*StressDevTrial(NDI+1:NTENS)
          coef2=compmod-SdevScaleFactor*shearmod2/3
          coef1=coef2+SdevScaleFactor*shearmod2
          !tangent stiffness matrix
          DDSDDE = 0.0_AbqRK
          forall(i=1:NDI) DDSDDE(i,i) = coef1
          forall(i=1:NDI, j=1:NDI, i.ne.j) DDSDDE(i,j) = coef2
          forall(i=NDI+1:NTENS) DDSDDE(i,i) = SdevScaleFactor*shearmod
          DDSDDE(1:NTENS,1:NTENS)=DDSDDE(1:NTENS,1:NTENS)-shearmod3**2/(shearmod3+HARD(1))*(SYIELD-HARD(1)*dPEEQ)/(StressMisesTrial**3)*matmul(reshape(StressDevTrial,(/NTENS,1/)),reshape(StressDevTrial,(/1,NTENS/)))
    end if

    !tangent correction (part associated with volumetric changes)
    if (JSTEP(3)==1) forall(i=1:NDI) DDSDDE(i,1:NDI)=DDSDDE(i,1:NDI)+STRESS(i)
        
    !Update State variables
    STATEV(1)=PEEQ
    if (size(STATEVHARD)>0) STATEV(NSTATVUMAT+1:NSTATV)=STATEVHARD(:)
    !write STRESS to STRESS_out
    STRESS_out(1:NTENS)=STRESS

    end subroutine
