********************************************************************************************
**  UMAT, FOR ABAQUS/STANDARD INCORPORATING ISOTROPIC LINEAR ELASTICITY, RATE INDEPENDENT **
**  VON MISES PLASTICTY AND LINEAR ISOTROPIC HARDENING. LARGE DEFORMATION FORMULATION     **
**  FOR 3D SOLID ELEMENTS. IMPLICIT INTEGRATION WITH CONSISTENT JACOBIAN                  **
********************************************************************************************
********************************************************************************************
**
*USER SUBROUTINE
      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,
     1                DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,
     2                PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,
     3                NPROPS,COORDS,DROT,PNEWDT,CELENT,DFGRD0,DFGRD1,
     4                NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
      IMPLICIT NONE
C
      CHARACTER*80 CMNAME    
C     
      DOUBLE PRECISION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),
     +        SSE,SPD,SCD,RPL,DDSDDT(NTENS),DRPLDE(NTENS),DRPLDT,
     +        STRAN(NTENS),DSTRAN(NTENS),TIME(2),DTIME,TEMP,DTEMP,
     +        PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3), 
     +        PNEWDT,CELENT,DFGRD0(3,3),DFGRD1(3,3)
C
      INTEGER NDI,NSHR,NTENS,NSTATV,NPROPS,NOEL,NPT,LAYER,KSPT,KSTEP,
     +        KINC
C
      DOUBLE PRECISION ZERO,ONE,TWO,THREE,SIX,NINE,TOLER,XIDEN(3,3),   
     +          XNV(6),DPSTRAN(6), STRESSOLD(6),DESTRAN(6),DSTRESS(6),
     +          XNDIR(3,3),STR(3,3),DPSTRN(3,3),DE,XNUE,SIGY0,h,p,r,
     +          EG,PJ,ZY,r0,dp,RES,XR,Q,XK,E,DSTR(3,3)
C
      INTEGER   I,J,K,L,M,N,ID
C
      PARAMETER (M=3,N=3,ID=3,ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,
     +          SIX=6.D0, NINE=9.D0, TOLER=1.D-5)     
C
C--------------------------------------------------------------------
C
C     SPECIFY MATERIAL PROPERTIES
C
      E = PROPS(1)
      XNUE = PROPS(2)
      SIGY0 = PROPS(3)
      h = PROPS(4)    
C
C      E = 21000.0
C      XNUE = 0.3
C      SIGY0 = 240.
C      h = 1206.0
C
C    ZERO MATRICES WITH UNSPECIFIED COMPONENTS
C
      DO I=1,M
      DO J=1,N
      STR(I,J)=0.0
      END DO
      END DO
      
C
C    RECOVER EFFECTIVE PLASTIC STRAIN, p, AND ISOTROPIC 
C    HARDENING VARIABLE, r,FROM PREVIOUS TIME STEP
C
      p = STATEV(1)
      r = STATEV(2)
C 
C    SET UP ELASTICITY MATRIX
C   
      DO I=1,6
      DO J=1,6
      DDSDDE(I,J)=ZERO
      END DO
      END DO
C
      DDSDDE(1,1)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*(ONE-XNUE)
      DDSDDE(1,2)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*XNUE
      DDSDDE(1,3)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*XNUE
      
      DDSDDE(2,1)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*XNUE
      DDSDDE(2,2)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*(ONE-XNUE)
      DDSDDE(2,3)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*XNUE
      
      DDSDDE(3,1)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*XNUE
      DDSDDE(3,2)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*XNUE
      DDSDDE(3,3)=(E/((ONE+XNUE)*(ONE-TWO*XNUE)))*(ONE-XNUE)
      
      DDSDDE(4,4)=(E/(ONE+XNUE))/TWO
      DDSDDE(5,5)=(E/(ONE+XNUE))/TWO
      DDSDDE(6,6)=(E/(ONE+XNUE))/TWO
C
      EG=(E/(ONE+XNUE))/TWO      
      XK = (E/(ONE-TWO*XNUE))/THREE
C
C     DEFINE IDENTITY MATRIX
C
      DO I=1,M
      DO J=1,N
        IF(I .EQ. J) THEN
          XIDEN(I,J)=1.0D0
        ELSE
          XIDEN(I,J)=0.0D0
        END IF
      END DO
      END DO
C
C     SAVE STRESS AT BEGINNING OF TIME STEP IN STRESSOLD
C
      DO K=1,6
      STRESSOLD(K) = STRESS(K)
      END DO
C
C     OBTAIN TRIAL (ELASTIC) STRESS
C
      CALL KMLT1(DDSDDE,DSTRAN,DSTRESS,NTENS)
      DO K=1,NTENS
      STRESS(K) = STRESS(K) + DSTRESS(K)
      END DO
