!DEC$ FREEFORM

! Rostyslav Skrypnyk, 28.02.2014
! Stephan Roth, 18.10.2017

!!!!!!!!!!!!!!!!!!!!!
! Utilities for UMAT
!!!!!!!!!!!!!!!!!!!!!
MODULE UMATUTILITIES
	USE ABQINTERFACE
	IMPLICIT NONE
	PRIVATE
	PUBLIC:: CheckMaterialParameters, GetNSTATV
	CONTAINS
		! Subroutine 'CheckMaterialParameters' checks correctness of material parameters in input file.
      SUBROUTINE CheckMaterialParameters(PROPS)
			  REAL(KIND=AbqRK),INTENT(IN):: PROPS(:)
			  ! Check if there is a correct number of properties specified in the input file.
        if ( size(PROPS)<2 ) then
            CALL STDB_ABQERR(-3_AbqIK, "CheckMaterialParameters reports: " // &
                 "Invalid number of properties in the input file, at least 2 required.", &
                 0, 0.0_AbqRK, " ")
        end if
        ! Check if Young's modulus is > zero
        if (PROPS(1)<=0) then
            CALL STDB_ABQERR(-3_AbqIK, "CheckMaterialParameters reports: " // &
                 "Young's modulus must be greater than zero, and it should be the first one in the list of properties.", &
                 0, 0.0_AbqRK, " ")
        end if
        if (PROPS(2)<=-1.0 .or. PROPS(2)>=0.5) then
            CALL STDB_ABQERR(-3_AbqIK, "CheckMaterialParameters reports: " // &
            "Poisson's ratio must belong interval ]-1, 0.5[, and it should be the second one in the list of properties.", &
            0, 0.0_AbqRK, " ")
        end if
      END SUBROUTINE CheckMaterialParameters
		! Function 'GetNSTATEV' provides number of state variables needed for UMAT
      FUNCTION GetNSTATV(NTENS,NDI)
			  INTEGER(KIND=AbqIK):: GetNSTATV
			  INTEGER(KIND=AbqIK), INTENT(IN):: NTENS,NDI

			  GetNSTATV=0
      END FUNCTION GetNSTATV
END MODULE UMATUTILITIES

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! UMAT for linear elastic 3D & 2D plane strain cases
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
MODULE UXMAT1
	USE ABQINTERFACE
	USE UMATUTILITIES ! To call routines that check input data
   IMPLICIT NONE
   PRIVATE
   PUBLIC::UXMAT, CheckMaterialParameters, GetNSTATV
   CONTAINS
		SUBROUTINE UXMAT(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,KSTEP,KINC)
			! Declare interface
			REAL(KIND=AbqRK),INTENT(IN):: TIME(:),PREDEF(:),STRAN(:),DSTRAN(:),DPRED(:),PROPS(:),COORDS(:),DROT(:,:),DFGRD0(:,:),DFGRD1(:,:), &
													DTIME,TEMP,DTEMP,CELENT
			INTEGER(KIND=AbqIK), INTENT(IN):: NTENS,NSTATV,NOEL,NPT,LAYER,KSPT,KINC,NPROPS,NDI,NSHR,KSTEP
			CHARACTER*80 CMNAME
			REAL(KIND=AbqRK):: STRESS(:), STATEV(:), DDSDDE(:,:), DDSDDT(:), DRPLDE(:), &
									 SSE,SPD,SCD,RPL,DRPLDT,PNEWDT
			! Declare local arguments
			REAL(KIND=AbqRK):: E, & ! Young's modulus
									 Nu, shearmod, & ! Poisson's ratio and shear modulus
									 coef0, coef1, coef2, coef3, coeftemp, & ! coefficients for stiffness tensor
									 StrainUpd(NTENS), StrainTrace ! store trace of the strain at the end of the time increment
			INTEGER(KIND=AbqIK):: i, j

      E = PROPS(1) ! Make sure that Young's modulus is the 1st property in input file.
			Nu = PROPS(2) ! Make sure that Poisson's ratio is the 2nd property in input file.

      ! Strain values at end of increment
			StrainUpd(:) =  STRAN(1:NTENS)+DSTRAN(1:NTENS)
			StrainTrace = SUM( StrainUpd(1:NDI) )

			! Assemble material stiffness matrix and compute stress
      coef0=E/(1.0_AbqRK+Nu)
      shearmod=coef0/2.0_AbqRK !=E/[2*(1+Nu)]
			coeftemp = shearmod/(0.5_AbqRK-Nu)!=E/[(1+Nu)*(1-2*Nu)]
			coef2 = Nu*coeftemp !=Nu*E/[(1+Nu)*(1-2*Nu)]
			coef1 = coeftemp - coef2 != (1-Nu)*E/[(1+Nu)*(1-2*Nu)]

      DDSDDE=0.0_AbqRK

			FORALL(i=1:NDI)
			   DDSDDE(i,i) = coef1
			   STRESS(i) = coef0*StrainUpd(i) + coef2*StrainTrace
			END FORALL

			FORALL(i=1:NDI, j=1:NDI, i.NE.j) DDSDDE(i,j) = coef2

               IF (NTENS .GT. NDI) THEN
    			  FORALL(i=NDI+1:NTENS)
			    DDSDDE(i,i) = shearmod
			    STRESS(i) = shearmod*StrainUpd(i)
			  END FORALL
               END IF

		END SUBROUTINE UXMAT
