Sources Esope | Cast3M

coul3
C COUL3     SOURCE    CHAT      05/01/12    22:24:13     5004
      SUBROUTINE COUL3(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI,
     .     DEPST,IFOUR,XMAT,NMATT,IVAL,DD,SIGF,DEFP,VARF,KERRE)
C--------------------------------------------------------------
C
C         PLASTICITE MODELE 3D MOHR COULOMB ELEMENTS JOINTS
C
C ENTREES
C     IB              = NUMERO DE L'ELEMENT
C     IGAUSS          = NUMERO DU POINT DE GAUSS
C     NSTRS           = NOMBRE DE COMPOSANTES DE CONTRAINTES
C     SIG0            = CONTRAINTES INITIALES (AU PAS PRECEDENT)
C     EPIN0           = DEFORMATIONS INITIALES (AU PAS PRECEDENT)
C     VAR0(NVARI)     = VARIABLES INTERNES DEBUT
C      VAR0(1)        = DEFORMATION PLASTIQUE EQUIVALENTE
C      VAR0(2)        = EPOUN (DEFO. PLASTIQUE DUE A L'OUVERTURE SEULE)
C      VAR0(3)        = STAT (ETAT DU JOINT)
C      VAR0(4)        = LAM1 ( MULTIPLICATEUR PLASTIQUE EN TRACTION )
C     DSIGT(NSTRS)    = INCREMENT DE CONTRAINTES TOTALES
C     XMAT(NCOMAT)    = COMPOSANTES DE MATERIAU
C     NMATT           = NOMBRE DE COMPOSANTES DE MATERIAU
C     IVAL(NCOMAT)    = INDICE DES COMPOSANTES DE MATERIAU
C     DD(NSTRS,NSTRS) = MATRICE DE HOOK
C
C  SORTIES
C     SIGF(NSTRS)     = CONTRAINTES FINALES
C     DEFP(NSTRS)     = INCREMENTS DE DEFORMATIONS INELASTIQUES FINALES
C     VARF(NVARI)     = VARIABLES INTERNES FINALES
C     KERRE           = INDICE D'ERREUR
C                     = 60 PAS DE CONVERGENCE APRES 10 ITERATIONS
C                          INTERNES (VARIABLE ISSU DANS COULM)
C                     = 61
C
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
C
      DIMENSION SIG0(*),EPIN0(*),VAR0(*),DEPST(*),XMAT(*)
      DIMENSION DD(NSTRS,NSTRS)
      DIMENSION SIGF(*),DEFP(*),VARF(*)
C
C
C ----------PARAMETRES POUR LE TEST DE NULLITE DES CONTRAINTES ET
C ----------DEFORMATIONS
C
      PRESP = +1.D-5
      PRESM = -1.D-5
      PREEP = +1.D-10
      PREEM = -1.D-10
C
C
C ----------PARAMETRE D'IMPPRESSION
C
      IMPRI = 0
C
C ----------NOMBRE MAXIMUM D'ITERATIONS INTERNES AUTORISEES
C
      MAXITE = 1000
C
C
C------CHARGEMENT DU EPSILON DE PRECISION DES ITERATIONS INTERNES
C
      PRECIS = 1.D-6
      PRECIN = -1.D-6
C
C------CHARGEMENT DES CARACTERISTIQUES MECANIQUES
C
      G1     = XMAT(1)
      G2     = XMAT(1)
      EC     = XMAT(2)
      EF     = XMAT(5)
      ECN    = XMAT(6)
      C      = XMAT(7)
      PHI    = XMAT(8)
      PHI    = (PHI*XPI)/180.D0
ccc          ZMU    = XMAT(8)
      ZMU    = 0.0D0
      ZMU    = (ZMU*XPI)/180.D0
      FTRC  = XMAT(4)
C  AM 28/7/95
C  FTRC EST EN 4 ET NON EN 9 A CAUSE DE RHO ET ALPH
C  VOIR COMMENTAIRE DANS ECOUL2
C
C   LA LIMITE EN TRACTION EST MISE BRUTALEMENT A 0
C   DES QU'ELLE EST DEPASSEE.
C
C
C------CHARGEMENT DE LA MATRICE DE HOOK INVERSE
C------POUR LE CALCUL DE L'INCREMENT DE DEFORMATION
C------A PARTIR DE L'INCREMENT DE CONTRAINTE
C------         DEPSN = D**-1 DSIGT
C                 ( D**-1 = DD )
C
      DD(1,1) = 1.0D0/G1
      DD(1,2) = 0.0D0
      DD(1,3) = 0.0D0
      DD(2,1) = 0.0D0
      DD(2,2) = 1.0D0/G2
      DD(2,3) = 0.0D0
      DD(3,1) = 0.0D0
      DD(3,2) = 0.0D0
      DD(3,3) = 1.0D0/EC
C
C           contrainte en fin de pente ec
C
      sigsei = -ecn * ec
c
C
C-----------INCREMENT DE DEFORMATION
C
      DEPS1 = DEPST(1)
      DEPS2 = DEPST(2)
      DEPSN = DEPST(3)
C
C             TRAITEMENT DU PASSAGE DIRECT ENTRE UN ETAT OUVERT
C             (I.E. EXISTENCE D'UN JEU AU TEMPS N) A UN ETAT
C             COMPRIME AU TEMPS N+1
C
      DEFPN=0.0D0
      INDIC=0
      IFIN =0
      STOTA1 = 0.0D0
      STOTA2 = 0.0D0
      DEFP11 = 0.0D0
      DEFP12 = 0.0D0
      DEFP21 = 0.0D0
      DEFP22 = 0.0D0
      ALP    = 1.0D0
      S1INT  = 0.0D0
      S2INT  = 0.0D0
      ZJEU = EPIN0(3)
C
C  AM  CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
C                                CAS DU RAPPROCHEMENT DES LEVRES
C        IF (DEPSN.LT.-1.D-4) THEN
      IF (DEPSN.LT.-1.D-4.AND.SIG0(3).GT.PRESP) THEN
C                                IL EXISTE UN JEU
         IF (ZJEU.GT.1.D-3) THEN
C                                PASSAGE DIRECT D'UN ETAT OUVERT
C                                A UN ETAT COMPRIME
cc               IF (ABS(DEPSN).GT.ZJEU) THEN
            valdif=(ABS(DEPSN) - ZJEU)/ABS(DEPSN)
ccc               IF (valdif.gt.3.05D-6) THEN
            IF (valdif.gt.1.D-3) THEN
C
               ALP = ZJEU/(ABS(DEPSN))
               STOTA1= DEPS1
               DEPS1 = ALP * DEPS1
               STOTA2= DEPS2
               DEPS2 = ALP * DEPS2
               STOCK = DEPSN
               DEPSN = ALP * DEPSN
               INDIC = 1
               GO TO 99
            ENDIF
C
         ENDIF
C
      ENDIF
C
 98   CONTINUE
C
C                                IL RESTE A PARCOURIR LE TRAJET
C                                (1-ALP) * DEPSN
C
      IF (INDIC.EQ.1) THEN
         DEPS1=(1.0D0-ALP)*STOTA1
         DEPS2=(1.0D0-ALP)*STOTA2
         DEPSN=(1.0D0-ALP)*STOCK
         INDIC=0
         IFIN =1
      ENDIF
C
 99   CONTINUE
C
C
C ----CAS GENERAL SANS TRANSITION DIRECTE D'UN ETAT OUVERT A UN ETAT
C ----A UN ETAT COMPRIME
C
C
C
C
      IF ((INDIC.EQ.1).OR.(IFIN.EQ.0)) THEN
C
C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS
C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT
C
C------CHARGEMENT DES VALEURS AU PAS PRECEDENT
C
C
C----------DEFORMATIONS
C
C             DEFORMATIONS INELASTIQUES
C
         EPSS1 = EPIN0(1)
         EPSS2 = EPIN0(2)
         EPSN1 = EPIN0(3)
C
C             DEFORMATIONS ELASTIQUES
C
         DEFS10 = DD(1,1)*SIG0(1)
         DEFS20 = DD(2,2)*SIG0(2)
c
         if   (sig0(3).lt.sigsei) then
            DEFN0 = ( -1.0d0*ecn ) + ( (sig0(3) - sigsei) / ef )
         else if ( (sig0(3).ge.sigsei) .and. (sig0(3).le.presp) ) then
            DEFN0 = DD(3,3)*SIG0(3)
         else if (sig0(3).gt.presp) then
C
C  AM  CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
C
C           DEFN0 = 0.0D0
            DEFN0 = DD(3,3)*SIG0(3)
         endif
C
C             DEFORMATIONS TOTALES AU PAS PRECEDENT N
C
         EPS1N  = EPSS1 + DEFS10
         EPS2N  = EPSS2 + DEFS20
         EPSNN  = EPSN1 + DEFN0
C
      ENDIF
C
C
C
C
C ----CAS DU 2ND PASSAGE : D'UN ETAT FERME ON PASSE A UN ETAT COMPRIME
C
C
C
C
      IF (IFIN.EQ.1) THEN
C
C--------------------CALCUL DE L ETAT DE CONTRAINTES, DE DEFORMATIONS
C--------------------DES VARIABLES INTERNES AU PAS PRECEDENT
C
C------CHARGEMENT DES VALEURS AU PAS PRECEDENT
C
C----------DEFORMATIONS
C
C             DEFORMATIONS INELASTIQUES
C
         EPSS1 = DEFP11
         EPSS2 = DEFP21
         EPSN1 = 0.0D0
C
C             DEFORMATIONS ELASTIQUES
C
         DEFS10 = DD(1,1)*S1INT
         DEFS20 = DD(2,2)*S2INT
         DEFN0  = 0.0D0
C
C             DEFORMATIONS TOTALES AU PAS PRECEDENT N
C
         EPS1N  = EPSS1 + DEFS10
         EPS2N  = EPSS2 + DEFS20
         EPSNN  = EPSN1 + DEFN0
C
      ENDIF
C
C             ESTIMATION DE LA DEFO. TOTALE AU PAS N+1
C
      EP1NP1 = EPS1N + DEPS1
      EP2NP1 = EPS2N + DEPS2
      EPNNP1 = EPSNN + DEPSN
C
C----------EPSILON ELASTIQUE (1 2 N) AU PAS N+1 A L'ITERATION 0
C
      EE1P10 = DEFS10 + DEPS1
      EE2P10 = DEFS20 + DEPS2
      EENP10 = DEFN0  + DEPSN
C
C----------SIGMA (1 2 N) AU PAS N+1 A L'ITERATION 0
C
      S1NP10 = G1 * EE1P10
      S2NP10 = G2 * EE2P10
c
      x6 = -1.0d0 * ecn
c
C  AM  CORRECTION POUR JOINT AVEC RESISTANCE EN TRACTION
C
C        if  (depsn.lt.preem) then
      IF  (DEPSN.LT.PREEM.AND.SIG0(3).LE.PRESP) THEN
c
c -----------fermeture du joint
c
         if (epnnp1.le.preep) then
c
c --------------fermeture du joint avec compression des levres
c
            if ( (eenp10.ge.x6) .and. (eenp10.le.preep) ) then
c
c ----------------compression : pente EC
c
               SNNP10 = XMAT(2) * EENP10
            endif
c
            if   (eenp10.lt.x6) then
c
c ----------------compression : pente Ef
c
               snnp10= ( xmat(2) * x6 ) +
     .              ( (eenp10 + xmat(6)) * xmat(5) )
            endif
c
         endif
c
         if ((epnnp1.gt.preep).or.(indic.eq.1)) then
c
c --------------rapprochement des levres sans compression
c
            SNNP10 = 0.0D0
c
         endif
c
      ELSE
c
c -----------ouverture du joint OU FERMETURE SI JOINT PAS CASSE
c
         SNNP10 = XMAT(2) * EENP10
c
      endif
C
C-----------VARIABLES INTERNES
C
C          DEFORMATION PLASTIQUE EQUIVALENTE
C
      EPSEQ=VAR0(1)
C
C          DEFORMATION PLASTIQUE DUE A L'OUVERTURE SEULE
C
      EPOUN=VAR0(2)
C
C          ETAT DU JOINT
C
C            STAT   = 0 LEVRES EN CONTACT, AVEC OU SANS CISAILLEMENT
C            STAT   = 1 JOINT EN OUVERTURE PENTE EC
C            STAT   = 2 JOINT EN FERMETURE PENTE EC
C            STAT   = 3 JOINT EN FERMETURE PENTE EF
C
C
C
      STAT =VAR0(3)
C
C           MULTIPLICATEUR PLASTIQUE EN TRACTION
C
      XLAM1=VAR0(4)
C
C      MISE A 0 DE LA LIMITE EN TRACTION SI DEJA CASSE OU SI JEU
C
      RTRAC=FTRC
      IF(XLAM1.GT.0.D0) RTRAC=0.D0
 
C
C
C-----------ETUDE DE L'ECOULEMENT SELON LES DIFFERENTS CRITERES
C
      CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
      CRINOR = SNNP10 - RTRAC
      CRICIP = CALN(S1NP10,S2NP10,SNNP10,PHI,C)
      ABSTAU = SQRT( (S1NP10 ** 2) + (S2NP10 ** 2) )
C
C----------AIGUILLAGE SUIVANT LA ZONE 0, 1, 2, 3
C
      IF ((CRICIS.LE.PRECIS).AND.(CRINOR.LE.PRECIS)) THEN
         GO TO 10
      ENDIF
C
      IF ((CRICIS.GT.PRECIS).AND.(CRICIP.LE.PRECIS)) THEN
         GO TO 40
      ENDIF
C
      IF ((CRICIP.GT.PRECIS).AND.(ABSTAU.GT.C)) THEN
         GO TO 20
      ENDIF
C
      IF ((ABSTAU.LE.C).AND.(CRINOR.GT.PRECIN)) THEN
         GO TO 30
      ENDIF
C
C-----------CAS DE LA CONVERGENCE A L ITERATION 0 :
C
 10   CONTINUE
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
      IF (IMPRI.EQ.1) THEN
         WRITE(6,*) 'PASSAGE N. 1'
         WRITE(6,*) 'LES DEUX CRITERES SONT VERIFIES'
         WRITE(6,*) 'CONVERGENCE A L ITERATION 0'
      ENDIF
C
      IF (INDIC.EQ.1) THEN
C                              FERMETURE DU JOINT SANS COMPRESSION
         S1INT = S1NP10
         DEFP11= 0.0D0
         S2INT = S2NP10
         DEFP21= 0.0D0
         DEFPN = -epin0(3)
         GO TO 98
C
      ELSE
C
         IF (IFIN.EQ.0) THEN
C
C
c
            defp(1)=0.0d0
            defp(2)=0.0d0
C
            if ( abs(snnp10).le.presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
               defp(3) = epnnp1 - epsn1
c
            else
c
c            le joint se deforme toujours elastiquement en fermeture
c
               DEFP(3) = 0.0D0
c
            endif
c
         ELSE
c
c             cas du passage direct d'un etat ouvert a un etat
c             comprime
c
            DEFP12  = 0.0D0
            DEFP(1) = 0.0d0
            DEFP22  = 0.0D0
            DEFP(2) = 0.0d0
cc                DEFP(3) = 0.0d0
            DEFP(3) = defpn
            IFIN = 0
         ENDIF
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delep2 = defp(2)
         delepn = defp(3)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp10).le.presp ) then
            epoun = defp(3)
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint  stat = 0   cisaillement pur, levres en contact
c                    stat = 1   ouverture du joint pente EC
c                    stat = 2   compression du joint pente EC
c                    stat = 3   compression du joint pente EF
c
         if ( abs(snnp10).le.presp) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-3 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp10.gt.sigsei).and.(snnp10.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp10.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
C
         SIGF(1) = S1NP10
         SIGF(2) = S2NP10
         SIGF(3) = SNNP10
C
         RETURN
c
      ENDIF
C
C----------CRITERES DE L'EFFORT NORMAL ET DE CISAILLEMENT
C----------TOUS LES DEUX VIOLES EN MEME TEMPS
C
 20   CONTINUE
C
C            ECOULEMENT SELON LES DEUX CRITERES A L'ITERATION ITER
C
      ITER = 1
      DELAM1=0.D0
C
 2    CONTINUE
C
      IF (ITER.GT.MAXITE) THEN
         WRITE(6,*) ' LES DEUX CRITERES SONT VIOLES '
         WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
         STOP
      ENDIF
C
      IF (IMPRI.EQ.1) THEN
         write(6,*) 'PASSAGE N. 2'
         write(6,*) 'LES DEUX CRITERES SONT VIOLES'
         write(6,*) 'ITERATION',ITER
      ENDIF
C
C
C            CALCUL DES CORRECTEURS PLASTIQUES
C
C
C                DERIVEE DU CRITERE F SUIVANT SIGMA
C
      DFIDS1 = 0.0D0
      DFIDS2 = 0.0D0
      DFIDSN = 1.0D0
C
C                DERIVEE DU CRITERE G SUIVANT SIGMA
C
      DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
      DGIDS1 = S1NP10 / DENO
      DGIDS2 = S2NP10 / DENO
      DGIDSN = TAN(PHI)
      DZGIDS = TAN(ZMU)
C
C                MATRICE DE HOOK * DFI/DSIGMA
C
      HODF1 = G1 * DFIDS1
      HODF2 = G2 * DFIDS2
      HODFN = EC * DFIDSN
C
C                MATRICE DE HOOK * DGI/DSIGMA
C
      HODG1 = G1 * DGIDS1
      HODG2 = G2 * DGIDS2
      HODGN = EC * DZGIDS
C
C                FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      FIFI = (DFIDS1 * HODF1) +  (DFIDS2 * HODF2) + (DFIDSN * HODFN)
C
C                FIGI = DFI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
C
      FIGI = (DFIDS1 * HODG1) +  (DFIDS2 * HODG2) + (DFIDSN * HODGN)
C
C                GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      GIFI = (DGIDS1 * HODF1) +  (DGIDS2 * HODF2) + (DGIDSN * HODFN)
C
C                GIFI = DGI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      GIGI = (DGIDS1 * HODG1) +  (DGIDS2 * HODG2) + (DGIDSN * HODGN)
C
C                RESOLUTION DU SYSTEME POUR TROUVER LES CORRECTEURS
C                PLASTIQUES DLAM1 ET DLAM2
C
      DET = (FIFI * GIGI) - (FIGI * GIFI)
C
      IF (ABS(DET).LE.PRECIS) THEN
         WRITE(6,*) ' SYSTEME INDETERMINE (F G): ITERATION ',ITER
         STOP
      ENDIF
C
      DLAM1 = (CRINOR * GIGI) - (CRICIS * FIGI)
      DLAM1 = DLAM1 / DET
      DELAM1 = DELAM1 + DLAM1
      RTRAC=0.D0
C
      DLAM2 = (CRICIS * FIFI) - (CRINOR * GIFI)
      DLAM2 = DLAM2 / DET
C
C                CALCUL DU NOUVEL ETAT DE CONTRAINTES
C
      S1NP11=S1NP10 - (G1 * ( (DLAM1*DFIDS1) + (DLAM2*DGIDS1)))
      S2NP11=S2NP10 - (G2 * ( (DLAM1*DFIDS2) + (DLAM2*DGIDS2)))
      SNNP11=SNNP10 - (EC * ( (DLAM1*DFIDSN) + (DLAM2*DZGIDS)))
C
C                ALLER EN 52 POUR LE TEST DE LA CONVERGENCE
C
C
 52   CONTINUE
C
C                CONVERGENCE CHECK
C
      CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
      CRINOR = SNNP11 - RTRAC
C
      IF ((ABS(CRICIS).LE.PRECIS).AND.(ABS(CRINOR).LE.PRECIS)) THEN
C
C----------LES DEUX CRITERES SONT VERIFIES
C
C
         IF (IMPRI.EQ.1) THEN
            write(6,*) 'PASSAGE N. 5 - 2'
            write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
            write(6,*) 'CONVERGENCE A L ITERATION',ITER
         ENDIF
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
         SIGF(1) = S1NP11
         SIGF(2) = S2NP11
         SIGF(3) = SNNP11
C
         DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
         DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
c
         if ( abs(sigf(3)) .le. presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
            DEFP(3) = EPNNP1 - EPSN1
         else
c
c        le joint se deforme toujours elastiquement en fermeture
c
            DEFP(3) = 0.0D0
c
         endif
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delep2 = defp(2)
         delepn = defp(3)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1 + DELAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp11).le.presp ) then
            epoun = delepn
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint
c                stat = 0  levres en contact, avec ou sans cisaillement
c                stat = 1  ouverture du joint pente EC
c                stat = 2  compression du joint pente EC
c                stat = 3  compression du joint pente EF
c
         if ( abs(snnp11).le.presp ) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp11.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