C
C     WRITE TRIAL STRESSES IN MATRIX FORM
C      
      DO K = 1,3
      STR(K,K) = STRESS(K)
      END DO
      STR(1,2) = STRESS(4)
      STR(2,1) = STRESS(4)
      STR(1,3) = STRESS(5)
      STR(3,1) = STRESS(5)
      STR(2,3) = STRESS(6)
      STR(3,2) = STRESS(6)                
C
C     CALCULATE DEVIATORIC TRIAL STRESS: SIGTRPR
C 
      CALL KDEVIA(STR,XIDEN,DSTR) 
C
C     CALCULATE EFFECTIVE TRIAL STRESS: SIGETR
C
      CALL KEFFP(DSTR,PJ)   
C
C    DETERMINE TRIAL FLOW DIRECTION FOR USE WITH THE JACOBIAN: SIGETR/|SIGETR|
C
      IF (PJ .NE. ZERO) THEN
      DO I=1,3
      DO J=1,3
      XNDIR(I,J) = DSTR(I,J)/PJ
      END DO
      END DO
      ELSE
      XNDIR(I,J) = ZERO
      END IF
C
C    ...AND WRITE IN VOIGT NOTATION
C
      DO K=1,3
      XNV(K) = XNDIR(K,K)
      END DO
      XNV(4) = XNDIR(1,2)
      XNV(5) = XNDIR(1,3)      
      XNV(6) = XNDIR(2,3)    
C
C    DETERMINE IF THE YIELD CONDITION IS SATISFIED
C     
      ZY = PJ - r - SIGY0
C
      IF (ZY.GT.0.) THEN
C
C     USE NEWTON ITERATION TO DETERMINE EFFECTIVE PLASTIC STRAIN INCREMENT
C 
      r0 = r
      dp = 0.D0
      RES = PJ-3.D0*EG*dp-r-SIGY0
      DO WHILE (DABS(RES) .GT. TOLER) 
      RES = PJ-3.D0*EG*dp-r-SIGY0   
      dp = dp + RES/(3.D0*EG+h)
      r = r0 + h*dp 
      RES = PJ-3.D0*EG*dp-r-SIGY0   
      END DO 
C
C     DETERMINE THE INCREMENTS IN PLASTIC STRAIN
C
      DO I = 1,3
       DO J = 1,3
       IF (PJ .NE. ZERO) THEN
       DPSTRN(I,J) = (THREE/TWO)*dp*DSTR(I,J)/PJ
       ELSE
       DPSTRN(I,J) = ZERO
       END IF
       END DO
      END DO   
C
C     WRITE THE STRAIN INCREMENTS IN VOIGT NOTATION (WITH ENGG SHEARS)
C
      DO K=1,3
      DPSTRAN(K) = DPSTRN(K,K)
      END DO
      DPSTRAN(4) = 2*DPSTRN(1,2)
      DPSTRAN(5) = 2*DPSTRN(1,3)
      DPSTRAN(6) = 2*DPSTRN(2,3)
C      
C     CALCULATE THE ELASTIC STRAIN INCREMENTS      
C      
      DO K=1,NTENS
      DESTRAN(K)=DSTRAN(K)-DPSTRAN(K)
      END DO  
C
C     DETERMINE STRESS INCREMENT
C
       CALL KMLT1(DDSDDE,DESTRAN,DSTRESS,NTENS)     
C
C      UPDATE THE STRESS, EFFECTIVE PLASTIC STRAIN
C      (NOTE: ISOTROPIC HARDENING VARIABLE ALREADY UPDATED)
C
       DO K = 1,NTENS
       STRESS(K) = STRESSOLD(K) + DSTRESS(K)
       END DO
C
       p = p + dp
C       
C     STORE UPDATED STATE VARIABLES
C
       STATEV(1) = p
       STATEV(2) = r    
       
      IF (NPT .EQ. 1) THEN 
       
      PRINT *, 'STATEV(1)=',STATEV(1)
      PRINT *, 'STATEV(2)=',STATEV(2)
      
      END IF
       
C       print *, 'statev in the umat=', statev     
C
C    DETERMINE MATERIAL JACOBIAN
C
      IF (PJ .NE. ZERO) THEN
      XR = (PJ-THREE*EG*dp)/PJ
      ELSE
      XR = ZERO
      END IF       
      
      Q = THREE/TWO*(h/(h+THREE*EG)-XR)     
      
      DO I=1,3
      DO J=1,3
       DDSDDE(I,J) = TWO*EG*Q*XNV(I)*XNV(J) 
     +  + TWO*EG*XR*XIDEN(I,J) + (XK-TWO/THREE*EG*XR)
      END DO
      END DO           
      
      DO K=1,3
      DO L=4,6
      DDSDDE(K,L) = TWO*EG*Q*XNV(K)*XNV(L)
      DDSDDE(L,K) = DDSDDE(K,L)
      END DO
      END DO
