三维编织umat有点小问题，欢迎大家讨论

zhuhao3802

纤维束子程序，能进行计算，但是不能进行失效判断。

zhuhao3802

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1     RPL,DDSDDT,DRPLDE,DRPLDT,
     2     STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3     NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4     CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C     
      INCLUDE 'ABA_PARAM.INC'
C     
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1     DDSDDE(NTENS,NTENS),
     2     DDSDDT(NTENS),DRPLDE(NTENS),
     3     STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     4     PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)

      DIMENSION CFULL(6,6),CDFULL(6,6)
      DIMENSION OLD_STRESS(6),STRANT(6),
     1      DOLD_STRESS(6),D_STRESS(6)
      DIMENSION DCDDL(6,6),DCDDT(6,6),DCDDZ(6,6),
     1      DDLDDE(6),DDTDDE(6),DDZDDE(6),DCTDDE(6),DTTDDE(6),
     2      DCLDDE(6),DTLDDE(6),DCZDDE(6),DTZDDE(6),
     3      ATEMP1(6),ATEMP2(6),ATEMP3(6)

      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,HALF=.5D0)

       EL=PROPS(1)
       ET=PROPS(2)
       VLT=PROPS(3)
       VTZ=PROPS(4)
       GLT=PROPS(5)
       XT=PROPS(6)
       XC=PROPS(7)
       YT=PROPS(8)
       YC=PROPS(9)
       SLT=PROPS(10)
       SLZ=SLT
       STZ=PROPS(11)
       GF=PROPS(12)
C   
C********************************************************************
C  INITIATION VARIABLES
      DO K1=1,6
       DDLDDE(K1)=ZERO
       DDTDDE(K1)=ZERO
       DDZDDE(K1)=ZERO
      END DO
C
C     CALCULATE THE STRAIN AT THE END OF THE INCREMENT   
      DO K1=1,NTENS
       STRANT(K1)=STRAN(K1)+DSTRAN(K1)
      END DO
C      
C     FILL THE 6X6 FULL STIFFNESS MATRIX
      DO K1=1,6
       DO K2=1,6
        CFULL(K1,K2)=ZERO
       END DO
      END DO
      GTZ=ET/(TWO*(ONE+VTZ))
      GZL=GLT
      VZL=VLT/EL*ET
      DET=ONE-TWO*VLT*VZL-VTZ**TWO-TWO*VLT*VZL*VTZ
      CFULL(1,1)=EL*(ONE-VTZ**TWO)/DET
      CFULL(2,2)=ET*(ONE-VLT*VZL)/DET
      CFULL(3,3)=CFULL(2,2)
      CFULL(1,2)=ET*(VLT+VLT*VTZ)/DET
      CFULL(1,3)=CFULL(1,2)
      CFULL(2,3)=ET*(VTZ+VLT*VZL)/DET
      CFULL(4,4)=GLT
      CFULL(5,5)=GZL
      CFULL(6,6)=GTZ
      DO I = 2, 6
         DO J = 1, I-1
            CFULL(I,J) = CFULL(J,I)
         END DO
      END DO

      print *,  CFULL
C
C     ASSIGNING STATE VARIABLE
      DLVOLD=STATEV(1)
      DTVOLD=STATEV(2)
      DZVOLD=STATEV(3)
      DLOLD=STATEV(4)
      DTOLD=STATEV(5)
      DZOLD=STATEV(6)
      DCLOLD=STATEV(7)
      DTLOLD=STATEV(8)
      DCTOLD=STATEV(9)
      DTTOLD=STATEV(10)
      DCZOLD=STATEV(11)
      DTZOLD=STATEV(12)
      XTLI=STATEV(19)
      XTLF=STATEV(20)
      XCLI=STATEV(21)
      XCLF=STATEV(22)
      XTTI=STATEV(23)
      XTTF=STATEV(24)
      XCTI=STATEV(25)
      XCTF=STATEV(26)
      XTZI=STATEV(27)
      XTZF=STATEV(28)
      XCZI=STATEV(29)
      XCZF=STATEV(30)
      ESTATE=STATEV(31)