c
         RETURN
      ELSE
C
         S1NP10=S1NP11
         S2NP10=S2NP11
         SNNP10=SNNP11
         CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
         CRINOR = SNNP10 - RTRAC
         ITER = ITER + 1
C
         GO TO 2
      ENDIF
C
C
C____________________________________________________________________
C
C
C----------CRITERE DE L'EFFORT NORMAL VIOLE MAIS
C----------CRITERE DE CISAILLEMENT VERIFIE
C
 30   CONTINUE
C
C            ECOULEMENT SELON LE CRITERE F A L'ITERATION ITER
C
      ITER = 1
      DELAM1=0.D0
C
 3    CONTINUE
C
      IF (ITER.GT.MAXITE) THEN
         WRITE(6,*) ' LE CRITERE NORMAL EST VIOLE '
         WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
         STOP
      ENDIF
C
      IF (IMPRI.EQ.1) THEN
         write(6,*) 'PASSAGE N. 3'
         write(6,*) 'CRITERE DE L EFFORT NORMAL VIOLE'
         write(6,*) 'ITERATION',ITER
      ENDIF
C
C
C
C            CALCUL DU CORRECTEUR PLASTIQUE
C
C
C                DERIVEE DU CRITERE F SUIVANT SIGMA
C
      DFIDS1 = 0.0D0
      DFIDS2 = 0.0D0
      DFIDSN = 1.0D0
