C ======================================================================
C User Subroutine UMAT for Abaqus:
C  AT2 formulation using the temperature/displacement element and
C  the staggered time integration scheme.:
C ======================================================================
C
C Copyright(C) Gergely Molnar, CNRS, INSA of Lyon, LaMCoS, UMR5259
C
C gergely.molnar@insa-lyon.fr
C http://www.molnar-research.com
C
C
C This software is governed by the CeCILL-B license under French law and
C abiding by the rules of distribution of free software. You can  use, 
C modify and/or redistribute the software under the terms of the CeCILL-B
C license as circulated by CEA, CNRS and INRIA at the following URL
C "http://www.cecill.info". 
C
C As a counterpart to the access to the source code and  rights to copy,
C modify and redistribute granted by the license, users are provided only
C with a limited warranty  and the software's author,  the holder of the
C economic rights,  and the successive licensors  have only  limited
C liability. 
C
C In this respect, the user's attention is drawn to the risks associated
C with loading,  using,  modifying and/or developing or reproducing the
C software by the user in light of its specific status of free software,
C that may mean  that it is complicated to manipulate,  and  that  also
C therefore means  that it is reserved for developers  and  experienced
C professionals having in-depth computer knowledge. Users are therefore
C encouraged to load and test the software's suitability as regards their
C requirements in conditions enabling the security of their systems and/or 
C data to be ensured and,  more generally, to use and operate it in the 
C same conditions as regards security. 
C
C The fact that you are presently reading this means that you have had
C knowledge of the CeCILL-B
C
C ======================================================================
C
C This implementation solves the phase-field fracture problem using
C the Abaqus user subroutines UMAT and HEATVAL within a coupled
C temperature-displacement analysis.
C
C The displacement field represents the mechanical response, while
C the temperature field is used as an auxiliary degree of freedom
C to represent the phase-field damage variable.
C
C UMAT defines the mechanical constitutive response, including
C stiffness degradation due to damage, while HEATVAL introduces
C the source term and its derivative required to solve the
C phase-field evolution equation through the thermal problem.
C
C This approach allows the coupled fracture problem to be solved
C using standard Abaqus finite elements without implementing
C additional user-defined elements.
C
C This implementation is inspired by the work of Yousef Navidtehrani,
C Covadonga Betegon, and Emilio Martinez-Paneda:
C
C "A Unified Abaqus Implementation of the Phase Field Fracture Method
C  Using Only a User Material Subroutine"
C Materials, 14(8), 1913 (2021).
C
C Material properties are defined in the Abaqus input file (*.inp).
C
C UMAT material parameters:
C PROPS(1) = Young's modulus (E)
C PROPS(2) = Poisson's ratio (nu)
C PROPS(3) = Phase-field length scale parameter (lc)
C PROPS(4) = Critical fracture energy (gc)
C PROPS(5) = Time increment control (0 = automatic, 1 = fixed/direct)
C
C ==============================================================
C ==============================================================
C                       Subroutine UMAT 
C                      Elastic material
C ==============================================================
C ==============================================================
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,JSTEP,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),
     4 JSTEP(4)
C 
       PARAMETER (ONE=1.0,TWO=2.0,THREE=3.0,SIX=6.0,HALF=0.5,
     1 TEN=10.0,ZERO=0.0, 
     2 DENGMAX=0.5,DTMIN=1.0D-9) 
       DATA NEWTON,TOLER/40,1.D-6/ 
C       
       REAL*8 EPSC(6),CMAT(6,6),EIGV(3),ALPHAI(3)
C
       REAL*8 EMOD,ENU,CLPAR,GCPAR,ELAMEG,ELAMEL,HIST,
     1  ANISOSWT,ENGN,ENGP,PARK,PHASE,EPSVOL
     2  GFUN,DENG,DENGMAXV,DTIMENEXT,DTSWT
C 
       COMMON/KUSER/USRVAR(8)
