Holzapfel Gasser Ogden VUMAT with damage

Hi all, let me start by saying that I wish I'd found this forum a month or 3 earlier.

For my graduation at the Erasmus MC (Netherlands) I'm researching the role of fibre dispersion and orientation in the damage initiation and propagation of bio-engineered tissue (before moving to actual tissue one day). My predecessor made a UMAT which incorporates the HGO model. From experiments in our group I was able to find the local fibre orientation and dispersion factor.

Now to calculate/represent the damage propagation in the modelled tissue, I added a damage factor to the model, which said that the element should be deleted once damage reached 1. However, as soon as the first element is deleted, the simulation diverges, gets negative eigenvalues, and stops working. I guess that is because when the element gets damaged, it is unable to take stresses, leading to no stiffness, leading to the infamous and dreaded buckling/ negative eigenvalues. From what I've read, my conclusion now is that only an Explicit simulation with a VUMAT will allow for the deletion of elements during a simulation.

Now the options I consider are the following:

1. Use partitions in my CAE, and give local material properties to certain tiles (including fibre dispersion and orientation) using the built in HGO model and add the damage to the model using VUSDFLD. However the damage is based on the strain energy density, which I don't think can be found as an output using built in material models.

2. Convert my UMAT into VUMAT. This seems the better option as it leads to no black boxes.

My question (finally) is:

1. Could anyone help out with converting my UMAT to a VUMAT?

2. On a more helicopter view, do you think this method of simulating damage propagation in fibrous tissue seems logical? 

I've added the code below, so hopefully we don't have to reinvent the wheel here together.

C      
      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),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
C
C ----------------------------------------------------------------
C LOCAL ARRAYS
C ----------------------------------------------------------------
C
      DIMENSION XKIRCH(3,3), DUMSTRSS(3,3), XAMAT(3,2), XI4(2), R(3,3),
     1 U(3,3), RT(3,3), XADUM(3,1), XADUM1(3,1)
C           
C
      PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0,
     1 SIX=6.D0)
C
C ----------------------------------------------------------------
C     NUMBER OF FIBRE FAMILIES
C ----------------------------------------------------------------
C
      NANISO = 1        
C
C ----------------------------------------------------------------
C     MATERIAL PROPERTIES
C ----------------------------------------------------------------
      C10 = PROPS(1)
      D1 = PROPS(2)
      xk1 = PROPS(3)
      xk2 = PROPS(4)
      xkap = PROPS(5)   
      THETAD1 = PROPS(6)
      THETAD2 = PROPS(7)
	  T1 = PROPS(8)
	  T2 = PROPS(9)
        
C       Convert degrees to radians        
      XPI = 3.14159265359
      THETAR1 = THETAD1*XPI/180.0
      THETAR2 = THETAD2*XPI/180.0
C        
C     Fibre vectors in the reference configuration
      XAMAT(1,1)=COS(thetar1); XAMAT(2,1)= SIN(thetar1); XAMAT(3,1)=0.
      XAMAT(1,2)=COS(thetar2); XAMAT(2,2)= SIN(thetar2); XAMAT(3,2)=0.
C
C     This is related to the calculation of the numberical Jacobian
C     0 indicates that the stress calculation is not part of the 
C     stress calculation.        
      iter=0       
C
C---------------------------------------------------------------------------
C         POLAR DECOMPOSITION
C---------------------------------------------------------------------------            
      call kpolarDecomp(DFGRD1, U, R)
C     Calculate the transpose of the rotation matrix      
      RT = transpose(R)
      call kprinter(R,3,3)
C
C-----------------------------------------------------------------
C        CALCULATE LOCAL FIBRE VECTOR IN THE REFERENCE CONFIG.
C-----------------------------------------------------------------                 
C    
      DO IAN = 1,NANISO
C         Pick out fibre vector          
          XADUM(:,1) = XAMAT(:,IAN)
C          
C         Update to local basis system       
          XADUM1 = MATMUL(RT,XADUM)