END MODULE UXMAT1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! UMAT for linear elastic 2D plane stress case
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
MODULE UMATPS
	USE ABQINTERFACE
	USE UMATUTILITIES ! To call routines that check input data
	IMPLICIT NONE
	PRIVATE
	PUBLIC:: UXMAT, CheckMaterialParameters, GetNSTATV
	CONTAINS
		SUBROUTINE UXMAT(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,KSTEP,KINC)
			! Declare interface
			REAL(KIND=AbqRK),INTENT(IN):: TIME(:),PREDEF(:),STRAN(:),DSTRAN(:),DPRED(:),PROPS(:),COORDS(:),DROT(:,:),DFGRD0(:,:),DFGRD1(:,:), &
													DTIME,TEMP,DTEMP,CELENT
			INTEGER(KIND=AbqIK), INTENT(IN):: NTENS,NSTATV,NOEL,NPT,LAYER,KSPT,KINC,NPROPS,NDI,NSHR,KSTEP
			CHARACTER*80 CMNAME
			REAL(KIND=AbqRK):: STRESS(:), STATEV(:), DDSDDE(:,:), DDSDDT(:), DRPLDE(:), &
									 SSE,SPD,SCD,RPL,DRPLDT,PNEWDT
			! Declare local arguments
			REAL(KIND=AbqRK):: E, & ! Young's modulus
									 Nu, & ! Poisson's ratio
									 coef1, coef2, coef3, & ! coefficients for stiffness tensor
									 StrainTrace ! store trace of the strain at the end of the time increment
			INTEGER(KIND=AbqIK):: i, j

         E = PROPS(1) ! Make sure that Young's modulus is the 1st property in input file.
			Nu = PROPS(2) ! Make sure that Poisson's ratio is the 2nd property in input file.

			StrainTrace = SUM( STRAN(1:NDI) ) + SUM( DSTRAN(1:NDI) )

			! Assemble stiffness tensor and compute stress
			DDSDDE = 0.0_AbqRK

			coef1 = E/( 1.0_AbqRK-Nu*Nu )
			coef2 = Nu*coef1
			coef3 = 0.5_AbqRK*(1.0_AbqRK-Nu)*coef1

			FORALL(i=1:NDI)
			   DDSDDE(i,i) = coef1
			   STRESS(i) = (coef1-coef2)*( STRAN(i)+DSTRAN(i) ) + coef2*StrainTrace
			END FORALL

			FORALL(i=1:NDI, j=1:NDI, i.NE.j) DDSDDE(i,j) = coef2

               IF (NSHR .GE. 1) THEN
			  FORALL(i=NDI+1:NDI+NSHR)
			    DDSDDE(i,i) = coef3
			    STRESS(i) = coef3*( STRAN(i)+DSTRAN(i) )
			  END FORALL
               END IF

		END SUBROUTINE UXMAT
END MODULE UMATPS