C ----------------------------------------------------------- 
C          Material properties
C ----------------------------------------------------------- 
C PROPS(1) = Young's modulus (E)
C PROPS(2) = Poisson's ratio (nu)
C PROPS(3) = Length scale parameter (lc)
C PROPS(4) = Crack surface energy (gc)
C PROPS(5) = time integration switch (0: automatic, 1: fixed)
C
C ----------------------------------------------------------- 
C
C	Elastic properties
C
       EMOD=PROPS(1)
       ENU=PROPS(2)
       CLPAR=PROPS(3)
       GCPAR=PROPS(4)
       DTSWT=PROPS(5)
       ELAMEG = EMOD/(TWO*(ONE+ENU))
       ELAMEL = ELAMEG*TWO*ENU/(ONE-TWO*ENU)
       PARK = TOLER
       ANISOSWT = ONE
C
C Transfer material properties to HETVAL subroutine
       HIST = STATEV(1)
       STATEV(2) = GCPAR
       STATEV(3) = CLPAR
C
C   Damage phase-field
       PHASE = TEMP + DTEMP
       IF (PHASE.GT.ONE-TOLER) THEN
        PHASE = ONE-TOLER
       ENDIF
       IF (PHASE.LT.ZERO) THEN
        PHASE = ZERO
       ENDIF
       GFUN=(ONE-PHASE)**TWO+PARK
C
C ======================================================================
C   Strain components
C ======================================================================
       DO K1=1,6
        EPSC(K1)=ZERO
       END DO
       DO K1=1,NTENS
        STRAN(K1)=STRAN(K1)+DSTRAN(K1)
       END DO
       DO K1=1,NDI
        EPSC(K1)=STRAN(K1)
       END DO
       DO K1=1,NTENS-NDI
        EPSC(K1+3)=STRAN(NDI+K1)
       END DO