C
C         Store result          
          XAMAT(:,IAN) = XADUM1(:,1)         
      END DO        
C	  
C-----------------------------------------------------------------
C     ZERO THE TANGENT MATRIX
C-----------------------------------------------------------------        
C
      DO I=1,NTENS
      DO J=1,NTENS
        DDSDDE(I,J) = ZERO
      END DO
      END DO          
C      
C-----------------------------------------------------------------
C     CALCULATE THE KIRCHHOFF STRESS 
C-----------------------------------------------------------------
C     
      call kstress_calc(DFGRD1, C10, D1, xkap, xk1, xk2, T1,T2, XKIRCH, XJ,
     1 XAMAT, XI1, XI4, NPT, NANISO, iter, XWE, XWK, RFACTOR,TK, EA, D2)
C     
        STATEV(1) = XI1
        STATEV(2) = XJ
        STATEV(3) = EA
        STATEV(4) = XI4(1)
        STATEV(5) = XI4(2)
        STATEV(6) = XWE
        STATEV(7) = XWK	
        STATEV(8) = RFACTOR	
        STATEV(9) = TK
        STATEV(10) = D2
        
C             
C-----------------------------------------------------------------
C CONVERT KIRCHHOFF STRESS TO CAUCHY STRESS               
C-----------------------------------------------------------------      
      DUMSTRSS = XKIRCH / XJ           
C     
C     Convert the stress tensor to a Voigt vector      
      call kmatrix2vector(DUMSTRSS, STRESS, nshr)         
      
C      
C-----------------------------------------------------------------
C CALCULATE THE TANGENT MATRIX
C-----------------------------------------------------------------
C      
      call kCTM(XKIRCH,DFGRD1,NTENS,PROPS,XAMAT,NPT,ITER,
     + NANISO,DDSDDE,NPROPS,XJ,NSHR)      
C
      RETURN
C
C      
      ENDSUBROUTINE UMAT
C
C-----------------------------------------------------------------------------
C           SUBROUTINES
C-----------------------------------------------------------------------------   
C      
      SUBROUTINE kstress_calc(DGRAD, C10, D1, xkap, xk1, xk2, T1,T2,
     1  XKIRCH, XJ, XAMAT, XI1, XI4, NPT, NANISO, iter, XWE, XWK, RFACTOR,TK,EA,D2)
C     
      INCLUDE 'ABA_PARAM.INC'
C      
      INTENT(IN) :: DGRAD, C10, D1, xkap, xk1,xk2,XAMAT,NANISO,iter,T1,T2
      INTENT(OUT):: XKIRCH, XJ, XI1, XI4, XWE, XWK, RFACTOR, TK, EA, D2 
C      
      DIMENSION DGRAD(3,3), BMAT(3,3), XKIRCH(3,3), XAMAT(3,2),
     1 B2MAT(3,3), XA(3), a_curnt(3), aoa(3,3), XANISOK(3,3),
     2 XANISOTOT(3,3), XI4(2)
C      
      PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0)
C
C LOCAL ARRAYS
C ----------------------------------------------------------------    
C      
C     Determinant of the deformation gradient
      XJ=DGRAD(1, 1)*DGRAD(2, 2)*DGRAD(3, 3)
     1   -DGRAD(1, 2)*DGRAD(2, 1)*DGRAD(3, 3)
     2   +DGRAD(1, 2)*DGRAD(2, 3)*DGRAD(3, 1)
     3   +DGRAD(1, 3)*DGRAD(3, 2)*DGRAD(2, 1)
     4   -DGRAD(1, 3)*DGRAD(3, 1)*DGRAD(2, 2)
     5   -DGRAD(2, 3)*DGRAD(3, 2)*DGRAD(1, 1)
C     
C        
C     Calculate left Cauchy-Green deformation tensor     
      BMAT = MATMUL(DGRAD,TRANSPOSE(DGRAD))