C
C                MATRICE DE HOOK * DGI/DSIGMA
C
      HODF1 = G1 * DFIDS1
      HODF2 = G2 * DFIDS2
      HODFN = EC * DFIDSN
C
C                FIFI = DFI/DSIGMA * MATRICE DE HOOK * DFI/DSIGMA
C
      FIFI = (DFIDS1 * HODF1) +  (DFIDS2 * HODF2) + (DFIDSN * HODFN)
C
C                RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
C                PLASTIQUE DLAM1
C
      IF (ABS(FIFI).LE.PRECIS) THEN
         WRITE(6,*) ' SYSTEME INDETERMINE (F) : ITERATION ',ITER
         STOP
      ENDIF
C
      DLAM1 = CRINOR / FIFI
      DELAM1 = DELAM1 + DLAM1
      RTRAC = 0.D0
C
C                CALCUL DU NOUVEL ETAT DE CONTRAINTES
C
      S1NP11=S1NP10 - (XMAT(1) *  (DLAM1*DFIDS1) )
      S2NP11=S2NP10 - (XMAT(1) *  (DLAM1*DFIDS2) )
      SNNP11=SNNP10 - ( EC     *  (DLAM1*DFIDSN) )
C
C                ALLER EN 53 POUR LE TEST DE LA CONVERGENCE
C
C
C
 53   CONTINUE