C
C ======================================================================
C	Stiffness tensor
C ======================================================================
C
C     ==================================================================
C     Calculating degrated elastic stiffness matrix
C     ==================================================================
       DO I=1,6
        DO J=1,6
        CMAT(I,J)=ZERO
        END DO
       END DO
       IF (PHASE.LT.TOLER) THEN
        CMAT(1,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(2,2)=CMAT(1,1)
        CMAT(3,3)=CMAT(1,1)
        CMAT(1,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(2,1)=CMAT(1,2)
        CMAT(1,3)=CMAT(1,2)
        CMAT(3,1)=CMAT(1,3)
        CMAT(2,3)=CMAT(1,2)
        CMAT(3,2)=CMAT(2,3)
        CMAT(4,4)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
        CMAT(5,5)=CMAT(4,4)
        CMAT(6,6)=CMAT(4,4)
        DO I=1,6
         DO J=1,6
          CMAT(I,J)=CMAT(I,J)*GFUN
         END DO
        END DO        
       ELSE
        CALL HMAT(EPSC,CMAT,ENU,EMOD,ANISOSWT,PHASE,PARK)
       ENDIF      
C       
        DO K1=1, NTENS
         DO K2=1, NTENS
          DDSDDE(K2, K1)=ZERO
         END DO
        END DO
        DO I=1,NDI
         DO J=1,NDI
          DDSDDE(I,J)=CMAT(I,J)
         END DO
        END DO
        DO I=1,NTENS-NDI
         DO J=1,NTENS
          DDSDDE(J,I+NDI)=CMAT(J,I+3)
          DDSDDE(I+NDI,J)=CMAT(I+3,J)
         END DO
        END DO
C
C ======================================================================
C	Calculate Stresses
C ======================================================================
C
       DO K1=1, NTENS
        STRESS(K1)=ZERO
        DO K2=1, NTENS
         STRESS(K1)=STRESS(K1)+DDSDDE(K1, K2)*STRAN(K2)
        END DO
       END DO
C
C ======================================================================
C	Energies
C ======================================================================
       DO K1=1,6
        EPSC(K1)=ZERO
       END DO
       DO K1=1,NDI
        EPSC(K1)=STRAN(K1)
       END DO
       DO K1=1,NTENS-NDI
        EPSC(K1+3)=STRAN(NDI+K1)
       END DO       
       CALL SPRIND(EPSC,EIGV,EIGVEC,2,3,3)
C
       ALPHA=ZERO
       EPSVOL = EIGV(1)+EIGV(2)+EIGV(3)
       IF ((EPSVOL.GT.TOLER).OR.(ANISOSWT.LT.TOLER)) THEN
        ALPHA=ONE
       ENDIF
C
       DO K1=1,3
        ALPHAI(K1)=ZERO
        IF ((EIGV(K1).GT.TOLER).OR.(ANISOSWT.LT.TOLER)) THEN
         ALPHAI(K1)=ONE
        ENDIF
       END DO
C       
       ENGP=(ELAMEL*(ALPHA*(EIGV(1)+EIGV(2)+EIGV(3)))**TWO)/
     1  TWO+ELAMEG*((EIGV(1)*ALPHAI(1))**TWO+(EIGV(2)*
     2  ALPHAI(2))**TWO+(EIGV(3)*ALPHAI(3))**TWO)
       ENGN=(ELAMEL*((ONE-ALPHA)*(EIGV(1)+EIGV(2)+EIGV(3)))**
     1  TWO)/TWO+ELAMEG*((EIGV(1)*(ONE-ALPHAI(1)))**TWO+(EIGV(2)*
     2  (ONE-ALPHAI(2)))**TWO+(EIGV(3)*(ONE-ALPHAI(3)))**TWO)
C
       SSE=ENGP*GFUN+ENGN
       DENG=ENGP-STATEV(1)
       STATEV(1) = MAX(ENGP,HIST)
       STATEV(4) = HIST
C       
C ======================================================================
C	Time step control
C ======================================================================
C 
        IF (DTSWT.LT.HALF) THEN
          DENGMAXV=GCPAR/2.0/CLPAR*DENGMAX
C --------- If automatic time integration is used
         IF (DENG.GT.DENGMAXV)THEN
C --------- If energy criterion is violated (dE>dE_max)
          IF (DTIME.GT.DTMIN*(ONE+TOLER)) THEN
C --------- If time-step still can be reduced
           PNEWDT=DENGMAXV/DENG/TEN
           DTIMENEXT=PNEWDT*DTIME
           IF (DTIMENEXT.LT.(DTMIN*(ONE+TOLER))) THEN
            PNEWDT=PNEWDT*DTMIN/DTIMENEXT
           ENDIF
          ELSE
C --------- If time-step is already too small
          PNEWDT=ONE
          ENDIF
         ENDIF
        ELSE
         PNEWDT=ONE
        ENDIF
C       
       RETURN
       END
C
C
      SUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX,
     1 PREDEF,DPRED)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
C
      DIMENSION TEMP(2),STATEV(*),PREDEF(*),TIME(2),FLUX(2),
     1 DPRED(*)
       PARAMETER (ONE=1.0,TWO=2.0,THREE=3.0,SIX=6.0,HALF=0.5,
     1 TEN=10.0,ZERO=0.0,FOUR=4.0) 
       DATA NEWTON,TOLER/40,1.D-6/ 
C
       REAL*8 PHASE,HIST,CLPAR,GCPAR
C
       PHASE = TEMP(1)
       GCPAR = STATEV(2) 
       CLPAR = STATEV(3)
       HIST = STATEV(4) 
C
C Defined normalised conductivity
       FLUX(1) = TWO*(ONE-PHASE)*HIST/(GCPAR*CLPAR)-
     1           PHASE/(CLPAR*CLPAR) 
       FLUX(2) = -TWO*HIST/(GCPAR*CLPAR)-ONE/(CLPAR*CLPAR) 
C
      RETURN
      END       
C      
C ==============================================================
C       Stiffness matrix based on asymmetric energy degradation
C ==============================================================
C
       SUBROUTINE HMAT(EPS,CMAT,ENU,EMOD,ANISOSWT,PHASE,PARK)
       INCLUDE 'ABA_PARAM.INC'
       PARAMETER(ZERO=0.D0,ONE=1.D0,MONE=-1.D0,TWO=2.D0,FOUR=4.D0,
     1  TOLER=1.0D-12,SIX=6.D0,FT=50.D0,THREE=3.D0,HALF=0.5D0,
     2  TS=27.D0,TEN=10.D0)
       INTEGER I, J, K, NDIM
       REAL*8 EPS(6),CLMAT(6,6),EIGV(6),ALPHAI(3),
     1  CMAT(6,6),EIGVEC(3,3),EIGVAL(3,3),TEPS(6,6),
     2  TSIGINV(6,6)
C     
       REAL*8 VMAXE, ALPHA, ANISOSWT, ENU, EMOD,EMU
     1 PHASE, CNT
C
C	Initialising CMAT
        DO I=1,6
         EIGV(I)=ZERO
         DO J=1,6
          CMAT(I,J)=ZERO
          CLMAT(I,J)=ZERO
         END DO
        END DO
C
C Obtaining eigenvalues and eigenvectors
        CALL SPRIND(EPS,EIGVAL,EIGVEC,2,3,3)
        DO I=1,3
         EIGV(I)=EIGVAL(I,I)
        END DO
C        
C ************ Starting spectral decomposition ************************        
       IF ((MINVAL(EIGV).GE.-TOLER).OR.(ANISOSWT.LT.TOLER)) THEN
C All eigenvalues are in tension, therefore the stiffness matrix can be
C     degraded without the decomposition
        CMAT(1,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(2,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(3,3)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(1,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(2,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(1,3)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(3,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(2,3)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(3,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(4,4)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
        CMAT(5,5)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
        CMAT(6,6)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
        DO I=1,6
         DO J=1,6
          CMAT(I,J)=CMAT(I,J)*((ONE-PHASE)**TWO+PARK)
         END DO
        END DO
C
       ELSEIF (MAXVAL(EIGV).LT.-TOLER) THEN
C All eigenvalues are in compression, therefore the stiffness
C     matrix does not need to be degraded 
        CMAT(1,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(2,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(3,3)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(ONE-ENU)
        CMAT(1,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(2,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(1,3)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(3,1)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(2,3)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(3,2)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*ENU
        CMAT(4,4)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
        CMAT(5,5)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
        CMAT(6,6)=EMOD/((ONE+ENU)*(ONE-TWO*ENU))*(HALF-ENU)
C
       ELSE
C Degradation switches       
        ALPHA=ZERO
        IF ((EIGV(1)+EIGV(2)+EIGV(3)).GT.TOLER) THEN
         ALPHA=ONE
        ENDIF
        DO K1=1,3
         ALPHAI(K1)=ZERO
         IF (EIGV(K1).GT.TOLER) THEN
          ALPHAI(K1)=ONE
         ENDIF
        END DO
C
C Elementary stiffness matrix in the principal directions
        GAMMA=ENU/(ONE-TWO*ENU)
        CLMAT(1,1)=((ONE-ALPHAI(1)*PHASE)**TWO+PARK)+GAMMA*
     1  ((ONE-ALPHA*PHASE)**TWO+PARK)
        CLMAT(2,2)=((ONE-ALPHAI(2)*PHASE)**TWO+PARK)+GAMMA*
     1  ((ONE-ALPHA*PHASE)**TWO+PARK)
        CLMAT(3,3)=((ONE-ALPHAI(3)*PHASE)**TWO+PARK)+GAMMA*
     1  ((ONE-ALPHA*PHASE)**TWO+PARK)
        CLMAT(1,2)=GAMMA*((ONE-ALPHA*PHASE)**TWO+PARK)
        CLMAT(1,3)=GAMMA*((ONE-ALPHA*PHASE)**TWO+PARK)
        CLMAT(2,3)=GAMMA*((ONE-ALPHA*PHASE)**TWO+PARK)
        CLMAT(2,1)=CLMAT(1,2)
        CLMAT(3,1)=CLMAT(1,3)
        CLMAT(3,2)=CLMAT(2,3)
        DO I=1,3
         DO J=1,3
          CLMAT(I,J)=CLMAT(I,J)*EMOD/(ONE+ENU)
         END DO
        END DO
        EMU=EMOD/(TWO*(ONE+ENU))
        CLMAT(4,4)=(ABS(EIGV(1))*((ONE-ALPHAI(1)*PHASE)**TWO+PARK)+
     1    ABS(EIGV(2))*((ONE-ALPHAI(2)*PHASE)**TWO+PARK))/
     2   (ABS(EIGV(1))+ABS(EIGV(2)))*EMU
        CLMAT(5,5)=(ABS(EIGV(1))*((ONE-ALPHAI(1)*PHASE)**TWO+PARK)+
     1    ABS(EIGV(3))*((ONE-ALPHAI(3)*PHASE)**TWO+PARK))/
     2   (ABS(EIGV(1))+ABS(EIGV(3)))*EMU
        CLMAT(6,6)=(ABS(EIGV(3))*((ONE-ALPHAI(3)*PHASE)**TWO+PARK)+
     1    ABS(EIGV(2))*((ONE-ALPHAI(2)*PHASE)**TWO+PARK))/
     2   (ABS(EIGV(3))+ABS(EIGV(2)))*EMU
C
C   Rotating the stiffness matrix in the original system
C
        CALL ROTMATS(EIGVEC,TEPS,TSIGINV)
        DO I=1,6
            DO J=1,6
                DO K=1,6
                    DO L=1,6
                        CMAT(I,J)=CMAT(I,J)+TSIGINV(I,K)*CLMAT(K,L)*TEPS(L,J);
                    END DO
                END DO
            END DO
        END DO
       ENDIF
C       
       RETURN
       END
C
      SUBROUTINE ROTMATS(EIGVEC,TEPS,TSIGINV)
       INCLUDE 'ABA_PARAM.INC'
       PARAMETER(TWO=2.D0)       
       REAL*8 EIGVEC(3,3), TEPS(6,6), TSIGINV(6,6)
C
C Rotation matrix for the strain tensor in Voigt notation
       TEPS(1,1)=EIGVEC(1,1)**TWO
       TEPS(1,2)=EIGVEC(2,1)**TWO
       TEPS(1,3)=EIGVEC(3,1)**TWO
       TEPS(1,4)=EIGVEC(2,1)*EIGVEC(1,1)
       TEPS(1,5)=EIGVEC(3,1)*EIGVEC(1,1)
       TEPS(1,6)=EIGVEC(2,1)*EIGVEC(3,1)
       TEPS(2,1)=EIGVEC(1,2)**TWO
       TEPS(2,2)=EIGVEC(2,2)**TWO
       TEPS(2,3)=EIGVEC(3,2)**TWO
       TEPS(2,4)=EIGVEC(2,2)*EIGVEC(1,2)
       TEPS(2,5)=EIGVEC(3,2)*EIGVEC(1,2)
       TEPS(2,6)=EIGVEC(2,2)*EIGVEC(3,2)
       TEPS(3,1)=EIGVEC(1,3)**TWO
       TEPS(3,2)=EIGVEC(2,3)**TWO
       TEPS(3,3)=EIGVEC(3,3)**TWO
       TEPS(3,4)=EIGVEC(2,3)*EIGVEC(1,3)
       TEPS(3,5)=EIGVEC(3,3)*EIGVEC(1,3)
       TEPS(3,6)=EIGVEC(2,3)*EIGVEC(3,3)
       TEPS(4,1)=TWO*EIGVEC(1,2)*EIGVEC(1,1)
       TEPS(4,2)=TWO*EIGVEC(2,2)*EIGVEC(2,1)
       TEPS(4,3)=TWO*EIGVEC(3,2)*EIGVEC(3,1)
       TEPS(4,4)=EIGVEC(2,2)*EIGVEC(1,1)+EIGVEC(1,2)*EIGVEC(2,1)
       TEPS(4,5)=EIGVEC(3,1)*EIGVEC(1,2)+EIGVEC(1,1)*EIGVEC(3,2)
       TEPS(4,6)=EIGVEC(2,2)*EIGVEC(3,1)+EIGVEC(3,2)*EIGVEC(2,1)
       TEPS(5,1)=TWO*EIGVEC(1,3)*EIGVEC(1,1)
       TEPS(5,2)=TWO*EIGVEC(2,3)*EIGVEC(2,1)
       TEPS(5,3)=TWO*EIGVEC(3,3)*EIGVEC(3,1)
       TEPS(5,4)=EIGVEC(1,1)*EIGVEC(2,3)+EIGVEC(2,1)*EIGVEC(1,3)
       TEPS(5,5)=EIGVEC(3,1)*EIGVEC(1,3)+EIGVEC(1,1)*EIGVEC(3,3)
       TEPS(5,6)=EIGVEC(2,3)*EIGVEC(3,1)+EIGVEC(3,3)*EIGVEC(2,1)
       TEPS(6,1)=TWO*EIGVEC(1,3)*EIGVEC(1,2)
       TEPS(6,2)=TWO*EIGVEC(2,3)*EIGVEC(2,2)
       TEPS(6,3)=TWO*EIGVEC(3,3)*EIGVEC(3,2)
       TEPS(6,4)=EIGVEC(2,2)*EIGVEC(1,3)+EIGVEC(2,3)*EIGVEC(1,2)
       TEPS(6,5)=EIGVEC(1,3)*EIGVEC(3,2)+EIGVEC(1,2)*EIGVEC(3,3)
       TEPS(6,6)=EIGVEC(2,3)*EIGVEC(3,2)+EIGVEC(3,3)*EIGVEC(2,2)
C
C Inverse of the rotation matrix for the stress tensor in Voigt notation
       TSIGINV(1,1)=EIGVEC(1,1)**TWO
       TSIGINV(1,2)=EIGVEC(1,2)**TWO
       TSIGINV(1,3)=EIGVEC(1,3)**TWO
       TSIGINV(1,4)=TWO*EIGVEC(1,2)*EIGVEC(1,1)
       TSIGINV(1,5)=TWO*EIGVEC(1,3)*EIGVEC(1,1)
       TSIGINV(1,6)=TWO*EIGVEC(1,2)*EIGVEC(1,3)
       TSIGINV(2,1)=EIGVEC(2,1)**TWO
       TSIGINV(2,2)=EIGVEC(2,2)**TWO
       TSIGINV(2,3)=EIGVEC(2,3)**TWO
       TSIGINV(2,4)=TWO*EIGVEC(2,2)*EIGVEC(2,1)
       TSIGINV(2,5)=TWO*EIGVEC(2,3)*EIGVEC(2,1)
       TSIGINV(2,6)=TWO*EIGVEC(2,2)*EIGVEC(2,3)
       TSIGINV(3,1)=EIGVEC(3,1)**TWO
       TSIGINV(3,2)=EIGVEC(3,2)**TWO
       TSIGINV(3,3)=EIGVEC(3,3)**TWO
       TSIGINV(3,4)=TWO*EIGVEC(3,2)*EIGVEC(3,1)
       TSIGINV(3,5)=TWO*EIGVEC(3,3)*EIGVEC(3,1)
       TSIGINV(3,6)=TWO*EIGVEC(3,2)*EIGVEC(3,3)
       TSIGINV(4,1)=EIGVEC(2,1)*EIGVEC(1,1)
       TSIGINV(4,2)=EIGVEC(2,2)*EIGVEC(1,2)
       TSIGINV(4,3)=EIGVEC(2,3)*EIGVEC(1,3)
       TSIGINV(4,4)=EIGVEC(2,2)*EIGVEC(1,1)+EIGVEC(2,1)*EIGVEC(1,2)
       TSIGINV(4,5)=EIGVEC(2,1)*EIGVEC(1,3)+EIGVEC(2,3)*EIGVEC(1,1)
       TSIGINV(4,6)=EIGVEC(2,3)*EIGVEC(1,2)+EIGVEC(2,2)*EIGVEC(1,3)
       TSIGINV(5,1)=EIGVEC(3,1)*EIGVEC(1,1)
       TSIGINV(5,2)=EIGVEC(3,2)*EIGVEC(1,2)
       TSIGINV(5,3)=EIGVEC(3,3)*EIGVEC(1,3)
       TSIGINV(5,4)=EIGVEC(3,2)*EIGVEC(1,1)+EIGVEC(3,1)*EIGVEC(1,2)
       TSIGINV(5,5)=EIGVEC(3,1)*EIGVEC(1,3)+EIGVEC(3,3)*EIGVEC(1,1)
       TSIGINV(5,6)=EIGVEC(3,3)*EIGVEC(1,2)+EIGVEC(3,2)*EIGVEC(1,3)
       TSIGINV(6,1)=EIGVEC(3,1)*EIGVEC(2,1)
       TSIGINV(6,2)=EIGVEC(3,2)*EIGVEC(2,2)
       TSIGINV(6,3)=EIGVEC(3,3)*EIGVEC(2,3)
       TSIGINV(6,4)=EIGVEC(2,2)*EIGVEC(3,1)+EIGVEC(2,1)*EIGVEC(3,2)
       TSIGINV(6,5)=EIGVEC(3,1)*EIGVEC(2,3)+EIGVEC(3,3)*EIGVEC(2,1)
       TSIGINV(6,6)=EIGVEC(2,3)*EIGVEC(3,2)+EIGVEC(2,2)*EIGVEC(3,3)       
       RETURN
       END