C       
C-----------------------------------------------------------------------
C     CALCULCATE THE INVARIANTS                
C-----------------------------------------------------------------------
C     
C     I1 invariant      
      XI1 = BMAT(1,1)+BMAT(2,2)+BMAT(3,3)
C        
C     Square of the B matrix
      B2MAT = MATMUL(BMAT,BMAT)
C      
C--------------------------------------------------------------------
C     CALCULATE THE ANISOTROPIC PORTION OF THE KIRCH STRESS
C--------------------------------------------------------------------       
C  
C     Zero the matrix      
      DO I=1,3
        DO J=1,3
          XANISOTOT(I,J) = ZERO
        END DO
      END DO
C                  
C--------------------------------------------------------------------
C     Loop over the number of fibre families    
C--------------------------------------------------------------------        
      DO IANISO = 1,NANISO                                     
C         
C         Fibre vector in reference configuration  
          XA = XAMAT(1:3,IANISO);
C                                              
C         Calculate fibre vector in current configuration: a = F*A
          a_curnt = MATMUL(DGRAD,XA)
C
C         Calculate structural tensor: a \otimes a     
          DO i=1,3
              DO j=1,3
                aoa(i,j) = a_curnt(i) * a_curnt(j)
              ENDDO
          ENDDO            
C            
C         Calculate the fibre invariant
          XI4(IANISO) = aoa(1,1) + aoa(2,2) + aoa(3,3)
C
          EA = xkap*(XI1 - 3.D0) + (1.D0 - 3.D0*xkap)*(XI4(IANISO)-1.D0)
          EAA = xk2*(EA)**2
C
	          XWE = C10*(XI1 - 3.D0) 
     1          + (xk1/2.D0*xk2)*((exp(EAA)) - 1.D0)
     2          + (1.D0/D1)*((XJ**TWO - 1.D0)/2.D0 - LOG(XJ))
C
	          TK = SQRT(TWO * XWE)
C			  
	            IF (TK .LT. T1) THEN
	                  DAM = ZERO
		      ELSE IF (TK .LE. T2) THEN
		            DAM = ((TK - T1) / (T2 - T1))**TWO 
			ELSE
	                  DAM = ONE
	            END IF
C
	          XWK = (1.D0 - DAM)*(C10*(XI1 - 3.D0) 
     1           + (xk1/2.D0*xk2)*((exp(EAA)) - 1.D0)) 
     2           + (1.D0/D1)*((XJ**TWO - 1.D0)/2.D0 - LOG(XJ))     
C 
              RFACTOR = ONE - DAM 	
C     If damage is higher than 0.2, delete element.                     		 	   
C
            IF (DAM .LE. 0.2) THEN
                D2= 1
            ELSE
                D2=0
            END IF 
    
        

C     CALCULATE THE ISOTROPIC PORTION OF THE KIRCHHOFF STRESS
C
C           PART 1
            COEFF1 = TWO*(C10 + exp(EAA)*xkap*EA*xk1)/(XJ**(TWO/THREE))
C
C            KIRCHHOFF STRESS PART 2
C
            TRBMAT= COEFF1*(BMAT(1,1)+BMAT(2,2)+BMAT(3,3))/THREE
C
C            KIRCHHOFF STRESS PART 1
C      
              XKIRCH = COEFF1 * BMAT
C
C            SUBTRACT THE PART 2
C    
               DO I = 1,3
                 XKIRCH(I,I) =XKIRCH(I,I) - TRBMAT
               END DO
               XKIRCH = XKIRCH * RFACTOR
C
C         FORM VOLUMETRIC PART, 3 
C      
             COEFF3 = 2*(XJ-ONE)*XJ/D1
C
C         ADD TO THE PREVIOUS PARTS
C      
                DO I = 1,3
                   XKIRCH(I,I) = XKIRCH(I,I) + COEFF3
                END DO