C
C                CONVERGENCE CHECK
C
      CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
      CRINOR = SNNP11 - RTRAC
C
      IF (ABS(CRINOR).LE.PRECIS) THEN
C
C----------LES DEUX CRITERES SONT VERIFIES
C
C
         IF (IMPRI.EQ.1) THEN
            write(6,*) 'PASSAGE N. 5 - 3'
            write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
            write(6,*) 'CONVERGENCE A L ITERATION',ITER
         ENDIF
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
         SIGF(1) = S1NP11
         SIGF(2) = S2NP11
         SIGF(3) = SNNP11
C
         DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
         DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
c
         if ( abs(sigf(3)) .le. presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
            DEFP(3) = EPNNP1 - EPSN1
         else
c
c        le joint se deforme toujours elastiquement en fermeture
c
            DEFP(3) = 0.0D0
c
         endif
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delep2 = defp(2)
         delepn = defp(3)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1 + DELAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp11).le.presp ) then
            epoun = delepn
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint
c                stat = 0  levres en contact, avec ou sans cisaillement
c                stat = 1  ouverture du joint pente EC
c                stat = 2  compression du joint pente EC
c                stat = 3  compression du joint pente EF
c
         if ( abs(snnp11).le.presp ) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp11.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
c
         RETURN
      ELSE