C      
      DDSDDE(4,5) = TWO*EG*Q*XNV(4)*XNV(5)
      DDSDDE(5,4) = DDSDDE(4,5)
      DDSDDE(4,6) = TWO*EG*Q*XNV(4)*XNV(6)
      DDSDDE(6,4) = DDSDDE(4,6)
      DDSDDE(5,6) = TWO*EG*Q*XNV(5)*XNV(6)
      DDSDDE(6,5) = DDSDDE(5,6)      
C            
      DDSDDE(4,4) = TWO*EG*Q*XNV(4)*XNV(4) + EG*XR
      DDSDDE(5,5) = TWO*EG*Q*XNV(5)*XNV(5) + EG*XR
      DDSDDE(6,6) = TWO*EG*Q*XNV(6)*XNV(6) + EG*XR     
C
      END IF
C     
      RETURN
      END
**
***********************************************
**           UTILITY    SUBROUTINES           *
***********************************************
**
**
***************************************************
**         MULTIPLY 6X6 MATRIX WITH 6X1 VECTOR    *
***************************************************
*USER SUBROUTINE
      SUBROUTINE KMLT1(DM1,DM2,DM,NTENS)
C      
      IMPLICIT NONE
C
      INTEGER M, I, K, NTENS
      DOUBLE PRECISION X, Y, DM1(6,6),DM2(6),DM(6)
      PARAMETER (M=6)
C
      DO I=1,6
      X=0.D0
      DO K=1,6 
      Y=DM1(I,K)*DM2(K)
      X=X+Y
      END DO
      DM(I)=X
      END DO
      RETURN
      END
**
**
***************************************
**          EFFECTIVE STRESS          *
**   (CONTRACTED MATRIX CALCULATION)  *
***************************************
*USER SUBROUTINE
      SUBROUTINE KEFFP(EFF1,VAL1)
C      
      IMPLICIT NONE
C
      INTEGER M, N, I, J
      PARAMETER (M=3,N=3)
      DOUBLE PRECISION EFF1(3,3), X, VAL1
C
      X=0.D0
      DO I=1,M
      DO J=1,N
       X=X+EFF1(I,J)*EFF1(I,J)
      END DO
      END DO
      IF(X .GE. 0.D0) THEN 
      VAL1=DSQRT((3.D0/2.D0)*X)
      END IF
      RETURN
      END
**
********************************************
**         DOT PRODUCT OF TWO VECTORS      *
********************************************
*USER SUBROUTINE
      SUBROUTINE DOTPROD(DM1,DM2,DM,NTENS)
C      
      IMPLICIT NONE
C
C      PARAMETER (M=6)
      INTEGER NTENS, K
C
      DOUBLE PRECISION DM1(6),DM2(6),Y,X,DM
C
      Y=0.D0
      DO K=1,NTENS 
      X=DM1(K)*DM2(K)
      Y=X+Y
      END DO
      DM=Y
      RETURN
      END
**      
***********************************************
**         DYADIC PRODUCT OF TWO VECTORS      *
***********************************************
      SUBROUTINE DYADICPROD(DM1,DM2,DM3,NTENS)
C      
      IMPLICIT NONE
C
C      PARAMETER (M=4)
      INTEGER NTENS, I, J
C
      DOUBLE PRECISION DM1(6),DM2(6),DM3(6,6)
C
      DO I=1,NTENS
       DO J=1,NTENS
        DM3(I,J) = DM1(I)*DM2(J)
       END DO
      END DO
C  
      RETURN
      END
**
***********************************************
**         DEVIATORIC STRESS CALCULATION      *
***********************************************
*USER SUBROUTINE
      SUBROUTINE KDEVIA(STRSS,XIDENTY,DEVITO)
C      
      IMPLICIT NONE
C
      INTEGER M, N, I, J
      PARAMETER (M=3,N=3)
      DOUBLE PRECISION X, Y, STRSS(3,3),XIDENTY(3,3),DEVITO(3,3)
C
      X=0.D0
      DO I=1,M
      DO J=1,N
      IF(I .EQ. J) THEN
      X=X+STRSS(I,J)
      END IF
      END DO
      END DO
C
      DO I=1,M
      DO J=1,N
      IF(I .EQ. J) THEN
        DEVITO(I,J)=STRSS(I,J)-((1.D0/3.D0)*X*XIDENTY(I,J))
      ELSE
        DEVITO(I,J)=STRSS(I,J)
      END IF
      END DO
      END DO
      RETURN
      END
**