C           
C   CALCULATE THE ANISOTROPIC PORTION OF THE KIRCHHOFF STRESS   
C          
C             Calculate the various parts of the equation for kirch stress
              p1 = 2.0*xk1;                      
C
              p2 = EA*(1.D0 - 3.D0*xkap)                  
C
              p4 = exp(EAA)
C
              XANISOK = RFACTOR*p1*p2*p4*aoa
C               
C             Sum the stress over the fibre families              
              XANISOTOT = XANISOTOT + XANISOK     
C        
      END DO
C 
C     Add the isotropic and anisotropic parts      
                XKIRCH = XKIRCH + XANISOTOT
C	
      RETURN
      END SUBROUTINE kstress_calc
C      
      SUBROUTINE kpolarDecomp(G,S,R)
c
c	declarations
      INCLUDE 'ABA_PARAM.INC'
      dimension G(3,3),R(3,3),S(3,3),COG(3,3),XP(3,3),ADP(3,3),
     1 XPBI(3,3)
c
c	constants
      EPS=1.0E-5
c
c	cofactors and determinant of g
      COG(1,1)=G(2,2)*G(3,3)-G(2,3)*G(3,2)
      COG(2,1)=G(1,3)*G(3,2)-G(1,2)*G(3,3)
      COG(3,1)=G(1,2)*G(2,3)-G(1,3)*G(2,2)
      COG(1,2)=G(2,3)*G(3,1)-G(2,1)*G(3,3)
      COG(2,2)=G(1,1)*G(3,3)-G(1,3)*G(3,1)
      COG(3,2)=G(1,3)*G(2,1)-G(1,1)*G(2,3)
      COG(1,3)=G(2,1)*G(3,2)-G(2,2)*G(3,1)
      COG(2,3)=G(1,2)*G(3,1)-G(1,1)*G(3,2)
      COG(3,3)=G(1,1)*G(2,2)-G(1,2)*G(2,1)
      G3=G(1,1)*COG(1,1)+G(2,1)*COG(2,1)+G(3,1)*COG(3,1)
c
c	P = trans(G) * G = S * S
      DO 10000 I=1,3
      XP(I,1)=G(1,I)*G(1,1)+G(2,I)*G(2,1)+G(3,I)*G(3,1)
      XP(I,2)=G(1,I)*G(1,2)+G(2,I)*G(2,2)+G(3,I)*G(3,2)
      XP(I,3)=G(1,I)*G(1,3)+G(2,I)*G(2,3)+G(3,I)*G(3,3)
10000 CONTINUE
c
c	adjoint of P
      ADP(1,1)=XP(2,2)*XP(3,3)-XP(2,3)*XP(3,2)
      ADP(2,2)=XP(1,1)*XP(3,3)-XP(1,3)*XP(3,1)
      ADP(3,3)=XP(1,1)*XP(2,2)-XP(1,2)*XP(2,1)
c
c	G invariants
      G1SQ=XP(1,1)+XP(2,2)+XP(3,3)
      G1=SQRT(G1SQ)
      G2SQ=ADP(1,1)+ADP(2,2)+ADP(3,3)
      G2=SQRT(G2SQ)
c
c	initialize iteration
      H1=G2/G1SQ
      H2=G3*G1/G2SQ
      X=1.0
      Y=1.0
c
c	iteration loop
10001 CONTINUE
      DEN=2.0*(X*Y-H1*H2)
      RES1=1.0-X*X+2.0*H1*Y
      RES2=1.0-Y*Y+2.0*H2*X
      DX=(Y*RES1+H1*RES2)/DEN
      DY=(H2*RES1+X*RES2)/DEN
      X=X+DX
      Y=Y+DY
      IF(ABS(DX/X).GT.EPS.OR.ABS(DY/Y).GT.EPS)GO TO 10001
c
c	BETA invariants
      BETA1=X*G1
      BETA2=Y*G2