C
         S1NP10=S1NP11
         S2NP10=S2NP11
         SNNP10=SNNP11
         CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
         CRINOR = SNNP10 - RTRAC
         ITER = ITER + 1
C
         GO TO 3
      ENDIF
C
C
C____________________________________________________________________
C
C
C
C----------CRITERES DE CISAILLEMENT VIOLE MAIS
C----------CRITERE DE L'EFFORT NORMAL VERIFIE
C
 40   CONTINUE
C
C            ECOULEMENT SELON LE CRITERE G A L'ITERATION ITER
C
      ITER = 1
C
 4    CONTINUE
C
      IF (ITER.GT.MAXITE) THEN
         WRITE(6,*) ' LE CRITERE DE CISAILLEMENT EST VIOLE '
         WRITE(6,*) ' DIVERGENCE APRES 1000 ITERATIONS INTERNES'
         STOP
      ENDIF
C
      IF (IMPRI.EQ.1) THEN
         write(6,*) 'PASSAGE N. 4'
         write(6,*) 'CRITERE DE CISAILLEMENT VIOLE'
         write(6,*) 'ITERATION',ITER
      ENDIF
C
C
C            CALCUL DU CORRECTEUR PLASTIQUE
C
C
C                DERIVEE DU CRITERE G SUIVANT SIGMA
C
      DENO = SQRT ( (S1NP10**2) + (S2NP10**2) )
      DGIDS1 = S1NP10 / DENO
      DGIDS2 = S2NP10 / DENO
      DZGIDS = TAN(ZMU)
      DGIDSN = TAN(PHI)
