!DEC$ FREEFORM

    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)
    
    !Hyperelastic SUBROUTINE UMAT, cannot be used for plane stress
    !assumed constant stiffness in Material configuration

    USE ABQinterface
                          
    IMPLICIT NONE
    
    !---------------------variables from Abaqus interface----------------------
    CHARACTER(LEN=80):: CMNAME
    INTEGER(KIND=AbqIK)::NOEL,NPT,LAYER,KSPT,KINC,NDI,NSHR,NTENS,NSTATV,NPROPS
    REAL(KIND=AbqRK),DIMENSION(NTENS):: STRESS,STRAN,DSTRAN
    REAL(KIND=AbqRK),DIMENSION(NSTATV):: STATEV
    REAL(KIND=AbqRK),DIMENSION(NTENS,NTENS):: DDSDDE
    REAL(KIND=AbqRK),DIMENSION(NTENS)::DDSDDT
    REAL(KIND=AbqRK),DIMENSION(NTENS):: DRPLDE
    REAL(KIND=AbqRK),DIMENSION(2):: TIME
    REAL(KIND=AbqRK),DIMENSION(1):: PREDEF,DPRED
    REAL(KIND=AbqRK),DIMENSION(NPROPS):: PROPS
    REAL(KIND=AbqRK),DIMENSION(3):: COORDS
    REAL(KIND=AbqRK),DIMENSION(3,3):: DROT,DFGRD0,DFGRD1
    INTEGER(KIND=AbqIK),DIMENSION(4):: JSTEP
    REAL(KIND=AbqRK):: DTIME,TEMP,DTEMP,SSE,SPD,SCD,CELENT,DRPLDT,PNEWDT,RPL
	REAL(KIND=AbqrK),DIMENSION(3,3,3,3):: Stiffness_PK2,Stiffness
	REAL(KIND=AbqrK),DIMENSION(9,9):: DDSDDE_PK1
	REAL(KIND=AbqrK),DIMENSION(3,3):: Cauchy,Kirchhoff
	REAL(KIND=AbqrK),DIMENSION(3,3):: PK2
	REAL(KIND=AbqrK),DIMENSION(3,3):: Green_Lagrange
	REAL(KIND=AbqrK),DIMENSION(3,3),PARAMETER:: IDENT=reshape([1.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,1.0_AbqRK,0.0_AbqRK,0.0_AbqRK,0.0_AbqRK,1.0_AbqRK],[3,3])
	INTEGER(KIND=AbqIK),DIMENSION(6,2),PARAMETER::ZO_3D=reshape([1,2,3,1,1,2,1,2,3,2,3,3],[6,2])
	INTEGER(KIND=AbqIK),DIMENSION(4,2),PARAMETER::ZO_2D=reshape([1,2,3,1,1,2,3,2],[4,2])
	REAL(KIND=AbqrK)::lambda,mu,DET
	INTEGER(KIND=AbqIK)::i,j,k,l,m,n,o,p
	
	!PROPS(1)=> Joungs Modulus
	!PROPS(2)=> poissions ratio
	
	!Lame constants
	lambda=PROPS(2)/(1.0_AbqRK-2.0_AbqRK*PROPS(2))/(1.0_AbqRK+PROPS(2))*PROPS(1)
	mu=1.0_AbqRK/2.0_AbqRK/(1.0_AbqRK+PROPS(2))*PROPS(1)

	!get Stiffness_PK2
    Stiffness_PK2=0.0_AbqRK
	forall(i=1:3,j=1:3,k=1:3,l=1:3) Stiffness_PK2(i,j,k,l)=lambda*IDENT(i,j)*IDENT(k,l)+mu*(IDENT(i,k)*IDENT(j,l)+IDENT(i,l)*IDENT(j,k))
	
	!get green lagrange strain
	Green_Lagrange=1.0_AbqRK/2.0_AbqRK*(matmul(transpose(DFGRD1),DFGRD1)-IDENT)
    
	!get second Piola Kirchhoff stress
    PK2=0.0_AbqRK
	forall(i=1:3,j=1:3,k=1:3,l=1:3) PK2(i,j)=PK2(i,j)+Stiffness_PK2(i,j,k,l)*Green_Lagrange(k,l)
	
    !get Kirchhoff stress
    Kirchhoff=matmul(matmul(DFGRD1,PK2),transpose(DFGRD1))
    
    !get determinant of deformation gradient
    DET = DFGRD1(1,1)*(DFGRD1(2,2)*DFGRD1(3,3) - DFGRD1(2,3)*DFGRD1(3,2)) - DFGRD1(1,2)*(DFGRD1(2,1)*DFGRD1(3,3)-DFGRD1(2,3)*DFGRD1(3,1)) + DFGRD1(1,3)*(DFGRD1(2,1)*DFGRD1(3,2)- DFGRD1(3,1)*DFGRD1(2,2))
    
    !get DDSDDE
    !push forward stiffness
    Stiffness=0.0_AbqRK
    forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3,o=1:3,p=1:3) Stiffness(i,j,k,l)=Stiffness(i,j,k,l)+Stiffness_PK2(m,n,o,p)*DFGRD1(i,m)*DFGRD1(j,n)*DFGRD1(k,o)*DFGRD1(l,p)
    !add contribution
    forall(i=1:3,j=1:3,k=1:3,l=1:3) Stiffness(i,j,k,l)=Stiffness(i,j,k,l)+0.5_AbqRK*(Kirchhoff(i,k)*IDENT(j,l)+Kirchhoff(j,k)*IDENT(i,l)+Kirchhoff(i,l)*IDENT(j,k)+Kirchhoff(j,l)*IDENT(i,k))
    !get voigt presentation of stiffness and Kirchhoff stress
    DDSDDE=0.0_AbqRK
    if (NTENS==6) then
        forall(o=1:NTENS,p=1:NTENS) DDSDDE(o,p)=DDSDDE(o,p)+Stiffness(ZO_3D(o,1),ZO_3D(o,2),ZO_3D(p,1),ZO_3D(p,2))
        forall(i=1:NTENS) STRESS(i)=Kirchhoff(ZO_3D(i,1),ZO_3D(i,2))
    else
        forall(o=1:NTENS,p=1:NTENS) DDSDDE(o,p)=DDSDDE(o,p)+Stiffness(ZO_2D(o,1),ZO_2D(o,2),ZO_2D(p,1),ZO_2D(p,2))
        forall(i=1:NTENS) STRESS(i)=Kirchhoff(ZO_2D(i,1),ZO_2D(i,2))
    end if
    DDSDDE=DDSDDE/DET
    STRESS=STRESS/DET
    
    END SUBROUTINE UMAT