c
c	invert ( trans(G) * G + BETA2 * identity )
      XP(1,1)=XP(1,1)+BETA2
      XP(2,2)=XP(2,2)+BETA2
      XP(3,3)=XP(3,3)+BETA2
      XPBI(1,1)=XP(2,2)*XP(3,3)-XP(2,3)*XP(3,2)
      XPBI(1,2)=XP(1,3)*XP(3,2)-XP(1,2)*XP(3,3)
      XPBI(1,3)=XP(1,2)*XP(2,3)-XP(1,3)*XP(2,2)
      XPBI(2,1)=XP(2,3)*XP(3,1)-XP(2,1)*XP(3,3)
      XPBI(2,2)=XP(1,1)*XP(3,3)-XP(1,3)*XP(3,1)
      XPBI(2,3)=XP(1,3)*XP(2,1)-XP(1,1)*XP(2,3)
      XPBI(3,1)=XP(2,1)*XP(3,2)-XP(2,2)*XP(3,1)
      XPBI(3,2)=XP(1,2)*XP(3,1)-XP(1,1)*XP(3,2)
      XPBI(3,3)=XP(1,1)*XP(2,2)-XP(1,2)*XP(2,1)
      DETPBI=XP(1,1)*XPBI(1,1)+XP(2,1)*XPBI(1,2)+XP(3,1)*XPBI(1,3)
c
c	R = (cofac(G)+BETA1*G) * inv(trans(G)*G+BETA2*identity)
      DO 10002 I=1,3
      R(I,1)=((COG(I,1)+BETA1*G(I,1))*XPBI(1,1)
     1       +(COG(I,2)+BETA1*G(I,2))*XPBI(2,1)
     2       +(COG(I,3)+BETA1*G(I,3))*XPBI(3,1))/DETPBI
      R(I,2)=((COG(I,1)+BETA1*G(I,1))*XPBI(1,2)
     1       +(COG(I,2)+BETA1*G(I,2))*XPBI(2,2)
     2       +(COG(I,3)+BETA1*G(I,3))*XPBI(3,2))/DETPBI
      R(I,3)=((COG(I,1)+BETA1*G(I,1))*XPBI(1,3)
     1       +(COG(I,2)+BETA1*G(I,2))*XPBI(2,3)
     2       +(COG(I,3)+BETA1*G(I,3))*XPBI(3,3))/DETPBI
10002 CONTINUE
c
c	S = trans(R) * G
      DO 10003 I=1,3
      S(I,1)=R(1,I)*G(1,1)+R(2,I)*G(2,1)+R(3,I)*G(3,1)
      S(I,2)=R(1,I)*G(1,2)+R(2,I)*G(2,2)+R(3,I)*G(3,2)
      S(I,3)=R(1,I)*G(1,3)+R(2,I)*G(2,3)+R(3,I)*G(3,3)
10003 CONTINUE
c
c	done
      RETURN
      END
C      
C     Subroutine to calculate the tangent matrix 
C     Subroutine to calculate the consistent tangent matrix using the 
C     perturbation method outlined by Miehe and Sun
C      
C          
      subroutine kCTM(XKIRCH1,DFGRD1,NTENS,PROPS,XAMAT,NPT,ITER,
     + NANISO,DDSDDE,NPROPS,XJ1,NSHR)
C      
      INCLUDE 'ABA_PARAM.INC'  
C      
      DIMENSION XKIRCH1(3,3), DFP(3,3), DFGRD_PERT(3,3),XKIRCH_PERT(3,3)
     1 , CMJ(3,3), CMJVEC(NTENS), ilist(6), jlist(6), XAMAT(3,2), XI4(2)
     2 , DFGRD1(3,3), DDSDDE(NTENS,NTENS), PROPS(NPROPS)