C
C                MATRICE DE HOOK * DGI/DSIGMA
C
      HODG1 = G1 * DGIDS1
      HODG2 = G2 * DGIDS2
      HODGN = EC * DZGIDS
C
C                GIGI = DGI/DSIGMA * MATRICE DE HOOK * DGI/DSIGMA
C
      GIGI = (DGIDS1 * HODG1) +  (DGIDS2 * HODG2) + (DGIDSN * HODGN)
C
C                RESOLUTION DE L'EQUATION POUR TROUVER LE CORRECTEUR
C                PLASTIQUE DLAM2
C
      IF (ABS(GIGI).LE.PRECIS) THEN
         WRITE(6,*) ' SYSTEME INDETERMINE (G) : ITERATION ',ITER
         STOP
      ENDIF
C
      DLAM2 = CRICIS / GIGI
C
C                CALCUL DU NOUVEL ETAT DE CONTRAINTES
C
cc       SNNP11=SNNP10 - (XMAT(2) *  (DLAM2*DGIDSN) )
C
      S1NP11=S1NP10 - (G1 *  (DLAM2*DGIDS1) )
      S2NP11=S2NP10 - (G2 *  (DLAM2*DGIDS2) )
      SNNP11=SNNP10 - (EC *  (DLAM2*DZGIDS) )