C      
      CALL CheckFailureIni(STRANT,STRESS,DL,DT,DZ,DCL,DTL,DCT,DTT,DCZ,
     1     DTZ,DTLOLD,DCLOLD,DTTOLD,DCTOLD,DTZOLD,DCZOLD,CELENT,CFULL,
     2     XTL,XCL,XTT,XCT,XTZ,XCZ,XTLI,XTLF,XCLI,XCLF,XTTI,XTTF,XCTI,
     3     XCTF,XTZI,XTZF,XCZI,XCZF)

C     SAVE THE OLD STRESS TO OLD_STRESS
      DO I = 1, NTENS
         OLD_STRESS(I) = STRESS(I)
      END DO
C     
C     CALL ROUTINE TO CALCULATE THE STRESS
C     CALCULATE THE STRESS IF THERE'S NO VISCOUS REGULARIZATION
      CALL GetStress(CFULL,CDFULL,DL,DT,DZ,D_STRESS,STRANT,NDI,NTENS)
C     CALCULATE THE STRESS IF THERE'S VISCOUS REGULARIZATION
      CALL GetStress(CFULL,CDFULL,DLV,DTV,DZV,STRESS,STRANT,NDI,NTENS)
C     GET THE OLD STRESS IF THERE'S NO VISCOUS REGULARIZATION
      DO I=1,NTENS
         DOLD_STRESS(I)=STATEV(I+12)
      END DO
C     SAVE THE CURRENT STRESS IF THERE'S NO VISCOUS REGULARIZATION
      DO I=1,NTENS
         STATEV(I+12)=D_STRESS(I)
      END DO
C        
C     UPDATED JACOBI MATRIX
C        
      CALL ELASTICDERIVATIVE (CFULL,DLV,DTV,DZV,DCDDL,DCDDT,DCDDZ)
      DO I=1,NTENS
       ATEMP1(I)=ZERO
       DO J=1,NTENS
        ATEMP1(I)=ATEMP1(I)+DCDDL(I,J)*STRANT(J)
       END DO
      END DO
      DO I=1,NTENS
       ATEMP2(I)=ZERO
       DO J=1,NTENS
        ATEMP2(I)=ATEMP2(I)+DCDDZ(I,J)*STRANT(J)
       END DO
      END DO
      DO I=1,NTENS
       ATEMP3(I)=ZERO
       DO J=1,NTENS
        ATEMP3(I)=ATEMP3(I)+DCDDT(I,J)*STRANT(J)
       END DO
      END DO
C
C     计算雅可比矩阵
C      
      DO I=1,NTENS
       DO J=1,NTENS
        DDSDDE(I,J)=CDFULL(I,J)+(ATEMP1(I)*DDLDDE(J)
     1     +ATEMP2(I)*DDZDDE(J)+ATEMP3(I)*DDTDDE(J))
     2     *DTIME/(DTIME+ETA)
       END DO
      END DO
C
      IF(DL.GT.0.99.OR.DT.GT.0.99.OR.
     1   DZ.GT.0.99) THEN
       ESTATE=1.
      ENDIF
C     
C     TO UPDATE THE STATE VARIABLE
C      
      STATEV(1)=DLV
      STATEV(2)=DTV
      STATEV(3)=DZV
      STATEV(4)=DL
      STATEV(5)=DT
      STATEV(6)=DZ
      STATEV(7)=DCL
      STATEV(8)=DTL
      STATEV(9)=DCT
      STATEV(10)=DTT
      STATEV(11)=DCZ
      STATEV(12)=DTZ
      STATEV(19)=XTLI
      STATEV(20)=XTLF
      STATEV(21)=XCLI
      STATEV(22)=XCLF
      STATEV(23)=XTTI
      STATEV(24)=XTTF
      STATEV(25)=XCTI
      STATEV(26)=XCTF
      STATEV(27)=XTZI
      STATEV(28)=XTZF
      STATEV(29)=XCZI
      STATEV(30)=XCZF
      STATEV(31)=ESTATE