C
C      
C LOCAL ARRAYS
C ----------------------------------------------------------------
C XKIRCH1         - True Kirchhoff stress (i.e. not based on perturbed def. grad.)
C DFGRD1          - True deformation gradient for this increment      
C DFP             - Increment of the perturbed deformation gradient
C DFGRD_PERT      - Perturbed def. grad.
C XKIRCH_PERT     - Perturbed Kirchhoff stress
C CMJ             - (:,:,I,J) Components of material jacobian
C CMJVEC          - Above in vector form
C ILIST, JLIST    - Set of index components to be perturbed upon
C DUMSTRSS        - DUMMY STRESS TENSOR
C XAMAT           - array containing reference configuration fibre vectors
C XI4             - array storing the fibre invariant (will not be used in this)
C DDSDDE          - Voigt notation matrix to store the tangent matrix 
C PROPS           - Properties read in from Abaqus .inp file
C ----------------------------------------------------------------
C          
C     Perturbation parameter
      eps = 1.0e-08     
C
C     Material properties      
      C10 = PROPS(1)
      D1 = PROPS(2)
      xk1 = PROPS(3)
      xk2 = PROPS(4)
      xkap = PROPS(5)      
      THETAD1 = PROPS(6)
      THETAD2 = PROPS(7)
	 T1 = PROPS(8)
	 T2 = PROPS(9)    
C                
C     The Kirchhoff stress is perturbed six times over the
C     i, j'th components of the deformation gradient. This
C     array lists each of the componets to be perturbed over.
C         i        |      j      
      ilist(1) = 1; jlist(1) = 1;
      ilist(2) = 2; jlist(2) = 2;
      ilist(3) = 3; jlist(3) = 3;
      ilist(4) = 1; jlist(4) = 2;
      ilist(5) = 1; jlist(5) = 3;
      ilist(6) = 2; jlist(6) = 3;        
C      
C     Loop over each of the six components of the deformation gradient
C     which are to be perturbed, each time calculating a new column of
C     the tangent matrix.      
      Perturbation: do iter = 1,NTENS          
C      
C-----------------------------------------------------------------
C         CREATE THE PERTURBED DEFORMATION GRADIENT
C-----------------------------------------------------------------     
C         Pick out the (i,j) component of the deformation gradient
C         to be perturbed      
          ii = ilist(iter)
          jj = jlist(iter)     
C     
C         Calculate the increment of the perturbed deformation gradient          
          call kdelF(ii, jj, DFGRD1, eps, DFP)
C          
C         Add to the "true" deformation gradient to create the perturbed
C         deformation gradient          
          DFGRD_PERT = DFGRD1 + DFP
C
C-----------------------------------------------------------------
C         CALCULATE KIRCHHOFF STRESS BASED ON THE PERTURBED DEFORMATION GRADIENT
C-----------------------------------------------------------------
C                             N.B.
          call kstress_calc(DFGRD_PERT, C10, D1, xkap, xk1, xk2, T1,T2, 
     1    XKIRCH_PERT, XJP, XAMAT, XI1, XI4, NPT, NANISO,
     2    iter, XWE, XWK, RFACTOR,TK, EA, D2)                    
C
C-----------------------------------------------------------------
C          
C         Difference between the perturbed(i,j) and unpert. stress
          CMJ = XKIRCH_PERT - XKIRCH1          
C         
C         Normalize by perturbation parameter and return to Cauchy stress       
          CMJ = CMJ/XJ1/eps
C      
C      
C-----------------------------------------------------------------
C         FORM THE TANGENT MATRIX
C-----------------------------------------------------------------      
C
C         Convert tensor to Voigt vector          
          call kmatrix2vector(CMJ, CMJVEC, NSHR)
C
C         Insert this vector into the iter'th column of DDSDDE          
          do insert = 1,NTENS
C          
              DDSDDE(insert,iter) = CMJVEC(insert)
C            
          end do
C		  
C
       end do Perturbation
       
       end subroutine kCTM