c
      if  (SNNP11.ge.0.0d0)  then
         SNNP11 = 0.0d0
      endif
c
cccssszzzeeerrr       if (SNNP11.GE.0.0d0) THEN
cccssszzzeeerrr             SNNP11 = 0.0d0
cccssszzzeeerrr       ENDIF
cccssszzzeeerrrc
cccssszzzeeerrr       if (abs(zmu).le.presp) then
cccssszzzeeerrr          differ = abs( depsn ) - epsn1
cccssszzzeeerrr          differ = abs(differ)/abs( depsn )
cccssszzzeeerrr          if (depsn.lt.preem) then
cccssszzzeeerrr             if (epnnp1.gt.preem) then
cccssszzzeeerrr                SNNP11 = 0.0d0
cccssszzzeeerrr             endif
cccssszzzeeerrr          else
cccssszzzeeerrrc
cccssszzzeeerrr             SNNP11 = 0.0d0
cccssszzzeeerrr          endif
cccssszzzeeerrr       endif
C
C                ALLER EN 54 POUR LE TEST DE LA CONVERGENCE
C
C
C
 54   CONTINUE
C
C                CONVERGENCE CHECK
C
      CRICIS = CALG(S1NP11,S2NP11,SNNP11,PHI,C)
      CRINOR = SNNP11 - RTRAC
C
      IF (ABS(CRICIS).LE.PRECIS) THEN
C
         IF (INDIC.EQ.1) THEN