C      open(6, file='results.txt',STATUS='UNKNOWN')
C      WRITE (6,'(1x,f10.8)) DT
C      CLOSE(6)   
C      
C     TO COMPUTE THE ENERGY
C
      DO I=1,NDI
       SSE=SSE+HALF*(STRESS(I)+OLD_STRESS(I))*DSTRAN(I)
      END DO
      DO I=NDI+1,NTENS
       SSE=SSE+(STRESS(I)+OLD_STRESS(I))*DSTRAN(I)
      END DO
C     TO COMPUTE THE INTERNAL,ENERGY WITHOUT VISCOUS REGULARIZATION
      DO I=1,NDI
       SCD=SCD+HALF*(STRESS(I)+OLD_STRESS(I)
     1     -D_STRESS(I)-DOLD_STRESS(I))*DSTRAN(I)
      END DO
      DO I=NDI+1,NTENS
       SCD=SCD+(STRESS(I)+OLD_STRESS(I)
     1    -D_STRESS(I)-DOLD_STRESS(I))*DSTRAN(I)
      END DO     

      RETURN
      END

C*****************************************************************************************************
!   CALCULATE THE STRESS BASEHD ON THE DAMAGE VARAIES
C*****************************************************************************************************

      SUBROUTINE GETSTRESS(CFULL,CDFULL,DLV,DTV,DZV,
     1            STRESS,STRANT,NDI,NTENS)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION CFULL(6,6),CDFULL(6,6),STRESS(6),STRANT(6)
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)
      DO K1=1,6
       DO K2=1,6
        CDFULL(K1,K2)=CFULL(K1,K2)
       END DO
      END DO
      IF((DLV .GT.ZERO).OR.(DTV.GT.ZERO).OR.(DZV.GT.ZERO))THEN
       RNL=ONE-DLV
       RNT=ONE-DTV
       RNZ=ONE-DZV
       RNLT=(TWO*(ONE-DLV)*(ONE-DTV))/(TWO-DLV-DTV)
       RNZL=(TWO*(ONE-DZV)*(ONE-DLV))/(TWO-DZV-DLV)
       RNTZ=(TWO*(ONE-DTV)*(ONE-DZV))/(TWO-DTV-DZV)
C      材料损失后的刚度矩阵
       CDFULL(1,1)=CFULL(1,1)*RNL**TW0
       CDFULL(2,2)=CFULL(2,2)*RNT**TWO
       CDFULL(3,3)=CFULL(3,3)*RNZ**TWO
       CDFULL(1,2)=CFULL(1,2)*RNL*RNT
       CDFULL(2,1)=CFULL(1,2)
       CDFULL(1,3)=CFULL(1,3)*RNL*RNZ
       CDFULL(3,1)=CFULL(1,3)
       CDFULL(2,3)=CFULL(2,3)*RNT*RNZ
       CDFULL(3,2)=CFULL(2,3)
       CDFULL(4,4)=CFULL(4,4)*RNLT**TWO
       CDFULL(5,5)=CFULL(5,5)*RNZL**TWO
       CDFULL(6,6)=CFULL(6,6)*RNTZ**TWO
      END IF
      print *, CDFULL
C     UPDATE THE STRESS STATE IF 3D CASE
      DO I=1,NTENS
       STRESS(I)=ZERO
       DO J=1,NTENS
        STRESS(I)=STRESS(I)+CDFULL(I,J)*STRANT(J)
       END DO
      END DO
      RETURN
      END