C      
C     Helper subroutines      
C * KDELF        -   Calculate the increment of the deformation gradient for
C                   a given perturbation in (i,j), with epsilon
C
C * KPRINTER     -   Print out a matrix of any size
C
C * KMTMS     -   Multiply two matrices
C
C * KMATRIX2VECTOR     -   Convert a 3x3 matrix to a 6x1 vector
C                   
C * KDOTPROD           - Dot product of two vectors
C
C------------------------------------------------------------      
      subroutine kdotprod(A, B, dotp, n)
      
      INCLUDE 'ABA_PARAM.INC'
      
      intent(in) :: A, B, n
      intent(out):: dotp      
      
      dimension A(n), B(n)
      dotp = 0.0
      
      do i = 1,n
        dotp = dotp + A(i)*B(i)
      end do
      
      end subroutine kdotprod    
C      
C-------------------------------------------------------      
C
C-------------------------------------------------------
C      
      subroutine kmatrix2vector(XMAT, VEC, NSHR)
      
      INCLUDE 'ABA_PARAM.INC'
      
      intent(in) :: XMAT, NSHR
      intent(out):: VEC
      
      dimension xmat(3,3), vec(6)
  
        do i=1,3
            vec(i) = xmat(i,i);
        end do
               
        vec(4) = xmat(1,2);
        
        IF (NSHR==3) then
            vec(5) = xmat(1,3);
            vec(6) = xmat(2,3);
        END IF
      
      end subroutine kmatrix2vector      
C
C-------------------------------------------------------      
C
C-------------------------------------------------------
C      
      SUBROUTINE KMTMS (M, N, L, A, KA, B, KB, C, KC)
      
      INCLUDE 'ABA_PARAM.INC'
C      
      intent(in) :: M, N, L, A, KA, B, KB, KC
      intent(out):: C      
C      
C
C    PRODUCT OF REAL MATRICES
C
      DIMENSION A(KA,N), B(KB,L), C(KC,L)
      DOUBLE PRECISION W
C       
C
      DO 30 J = 1,L
         DO 20 I = 1,M
            W = 0.D0
            DO 10 K = 1,N
               W = W + A(I,K) * B(K,J)
   10       CONTINUE
            C(I,J) = W
   20    CONTINUE
   30 CONTINUE
      RETURN
      END SUBROUTINE KMTMS      
C
C-------------------------------------------------------      
C
C-------------------------------------------------------
C      
      subroutine kprinter(tens, m, n)
      
      INCLUDE 'ABA_PARAM.INC'
      
      intent(in):: tens, m, n      
      
      dimension tens(m,n)
        
        write(6,*)
        do i = 1,m
        do j = 1,n
            write(6,'(e19.9)',advance='no'),tens(i,j)
        end do
        write(6,*)
        end do
        write(6,*)
      return
      end subroutine kprinter      
C
C-------------------------------------------------------      
C
C-------------------------------------------------------
C  
      subroutine kdelF(m, n, DGRAD, eps, DF)
      
      INCLUDE 'ABA_PARAM.INC'
      
      intent (in) :: DGRAD, eps, m, n
      intent (out):: DF      
        
C Input: the index's i & j; The current deformation gradient (DGRAD). The perturbation
C        increment (eps)
C
C Output: The perturbed increment DF
C
      dimension dyad1(3,3), dyad2(3,3), DGRAD(3,3), DF(3,3), DFp1(3,3)
        
c Zero the dyad matrices
c
      do i = 1,3
        do j = 1,3
            dyad1(i,j) = zero
            dyad2(i,j) = zero
        end do
      end do
      
c Place the 1's in the correct location        
      dyad1(m,n) = 1.0;
c
      dyad2(n,m) = 1.0;
c      
c      KMTMS (M, N, L, A, KA, B, KB, C, KC)  
      call KMTMS(3, 3, 3, dyad1, 3, DGRAD, 3, DFp1, 3)
      DF = DFp1
      
      
      call KMTMS(3, 3, 3, dyad2, 3, DGRAD, 3, DFp1, 3)
      DF = DF + DFp1
           
      
      DF = 0.5*DF*eps            
      
      end subroutine kdelF
C