C
C----------1ER TRAJET : ON PARCOURT ALP * DEPSN
C
            DEFPN = -epin0(3)
            S1INT = S1NP11
            S2INT = S2NP11
            DEFP11 = EP1NP1 - ( DD(1,1) * S1NP11 )
            DEFP21 = EP2NP1 - ( DD(2,2) * S2NP11 )
            GO TO 98
         ELSE
            IF (IFIN.EQ.0) THEN
C
C
C----------LES DEUX CRITERES SONT VERIFIES
C
C
               IF (IMPRI.EQ.1) THEN
                  write(6,*) 'PASSAGE N. 5 - 4'
                  write(6,*) 'LES DEUX CRITERES SONT VERIFIES'
                  write(6,*) 'CONVERGENCE A L ITERATION',ITER
               ENDIF
C
C----------CHARGEMENT DES VALEURS DE CONTRAINTES, DEFORMATIONS,
C----------VARIABLES INTERNES APRES CONVERGENCE INTERNE
C
C
               SIGF(1) = S1NP11
               SIGF(2) = S2NP11
               SIGF(3) = SNNP11
C
               DEFP(1) = ( EP1NP1 - ( DD(1,1) * SIGF(1) ) ) - EPSS1
               DEFP(2) = ( EP2NP1 - ( DD(2,2) * SIGF(2) ) ) - EPSS2
c
               if ( abs(sigf(3)) .le. presp ) then
c
c            deformation plastique du pas (n+1) moins
c            deformation plastique du pas (n)
c
                  DEFP(3) = EPNNP1 - EPSN1
               else
c
c        le joint se deforme toujours elastiquement en fermeture
c
                  DEFP(3) = 0.0D0
c
               endif
            ELSE
C
C----------CAS DU PASSAGE DIRECT D'UN ETAT OUVERT A UN ETAT COMPRIME
C
               DEFP12  = EP1NP1 - (DD(1,1)*S1NP11)
               DEFP22  = EP2NP1 - (DD(2,2)*S2NP11)
               DEFP(1) =  DEFP12 - epin0(1)
               DEFP(2) =  DEFP22 - epin0(2)
               DEFP(3) = -epin0(3)
               SIGF(1) = S1NP11
               SIGF(2) = S2NP11
               SIGF(3) = SNNP11
               IFIN    = 0
            ENDIF
C
         ENDIF
C
C     CALCUL DES VARIABLES INTERNES : DEPEQ  STAT  EPOUN
C
c
c      calcul du delta epsilon plastique entre
c      le pas actuel et le pas precedent
c
         delep1 = defp(1)
         delep2 = defp(2)
         delepn = defp(3)
c
c      produit contracte pour le calcul de la defo. plasti. cumulee
c
         pro1=(delep1 * delep1) + (delep2 * delep2) + (delepn * delepn)
         depeq = sqrt( (2.0D0/3.0D0) * pro1 )
 
         VARF(1) = EPSEQ + DEPEQ
         VARF(4) = XLAM1
c
c      calcul de la deformation normale due a l'ouverture seule
c
         if ( abs(snnp11).le.presp ) then
            epoun = delepn
            VARF(2) = EPOUN
         else
            epoun = 0.0d0
            VARF(2) = EPOUN
         endif
c
c     etat du joint
c                stat = 0  levres en contact, avec ou sans cisaillement
c                stat = 1  ouverture du joint pente EC
c                stat = 2  compression du joint pente EC
c                stat = 3  compression du joint pente EF
c
         if ( abs(snnp11).le.presp ) then
            stat = 1.D0
            if ( abs(epnnp1).le.1.D-8 ) then
               stat = 0.D0
            endif
            VARF(3) = STAT
         endif
         if ( (snnp11.gt.sigsei).and.(snnp11.lt.presm) ) then
            stat = 2.D0
            VARF(3) = STAT
         endif
         if (snnp11.le.sigsei) then
            stat = 3.D0
            VARF(3) = STAT
         endif
c
         RETURN
      ELSE
C
         S1NP10=S1NP11
         S2NP10=S2NP11
         SNNP10=SNNP11
         CRICIS = CALG(S1NP10,S2NP10,SNNP10,PHI,C)
         CRINOR = SNNP10 - RTRAC
         ITER = ITER + 1
C
         GO TO 4
      ENDIF
C
C
C____________________________________________________________________
C
      RETURN
      END