C*****************************************************************************************************
C   CALCULATE THE STRESS BASEHD ON THE DAMAGE VARAIES
C*****************************************************************************************************   
      SUBROUTINE CheckFailureIni(STRANT,STRESS,DL,DT,DZ,DCL,DTL,DCT,DTT,
     1     DTZ,DTLOLD,DCLOLD,DTTOLD,DCTOLD,DTZOLD,DCZOLD,CELENT,CFULL,
     2     XTL,XCL,XTT,XCT,XTZ,XCZ,XTLI,XTLF,XCLI,XCLF,XTTI,XTTF,XCTI,
     3     XCTF,XTZI,XTZF,XCZI,XCZF,DCZ)  
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION DTLDDE(6),DCLDDE(6),DTTDDE(6),DCTDDE(6),
     1        DTZDDE(6),DCZDDE(6),STRANT(6),DDLDDE(6),
     2        DDTDDE(6),DDZDDE(6),STRESS(6)
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)
      IF(STRESS(1).GT.ZERO)THEN
          FTL=(STRESS(1)/XT)**2+(STRESS(4)/SLT)**2+(STRESS(5)/SLZ)**2   
          IF(STRANT(1).LE.ZERO)THEN
             FTL=ZERO
          END IF
      ENDIF
c      print *, FTL
      IF (FTL.GT.ONE.AND.XTLF.EQ.ZERO) THEN
          XTL=CELENT*SQRT(STRANT(1)**2+STRANT(4)**2+STRANT(5)**2)
          SQTL=CELENT*(STRESS(1)*STRANT(1)+STRESS(4)
     1         *STRANT(4)+STRESS(5)*STRANT(5))/XTL
      END IF   
      CALL DAMAGECHECK(FTL,STRESS,STRANT,GF,1,4,5,CMNAME,CELENT,
     1    SQTL,XTLF,XTLI,XTL,DTL,DTLOLD,DTLDDE)
      IF (DTL.GT.DTLOLD) THEN
              DTLDDE(1)=(XTLF/(XTLF-XTLI)-DTL)*STRANT(1)/
     1         (STRANT(1)**2+STRANT(4)**2+STRANT(5)**2)/(CELENT**2)
              DTLDDE(4)=(XTLF/(XTLF-XTLI)-DTL)*STRANT(4)/
     1         (STRANT(1)**2+STRANT(4)**2+STRANT(5)**2)/(CELENT**2)
              DTLDDE(5)=(XTLF/(XTLF-XTLI)-DTL)*STRANT(5)/
     1         (STRANT(1)**2+STRANT(4)**2+STRANT(5)**2)/(CELENT**2)
      END IF
         DTL=MAX(DTL,DTLOLD)      
C*****************************************************************************************************
       IF(STRESS(1).LE.ZERO)THEN
          FCL=(STRESS(1)/XC)**2
          IF(STRANT(1).GT.ZERO)THEN
            FCL=ZERO
          END IF
       ENDIF
       IF (FCL.GT.ONE.AND.XCLF.EQ.ZERO) THEN
           XCL=CELENT*(-STRANT(1))
           SQCL=CELENT*(-STRESS(1))*(-STRANT(1))/XCL
       END IF
       CALL DAMAGECHECK(FCL,STRESS,STRANT,GF,1,CMNAME,CELENT,
     1         SQCL,XCLF,XCLI,XCL,DCL,DCLOLD,DCLDDE)
       IF (DCL.GT.DCLOLD) THEN
           DCLDDE(1)=(XTLF/(XTLF-XTLI)-DCL)/(-STRANT(1))
       END IF
         DCL=MAX(DCL,DCLOLD)
C*****************************************************************************************************
       IF (XCLF.GT.ZERO)THEN
        DL=DCL
        DO I=1,6
         DDLDDE(I)=DCLDDE(I)
        END DO
       ELSE IF (XTLF.GT.ZERO)THEN
        DL=DTL
        DO I=1,6
         DDLDDE(I)=DTLDDE(I)
        END DO
       ELSE
        DL=ZERO
       ENDIF
        DL=MAX(DCL,DTL)      
       IF (XCLF.GT.ZERO.OR.XTLF.GT.ZERO)THEN
        DLV=ETA/(ETA+DTIME)*DLVOLD+DTIME/(ETA+DTIME)*DL
       END IF

