!DEC$ FREEFORM

!=====================================================================
! - general empty Interface for UMAT to be used inside UelXCU
! - contains general plane stress algorithm
!=====================================================================

MODULE UXMATModule
	USE ABQINTERFACE
	IMPLICIT NONE
	PRIVATE
	PUBLIC::UXMAT, CheckMaterialParameters, GetNSTATV
	CONTAINS
		                                                                                                                                                           
		INCLUDE 'UMAT.f'	!include the actual UMAT
		INCLUDE 'CheckMaterialParameters.f' !CheckMaterialParameters
		INCLUDE 'GetNSTATV.f' !GetNSTATV

		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,JSTEP,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,JSTEP(:)
			CHARACTER*80 CMNAME
			REAL(KIND=AbqRK)::STRESS(:),STATEV(:),DDSDDE(:,:),DDSDDT(:),DRPLDE(:),&
							SSE,SPD,SCD,RPL,DRPLDT,PNEWDT,STRESS_PS(4),DSTRAN_PS(4),&
							DDSDDE_PS(4,4),STRAN_PS(4),PNEWDT_PS,DDSDDT_PS(4),&
							DRPLDE_PS(4),STATEV_PS(NSTATV-1),DFGRD0_PS(3,3),DFGRD1_PS(3,3)
			INTEGER(KIND=AbqIK),PARAMETER:: NDI_PS=3,NTENS_PS=4,NSHR_PS=1
			INTEGER(KIND=AbqIK):: I
			REAL(KIND=AbqRK),PARAMETER:: TOLER=0.00001
			INTEGER:: NEWTON=15
			
			IF (NTENS==3) THEN !plane stress case
				!STRAN of last increment
				STRAN_PS(1:2)=STRAN(1:2)
				STRAN_PS(3)=STATEV(1)
				STRAN_PS(4)=STRAN(3)
				!deformation gradient of last increment
				DFGRD0_PS=DFGRD0
				IF (JSTEP(3)==1) DFGRD0_PS(3,3)=EXP(STRAN_PS(3))
				!deformation gradient of current increment
				DFGRD1_PS=DFGRD1
				IF (JSTEP(3)==1) DFGRD1_PS(3,3)=EXP(STRAN_PS(3)) !initial guess
				!DSTRAN
				DSTRAN_PS(1:2)=DSTRAN(1:2)
				DSTRAN_PS(4)=DSTRAN(3)
				DSTRAN_PS(3)=0.0_AbqRK !initial guess

				DO I=1,NEWTON !iterate for DSTRAN_PS(3) to unti STRESS_PS(3)==0.0
				   !STATEV of last increment
				   STATEV_PS=STATEV(2:NSTATV)
				   !STRESS of last increment
				   STRESS_PS(1:2)=STRESS(1:2)
				   STRESS_PS(4)=STRESS(3)
				   STRESS_PS(3)=0.0_AbqRK
				   
				   PNEWDT_PS=1.0_AbqRK
				   ! Call UMAT in FORTRAN77 style  
				   CALL UMAT(STRESS_PS,STATEV_PS,DDSDDE_PS,SSE,SPD,SCD, RPL,DDSDDT_PS,&
							DRPLDE_PS,DRPLDT, STRAN_PS,DSTRAN_PS,TIME,DTIME,TEMP,DTEMP,&
							PREDEF,DPRED,CMNAME,NDI_PS,NSHR_PS,NTENS_PS,NSTATV-1,&
							PROPS,NPROPS,COORDS,DROT,PNEWDT_PS, CELENT,DFGRD0_PS,&
							DFGRD1_PS,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
							 
					IF (ABS(STRESS_PS(3)/MAXVAL(ABS(STRESS_PS)))<=TOLER .OR.&
					    MAXVAL(ABS(STRESS_PS))==0.0_AbqRK) THEN
					
						IF (I<=NEWTON/2 .AND. PNEWDT_PS>=1.0_AbqRK) THEN
							PNEWDT=1.5_AbqRK
						ELSE
							PNEWDT=PNEWDT_PS
						END IF
						
						EXIT
						
					ELSE IF (I==NEWTON) THEN
					
						PNEWDT=0.5_AbqRK
					END IF
                    
                    !update strain in x3 direction
                    DSTRAN_PS(3)=DSTRAN_PS(3)-STRESS_PS(3)/DDSDDE_PS(3,3)
                    !update deformation gradient
                    IF (JSTEP(3)==1) DFGRD1_PS(3,3)=EXP(STRAN_PS(3)+DSTRAN_PS(3))
                    
                  END DO

                  !update STATEV
                  STATEV(1)=STATEV(1)+DSTRAN_PS(3) !updated stran in x3 direction
                  STATEV(2:NSTATV)=STATEV_PS
                  
                  !output reduced STRESS and reduced DDSDDE
                  
                  STRESS(1:2)=STRESS_PS(1:2)
                  STRESS(3)=STRESS_PS(4)
                  
                  DDSDDE(1,1)=DDSDDE_PS(1,1)-1/DDSDDE_PS(3,3)*DDSDDE_PS(1,3)**2
                  DDSDDE(1,2)=DDSDDE_PS(1,2)-1/DDSDDE_PS(3,3)*DDSDDE_PS(1,3)*DDSDDE_PS(2,3)
                  DDSDDE(2,1)=DDSDDE(1,2)
                  DDSDDE(1,3)=DDSDDE_PS(1,4)-1/DDSDDE_PS(3,3)*DDSDDE_PS(1,3)*DDSDDE_PS(4,3)
                  DDSDDE(3,1)=DDSDDE(1,3)
                  DDSDDE(2,2)=DDSDDE_PS(2,2)-1/DDSDDE_PS(3,3)*DDSDDE_PS(2,3)**2
                  DDSDDE(2,3)=DDSDDE_PS(2,4)-1/DDSDDE_PS(3,3)*DDSDDE_PS(2,3)*DDSDDE_PS(4,3)
                  DDSDDE(3,2)=DDSDDE(2,3)
                  DDSDDE(3,3)=DDSDDE_PS(4,4)-1/DDSDDE_PS(3,3)*DDSDDE_PS(4,3)**2
                     
            ELSE !3D, plane strain, axissymmetric
				! Call UMAT in FORTRAN77 style                                                                            
				CALL 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)
            END IF
            
      END SUBROUTINE UXMAT
      
END MODULE UXMATModule
