C     Authors               : Jean-Baptiste LEBLOND and Leo MORIN
C     Contact (J.B. Leblond): jbl@lmm.jussieu.fr
C 
C     Copyright 2015 Sorbonne Universites, Universite Pierre et Marie Curie
C     (Paris 6), Institut Jean Le Rond d'Alembert, Tour 55-65, 4 place Jussieu,
C     75252 Paris Cedex 05, France
C
C     This program is free software: you can redistribute it and/or modify
C     it under the terms of the GNU Lesser General Public License as 
C     published by the Free Software Foundation, either version 3 of the
C     License, or (at your option) any later version.
C
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C     GNU Lesser General Public License for more details.
C
C     You should have received a copy of the GNU Lesser General Public License
C     along with this program. If not, see <http://www.gnu.org/licenses/>.
C	
C 
C 
C
C	 
C 	  
	  SUBROUTINE 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,KSTEP,KINC)
C
C-----------------------------------------------------------------------
C
C     Purpose : Projection of the elastic predictor onto Madou-Leblond's 
C               yield surface - interface with subroutine PROJML
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I/O  REAL      STRESS    Input : stress tensor at time T
C                              Output: stress tensor at time T+DT
C                              (Voigt's notations)
C     I/O  REAL      STATEV    Input : Vector of internal variables at 
C                                      time T
C                              Output: Vector of internal variables at 
C                                      time T+DT
C     O    REAL      DDSDDE    Tangent matrix (Voigt's notations)
C     I/O  REAL      SSE       Unused
C     I/O  REAL      SPD       Unused
C     I/O  REAL      SCD       Unused
C     O    REAL      RPL       Unused
C     O    REAL      DDSDDT    Unused
C     O    REAL      DRPLDE    Unused
C     O    REAL      DRPLDT    Unused
C     I    REAL      STRAN     Unused
C     I    REAL      DSTRAN    Increment of strain between times T 
C                              and T+DT (Voigt's notations)
C     I    REAL      TIME      Unused
C     I    REAL      DTIME     Unused
C     I    REAL      TEMP      Unused
C     I    REAL      DTEMP     Unused
C     I    REAL      PREDEF    Unused
C     I    REAL      DPRED     Unused
C     I    CHARACTER CMNAME    Name of material (MADLEB)
C     I    INTEGER   NDI       Number of diagonal stress or strain   
C                              components (3)
C     I    INTEGER   NSHR      Number of off-diagonal stress or strain    
C                              components (3)
C     I    INTEGER   NTENS     Total number of stress or strain components (6)  
C     I    INTEGER   NSTATV    Number of internal variables (9) 
C     I    REAL      PROPS     Vector containing material constants [elastic    
C                              constants, data pertaining to voids, 
C                              (equivalent strain)-(equivalent stress) curve
C                              in discrete form]
C     I    INTEGER   NPROPS    Number of material constants [(2 + 16     
C                              + 2*(nb. of points on the stress-strain curve)]
C     I    REAL      COORDS    Unused
C     I    REAL      DROT      Rotation of the material basis between times
C                              T and T+DT
C     I    REAL      CELENT    Unused
C     I    REAL      DFGRD0    Unused
C     I    REAL      DFGRD1    Unused
C     I    INTEGER   NOEL      Element number
C     I    INTEGER   NPT       Integration point number
C     I    INTEGER   LAYER     Unused
C     I    INTEGER   KSPR      Unused
C     I    INTEGER   KSTEP     Unused
C     I    INTEGER   KINC      Unused
C
C  ----------------------------------------------------------------------
C
      INCLUDE 'ABA_PARAM.inc'  
C
C     Arguments
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     .          DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     .          STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     .          PROPS(NPROPS),COORDS(3),DROT(3,3),
     .          DFGRD0(3,3),DFGRD1(3,3)
      CHARACTER*80 CMNAME
C
C     Local variables
C     Arguments of subroutine STDB_ABQERR
      DIMENSION REALV(1),INTV(4)
      CHARACTER*8 CHARV(1)
C     Other variables
      REAL*8    ELDAT(2),PLDAT(201),VOIDAT(16),SIG0(6),INTPAR0(8),
     .          GRTRSF(3,3),LAMBDA,MU,SIG1(6),TRDDEF,INTPAR1(8),QFG1G,
     .          TGTMAT(6,6)
      INTEGER   NPLINFO,NPTCRV,I,ICONV,ITGTMT,J,ERRINFO(3)
C
      INCLUDE 'STRINGS.inc'
C	  
C     Preliminaries
C     *************
C
C     If there is an error the programme is stopped  
      LOP = - 3
C     INTV(1) contains the element number and INTV(2) the integration
C     point number
      INTV(1) = NOEL
      INTV(2) = NPT
C 
      IF (CMNAME.NE.'MADLEB') 
     .                 CALL STDB_ABQERR(LOP,STRINGA1,INTV,REALV,CHARV)      
      IF (NDI.NE.3)    CALL STDB_ABQERR(LOP,STRINGA2,INTV,REALV,CHARV)
      IF (NSHR.NE.3)   CALL STDB_ABQERR(LOP,STRINGA3,INTV,REALV,CHARV)
      IF (NTENS.NE.6)  CALL STDB_ABQERR(LOP,STRINGA4,INTV,REALV,CHARV)
      IF (NSTATV.NE.9) CALL STDB_ABQERR(LOP,STRINGA5,INTV,REALV,CHARV)
C     There must be 2 elastic constants, 16 constants pertaining to 
C     voids and at least 4 plastic constants (2 points on the 
C     stress-strain curve) 
      IF (NPROPS.LT.22) CALL STDB_ABQERR(LOP,STRINGA6,INTV,REALV,CHARV)
C 
C     Preparation of inputs of subroutine PROJML
C     ******************************************
C 
C     Elastic data: Young's modulus, Poisson's ratio
C     ----------------------------------------------
      ELDAT(1) = PROPS(1)
      ELDAT(2) = PROPS(2)
C
C     Plastic data: stress-strain curve 
C     ---------------------------------
C     Number of plastic constants
      NPLINFO = NPROPS - 18 
C     Number of points on the stress-strain curve
      NPTCRV  = NPLINFO/2 
      IF (NPTCRV.GT.100) 
     .  CALL STDB_ABQERR(LOP,STRINGA7,INTV,REALV,CHARV)
      PLDAT(1) = DFLOAT(NPTCRV)
C     Points on the stress-strain curve
      DO I=1,NPLINFO
        PLDAT(I+1) = PROPS(I+18)
      ENDDO
C 
C     Data pertaining to voids
C     ------------------------
      DO I=1,16
        VOIDAT(I) = PROPS(I+2)	  
      ENDDO
C 
C     Stress tensor at time T
C     -----------------------
C     Rem. SIG0 as defined below does not really contain the stress 
C          tensor at time T, as required by subroutine PROJML, since 
C          the rotation DROT has been applied to this tensor prior 
C          to the call to subroutine UMAT. This is unimportant since
C          SIG0 is used in PROJML only to test whether the stress
C          tensor at time T is zero, which is independent of whether
C          it has been multiplied by DROT or not.  
      DO I=1,NTENS
        SIG0(I) = STRESS(I)
      ENDDO 
C 
C     Internal parameters at time T
C     -----------------------------
      DO I=1,NSTATV-1
        INTPAR0(I) = STATEV(I)
      ENDDO
C
C     Gradient of transformation between times T and T+DT
C     ---------------------------------------------------
C     Rem. The formulae below account for both the deformation and 
C          rotation parts of the transformation gradient, for 
C          consistency with the arguments of subroutine PROJML. 
C          They could however be simplified by omitting the 
C          deformation part since PROJML uses only the 
C          antisymmetric part of the transformation gradient.
      GRTRSF(1,1) = DSTRAN(1) + DROT(1,1) 
      GRTRSF(2,2) = DSTRAN(2) + DROT(2,2) 
      GRTRSF(3,3) = DSTRAN(3) + DROT(3,3)
      GRTRSF(1,2) = 0.5D0*DSTRAN(4) + DROT(1,2) 
      GRTRSF(2,1) = 0.5D0*DSTRAN(4) + DROT(2,1) 
      GRTRSF(1,3) = 0.5D0*DSTRAN(5) + DROT(1,3) 
      GRTRSF(3,1) = 0.5D0*DSTRAN(5) + DROT(3,1) 
      GRTRSF(2,3) = 0.5D0*DSTRAN(6) + DROT(2,3) 
      GRTRSF(3,2) = 0.5D0*DSTRAN(6) + DROT(3,2) 
C
C     Codes for convergence of the global elastoplastic iterations
C     and request of the tangent matrix     
C     ------------------------------------------------------------
C     These codes being not provided in the arguments of UMAT 
C     are set to unity 
      ICONV  = 1
      ITGTMT = 1
C 
C     Elastic stress predictor
C     ------------------------
C     The terms pertaining to the objective derivative of the 
C     stress tensor have already been accounted for since the 
C     rotation DROT has been applied to STRESS.
      LAMBDA = ELDAT(1)*ELDAT(2)/(1.D0+ELDAT(2))/(1.D0-2.D0*ELDAT(2))
      MU     = 0.5D0*ELDAT(1)/(1.D0+ELDAT(2))
      DO I=1,NTENS
        SIG1(I) = STRESS(I)
      ENDDO
      TRDDEF = DSTRAN(1) + DSTRAN(2) + DSTRAN(3)
      DO I=1,NDI
        SIG1(I) = SIG1(I) + LAMBDA*TRDDEF + 2.D0*MU*DSTRAN(I)
      ENDDO
      DO I=1,NSHR
        J = I+NDI
        SIG1(J) = SIG1(J) + MU*DSTRAN(J)
      ENDDO
C 
C     Call to subroutine PROJML (projection onto Madou-Leblond's 
C     yield locus)
C     **********************************************************
C 
      CALL PROJML(ELDAT,PLDAT,VOIDAT,SIG0,INTPAR0,GRTRSF,
     .            ICONV,ITGTMT,SIG1,INTPAR1,QFG1G,TGTMAT,ERRINFO)
C 
C     Error messages
C     **************
C
C     INTV(3) contains the number of the current iteration on the 
C     strain hardening parameter and INTV(4) the number of the
C     current iteration on the plastic multiplier 
      INTV(3) = ERRINFO(2)
	  INTV(4) = ERRINFO(3)
C
      IF (ERRINFO(1).EQ.1) THEN
       CALL STDB_ABQERR(LOP,STRINGB1,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.2) THEN
       CALL STDB_ABQERR(LOP,STRINGB2,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.3) THEN
       CALL STDB_ABQERR(LOP,STRINGB3,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.4) THEN
       CALL STDB_ABQERR(LOP,STRINGB4,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.5) THEN
       CALL STDB_ABQERR(LOP,STRINGB5,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.6) THEN
       CALL STDB_ABQERR(LOP,STRINGB6,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.7) THEN
       CALL STDB_ABQERR(LOP,STRINGB7,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.8) THEN
       CALL STDB_ABQERR(LOP,STRINGB8,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.9) THEN
       CALL STDB_ABQERR(LOP,STRINGB9,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.10) THEN
       CALL STDB_ABQERR(LOP,STRINGB10,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.11) THEN
       CALL STDB_ABQERR(LOP,STRINGB11,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.12) THEN
       CALL STDB_ABQERR(LOP,STRINGB12,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.13) THEN
       CALL STDB_ABQERR(LOP,STRINGB13,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.14) THEN
       CALL STDB_ABQERR(LOP,STRINGB14,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.15) THEN
       CALL STDB_ABQERR(LOP,STRINGB15,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.16) THEN
       CALL STDB_ABQERR(LOP,STRINGB16,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.17) THEN
       CALL STDB_ABQERR(LOP,STRINGB17,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.18) THEN
       CALL STDB_ABQERR(LOP,STRINGB18,INTV,REALV,CHARV)
	  ELSEIF (ERRINFO(1).EQ.19) THEN
       CALL STDB_ABQERR(LOP,STRINGB19,INTV,REALV,CHARV)
      ENDIF
C 
C     Treatment of outputs of subroutine PROJML
C     *****************************************
C 
C     Stress tensor at time T+DT
C     --------------------------
      DO I=1,NTENS
        STRESS(I) = SIG1(I) 
      ENDDO 
C 
C     Internal parameters at time T+DT 
C     --------------------------------
C     STATEV(1:NSTATV-1) contains the true internal variables of the material
C     STATEV(NSTATV) contains the "damage parameter" QFG1G of the material 
C     which is not a true internal variable but may be used for postprocessing
      DO I=1,NSTATV-1
        STATEV(I) = INTPAR1(I) 
      ENDDO
      STATEV(NSTATV) = QFG1G
C 
C     Tangent matrix
C     --------------
      DO I=1,NTENS
        DO J=1,NTENS
          DDSDDE(I,J) = TGTMAT(I,J)
        ENDDO 
      ENDDO		
C
      RETURN
      END	
C
C
C
C
C
C
      SUBROUTINE CHGBAS(MATG,MATL,MATR,ICODE)
C
C-----------------------------------------------------------------------
C
C     Purpose : Transformation of an arbitrary (not necessarily symmetric)
C               3*3 matrix from the global basis to the local one or vice-versa
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I/O  REAL*8  MATG    3X3 matrix in the global basis
C     I/O  REAL*8  MATL    3X3 matrix in the local basis
C     I    REAL*8  MATR    Transformation matrix (rotation) from the global basis 
C                          to the local one (the columns of this matrix provide 
C                          the vectors of the local basis, expressed in the global 
C                          one)
C     I    INTEGER ICODE   Code:
C                          =1 : transformation from the global basis
C                                to the local one
C                          =2 : transformation from the local basis
C                                to the global one
C
C  ----------------------------------------------------------------------
C
      IMPLICIT NONE
C
C     Arguments
      REAL*8 MATG(3,3),MATL(3,3),MATR(3,3)
      INTEGER ICODE
C
C     Local variables
      REAL*8 A(3,3),B(3,3),ROT(3,3)  
      INTEGER I,J,K,L
C
      DO I=1,3
        DO J=1,3
          IF (ICODE.EQ.1) THEN		  
            A(I,J)   = MATG(I,J)
            ROT(I,J) = MATR(I,J)
          ELSEIF (ICODE.EQ.2) THEN
            A(I,J)   = MATL(I,J)
            ROT(I,J) = MATR(J,I)	
          ENDIF			
          B(I,J) = 0.D0
        ENDDO		
      ENDDO		  
C		  
      DO I=1,3
        DO J=1,3
          DO K=1,3
		    DO L=1,3
		      B(I,J) = B(I,J) + ROT(K,I)*A(K,L)*ROT(L,J)
            ENDDO		
          ENDDO
          IF(ICODE.EQ.1) THEN
            MATL(I,J) = B(I,J)
          ELSEIF (ICODE.EQ.2) THEN
            MATG(I,J) = B(I,J)
          ENDIF
        ENDDO		
      ENDDO
C	  
      RETURN
      END	
C
C
C
C
C
C
      SUBROUTINE CONTML(AMINI,BMINI,CMINI,F,P,IERR)
C
C-----------------------------------------------------------------------
C
C     Purpose  : Account for the possibility of minimum values of the semi-axes 
C                A,B,C of the void in Madou-Leblond's model. These minimum values 
C                may be due to either contact of the void's boundary with an 
C                enclosed inclusion, or autocontact between opposite parts of 
C                the void's boundary (in this case the minimum values must be 
C                slightly positive, zero values are not allowed).
C     Remark:    The modeling of contact is approximate, in that the evolution laws
C                of the porosity and the shape of the voids account for contact,
C                but neither the yield criterion nor the plastic flow rule do.
C     Reference: In construction
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I    REAL*8   AMINI     Minimum value of the major semi-axis of the void
C     I    REAL*8   BMINI     Minimum value of the intermediate semi-axis of the 
C                             void
C     I    REAL*8   CMINI     Minimum value of the minor semi-axis of the void
C     I/O  REAL*8   F         Porosity at time T+DT
C                             Input : porosity prior to correction due to minimum
C                                     values of the void axes
C                             Output: porosity after correction due to minimum
C                                     values of the void axes
C     I/O  REAL*8   P         Quadratic form characterizing the orientation and 
C                             semi-axes of the void at time T+DT (global basis,
C                             tensorial form)
C                             Input : quadratic form prior to correction due to 
C                                     minimum values of the void axes
C                             Output: quadratic form after correction due to 
C                                     minimum values of the void axes
C     O    INTEGER  IERR      Error code:
C                             =1, error in the diagonalization of P
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C
C     Arguments
      REAL*8  AMINI,BMINI,CMINI,F,P(3,3)
      INTEGER IERR
C
C     Local variables
      REAL*8  EIGNVL(3),WORK(8),A,B,C,MATR(3,3)
      INTEGER LWORK,INFO,I,J
C
C     Initialization of error code
C     ****************************
C
      IERR = 0
C
C     Diagonalization of the quadratic form P
C     ***************************************
C     Input : P(3,3)       symmetric matrix to be diagonalized
C     Output: EIGNVL(3)    eigenvalues in ascending order
C             P(3,3)       matrix whose columns contain the corresponding 
C                          orthonormal eigenvectors
      LWORK = 8
      CALL DSYEV('V','L',3,P,3,EIGNVL,WORK,LWORK,INFO)
      IF ( (EIGNVL(1).LE.0.D0) .OR. (EIGNVL(2).LE.0.D0) 
     .     .OR. (EIGNVL(3).LE.0.D0) .OR. (INFO.NE.0) ) THEN
        IERR = 1
        RETURN
      ENDIF
C 
C     Semi-axes of the void at time T+DT
      A = 1.D0/SQRT(EIGNVL(1))
      B = 1.D0/SQRT(EIGNVL(2))
      C = 1.D0/SQRT(EIGNVL(3))
C
C     Tests and possibly corrections on the semi-axes and the porosity 
C     at time T+DT
C     ****************************************************************
C
      IF (A.LT.AMINI) THEN
        F = F*AMINI/A
        A = AMINI
      ENDIF
C 
      IF (B.LT.BMINI) THEN
        F = F*BMINI/B
        B = BMINI
      ENDIF
C 
      IF (C.LT.CMINI) THEN
        F = F*CMINI/C
        C = CMINI
      ENDIF
C
C     Corrected quadratic form P at time T+DT     
C     ***************************************
C
      CALL TRNSFR(P,MATR,9)
      DO I=1,3
        DO J=1,3
          P(I,J) = MATR(I,1)*MATR(J,1)/A**2
     .             + MATR(I,2)*MATR(J,2)/B**2
     .             + MATR(I,3)*MATR(J,3)/C**2
        ENDDO
      ENDDO	  
C
      RETURN
      END 
C
C
C
C
C
C
      SUBROUTINE CRITML(F,A,B,C,
     .                  K,Q,G,KAPPA,HX,HY,HZ,
     .                  IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC,
     .                  IERR)
C
C-----------------------------------------------------------------------
C
C     Purpose  : Calculate the coefficients of Madou and Leblond's 
C                criterion for plastic porous materials containing arbitrary 
C                ellipsoidal voids
C     Reference: K. Madou, J.B. Leblond, A Gurson-type criterion for porous 
C                ductile solids containing arbitrary ellipsoidal voids -
C                II: Determination of yield criterion parameters, 
C                Journal of the Mechanics and Physics of Solids, vol. 60, 
C                pp. 1037-1058 (2012) [especially Section 7, pp. 1051-1052, 
C                Appendix B, pp. 1056-1057, and Appendix C, p. 1057]
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I    REAL*8   F           Porosity
C     I    REAL*8   A           Major semi-axis of the void
C     I    REAL*8   B           Intermediate semi-axis of the void
C     I    REAL*8   C           Minor semi-axis of the void
C
C     O    REAL*8   K           Parameter related to the void shape:
C                               =0, spherical or prolate spheroidal void
C                               =1, oblate spheroidal void 
C                               0<K<1, general ellipsoidal void
C     O    REAL*8   Q           Quadratic form of the stress components 
C                               (first term of the criterion), stored in the 
C                               order Q11,Q12,Q13,Q22,Q23,Q33,Q44,Q55,Q66 
C                               (Voigt's notation)
C     O    REAL*8   G           "Second porosity"
C     O    REAL*8   KAPPA       Coefficient of SIGMA_H in the hyperbolic cosine
C                               (second term of the criterion)
C     O    REAL*8   HX,HY,HZ    Coefficients defining SIGMA_H 
C                               (second term of the criterion)
C
C     O    REAL*8   IA,IB,IC    Integrals IA,IB,IC
C     O    REAL*8   IAA,IBB,ICC Integrals IAA,IBB,ICC
C     O    REAL*8   IAB,IAC,IBC Integrals IAB,IAC,IBC
C                               (All of these integrals are used again in 
C                               subroutine EVOLML)
C
C     O    INTEGER  IERR        Error code:
C                               =1, error on the semi-axes of the void
C                               =2, the Newton iterations on LAMBDA do not 
C                                   converge
C                               =3, the arguments of subroutine FUNCRD for the 
C                                   calculation of IA are invalid
C                               =4, the iterations needed to calculate IA in FUNCRD 
C                                   do not converge
C                               =5, the arguments of subroutine FUNCRD for the 
C                                   calculation of IB are invalid
C                               =6, the iterations needed to calculate IB in FUNCRD 
C                                   do not converge
C                               =7, non-invertible matrix T (zero determinant) 
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C
C     Arguments
      REAL*8  F,A,B,C
      REAL*8  K,Q(9),G,KAPPA,HX,HY,HZ,
     .        IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC
      INTEGER IERR
C
C     Local variables
      REAL*8  A2,B2,C2,F2,COEF0,COEF1,COEF2,SMALL,EPS,LAMBDA,
     .        ERROR,SLOPE,AEXT,BEXT,CEXT,ABAR,BBAR,K2,
     .        G1,GF,EXZ,EEXZ,FBAR,ALPHA,BETA,EEXZ2,EEXZ4,
     .        HYPROL,HYOBL,HXPROL,HXOBL,PI,
     .        T(3,3),TXYXY,TYZYZ,TXZXZ,TINV(3,3)
      INTEGER NITLMB,ITE,IERR1
C
C     Initialization of error code
C     ****************************
C
      IERR = 0
C
C     Tests on the semi-axes of the void
C     **********************************
C
      IF ( (A.LT.B) .OR. (B.LT.C) ) THEN
        IERR = 1
        RETURN
      ENDIF
C
C     Geometrical parameters
C     **********************
C
C     Calculation of LAMBDA: Newton method
C
      A2 = A*A
      B2 = B*B
      C2 = C*C
      F2 = F*F
      COEF0 = A2*B2*C2*(1.D0 - 1.D0/F2) 
      COEF1 = A2*B2 + B2*C2 + A2*C2 
      COEF2 = A2 + B2 + C2 
C
      NITLMB = 20 
C     The accuracy requested EPS is a small fraction of COEF0
      SMALL = 1.D-9
      EPS   = A2*B2*C2/F2*SMALL
C     Initialization of LAMBDA assuming equal semi-axes
      LAMBDA = (A2*B2*C2/F2)**(1.D0/3.D0) - C2 
C
      DO ITE=1,NITLMB
        ERROR = LAMBDA**3 + COEF2*LAMBDA**2
     .		    + COEF1*LAMBDA + COEF0 
        IF (ABS(ERROR).LE.EPS) GOTO 10
        SLOPE  = 3.D0*LAMBDA**2 + 2.D0*COEF2*LAMBDA + COEF1 
        LAMBDA = LAMBDA - ERROR/SLOPE	  
      ENDDO
C
C     No convergence on LAMBDA
      IERR = 2
      RETURN
C
C     Convergence on LAMBDA
   10 CONTINUE
C
C     Exterior semi-axes
C
      AEXT = SQRT(A2 + LAMBDA)	 
	  BEXT = SQRT(B2 + LAMBDA) 
	  CEXT = SQRT(C2 + LAMBDA) 
C
C     Semi-axes ABAR, BBAR (CBAR=0) of the completely flat confocal ellipsoid
C     Dimensionless parameter K
C
      ABAR = SQRT(A2 - C2)
	  BBAR = SQRT(B2 - C2)

      IF(A.EQ.C) THEN
C       Special case: spherical void - conventionally K = 0
        K = 0.D0
      ELSE
C       Generic case
        K = BBAR/ABAR
      ENDIF	  
      K2 = K*K
C
C     Second porosity G and related parameters G1, GF
C
      G  = ABAR*BBAR*BBAR/AEXT/BEXT/CEXT 
      G1 = G/(1.D0+G)
	  GF = G/(F+G)
C
C     Eccentricities of the inner and outer ellipsoids in the XZ plane
C
      EXZ  = ABAR/A
      EEXZ = ABAR/AEXT
C
C     Coefficient KAPPA
C     *****************
C
C     Mean value FBAR of the function F (prolate/generic cases)
      IF (B.EQ.C) THEN
C       Prolate case
        FBAR = 1.D0+(1.D0-SQRT(3.D0)/2.D0)/LOG(F)*LOG((11.D0
     .	       +5.D0*EXZ**3/(1.D0-EXZ**2))/(11.D0
     .         +5.D0*F*EXZ**3/(1.D0-EXZ**2)))
      ELSE
C       Other cases      
        FBAR = 1.D0+1.D0/LOG(GF/G1)*( -(1.D0-K)*(1.D0
     .         -SQRT(3.D0)/2.D0)*LOG((11.D0*K2+5.D0*GF)/(11.D0*K2
     .         +5.D0*G1))+3.D0/5.D0*(1.D0-K)**2*LOG((8.D0
     .         -5.D0*G1)/(8.D0-5.D0*GF))
     .         +13.D0/10.D0*K*(GF-G1)-3.D0/10.D0*K*(GF**5-G1**5) )
      ENDIF
      KAPPA = 3.D0/2.D0/FBAR
C
C     Coefficients HX,HY,HZ
C     *********************
C
C     Ancillary coefficients: ALPHA,BETA,HXPROL,HYPROL,HXOBL,HYOBL
C
      ALPHA = 4.D0*K2/(1.D0+9.D0*K2)
      BETA  = 3.D0*K2/(1.D0+30.D0*K2)
C
      EEXZ2  = EEXZ*EEXZ
      EEXZ4  = EEXZ2*EEXZ2	  
      HYPROL = (1.D0+EEXZ2-EEXZ4/2.D0)/3.D0
      HYOBL  = (2.D0-7.D0*EEXZ2+5.D0*EEXZ4)/
     .         (2.D0-7.D0*EEXZ2+10.D0*EEXZ4)/3.D0
      HXPROL = 1.D0-2.D0*HYPROL
      HXOBL  = HYOBL
C
C     Calculation of coefficients HX,HY,HZ (prolate/generic cases)
C
      HX = (1.D0-K2)*HXPROL + K2*HXOBL
      IF (B.EQ.C) THEN
C       Prolate case
        HY = HYPROL
      ELSE
C       Other cases      
        HY = (1.D0-K)*HYPROL + K*HYOBL +(1.D0-K)/2.D0*(ALPHA**2
     .       +BETA**2)/ALPHA*EEXZ**(3.D0/2.D0)*(1.D0-ALPHA-EEXZ)
     .       /((1.D0-ALPHA-EEXZ)**2+BETA**2)
      ENDIF
      HZ = 1.D0-HX-HY
C
C     Elliptic integrals IA, IB, etc.
C     *******************************
C
      CALL FUNCRD(B2,C2,A2,IA,IERR1)
      IF (IERR1.NE.0) THEN
        IERR = IERR1 + 2
        RETURN
      ENDIF
      CALL FUNCRD(C2,A2,B2,IB,IERR1)
      IF (IERR1.NE.0) THEN
        IERR = IERR1 + 4
        RETURN
      ENDIF
C
C     Integrals IA,IB,IC
      PI    = 4.D0*ATAN(1.D0)
      COEF0 = 4.D0*PI/3.D0*A*B*C
      IA    = IA*COEF0
      IB    = IB*COEF0
      IC    = 4.D0*PI - IA - IB
C
C     Integrals IAB,IAC,IBC
      IF ( (A.EQ.C) ) THEN
C       Spherical case
        IAB = 4.D0*PI/15.D0/A2
        IAC = IAB
        IBC = IAB	  		 
      ELSEIF( (A.EQ.B) .AND. (B.NE.C) ) THEN
C       Oblate case
        IAC = (IC - IA)/3.D0/(A2 - C2)
        IBC = (IC - IB)/3.D0/(B2 - C2)
        IAB = PI/3.D0/A2 - IAC/4.D0 
      ELSEIF( (B.EQ.C) .AND. (A.NE.B) ) THEN
C       Prolate case
        IAB = (IB - IA)/3.D0/(A2 - B2)
        IAC = (IC - IA)/3.D0/(A2 - C2)
        IBC = PI/3.D0/B2 - IAB/4.D0
      ELSE	   
C       Generic case
        IAB = (IB - IA)/3.D0/(A2 - B2)
        IAC = (IC - IA)/3.D0/(A2 - C2)
        IBC = (IC - IB)/3.D0/(B2 - C2)	
      ENDIF
C
C     Integrals IAA,IBB,ICC 
      COEF0 = 4.D0*PI/3.D0
      IAA   = COEF0/A2 - IAB - IAC
      IBB   = COEF0/B2 - IAB - IBC
      ICC   = COEF0/C2 - IAC - IBC
C	   
C     Quadratic form Q of the stress components Q(1-9)
C     ************************************************
C
C     Components T(1,1)=TXXXX, T(1,2)=TXXYY, T(1,3)=TXXZZ, 
C                T(2,2)=TYYYY, T(2,3)=TYYZZ, T(3,3)=TZZZZ
C                TXYXY, TXZXZ, TYZYZ 
C     of the tensor T
      COEF0  = 8.D0*PI
      T(1,1) = 2.D0 - (3.D0*A2*IAA + IA)*2.D0/COEF0
      T(2,2) = 2.D0 - (3.D0*B2*IBB + IB)*2.D0/COEF0
      T(3,3) = 2.D0 - (3.D0*C2*ICC + IC)*2.D0/COEF0	  
      T(1,2) = 1.D0 - (IA + IB + 3.D0*(A2+B2)*IAB)/COEF0
      T(2,1) = T(1,2)	  
      T(1,3) = 1.D0 - (IA + IC + 3.D0*(A2+C2)*IAC)/COEF0	
      T(3,1) = T(1,3)	  
      T(2,3) = 1.D0 - (IB + IC + 3.D0*(B2+C2)*IBC)/COEF0
      T(3,2) = T(2,3)	  
      TXYXY  = 1.D0/2.D0 - 3.D0*(A2 + B2)*IAB/COEF0
      TXZXZ  = 1.D0/2.D0 - 3.D0*(A2 + C2)*IAC/COEF0
      TYZYZ  = 1.D0/2.D0 - 3.D0*(B2 + C2)*IBC/COEF0 
C
C     Inverse of T
      CALL INVMAT(T,TINV,IERR1)
      IF(IERR1.EQ.1) THEN
        IERR = 7
        RETURN
      ENDIF
C
C     Willis's quadratic form (nonlinear bound for ellipsoidal voids)
C     Q(1-6) ("diagonal" terms)
      COEF0 = 1.D0-F
      COEF1 = 3.D0*F/2.D0
      Q(1) = COEF0 + COEF1*TINV(1,1)
      Q(2) = - COEF0/2.D0 + COEF1*TINV(1,2)    
      Q(3) = - COEF0/2.D0 + COEF1*TINV(1,3)	 
      Q(4) = COEF0 + COEF1*TINV(2,2) 
      Q(5) = - COEF0/2.D0 + COEF1*TINV(2,3)
      Q(6) = COEF0 + COEF1*TINV(3,3)
C     Q(7-9) ("off-diagonal terms)
      Q(7) = 3.D0*COEF0 + COEF1/TXYXY
      Q(8) = 3.D0*COEF0 + COEF1/TXZXZ
      Q(9) = 3.D0*COEF0 + COEF1/TYZYZ
C
C     Madou-Leblond's quadratic form: corrections to Willis's quadratic form
      COEF0 = (1.D0+G)*(F+G)*KAPPA**2
      Q(1)  = Q(1) - COEF0*HX*HX
      Q(2)  = Q(2) - COEF0*HX*HY
      Q(3)  = Q(3) - COEF0*HX*HZ	  
      Q(4)  = Q(4) - COEF0*HY*HY
      Q(5)  = Q(5) - COEF0*HY*HZ
      Q(6)  = Q(6) - COEF0*HZ*HZ
C
      RETURN
      END 
C
C
C
C
C
C
      SUBROUTINE EVOLML(F,A,B,C,K,SIGLV,
     .                  IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC,
     .                  L,R,IERR)
C
C-----------------------------------------------------------------------
C
C     Purpose  : Calculate the coefficients in Madou and Leblond's expressions 
C                of the void's strain and rotation rates
C     Reference: K. Madou, J.B. Leblond, L. Morin, Numerical studies of 
C                porous ductile materials containing arbitrary ellipsoidal voids 
C                - II: Evolution of the length and orientation of the void axes,
C                European Journal of Mechanics A/Solids, vol. 42, pp. 490-507 
C                (2013) [especially Section 6, p. 505, and Appendix A, pp. 506-507]
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I    REAL*8   F           Porosity
C     I    REAL*8   A           Major semi-axis of the void
C     I    REAL*8   B           Intermediate semi-axis of the void
C     I    REAL*8   C           Minor semi-axis of the void
C     I    REAL*8   K           Parameter related to the void shape:
C                               =0, spherical or prolate spheroidal void
C                               =1, oblate spheroidal void 
C                               0<K<1, general ellipsoidal void
C     I    REAL*8   SIGLV       Stress tensor (SIG) in the local (L) principal 
C                               basis of the ellipsoidal void, in Voigt's (V) 
C                               "stress vector" form 
C
C     I    REAL*8   IA,IB,IC    Integrals IA,IB,IC
C     I    REAL*8   IAA,IBB,ICC Integrals IAA,IBB,ICC
C     I    REAL*8   IAB,IAC,IBC Integrals IAB,IAC,IBC
C                               (All of these integrals have been previously 
C                               calculated in subroutine CRITML)
C 
C     O    REAL*8   L           Localization tensor of the void's strain rate, 
C                               stored in the order L11,L21,L31,L12,L22,L32,
C                               L13,L23,L33,L44,L55,L66 (Voigt's notation)
C     O    REAL*8   R           Localization tensor of the void's rotation rate,
C                               stored in the order RXYXY,RXZXZ,RYZYZ
C     O    INTEGER  IERR        Error code:
C                               =1, error in the inversion of Id-(1-f)S (zero 
C                               determinant)
C                               =2, the equivalent stress is zero
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C
C     Arguments
      REAL*8  F,A,B,C,K,SIGLV(6),IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC,
     .        L(12),R(3)
      INTEGER IERR
C
C     Local variables
      REAL*8  A2,B2,C2,ABAR,EXZ,PI,COEF,
     .        SXXXX,SYYYY,SZZZZ,SXXYY,SXXZZ,SYYXX,SYYZZ,SZZXX,SZZYY,
     .        SXYXY,SXZXZ,SYZYZ,
     .        MAT(3,3),MAT44,MAT55,MAT66,LEL(3,3),LEL44,LEL55,LEL66,
     .        LELBAR11,LELBAR21,LELBAR22,LELBAR23,LELBAR31,LELBAR32,
     .        LELBAR33,LELBAR44,LELBAR55,LELBAR66
      REAL*8  SIGLVM,SIGLVD(6),SIGLVE,T,DETSIG,CS3PHI,TH,X,SGN,SQRTF,
     .        HSPH,HCYL33,HCYL44,EXZM,EXZN,HPROL22,HPROL33,HPROL44,
     .        HOBL,H22,H23,H32,H33,H44,H55,H66,
     .        LBAR11,LBAR21,LBAR22,LBAR23,LBAR31,LBAR32,LBAR33,
     .        LBAR44,LBAR55,LBAR66,PIXYXY,PIXZXZ,PIYZYZ,
     .        RELXYXY,RELXZXZ,RELYZYZ,KCYLXY,KPROLXY,KXY,KXZ,KYZ
      INTEGER M,N,IERR1
C
C     Initialization of error code
C     ****************************
C
      IERR = 0
C
C     Strain rate of the void
C     ***********************
C
C     Geometric parameters
C
      A2 = A**2
      B2 = B**2
      C2 = C**2
C     Major semi-axis ABAR of the completely flat confocal ellipsoid
      ABAR = SQRT(A2-C2)
C     Eccentricity of the void in the XZ plane
      EXZ = ABAR/A
C
C     Eshelby's (non-symmetric) first tensor S
C
      PI    = 4.D0*ATAN(1.D0)
      COEF  = 3.D0/4.D0/PI
      SXXXX = COEF*A2*IAA
      SYYYY = COEF*B2*IBB
      SZZZZ = COEF*C2*ICC
      SXXYY = COEF*B2*IAB
      SYYXX = COEF*A2*IAB
      SXXZZ = COEF*C2*IAC
      SZZXX = COEF*A2*IAC
      SYYZZ = COEF*C2*IBC
      SZZYY = COEF*B2*IBC
      SXYXY = (SXXYY+SYYXX)/2.D0
      SXZXZ = (SXXZZ+SZZXX)/2.D0
      SYZYZ = (SYYZZ+SZZYY)/2.D0
C
C     (Non-symmetric) tensor MAT = Id-(1-f)S in Voigt's form
C
      COEF     = 1.D0-F
      MAT(1,1) = 1.D0 - COEF*SXXXX
      MAT(2,2) = 1.D0 - COEF*SYYYY
      MAT(3,3) = 1.D0 - COEF*SZZZZ
      MAT(1,2) = - COEF*SXXYY
      MAT(2,1) = - COEF*SYYXX
      MAT(1,3) = - COEF*SXXZZ
      MAT(3,1) = - COEF*SZZXX
      MAT(2,3) = - COEF*SYYZZ
      MAT(3,2) = - COEF*SZZYY
C     Remark: In Voigt's form MAT44 = 1 - 2*(1-F)*SXYXY, similarly for MAT55 
C             and MAT66  
      MAT44    = 1.D0 - 2.D0*COEF*SXYXY
      MAT55    = 1.D0 - 2.D0*COEF*SXZXZ
      MAT66    = 1.D0 - 2.D0*COEF*SYZYZ
C
C     (Non-symmetric) elastic localization tensor LEL 
C     of the void's strain rate (inverse of MAT)
C
      CALL INVMAT(MAT,LEL,IERR1)
      IF(IERR1.EQ.1) THEN
        IERR = 1
        RETURN
      ENDIF
C     Remark: LEL44=2*LELXYXY, similarly for LEL55 and LEL66  
      LEL44 = 1.D0/MAT44	  
      LEL55 = 1.D0/MAT55
      LEL66 = 1.D0/MAT66
C
C     Transformation of LEL into LELBAR (for the plastic correction of LEL)
C
      LELBAR11 = 1.D0/F
      LELBAR21 = LEL(1,1) + LEL(1,2) + LEL(1,3) - 1.D0/F
      LELBAR22 = LEL(1,1) - (LEL(1,2)+LEL(1,3))/2.D0
      LELBAR23 = (LEL(1,2)-LEL(1,3))/2.D0
      LELBAR31 = LEL(2,1) + LEL(2,2) + LEL(2,3) - LEL(3,1) 
     .           - LEL(3,2) - LEL(3,3)
      LELBAR32 = LEL(2,1) - LEL(3,1) + (-LEL(2,2)-LEL(2,3)
     .           +LEL(3,2)+LEL(3,3))/2.D0
      LELBAR33 = (LEL(2,2)-LEL(2,3)-LEL(3,2)+LEL(3,3))/2.D0
      LELBAR44 = LEL44
      LELBAR55 = LEL55
      LELBAR66 = LEL66
C
C     Triaxiality and Lode's angle
C
      CALL STRMDE(SIGLV,SIGLVM,SIGLVD,SIGLVE)
      IF(SIGLVE.EQ.0.) THEN
        IERR = 2	  
        RETURN
      ENDIF
      T = SIGLVM/SIGLVE
      DETSIG = SIGLVD(1)*SIGLVD(2)*SIGLVD(3)
     .         + 2.D0*SIGLVD(4)*SIGLVD(5)*SIGLVD(6)
     .         - SIGLVD(4)*SIGLVD(4)*SIGLVD(3)
     .         - SIGLVD(5)*SIGLVD(5)*SIGLVD(2)
     .         - SIGLVD(6)*SIGLVD(6)*SIGLVD(1)
      CS3PHI = 27.D0/2.D0*DETSIG/SIGLVE**3
C
C     Transverse invariants
C
      TH = (SIGLV(2)+SIGLV(3))/2.D0/SIGLVE
      X  = (SIGLV(1)-(SIGLV(2)+SIGLV(3))/2.D0)/SIGLVE
C
C     "HIJ" coefficients
C
      SGN = 1.D0
      IF (T.LT.0.) SGN = -1.D0
      SQRTF  = SQRT(F)
C
      HSPH   = (9.D0-10.D0*SQRTF+8.*F)/7.D0-(1.D0-SQRTF)**8/5.D0*
     .         11.D0*T**4/(11.D0+ABS(T)**3)*(1.D0+SGN*CS3PHI/3.D0)	  	  
      HCYL33 = (3.D0-5.D0*SQRTF+4.D0*F)/2.D0-(1.D0-SQRTF)**5*9.D0*TH**4/
     .         (20.D0+TH**4)
      HCYL44 = 1.D0+11.D0*(1.D0-SQRTF)**11/2.D0
     .         -(1.D0-SQRTF)**11*84.D0*X**2/(1.D0+20.D0*X**2)
C
      M    = 30
      EXZM = EXZ**M
      N    = 50
      EXZN = EXZ**N
C
      HPROL22 = (1.D0-EXZM)*HSPH + EXZM
      HPROL33 = (1.D0-EXZM)*HSPH + EXZM*HCYL33	   
      HPROL44 = (1.D0-EXZM)*HSPH + EXZM*HCYL44	  
      HOBL    = (1.D0-EXZN)*HSPH + EXZN
C	  
      H22 = (1.D0-K)*HPROL22 + K*HOBL
      H23 = H22
      H32 = H22
      H33 = (1.D0-K)*HPROL33 + K*HOBL
      H44 = (1.D0-K)*HPROL44 + K*HOBL
      H55 = H44
      H66 = H33	  
C    
C     Plastic correction of LELBAR
C	  
      LBAR11 = LELBAR11
      LBAR21 = LELBAR21
      LBAR22 = H22*LELBAR22
      LBAR23 = H23*LELBAR23
      LBAR31 = LELBAR31
      LBAR32 = H32*LELBAR32
      LBAR33 = H33*LELBAR33
      LBAR44 = H44*LELBAR44	  
      LBAR55 = H55*LELBAR55
      LBAR66 = H66*LELBAR66
C
C     Inverse transformation of LBAR into L (localization tensor 
C     for the void's strain rate)
C
      L(1) = (LBAR21+2.D0*LBAR22+LBAR11)/3.D0
      L(2) = (-LBAR21-2.D0*LBAR22+LBAR31+2.D0*LBAR32+2.D0*LBAR11)/6.D0
      L(3) = (-LBAR21-2.D0*LBAR22-LBAR31-2.D0*LBAR32+2.D0*LBAR11)/6.D0
      L(4) = (LBAR21-LBAR22+3.D0*LBAR23+LBAR11)/3.D0
      L(5) = (-LBAR21+LBAR22-3.D0*LBAR23+LBAR31-LBAR32+3.D0*LBAR33
     .       +2.D0*LBAR11)/6.D0
      L(6) = (-LBAR21+LBAR22-3.D0*LBAR23-LBAR31+LBAR32-3.D0*LBAR33
     .       +2.D0*LBAR11)/6.D0
      L(7) = (LBAR21-LBAR22-3.D0*LBAR23+LBAR11)/3.D0
      L(8) = (-LBAR21+LBAR22+3.D0*LBAR23+LBAR31-LBAR32-3.D0*LBAR33
     .       +2.D0*LBAR11)/6.D0
      L(9) = (-LBAR21+LBAR22+3.D0*LBAR23-LBAR31+LBAR32+3.D0*LBAR33
     .       +2.D0*LBAR11)/6.D0 
      L(10) = LBAR44
      L(11) = LBAR55
      L(12) = LBAR66
C
C     Rotation rate of the voids
C     **************************
C
C     Eshelby's second tensor PI
C
      COEF   = 8.D0*PI
      PIXYXY = (IB-IA)/COEF
      PIXZXZ = (IC-IA)/COEF	  
      PIYZYZ = (IC-IB)/COEF
C
C     Elastic localization tensor REL for the void's rotation rate
C
      COEF    = 1.D0-F
      RELXYXY = COEF*PIXYXY*LEL44
      RELXZXZ = COEF*PIXZXZ*LEL55
      RELYZYZ = COEF*PIYZYZ*LEL66
C
C     "KIJ" coefficients
C
      KCYLXY  = 1.D0+11.D0*(1.D0-SQRTF)**11-(1.D0-SQRTF)**11*170.D0*
     .          X**2/(1.D0+20.D0*X**2)
      KPROLXY = 1.D0 - EXZM + EXZM*KCYLXY
      KXY     = (1.D0-K)*KPROLXY + K	  
      KXZ     = KXY		  
      KYZ     = 1.D0
C
C     Localization tensor R for the void's rotation rate
C
      R(1) = KXY*RELXYXY
      R(2) = KXZ*RELXZXZ
      R(3) = KYZ*RELYZYZ	  
C
      RETURN
      END	  
C
C
C
C
C
C
      SUBROUTINE FUNCRD(X,Y,Z,RD,IERR)
C
C-----------------------------------------------------------------------
C
C     Purpose  : Calculation of Carlson's function RD(X,Y,Z) 
C     Reference: W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery,
C                Numerical Recipes in Fortran 77 - The Art of Scientific 
C                Computing, Second Edition, Cambridge University Press (1992), 
C                pp. 257-258
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I    REAL*8   X,Y,Z     Arguments of RD
C     O    REAL*8   RD        Carlson's function
C     O    INTEGER  IERR      Error code: 
C                             = 1 if the arguments are invalid
C                             = 2 if the iterations do not converge
C
C-----------------------------------------------------------------------

      IMPLICIT NONE
C
C     Arguments
      REAL*8 X,Y,Z,RD
      INTEGER IERR
C
C     Local variables
      REAL*8 ERRTOL,SMALL,BIG,C1,C2,C3,C4,C5,C6
      REAL*8 ALAMB,AVE,DELX,DELY,DELZ,EA,EB,EC,ED,EE,FAC,
     .       SQRTX,SQRTY,SQRTZ,SOM,XT,YT,ZT
      INTEGER ITE
C	  
      PARAMETER (ERRTOL=.05D0,SMALL=1.D-25,BIG=4.5D21,C1=3.D0/14.D0,
     .           C2=1.D0/6.D0,C3=9.D0/22.D0,C4=3.D0/26.D0,C5=.25D0*C3,
     .           C6=1.5D0*C4)
C
      IERR = 0
C
      IF( (MIN(X,Y).LT.0.D0) .OR. (MIN(X+Y,Z).LT.SMALL) .OR.
     .    (MAX(X,Y,Z).GT.BIG) ) THEN
        IERR = 1
        RETURN	  
      ENDIF
C
      XT  = X
      YT  = Y
      ZT  = Z
      SOM = 0.D0
      FAC = 1.D0
C	  
      DO ITE=1,20	  
        SQRTX = SQRT(XT)
        SQRTY = SQRT(YT)
        SQRTZ = SQRT(ZT)
        ALAMB = SQRTX*(SQRTY+SQRTZ)+SQRTY*SQRTZ
        SOM   = SOM + FAC/(SQRTZ*(ZT+ALAMB))
		FAC   = .25D0*FAC
        XT    = .25D0*(XT+ALAMB)
        YT    = .25D0*(YT+ALAMB)
        ZT    = .25D0*(ZT+ALAMB)
        AVE   = .2D0*(XT+YT+3.D0*ZT)
        DELX  = (AVE-XT)/AVE
        DELY  = (AVE-YT)/AVE
        DELZ  = (AVE-ZT)/AVE
        IF ( MAX(ABS(DELX),ABS(DELY),ABS(DELZ)).LT.ERRTOL ) GOTO 10
      ENDDO
C
      IERR = 2
      RETURN	  
C
   10 CONTINUE	  
      EA = DELX*DELY
      EB = DELZ*DELZ
      EC = EA - EB
      ED = EA - 6.D0*EB
      EE = ED + EC + EC
      RD = 3.D0*SOM + FAC*(1.D0+ED*(-C1+C5*ED-C6*DELZ*EE)
     .     + DELZ*(C2*EE+DELZ*(-C3*EC+DELZ*C4*EA)))/(AVE*SQRT(AVE))
C	  
      RETURN
      END   
C
C
C
C
C
C	  
      SUBROUTINE INVMAT(MAT,MATINV,IERR)
C 
C-----------------------------------------------------------------------
C
C     Purpose: Inversion of an arbitrary (not necessarily symmetric) 
C              3*3 matrix 
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I   REAL*8    MAT      Matrix to be inverted
C     O   REAL*8    MATINV   Inverse of MAT
C     O   INTEGER   IERR     Error code: = 1 if the determinant is zero
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C  
C     Arguments
      REAL*8 MAT(3,3),MATINV(3,3)
      INTEGER IERR
C
C     Local variables
      REAL*8 DET
C	  
      IERR = 0
C
      DET = MAT(1,1)*MAT(2,2)*MAT(3,3)
     .      + MAT(1,2)*MAT(2,3)*MAT(3,1)
     .	    + MAT(1,3)*MAT(2,1)*MAT(3,2)
     .      - MAT(1,3)*MAT(2,2)*MAT(3,1)
     .      - MAT(1,1)*MAT(2,3)*MAT(3,2)	
     .      - MAT(1,2)*MAT(2,1)*MAT(3,3) 
      IF (DET.EQ.0.D0) THEN
        IERR = 1
        RETURN
      ENDIF
C			
      MATINV(1,1) = (MAT(2,2)*MAT(3,3)-MAT(2,3)*MAT(3,2))/DET
      MATINV(1,2) = (MAT(1,3)*MAT(3,2)-MAT(1,2)*MAT(3,3))/DET
      MATINV(1,3) = (MAT(1,2)*MAT(2,3)-MAT(1,3)*MAT(2,2))/DET
      MATINV(2,1) = (MAT(2,3)*MAT(3,1)-MAT(2,1)*MAT(3,3))/DET
      MATINV(2,2) = (MAT(1,1)*MAT(3,3)-MAT(1,3)*MAT(3,1))/DET
      MATINV(2,3) = (MAT(1,3)*MAT(2,1)-MAT(1,1)*MAT(2,3))/DET
      MATINV(3,1) = (MAT(2,1)*MAT(3,2)-MAT(2,2)*MAT(3,1))/DET
      MATINV(3,2) = (MAT(1,2)*MAT(3,1)-MAT(1,1)*MAT(3,2))/DET
      MATINV(3,3) = (MAT(1,1)*MAT(2,2)-MAT(1,2)*MAT(2,1))/DET
C	  	  
      RETURN
      END	
C
C
C
C
C
C	
      SUBROUTINE PROJML(ELDAT,PLDAT,VOIDAT,SIG0,INTPAR0,GRTRSF,
     .                  ICONV,ITGTMT,SIG1,INTPAR1,QFG1G,TGTMAT,ERRINFO)
C
C-----------------------------------------------------------------------
C
C     Purpose  : Plastic correction (projection onto the yield locus) 
C                of the elastic stress predictor in Leblond-Madou's model  
C     Reference: In construction
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I    REAL*8   ELDAT        Vector containing elastic data
C                                ELDAT(1) = Young's modulus
C                                ELDAT(2) = Poisson's ratio
C     I    REAL*8   PLDAT        Vector containing plastic data 
C                                PLDAT(1)     = Number of points (EPSILON_EQ,SIGMA)
C                                               on the uniaxial stress-strain curve 
C                                               of the sound material
C                                PLDAT(2*I)   = I-th value of EPSILON_EQ on the 
C                                               stress-strain curve
C                                PLDAT(2*I+1) = I-th value of SIGMA on the 
C                                               stress-strain curve
C     I    REAL*8   VOIDAT       Vector containing data pertaining to voids
C                                VOIDAT(1)    = Initial porosity
C                                VOIDAT(2-7)  = Initial quadratic form characterizing
C                                               the orientation and semi-axes of the 
C                                               void, stored in "Voigt's stress vector
C                                               form"
C                                Remark: this quadratic form is defined by 
C                                        P(U) = (U.E_X)**2/A**2 + (U.E_Y)**2/B**2 
C                                               + (U.E_Z)**2/C**2
C                                        where A,B,C are the semi-axes of the void 
C                                        and E_X,E_Y,E_Z the unit vectors collinear 
C                                        to these axes 
C                                VOIDAT(8)    = Tvergaard's "Q"-parameter
C                                Remark: this parameter is modified according to 
C                                        the void shape
C                                VOIDAT(9,10) = Parameters pertaining to coalescence
C                                  VOIDAT(9)  : "Critical" porosity [F+G]_C 
C                                  VOIDAT(10) : "Accelerating factor" DLT (equivalent
C                                               of parameter (1/Q-F_C)/(F_U-F_C) in
C                                               the GTN model)
C                                VOIDAT(11-13)= Minimum values of the semi-axes of 
C                                               the void 
C                                  VOIDAT(11) : Minimum value of major semi-axis
C                                  VOIDAT(12) : Minimum value of intermediate
C                                               semi-axis                                             
C                                  VOIDAT(13) : Minimum value of minor semi-axis
C                                VOIDAT(14-16)= Parameters pertaining to continuous
C                                               nucleation - Chu and Needleman's
C                                               modelling
C                                  VOIDAT(14) : Maximum porosity which can be generated
C                                  VOIDAT(15) : Mean strain 
C                                  VOIDAT(16) : Standard deviation
C     I    REAL*8   SIG0         Stress tensor at time T, stored in "Voigt's 
C                                stress vector form"
C     I    REAL*8   INTPAR0      Vector containing internal parameters at time T
C                                INTPAR0(1)   = Porosity
C                                INTPAR0(2-7) = Quadratic form characterizing the 
C                                               orientation and semi-axes of the 
C                                               void, stored in "Voigt's stress 
C                                               vector form"
C                                INTPAR0(8)   = (Isotropic) strain hardening parameter
C     I    REAL*8   GRTRSF       Gradient of transformation between times T and T+DT 
C     I    INTEGER  ICONV        Convergence code of the global elastoplastic
C                                iterations:
C                                ICONV=0: Convergence is not achieved, internal
C                                         parameters at time T+DT are not required
C                                ICONV=1: Convergence is achieved, internal parameters 
C                                         at time T+DT are required
C                                Remark:  if not available in the calling 
C                                         programme, ICONV must be set to 1 there
C     I    INTEGER  ITGTMT       Code for calculation of the consistent tangent matrix:
C                                ITGTMT=0: Calculation of the consistent tangent 
C                                          matrix is not required 
C                                ITGTMT=1: Calculation of the consistent tangent 
C                                          matrix is required 
C                                Remarks: 1) If not available in the calling 
C                                            programme, ITGTMT must be set to 1 
C                                            there
C                                         2) If ITGTMT=0, the tangent matrix 
C                                            TGTMAT is equated to the elastic 
C                                            stiffness matrix
C
C     I/O  REAL*8   SIG1         Stress tensor at time T+DT, stored in "Voigt's 
C                                stress vector form"
C                                Input:  elastic stress predictor
C                                Output: plastically corrected stress tensor
C
C     O    REAL*8   INTPAR1      Vector containing internal parameters at time 
C                                T+DT (same conventions as for INTPAR0)
C     O    REAL*8   QFG1G        Q(F+G)*/(1+G): parameter characterizing damage -
C                                value 0 when the material is sound, 1 when it is 
C                                completely destroyed (output for postprocessing only)
C     O    REAL*8   TGTMAT       Consistent tangent matrix (not programmed yet 
C                                - in practice a BFGS method not using this 
C                                matrix works well)
C     O    INTEGER               Vector containing information about errors:
C                                ERRINFO(1): error code:
C                                  =1, the number of points on the stress-strain 
C                                      curve is too large
C                                  =2, error in subroutine DSYEV
C                                  =3 to 9, error in subroutine CRITML:
C                                    =3, error on the semi-axes of the void
C                                    =4, the iterations needed to calculate 
C                                        LAMBDA do not converge
C                                    =5, the arguments of subroutine FUNCRD for 
C                                        the calculation of IA are invalid
C                                    =6, the iterations needed to calculate IA 
C                                        in FUNCRD do not converge
C                                    =7, the arguments of subroutine FUNCRD for 
C                                        the calculation of IB are invalid
C                                    =8, the iterations needed to calculate IB 
C                                        in FUNCRD do not converge
C                                    =9, non-invertible matrix T (zero determinant) 
C                                  =10 or 11, error in subroutine YLDSTR:
C                                    =10, the first value of the equivalent strain 
C                                         in the vector STRAIN is nonzero
C                                    =11, the first value of the yield stress in 
C                                         the vector STRESS is zero
C                                  =12, non-invertible matrix MAT (zero 
C                                       determinant)
C                                  =13, the Newton iterations on SIGH do not 
C                                       converge
C                                  =14, unsuccessful relaxation on DLAMB 
C                                  =15, the Newton iterations on the plastic 
C                                       multiplier do not converge 
C                                  =16, the fixed point iterations on the strain 
C                                       hardening parameter do not converge 
C                                  =17 or 18, error in subroutine EVOLML:
C                                    =17, error in the inversion of Id-(1-f)S 
C                                         (zero determinant)
C                                    =18, the equivalent stress is zero 
C                                  =19, error in subroutine CONTML
C                                ERRINFO(2): number of current iteration on the 
C                                            strain hardening parameter
C                                ERRINFO(3): number of current iteration on the 
C                                            plastic multiplier 
C 
C-----------------------------------------------------------------------
C
C     Remarks:
C     1) The projection algorithm is implicit with respect to the components of 
C        the plastic strain increment and the hardening parameter, but explicit 
C        with respect to the porosity and the quadratic form characterizing the 
C        orientation and semi-axes of the void.
C     2) Stress-type 2nd-rank 3*3 tensors SIG are stored in "Voigt's stress vector
C        form", that is as 6-vectors of components SIGXX,SIGYY,SIGZZ,SIGXY,SIGXZ,
C        SIGYZ. Deformation-type 2nd-rank 3*3 tensors DEF are stored in "Voigt's 
C        deformation vector form", that is as 6-vectors of components DEFXX,DEFYY,
C        DEFZZ,2*DEFXY,2*DEFXZ,2*DEFYZ.
C     3) Often but not systematically:
C        * Tensors expressed in the global basis are denoted with an extra symbol 
C          "G"
C        * Tensors expressed in the local basis of the ellipsoidal void are denoted 
C          with an extra symbol "L"  
C        * Tensors expressed in tensorial form are denoted with an extra symbol 
C          "T"    
C        * Tensors expressed in Voigt's form are denoted with an extra symbol "V"  
C        For instance AGT is a tensor expressed in the global basis in tensorial 
C                     form;
C                     ALV is the same tensor but expressed in the local basis and 
C                     in Voigt's form.
C
C-----------------------------------------------------------------------
C     
      IMPLICIT NONE
C
C     Maximum number of points on the experimental stress-strain curve
      INTEGER NCRVM
      PARAMETER (NCRVM=100)
C     Maximum dimension of the vector containing plastic data (=2*NCRVM+1)
      INTEGER NPLM
      PARAMETER (NPLM=201)
C
C     Arguments
      REAL*8  ELDAT(2),PLDAT(NPLM),VOIDAT(16),SIG0(6),INTPAR0(8),
     .        GRTRSF(3,3),SIG1(6),INTPAR1(8),QFG1G,TGTMAT(6,6)
      INTEGER ICONV,ITGTMT,ERRINFO(3)
C
C     Local variables
      REAL*8  TEST,E,NU,LAMBDA,MU,STRAIN(NCRVM),STRESS(NCRVM),
     .        QTV,FGC,DLT,AMINI,BMINI,CMINI,FMAX,EN,SN,
     .        F0,P0GV(6),EPSB0,P0GT(3,3),P0GTS(3,3),EIGNVL(3),WORK(8),
     .        A,B,C,ROT(3,3),DET,K,Q(9),G,KAPPA,HX,HY,HZ,
     .        IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC
      REAL*8  COEFA1,QPROL,QOBL,SIG1GT(3,3),SIG1LT(3,3),
     .        SIG1LV(6),SIG1LVS(6),QLT(3,3),QLT44,QLT55,QLT66,SIGB,
     .        QFG,COEFA2,SIGH,QUAD,CRIT,F1,DTRVOIDG(3,3),P1GT(3,3),
     .        P1GV(6),EPSB1,DEPSB,DDEPSB,SIGB2,SIG1M,SIG1D(6),SIG1Q
      REAL*8  DLAMB,PRECLAMB,COEFB1,MAT(3,3),DMATDL(3,3),MATINV(3,3),
     .        MATINV44,DMATDL44,MATINV55,DMATDL55,MATINV66,DMATDL66,
     .        SIG1TRL(3),COEFB2,COEFB3(3),PRECSIGH,SH,CH,COEFC,
     .        SIG1AUX(6),DS1XDSH(3),DS1DSH(3),FUNCG,DFUNCG,FUNCF
      REAL*8  DS1XDL(3),DMDLS1(3),DS1DL(6),DFDL,DFDSH,DGDL,DFUNCF,
     .        DDLAMB,DEPLV(6),PSCAL,DDEPSB1,PRECEB1,C1,C2,PRECEB2,DF,
     .        PI,ARG,L(12),R(3),DEVOIDLV(6),GRTRSFL(3,3),
     .        DRMATLXY,DRMATLXZ,DRMATLYZ,DRVOIDLXY,DRVOIDLXZ,DRVOIDLYZ,
     .        DTRVOIDL(3,3)
      INTEGER IRUPT,I,J,NCRV,LWORK,INFO,IERR,M,NITEPSB,ITEPSB,
     .        NITLAMB,ITLAMB,NITSIGH,ITSIGH,NITEREL,ITEREL
C
C
C     Initialization of error code
C     ****************************
C
      ERRINFO(1) = 0
C
C     Initialization of rupture code
C     ******************************
C 
      IRUPT = 0
      TEST  = 0.D0
      DO I=1,6 
        TEST = TEST + SIG0(I)**2
      ENDDO
C     If, at time T, the porosity is nonzero (this excludes the initial instant) 
C     but all stress components vanish, the material is considered to be destroyed 
C     at time T+DT (irreversible destruction)
      IF ( (INTPAR0(1).NE.0.D0) . AND. (TEST.EQ.0.D0) ) IRUPT = 1
C 	  
C     Elastic data
C     ************
C
C     E,NU
      E  = ELDAT(1)
      NU = ELDAT(2)
C
C     LAMBDA,MU
      LAMBDA = E*NU/(1.D0+NU)/(1.D0-2.D0*NU)
      MU     = E/2.D0/(1.D0+NU)
C 
C     Initialization of the tangent matrix TGTMAT
C     *******************************************
C     TGTMAT is equated to the elastic stiffness matrix
C 
      DO I=1,3
        DO J=1,3 	  
          TGTMAT(I+3,J) = 0.D0
          TGTMAT(I,J+3) = 0.D0
          IF (I.EQ.J) THEN
C           Diagonal components		  
            TGTMAT(I,I)     = LAMBDA + 2.D0*MU	
            TGTMAT(I+3,I+3) = MU
          ELSE
C           Off-diagonal components 
            TGTMAT(I,J)     = LAMBDA	
            TGTMAT(I+3,J+3) = 0.D0
          ENDIF 
        ENDDO
      ENDDO
C 
C     Plastic data
C     ************
C  
C     Number of points on the stress-strain curve 
      NCRV = INT(PLDAT(1))
      IF (NCRV.GT.NCRVM) THEN
        ERRINFO(1) = 1
        RETURN	
      ENDIF
C 
C     Points on the stress-strain curve 
      DO I=1,NCRV
        STRAIN(I) = PLDAT(2*I)
        STRESS(I) = PLDAT(2*I+1)
      ENDDO	  
C
C     Data pertaining to voids
C     ************************
C 
C     Initialization of the porosity and the quadratic form
C     characterizing the orientation and semi-axes of the void
      IF (INTPAR0(1).EQ.0.D0) THEN
        CALL TRNSFR(VOIDAT,INTPAR0,7)
C       Precaution
        INTPAR0(8) = 0.D0
      ENDIF
C 
C     Tvergaard's Q-parameter
C     Remark: Q comes as a multiplicative factor of F+G in the yield criterion
      QTV = VOIDAT(8)
C
C     Parameters pertaining to coalescence
      FGC = VOIDAT(9)
      DLT = VOIDAT(10)
C 
C     Minimum values of the semi-axes of the voids
      AMINI = VOIDAT(11) 
      BMINI = VOIDAT(12) 
      CMINI = VOIDAT(13) 
C 
C     Parameters pertaining to continuous nucleation
      FMAX = VOIDAT(14)
      EN   = VOIDAT(15)
      SN   = VOIDAT(16)
C 
C     Internal parameters at time T
C     *****************************
C 
C     Porosity
      F0 = INTPAR0(1)
C     Quadratic form P0 characterizing the orientation and semi-axes of the void
C     - global basis, Voigt's form
      CALL TRNSFR(INTPAR0(2),P0GV,6)
C     Strain hardening parameter - the absolute value is introduced because this 
C     parameter is conventionally considered as negative in the case of elastic 
C     unloading 
      EPSB0 = ABS(INTPAR0(8))
C 
C     Semi-axes of the void and transformation matrix (rotation) 
C     from the global basis to the local one 
C     **********************************************************
C
C     Transformation of P0GV into tensorial form
      CALL VGTTNS(P0GV,P0GT,1)
C     P0GT must be saved because it is destroyed by subroutine DSYEV
      CALL TRNSFR(P0GT,P0GTS,9)
C
C     Diagonalization of the quadratic form P0GT  
C     Input : P0GT(3,3)    symmetric matrix to be diagonalized
C     Output: EIGNVL(3)    eigenvalues in ascending order
C             P0GT(3,3)    matrix whose columns contain the corresponding 
C                          orthonormal eigenvectors
      LWORK = 8
      CALL DSYEV('V','L',3,P0GT,3,EIGNVL,WORK,LWORK,INFO)
      IF ( (EIGNVL(1).LE.0.D0) .OR. (EIGNVL(2).LE.0.D0) 
     .     .OR. (EIGNVL(3).LE.0.D0) .OR. (INFO.NE.0) ) THEN
        ERRINFO(1) = 2
        RETURN
      ENDIF
C 
C     Semi-axes of the void
      A = 1.D0/SQRT(EIGNVL(1))
      B = 1.D0/SQRT(EIGNVL(2))
      C = 1.D0/SQRT(EIGNVL(3))
C 
C     Transformation matrix (rotation) from the global basis to the local one 
      CALL TRNSFR(P0GT,ROT,9)	  
C     Change of sign of the first column of ROT if its determinant is negative
C     (subroutine DSYEV does not warrant that the matrix giving the orthonormal 
C     eigenvectors is a rotation)
      DET = ROT(1,1)*ROT(2,2)*ROT(3,3) + ROT(1,2)*ROT(2,3)*ROT(3,1)
     .	    + ROT(1,3)*ROT(2,1)*ROT(3,2) - ROT(1,3)*ROT(2,2)*ROT(3,1)
     .      - ROT(1,2)*ROT(2,1)*ROT(3,3) - ROT(1,1)*ROT(2,3)*ROT(3,2)
      IF(DET.LT.0.) THEN
        DO I=1,3
          ROT(I,1) = - ROT(I,1)
        ENDDO
      ENDIF
C
C     Restoration of P0GT
      CALL TRNSFR(P0GTS,P0GT,9)
C
C     Parameters of the yield criterion 
C     *********************************
C
      CALL CRITML(F0,A,B,C,K,Q,G,KAPPA,HX,HY,HZ,
     .            IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC,IERR)
      IF (IERR.NE.0) THEN
        ERRINFO(1) = IERR+2
        RETURN
      ENDIF
C
C     Modification of Tvergaard's Q-parameter according to the void shape
C     *******************************************************************
C     Reference: In construction 
C
C     (Delete the following 4 lines if the modification of QTV is not desired) 
      COEFA1 = 2.D0*A/C/(1.D0+(A/C)**2)
      QPROL  = 1.D0+(QTV-1.D0)*COEFA1**(3.D0/2.D0)
      QOBL   = 1.D0+(QTV-1.D0)*COEFA1
      QTV    = (1.D0-K)*QPROL + K*QOBL
C
C     Heuristic modeling of coalescence: if F+G is larger than [F+G]_C, 
C     it is replaced by [F+G]_C + DLT*(F+G - [F+G]_C)
C     *****************************************************************
C     This is equivalent to a suitable modification of Tvergaard's Q-parameter
C
      IF ((F0+G).GT.FGC) QTV = DLT*QTV + (1.D0-DLT)*QTV*FGC/(F0+G)
C
C     Case of total destruction
C     *************************
C
C     When QFG1G becomes unity the stress tensor is necessarily zero
      QFG1G = QTV*(F0+G)/(1.D0+G)
      IF (QFG1G.GE.1.D0) IRUPT = 1
C 
      IF (IRUPT.EQ.1) THEN
C       Total damage - the stress components must be set to zero, and the tangent 
C       matrix possibly calculated  
C       Remark: activate the instructions with "CC" below when the tangent matrix 
C       is programmed
CC        IF (ITGTMT.EQ.1) THEN
C         The tangent matrix is calculated with QFG1G = 0.99 (large but not total 
C         damage) in order to avoid getting a zero result - the stress components 
C         will be set to zero at the end
CC          QFG1G = 0.99D0
CC          QTV = QFG1G*(1.D0+G)/(F0+G)
CC        ELSE
C         The damage parameter is set to unity
          QFG1G = 1.D0
C         All components of the stress tensor are set to zero		
          DO I=1,6 
            SIG1(I) = 0.D0
          ENDDO 
C         The evolutions of internal parameters are disregarded
          CALL TRNSFR(INTPAR0,INTPAR1,8)
          RETURN
CC        ENDIF		  
      ENDIF		
C
C     Elastic stress predictor in the local basis, Voigt's form
C     *********************************************************
C     Starting point: global basis, Voigt's form
C
C     Global basis, Voigt's form -> Global basis, tensorial form
      CALL VGTTNS(SIG1,SIG1GT,1)
C     Global basis, tensorial form -> Local basis, tensorial form
      CALL CHGBAS(SIG1GT,SIG1LT,ROT,1)
C     Local basis, tensorial form -> Local basis, Voigt's form
      CALL VGTTNS(SIG1LV,SIG1LT,2)
C     SIG1LV must be saved because it will be used later - the elastic stress 
C     predictor is in SIG1LVS from now on
      CALL TRNSFR(SIG1LV,SIG1LVS,6)
C   
C     Matrix form of the quadratic form Q (first term of the criterion)
C     of the components of the stress tensor - local basis
C     *****************************************************************
C
      QLT(1,1) = Q(1)
      QLT(1,2) = Q(2)
      QLT(1,3) = Q(3)
      QLT(2,1) = QLT(1,2)
      QLT(2,2) = Q(4)
      QLT(2,3) = Q(5)
      QLT(3,1) = QLT(1,3)
      QLT(3,2) = QLT(2,3)
      QLT(3,3) = Q(6)
      QLT44    = Q(7)
      QLT55    = Q(8)
      QLT66    = Q(9)
C 
C     Current yield stress for the strain hardening parameter at time T
C     *****************************************************************
C 
      CALL YLDSTR(STRAIN,STRESS,NCRV,EPSB0,SIGB,IERR)
      IF (IERR.NE.0) THEN
        ERRINFO(1) = IERR+9
        ERRINFO(2) = 0
        RETURN
      ENDIF
C 
C     Yield criterion
C     ***************
C  
      QFG     = QTV*(F0+G)
      COEFA1  = 2.D0*(1.D0+G)*QFG
      COEFA2  = (1.D0+G)**2 + QFG**2
      SIGH    = HX*SIG1LV(1)+HY*SIG1LV(2)+HZ*SIG1LV(3)
      QUAD    = 0.D0
      DO I=1,3 
        DO J=1,3
          QUAD = QUAD + QLT(I,J)*SIG1LV(I)*SIG1LV(J)
        ENDDO 
      ENDDO	  
      QUAD = QUAD + QLT44*SIG1LV(4)**2 + QLT55*SIG1LV(5)**2
     .       + QLT66*SIG1LV(6)**2
      CRIT = QUAD/SIGB**2 + COEFA1*COSH(KAPPA*SIGH/SIGB) - COEFA2	  	
C
C     Elastic unloading
C     *****************
C 
      IF (CRIT.LE.0.D0) THEN
        IF (ICONV.EQ.1) THEN
C         Convergence of the global elastoplastic iterations achieved - 
C         the internal parameters at time T+DT must be calculated and stored 
C 
C         Porosity at time T+DT
          F1 = F0
C 
C         Quadratic form P1 characterizing the shape and orientation of the void 
C         at time T+DT - the void's increment of rotation is (1/2)[GRTRSF-(GRTRSF)^T]
          DO I=1,3 
            DO J=1,3
              DTRVOIDG(I,J)=(GRTRSF(I,J)-GRTRSF(J,I))/2.D0
            ENDDO				
          ENDDO
C         P1, global basis, tensorial form 
          DO I=1,3
            DO J=1,3
              P1GT(I,J) = P0GT(I,J)
              DO M=1,3
                P1GT(I,J) = P1GT(I,J) - P0GT(I,M)*DTRVOIDG(M,J)
     .                      - DTRVOIDG(M,I)*P0GT(M,J)
              ENDDO
            ENDDO
          ENDDO
C         Transformation of P1GT into Voigt's form
          CALL VGTTNS(P1GV,P1GT,2)
C 
C         Strain hardening parameter at time T+DT (conventionally negative in the 
C         case of elastic unloading)
          EPSB1 = - EPSB0
C 		  
C         Storage of internal parameters at time T+DT
          INTPAR1(1) = F1
          CALL TRNSFR(P1GV,INTPAR1(2),6)
          INTPAR1(8) = EPSB1
        ENDIF	   
        RETURN
      ENDIF
C 
C     Iterations on the strain hardening parameter (fixed point method)
C     *****************************************************************
C
C     Initialization
      DEPSB   = 0.D0
      DDEPSB  = 0.D0
      NITEPSB = 50
C 
C     Iterations
      DO ITEPSB=1,NITEPSB
C     
C       Current strain hardening parameter
        EPSB1 = EPSB0 + DEPSB
C       Current yield stress
        CALL YLDSTR(STRAIN,STRESS,NCRV,EPSB1,SIGB,IERR)
        IF (IERR.NE.0) THEN
          ERRINFO(1) = IERR+9
          ERRINFO(2) = ITEPSB
          RETURN
        ENDIF
        SIGB2 = SIGB**2
C
C       Iterations on the plastic multiplier (Newton method)
C       ----------------------------------------------------
C
C       Initialization - the value corresponding to von Mises's criterion is 
C       adopted
        CALL STRMDE(SIG1,SIG1M,SIG1D,SIG1Q)
C       The absolute value is a precaution 
        DLAMB    = ABS((SIG1Q-SIGB)*SIGB/6.D0/MU)
        NITLAMB  = 50
        PRECLAMB = 1.D-9
C 
C       Iterations
        DO ITLAMB =1,NITLAMB
C
C         Matrix MAT = (Id + 2*DLAMB/SIGB**2*TGTMAT*QLT) (TGTMAT = elastic 
C         stiffness matrix at this stage) - upper left submatrix corresponding to 
C         diagonal components of SIG1
          COEFB1 = 2.D0*DLAMB/SIGB2
          DO I=1,3
            DO J=1,3
              MAT(I,J) = COEFB1*(TGTMAT(I,1)*QLT(1,J)
     .                   + TGTMAT(I,2)*QLT(2,J) + TGTMAT(I,3)*QLT(3,J))
C             Derivative of submatrix with respect to DLAMB
              DMATDL(I,J) = MAT(I,J)/DLAMB 
            ENDDO
            MAT(I,I) = MAT(I,I) + 1.D0
          ENDDO			 
C         Inverse of submatrix
          CALL INVMAT(MAT,MATINV,IERR)
          IF (IERR.NE.0) THEN
            ERRINFO(1) = 12
            ERRINFO(2) = ITEPSB
            ERRINFO(3) = ITLAMB
            RETURN
          ENDIF
C 
C         Matrix MAT = (Id + 2*DLAMB/SIGB**2*TGTMAT*QLT) - lower right submatrix 
C         corresponding to off-diagonal components of SIG1
C         Derivative with respect to DLAMB and inverse
          MATINV44 = COEFB1*TGTMAT(4,4)*QLT44
          DMATDL44 = MATINV44/DLAMB
          MATINV44 = 1.D0/(1.D0 + MATINV44)
          MATINV55 = COEFB1*TGTMAT(5,5)*QLT55
          DMATDL55 = MATINV55/DLAMB
          MATINV55 = 1.D0/(1.D0 + MATINV55)
          MATINV66 = COEFB1*TGTMAT(6,6)*QLT66
          DMATDL66 = MATINV66/DLAMB
          MATINV66 = 1.D0/(1.D0 + MATINV66)
C
C         Iterations on SIGH (Newton method)
C         ..................................
C 
C         Initialization
C         Trial value SIG1TRL=MAT**(-1)*SIG1LVS of the diagonal part of SIG1LV
          DO I=1,3
            SIG1TRL(I) = MATINV(I,1)*SIG1LVS(1) + MATINV(I,2)*SIG1LVS(2)
     .                   + MATINV(I,3)*SIG1LVS(3) 
          ENDDO
C         The trial value H*SIG1TRL is used for SIGH
          SIGH    = HX*SIG1TRL(1) + HY*SIG1TRL(2) + HZ*SIG1TRL(3)
          COEFB2  = COEFA1*KAPPA*DLAMB/SIGB
          DO I=1,3
            COEFB3(I) = COEFB2*(TGTMAT(I,1)*HX + TGTMAT(I,2)*HY
     .                  + TGTMAT(I,3)*HZ)
          ENDDO
          NITSIGH  = 50
          PRECSIGH = 1.D-9
C
C         Iterations
          DO ITSIGH=1,NITSIGH
            SH = SINH(KAPPA*SIGH/SIGB)
            CH = SQRT(SH**2 + 1.D0)			   
C
C           Auxiliary stress tensor 
C           SIG1AUX = SIG1LVS - 2*(1+G)*QTV*(F0+G)*KAPPA*DLAMB/SIGB*SH*TGTMAT*H

            COEFC = KAPPA/SIGB*CH
            DO I=1,3
              SIG1AUX(I) = SIG1LVS(I) - COEFB3(I)*SH
C             Derivative of SIG1AUX with respect to SIGH
              DS1XDSH(I)   = - COEFC*COEFB3(I)
              SIG1AUX(I+3) = SIG1LVS(I+3)
            ENDDO  
C
C           Stress tensor at time T+DT - local basis, Voigt's form
            DO I=1,3
              SIG1LV(I) = MATINV(I,1)*SIG1AUX(1)
     .                    + MATINV(I,2)*SIG1AUX(2)
     .                    + MATINV(I,3)*SIG1AUX(3)
            ENDDO
            SIG1LV(4) = MATINV44*SIG1AUX(4)
            SIG1LV(5) = MATINV55*SIG1AUX(5)
            SIG1LV(6) = MATINV66*SIG1AUX(6)
C
C           Derivative of SIG1LV with respect to SIGH
            DO I=1,3
              DS1DSH(I) = MATINV(I,1)*DS1XDSH(1)+MATINV(I,2)*DS1XDSH(2)
     .                    +MATINV(I,3)*DS1XDSH(3) 
            ENDDO	 
C
C           Function G(SIGH) to be equated to zero
            FUNCG = (SIGH-HX*SIG1LV(1)-HY*SIG1LV(2)-HZ*SIG1LV(3))/SIGB
C           Derivative of FUNCG with respect to SIGH			
            DFUNCG = (1.D0-HX*DS1DSH(1)-HY*DS1DSH(2)-HZ*DS1DSH(3))/SIGB
C
C           Convergence test on SIGH
            IF (ABS(FUNCG).LE.PRECSIGH) GOTO 10
C
C           Update of SIGH (Newton method)
            SIGH = SIGH - FUNCG/DFUNCG	
C
C         End of loop on iterations on SIGH  
C         .................................
          ENDDO
C 
C         No convergence on SIGH 
          ERRINFO(1) = 13 
          ERRINFO(2) = ITEPSB
          ERRINFO(3) = ITLAMB
          RETURN 
C
C         Convergence on SIGH 
  10      CONTINUE
C 
C         Yield criterion - function F(DLAMB) to be equated to zero
          QUAD = 0.D0
          DO I=1,3 
            DO J=1,3
              QUAD = QUAD + QLT(I,J)*SIG1LV(I)*SIG1LV(J)
            ENDDO 
          ENDDO	  
          QUAD = QUAD + QLT44*SIG1LV(4)**2 + QLT55*SIG1LV(5)**2
     .           + QLT66*SIG1LV(6)**2
          FUNCF = QUAD/SIGB2 + COEFA1*CH - COEFA2	  	         
C
C         Convergence test on DLAMB
          IF (ABS(FUNCF).LE.PRECLAMB) GOTO 30
C
C         Derivative of FUNCF with respect to DLAMB
C         Beginning................................
C
C         Step 1: Derivative of FUNCF with respect to DLAMB at fixed SIGH
C         Derivative of SIG1AUX with respect to DLAMB
          DO I=1,3
            DS1XDL(I) = - COEFB3(I)*SH/DLAMB
          ENDDO			
C
C         (Derivative of MAT with respect to DLAMB)*SIG1LV
          DO I=1,3
		    DMDLS1(I) = DMATDL(I,1)*SIG1LV(1) + DMATDL(I,2)*SIG1LV(2) 
     .                  + DMATDL(I,3)*SIG1LV(3)
          ENDDO 
C
C         Derivative of SIG1LV with respect to DLAMB
          DO I=1,3
            DS1DL(I) = MATINV(I,1)*(DS1XDL(1) - DMDLS1(1))
     .                 + MATINV(I,2)*(DS1XDL(2) - DMDLS1(2))
     .                 + MATINV(I,3)*(DS1XDL(3) - DMDLS1(3))
          ENDDO
          DS1DL(4) = - MATINV44*DMATDL44*SIG1LV(4)
          DS1DL(5) = - MATINV55*DMATDL55*SIG1LV(5)
          DS1DL(6) = - MATINV66*DMATDL66*SIG1LV(6)
C
C         Derivative of FUNCF with respect to DLAMB (at fixed SIGH)
          DFDL = 0.D0
          DO I=1,3
            DO J=1,3
              DFDL = DFDL + QLT(I,J)*SIG1LV(I)*DS1DL(J)
            ENDDO
          ENDDO	
          DFDL = DFDL + QLT44*SIG1LV(4)*DS1DL(4) 
     .           + QLT55*SIG1LV(5)*DS1DL(5) + QLT66*SIG1LV(6)*DS1DL(6)
          DFDL = 2.D0*DFDL/SIGB2 
C
C         Step 2 : Derivative of FUNCF with respect to SIGH at fixed DLAMB
          DFDSH = COEFA1*KAPPA/SIGB*SH
          DO I=1,3
            DO J=1,3
              DFDSH = DFDSH + 2.D0*QLT(I,J)*SIG1LV(I)*DS1DSH(J)/SIGB2
            ENDDO
          ENDDO 
C
C         Step 3 : Derivative of FUNCG with respect to DLAMB at fixed SIGH
          DGDL = - (HX*DS1DL(1)+HY*DS1DL(2)+HZ*DS1DL(3))/SIGB
C
C         Step 4 : Total derivative of FUNCF with respect to DLAMB (with SIGH 
C         depending on DLAMB)
          DFUNCF = DFDL - DFDSH*DGDL/DFUNCG
C
C         End......................................
C
C         Update of DLAMB (Newton method)
C         Beginning......................
C 
          DDLAMB = - FUNCF/DFUNCF
C
C         Relaxation on DLAMB (must remain positive)
          NITEREL = 10
          DO ITEREL=1,NITEREL
            IF((DLAMB+DDLAMB).GE.0.D0) THEN
              GOTO 20
            ELSE
              DDLAMB = DDLAMB/2.D0
            ENDIF
          ENDDO
C
C         Unsuccessful relaxation on DLAMB 
          ERRINFO(1) = 14 
          ERRINFO(2) = ITEPSB
          ERRINFO(3) = ITLAMB
          RETURN 
C				 
C         Successful relaxation on DLAMB 
  20      CONTINUE
C 
C         New value of DLAMB 
          DLAMB = DLAMB + DDLAMB
C
C         End......................................
C 
C       End of loop on iterations on the plastic multiplier  
C       ---------------------------------------------------
        ENDDO   
C 		
C       No convergence on the plastic multiplier 
        ERRINFO(1) = 15
        ERRINFO(2) = ITEPSB
        RETURN 
C 
C       Convergence on the plastic multiplier 
  30    CONTINUE
C
C       Update of the strain hardening parameter (fixed point method)
C       Beginning----------------------------------------------------
C 
C       Increment of plastic strain - local basis, Voigt's form
        DO I=1,3
          DEPLV(I) = 2.D0*(QLT(I,1)*SIG1LV(1) + QLT(I,2)*SIG1LV(2)
     .               + QLT(I,3)*SIG1LV(3))/SIGB2
        ENDDO
        DEPLV(1) = DLAMB*DEPLV(1) + COEFB2*HX*SH
        DEPLV(2) = DLAMB*DEPLV(2) + COEFB2*HY*SH
        DEPLV(3) = DLAMB*DEPLV(3) + COEFB2*HZ*SH
        DEPLV(4) = COEFB1*QLT44*SIG1LV(4)
        DEPLV(5) = COEFB1*QLT55*SIG1LV(5)
        DEPLV(6) = COEFB1*QLT66*SIG1LV(6)
C 
C       Scalar product of the stress and the increment of plastic strain 
C       (cannot be negative)
        PSCAL = 0.D0
        DO I=1,6 
          PSCAL = PSCAL + SIG1LV(I)*DEPLV(I)
        ENDDO 
C 
C       New value of DEPSB (cannot be negative)
        DDEPSB1 = DDEPSB
        DDEPSB  = PSCAL/SIGB/(1.D0-F0) - DEPSB
C       Precaution
        IF ((DDEPSB1*DDEPSB).LT.0.D0) DDEPSB = DDEPSB/2.D0
        DEPSB = DEPSB + DDEPSB
C 
C       End----------------------------------------------------------
C
C       Convergence test on the strain hardening parameter
        PRECEB1 = 1.D-6
C       The values of C1 and C2 are adjusted in such a way that PRECEB2=1.D-5 for 
C       DEPSB=1.D-2 and PRECEB2=1.D-6 for DEPSB=1.D-6
        C1 = 10.D0**(-4.5D0)
        C2 = 0.25D0
        PRECEB2 = C1*DEPSB**C2
        IF ( (DEPSB.LE.PRECEB1) .OR. (ABS(DDEPSB).LE.PRECEB2) ) GOTO 40
C 
C     End of loop on iterations on the strain hardening parameter	  
C     ***********************************************************
      ENDDO
C 
C     No convergence on the strain hardening parameter
C     ************************************************
      ERRINFO(1) = 16
      RETURN 
C 
C     Convergence on the strain hardening parameter
C     ************************************************
  40  CONTINUE 
C 
C     Stress tensor at time T+DT in the global basis, Voigt's form
C     ************************************************************
C     Starting point: local basis, Voigt's form
C
C     Local basis, Voigt's form -> Local basis, tensorial form
      CALL VGTTNS(SIG1LV,SIG1LT,1)
C     Local basis, tensorial form -> Global basis, tensorial form
      CALL CHGBAS(SIG1GT,SIG1LT,ROT,2)
C     Global basis, tensorial form -> Global basis, Voigt's form
      CALL VGTTNS(SIG1,SIG1GT,2)
C 
C     Internal parameters at time T+DT (if convergence of the global
C     elastoplastic iterations is achieved) 
C     **************************************************************
C 
      IF (ICONV.EQ.1) THEN
C 
C       Porosity at time T+DT (growth term only)
        DF = (1.D0-F0)*(DEPLV(1)+DEPLV(2)+DEPLV(3))
        F1 = F0 + DF
C 
C       Localization tensors for the void's strain and rotation rates
        CALL EVOLML(F0,A,B,C,K,SIG1LV,
     .              IA,IB,IC,IAA,IBB,ICC,IAB,IAC,IBC,
     .              L,R,IERR)
        IF (IERR.NE.0) THEN
          ERRINFO(1) = IERR+16		
          RETURN
        ENDIF
C 
C       Increment of deformation of the void - local basis, Voigt's form
        DO I=1,3
          DEVOIDLV(I)   = L(I)*DEPLV(1)+L(I+3)*DEPLV(2)+L(I+6)*DEPLV(3)
          DEVOIDLV(I+3) = L(I+9)*DEPLV(I+3)
        ENDDO 
C 
C       Gradient of the increment of displacement in the local basis
        CALL CHGBAS(GRTRSF,GRTRSFL,ROT,1) 
C 
C       Increment of rotation of the material -local basis
        DRMATLXY = (GRTRSFL(1,2)-GRTRSFL(2,1))/2.D0
        DRMATLXZ = (GRTRSFL(1,3)-GRTRSFL(3,1))/2.D0
        DRMATLYZ = (GRTRSFL(2,3)-GRTRSFL(3,2))/2.D0
C 
C       Increment of rotation of the void - local basis
C       Remark: this increment of rotation is calculated using the material's  
C       total increment of rotation since its plastic increment of rotation
C       is unknown                            
        DRVOIDLXY = DRMATLXY + R(1)*DEPLV(4)
        DRVOIDLXZ = DRMATLXZ + R(2)*DEPLV(5)
        DRVOIDLYZ = DRMATLYZ + R(3)*DEPLV(6)
C 
C       Increment of transformation of the void - local basis
        DO I=1,3
          DTRVOIDL(I,I) = DEVOIDLV(I)
        ENDDO
        DTRVOIDL(1,2) = DEVOIDLV(4)/2.D0 + DRVOIDLXY 
        DTRVOIDL(2,1) = DEVOIDLV(4)/2.D0 - DRVOIDLXY 
        DTRVOIDL(1,3) = DEVOIDLV(5)/2.D0 + DRVOIDLXZ 
        DTRVOIDL(3,1) = DEVOIDLV(5)/2.D0 - DRVOIDLXZ
        DTRVOIDL(2,3) = DEVOIDLV(6)/2.D0 + DRVOIDLYZ 
        DTRVOIDL(3,2) = DEVOIDLV(6)/2.D0 - DRVOIDLYZ 		
C 
C       Increment of transformation of the void in the global basis
        CALL CHGBAS(DTRVOIDG,DTRVOIDL,ROT,2)
C 
C       Quadratic form P1 characterizing the orientation and semi-axes of the void 
C       at time T+DT - global basis, tensorial form
        DO I=1,3
          DO J=1,3
            P1GT(I,J) = P0GT(I,J)
            DO M=1,3
              P1GT(I,J) = P1GT(I,J) - P0GT(I,M)*DTRVOIDG(M,J)
     .                    - DTRVOIDG(M,I)*P0GT(M,J)
            ENDDO
          ENDDO
        ENDDO
C 
C       Possible modification of F1 and P1 due to the minimum values of the semi-axes
C       of the void (the correction of the porosity is done on that part arising
C       from growth only - nucleation is considered independent of contact)
        CALL CONTML(AMINI,BMINI,CMINI,F1,P1GT,IERR)
        IF (IERR.NE.0) THEN
          ERRINFO(1) = 19		
          RETURN
        ENDIF
C 		  
C       Correction of F1 due to nucleation (the newly created voids are assumed 
C       to appear with axes identical in size and orientation to those of 
C       pre-existing voids)
        IF ( (FMAX.GT.0.D0) .AND. (SN.GT.0.D0) ) THEN
          PI  = 4.D0*ATAN(1.D0)
          ARG = - 0.5D0*((EPSB1-EN)/SN)**2
          DF  = FMAX/SN/SQRT(2.D0*PI)*EXP(ARG)*DEPSB
          F1  = F1 + DF
        ENDIF
C 
C       Transformation of P1GT into Voigt's form
        CALL VGTTNS(P1GV,P1GT,2)
C
C       Storage of internal parameters at time T+DT
        INTPAR1(1) = F1
        CALL TRNSFR(P1GV,INTPAR1(2),6)
        INTPAR1(8) = EPSB1
C 
C     End of calculation of internal parameters at time T+DT
C     ******************************************************
      ENDIF
C 
C     Tangent matrix
C     **************
C
      IF (ITGTMT.EQ.1) THEN
C
C       Write instructions to compute the tangent matrix TGTMAT
C       In the present stage of the programme TGTMAT remains 
C       identical to the elastic stiffness matrix 
C
      ENDIF
C
C     Case of total destruction - all stress components are finally set to zero
C     *************************************************************************
C
C     Remark: activate the instructions with "CC" below when the tangent matrix 
C     is programmed
CC      IF (IRUPT.EQ.1) THEN
C       The damage parameter is set to unity
CC        QFG1G = 1.D0
C       All components of the stress tensor are set to zero		
CC        DO I=1,6 
CC          SIG1(I) = 0.D0
CC        ENDDO 
C       The evolutions of internal parameters are disregarded
CC        CALL TRNSFR(INTPAR0,INTPAR1,8)
CC      ENDIF
C
      RETURN
      END
C
C
C
C
C
C	
      SUBROUTINE STRMDE(SIG,SIGM,SIGDEV,SIGEQ)
C
C-----------------------------------------------------------------------
C
C     Purpose : Calculation of the mean part, deviatoric part, and  
C               "von Mises norm" (equivalent stress) of a given stress tensor
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I  REAL*8  SIG      Stress tensor in Voigt's "stress vector" form
C     O  REAL*8  SIGM     Mean part (1/3 of the trace) of SIG
C     O  REAL*8  SIGDEV   Deviatoric part of SIG
C     O  REAL*8  SIGEQ    "Von Mises norm" of SIG (equivalent stress)
C
C  ---------------------------------------------------------------------
C
      IMPLICIT NONE
C
C     Arguments
      REAL*8 SIG(6),SIGM,SIGDEV(6),SIGEQ
C
C     Local variables
      INTEGER I
C
      SIGM = (SIG(1)+SIG(2)+SIG(3))/3.D0
C
      DO I=1,3
        SIGDEV(I)   = SIG(I) - SIGM
        SIGDEV(I+3) = SIG(I+3)		
      ENDDO
C
      SIGEQ = SIGDEV(1)**2 + SIGDEV(2)**2 + SIGDEV(3)**2 
     .        + 2.D0*(SIGDEV(4)**2+SIGDEV(5)**2+SIGDEV(6)**2)
      SIGEQ = SQRT(3.D0/2.D0*SIGEQ)    
C
      RETURN
      END	
C
C
C
C
C
C	
      SUBROUTINE TRNSFR(A,B,N)
C 
C-----------------------------------------------------------------------
C
C     Purpose: Transfer of a vector into another one 
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I   REAL*8    A        Vector to be transferred
C     I   INTEGER   N        Dimension of A and B
C     O   REAL*8    B        Vector into which A is to be transferred
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C  
C     Arguments
      REAL*8 A(N),B(N)
      INTEGER N
C
C     Local variables
      INTEGER I 
C	  
      DO I=1,N
        B(I) = A(I)
      ENDDO
C	  	  
      RETURN
      END	
C
C
C
C
C
C	
      SUBROUTINE VGTTNS(AV,AT,ICODE)
C
C-----------------------------------------------------------------------
C
C     Purpose : Transformation of a 3X3 symmetric matrix stored in Voigt's 
C               "stress vector" form to tensorial form or vice-versa
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I/O  REAL*8  AV    Symmetric matrix in Voigt's "stress vector" form
C     I/O  REAL*8  AT    Symmetric matrix in tensorial form
C     I    INTEGER ICODE Code:
C                          = 1 : transformation from Voigt's form to tensorial 
C                                form
C                          = 2 : transformation from tensorial form to Voigt's 
C                                form
C
C  ----------------------------------------------------------------------
C
      IMPLICIT NONE
C
C     Arguments
      REAL*8 AV(6),AT(3,3)
      INTEGER ICODE
C
      IF (ICODE.EQ.1) THEN
        AT(1,1) = AV(1)
        AT(2,2) = AV(2)
        AT(3,3) = AV(3)
        AT(1,2) = AV(4)
        AT(2,1) = AT(1,2)
        AT(1,3) = AV(5)
        AT(3,1) = AT(1,3)
        AT(2,3) = AV(6)
        AT(3,2) = AT(2,3)
      ELSEIF (ICODE.EQ.2) THEN
        AV(1) = AT(1,1)
        AV(2) = AT(2,2)
        AV(3) = AT(3,3)
        AV(4) = AT(1,2)
        AV(5) = AT(1,3)
        AV(6) = AT(2,3)		
      ENDIF
C	  
      RETURN
      END	
C
C
C
C
C
C		  
      SUBROUTINE YLDSTR(STRAIN,STRESS,N,EPSB,SIGB,IERR)
C 
C-----------------------------------------------------------------------
C
C     Purpose: Calculation of the current yield stress (isotropic hardening)
C              by linear interpolation of the known stress-strain curve
C              in simple tension 
C
C-----------------------------------------------------------------------
C
C     Arguments (I=Input, O=Output):
C
C     I   REAL*8    STRAIN   Vector of values of the equivalent strain 
C                            (material data)
C     I   REAL*8    STRESS   Vector of corresponding values of the yield stress
C                            (material data)
C     I   INTEGER   N        Common dimension of the vectors STRAIN and STRESS
C     I   REAL*8    EPSB     Value of EPSILON_BAR
C     O   REAL*8    SIGB     Corresponding value of SIGMA_BAR      
C     O   INTEGER   IERR     Error code: 
C                            =1 if the first value of the equivalent strain
C                            in the vector STRAIN is nonzero
C                            =2 if the first value of the yield stress in the
C                            vector STRESS is zero
C
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C  
C     Arguments
      REAL*8 STRAIN(N),STRESS(N),EPSB,SIGB
      INTEGER N,IERR
C
C     Local variables
      REAL*8 COEF
      INTEGER I
C	  
      IERR = 0
C
      IF (STRAIN(1).NE.0.D0) THEN
        IERR = 1
        RETURN		
      ENDIF
      IF (STRESS(1).EQ.0.D0) THEN
        IERR = 2
        RETURN		
      ENDIF
C
C     The yield stress is considered as constant outside the 
C     interval of definition [STRAIN(1),STRAIN(N)]
      IF (EPSB.LE.STRAIN(1)) THEN
        SIGB = STRESS(1)
      ELSEIF (EPSB.GE.STRAIN(N)) THEN
        SIGB = STRESS(N)
      ELSE
        DO I=1,N-1
          IF (EPSB.LT.STRAIN(I+1)) GOTO 10
        ENDDO
   10  	COEF = ( EPSB-STRAIN(I) )/( STRAIN(I+1)-STRAIN(I) )
        SIGB = (1.D0-COEF)*STRESS(I) + COEF*STRESS(I+1)   
      ENDIF
C
      RETURN
      END