C*****************************************************************************************************      
C*****************************************************************************************************
      IF(STRESS(2).GT.ZERO)THEN
        FTT=((STRESS(2)+STRESS(3))/YT)**2
     1     +(STRESS(6)**2-STRESS(2)*STRESS(3))/(STZ**2)
     2     +(STRESS(4)/SLT)**2+(STRESS(5)/SLZ)**2
        IF(STRANT(2).LE.ZERO)THEN
         FTT=ZERO
        END IF
      END IF
      IF (FTT.GT.ONE.AND.XTTF.EQ.ZERO) THEN
          XTT=CELENT*SQRT(STRANT(2)**2+STRANT(4)**2+STRANT(6)**2)
          SQTT=CELENT*(STRESS(2)*STRANT(2)+STRESS(4)
     1         *STRANT(4)+STRESS(6)*STRANT(6))/XTT        
      END IF
      CALL DAMAGECHECK(FTT,STRESS,STRANT,GF,2,4,6,CMNAME,CELENT,
     1           SQTT,XTTF,XTTI,XTT,DTT,DTTOLD,DTTDDE)
      IF (DTT.GT.DTTOLD) THEN
              DTTDDE(2)=(XTTF/(XTTF-XTTI)-DTT)*STRANT(2)/
     1         (STRANT(2)**2+STRANT(4)**2+STRANT(6)**2)/(CELENT**2)
              DTTDDE(4)=(XTTF/(XTTF-XTTI)-DTT)*STRANT(4)/
     1         (STRANT(2)**2+STRANT(4)**2+STRANT(6)**2)/(CELENT**2)
              DTTDDE(6)=(XTTF/(XTTF-XTTI)-DTT)*STRANT(6)/
     1         (STRANT(2)**2+STRANT(4)**2+STRANT(6)**2)/(CELENT**2)
      END IF
         DTT=MAX(DTT,DTTOLD)
C*****************************************************************************************************
      IF(STRESS(2).LE.ZERO)THEN
        FCT=YC/((YC/(2*STZ))**2-1)*(STRESS(2)+STRESS(3))+((STRESS(2)
     1      +STRESS(3))/(2*STZ))**2+(STRESS(6)**2-STRESS(2)*STRESS(3))
     2      /(STZ**2)+(STRESS(4)/SLT)**2+(STRESS(5)/SLZ)**2

        IF(STRANT(2).GT.ZERO)THEN
         FCT=ZERO
        END IF
      ENDIF
      IF (FCT.GT.ONE.AND.XCTF.EQ.ZERO) THEN
           XCT=CELENT*(-STRANT(2))
           SQCT=CELENT*(-STRESS(2))*(-STRANT(2))/XCT
      END IF
      CALL DAMAGECHECK(FCT,STRESS,STRANT,GF,2,CMNAME,CELENT,
     1        SQCT,XCTF,XCTI,XCT,DCT,DCTOLD,DCTDDE)
      IF (DCT.GT.DCTOLD) THEN
              DCTDDE(2)=(XCTF/(XCTF-XCTI)-DCT)/(-STRANT(2))
      END IF
         DCT=MAX(DCT,DCTOLD)
C*****************************************************************************************************
       IF (XCTF.GT.ZERO)THEN
        DT=DCT
        DO I=1,6
         DDTDDE(I)=DCTDDE(I)
        END DO
       ELSE IF (XTTF.GT.ZERO)THEN
        DT=DTT
        DO I=1,6
         DDTDDE(I)=DTTDDE(I)
        END DO
       ELSE
        DT=ZERO
       ENDIF
         DT=MAX(DCT,DTT)      
       IF (XCTF.GT.ZERO.OR.XTTF.GT.ZERO)THEN
        DTV=ETA/(ETA+DTIME)*DTVOLD+DTIME/(ETA+DTIME)*DT
       END IF

C*****************************************************************************************************
C*****************************************************************************************************
      IF(STRESS(3).GT.ZERO)THEN
        FTZ=((STRESS(2)+STRESS(3))/YT)**2+(STRESS(6)**2-STRESS(2)
     1   *STRESS(3))/(STZ**2)+(STRESS(4)/SLT)**2+(STRESS(5)/SLZ)**2
        IF(STRANT(3).LE.ZERO)THEN
         FTZ=ZERO
        END IF
      ENDIF
      IF (FTZ.GT.ONE.AND.XTZF.EQ.ZERO) THEN
          XTZ=CELENT*SQRT(STRANT(3)*2+STRANT(6)**2+E31**2)
          SQTZ=CELENT*(STRESS(3)*STRANT(3)+STRESS(6)
     1         *STRANT(6)+S31*E31)/XTZ
      END IF
      CALL DAMAGECHECK(FTZ,STRESS,STRANT,GF,3,5,6,CMNAME,CELENT,
     1               SQTZ,XTZF,XTZI,XTZ,DTZ,DTZOLD,DTZDDE)
      IF (DTZ.GT.DTZOLD) THEN
               DTZDDE(3)=(XTZF/(XTZF-XTZI)-DTZ)*STRANT(3)/
     1         (STRANT(3)**2+E31**2+STRANT(6)**2)/(CELENT**2)
               DTZDDE(5)=(XTZF/(XTZF-XTZI)-DTZ)*E31/
     1         (STRANT(3)**2+E31**2+STRANT(6)**2)/(CELENT**2)
               DTZDDE(6)=(XTZF/(XTZF-XTZI)-DTZ)*STRANT(6)/
     1         (STRANT(3)**2+E31**2+STRANT(6)**2)/(CELENT**2)
      END IF
         DTZ=MAX(DTZ,DTZOLD)
C*****************************************************************************************************         
       IF(STRESS(3).LE.ZERO)THEN
        FCZ=YC/((YC/(2*STZ))**2-1)*(STRESS(2)+STRESS(3))+((STRESS(2)
     1    +STRESS(3))/(2*STZ))**2+(STRESS(6)**2-STRESS(2)*STRESS(3))
     2    /(STZ**2)+(STRESS(4)/SLT)**2+(STRESS(5)/SLZ)**2
        IF(STRANT(3).GT.ZERO)THEN
         FCZ=ZERO
        END IF
       ENDIF
       IF (FCZ.GT.ONE.AND.XCZF.EQ.ZERO) THEN
           XCZ=CLENT*(-STRANT(3))
           SQCZ=CELENT*(-STRESS(3))*(-STRANT(3))/XCZ
       END IF
       CALL DAMAGECHECK(FCZ,STRESS,STRANT,GF,3,CMNAME,CELENT,
     1               SQCZ,XCZF,XCZI,XCZ,DCZ,DCZOLD,DCZDDE)      
       IF (DCZ.GT.DCZOLD) THEN
            DCZDDE(3)=(XCZF/(XCZF-XCZI)-DCZ)/(-STRANT(3))
       END IF
           DCZ=MAX(DCZ,DCZOLD)
C*****************************************************************************************************
       IF (XCZF.GT.ZERO)THEN
        DZ=DCZ
        DO I=1,6
         DDZDDE(I)=DCZDDE(I)
        END DO
       ELSE IF (XTZF.GT.ZERO)THEN
        DZ=DTZ
        DO I=1,6
         DDZDDE(I)=DTZDDE(I)
        END DO
       ELSE
        DZ=ZERO
       ENDIF
       DZ=MAX(DCZ,DTZ)
       IF (XCZF.GT.ZERO.OR.XTZF.GT.ZERO)THEN
        DZV=ETA/(ETA+DTIME)*DZVOLD+DTIME/(ETA+DTIME)*DZ
       END IF   
      RETURN
      END

C*****************************************************************************************************
C              INITIATION DAMAGE CHECK
C*****************************************************************************************************

      SUBROUTINE  DAMAGECHECK(FX,CMNAME,CELENT,GF,  
     1                  SQX,XXF,XXI,XX,DD,DDOLD,DDXDDE)
      INCLUDE 'ABA_PARAM.INC'
      CHARACTER*80 CMNAME
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,HALF=0.5D0)
      IF(FX.GT.ONE.AND.XXF.EQ.ZERO)THEN
        XXF=TWO*GF*sqrt(FX)/SQX
        XXI=XX/sqrt(FX)
      END IF
      IF(XXF.GT.ZERO)THEN
         IF(XX.GE.XXF)THEN
          DD=0.99
         ELSE IF(XX.LE.XXI) THEN
          DD=0.
         ELSE
          DD=XXF*(XX-XXI)/(XX*(XXF-XXI))
         ENDIF
         IF(DD.LT.DDOLD)THEN
          DD=DDOLD
         ENDIF
      END IF
      RETURN
      END
C*****************************************************************************************************
C     GET THE DERIVATIVE MATRIX OF DAMAGEID MATRIX OVER DAMAGE VARIABLE
C*****************************************************************************************************

      SUBROUTINE ELASTICDERIVATIVE(CFULL,DLV,DTV,DZV,
     1       DCDDL,DCDDT,DCDDZ)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION CFULL(6,6),DCDDT(6,6),DCDDZ(6,6),DCDDL(6,6)
      PARAMETER (ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,EIGHT=8.D0)
C     INITIALIZE THE DATA TO ZERO
      DO I=1,6
       DO J=1,6
        DCDDT(I,J)=ZERO
        DCDDL(I,J)=ZERO
        DCDDZ(I,J)=ZERO
       END DO
      END DO
C     bL,BT,BZ
      RNL=ONE-DLV
      RNT=ONE-DTV
      RNZ=ONE-DZV
C     BLT,BTZ,BZL
      RNTZ=(TWO*(ONE-DTV)*(ONE-DZV))/(TWO-DTV-DZV)
      RNZL=(TWO*(ONE-DZV)*(ONE-DLV))/(TWO-DZV-DLV)
      RNLT=(TWO*(ONE-DLV)*(ONE-DTV))/(TWO-DLV-DTV)
C     CAL CULCULATE DC/DDL
      DCDDL(1,1)=-TWO*RNL*CFULL(1,1)
      DCDDL(1,2)=-RNT*CFULL(1,2)
      DCDDL(2,1)= DCDDL(1,2)
      DCDDL(1,3)=-RNZ*CFULL(1,3)
      DCDDL(3,1)= DCDDL(1,3)
      DCDDL(4,4)=-EIGHT*RNT**THREE*RNL/((RNT+RNL)**THREE)
     1         *CFULL(4,4)
      DCDDL(5,5)=-EIGHT*RNZ**THREE*RNL/((RNZ+RNL)**THREE)
     1         *CFULL(5,5)
C     CALCULATE DC/DDT
      DCDDT(1,2)=-RNL*CFUlL(1,2)
      DCDDT(2,1)= DCDDT(1,2)
      DCDDT(2,2)=-TWO*RNT*CFULL(2,2)
      DCDDT(2,3)=-RNZ*CFULL(2,3)
      DCDDT(3,2)= DCDDT(2,3)
      DCDDT(4,4)=-EIGHT*RNL**THREE*RNT/((RNT+RNL)**THREE)
     1        *CFULL(4,4)
      DCDDT(6,6)=-EIGHT*RNZ**THREE*RNT/((RNZ+RNT)**THREE)
     1        *CFULL(6,6)
C     CALCULATE DC/DDZ
      DCDDZ(1,3)=-RNL*CFULL(1,3)
      DCDDZ(3,1)=-RNL*CFULL(1,3)
      DCDDZ(2,3)=-RNT*CFUlL(2,3)
      DCDDZ(3,2)=DCDDZ(2,3)
      DCDDZ(3,3)=-TWO*RNZ*CFULL(3,3)
      DCDDZ(5,5)=-EIGHT*RNL**THREE*RNZ/((RNL+RNZ)**THREE)
     1       *CFULL(5,5)
      DCDDZ(6,6)=-EIGHT*RNT**THREE*RNZ/((RNZ+RNT)**THREE)
     1       *CFULL(6,6)
      RETURN
      END