C TUFIEN    SOURCE    CHAT      05/01/13    03:55:36     5004           
      SUBROUTINE TUFIEN(XM,XP,DELTAM,DELTAP,THETA,XJP,PRECIS,THETA0,
     &     XM0,XP0,XJ1C,T,RAYOM,EPAI,YOUN)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C=======================================================================
C          CE SOUS PROGRAMME SERT A CALCULER LE NOUVEL ANGLE DE LA
C          FISSURE DANS LE CAS D UN INCREMENT DE CHARGE POUR
C          LEQUEL IL Y A PROPAGATION EN ELASTICITE
C          SI ON PLASTIFIE ON RENVOIT LA PARTIE DE L INCREMENT
C          QU IL VA FALLOIR ECOULER
C          IL EST APPELE PAR TUFPLA
C     
C          ENTREES:XM     VALEUR DU MOMENT AU DEBUT DU PAS
C                  XP     VALEUR DE L EFFORT NORMAL AU DEBUT DU PAS
C                  DELTAM  INCREMENT DE MOMENT
C                  DELTAP  INCREMENT D EFFORT NORMAL
C                  THETA   ANGLE DE LA FISSURE AU DEBUT DU PAS
C                  XJP    VALEUR DE LA PARTIE PLASTIQUE DE J
C     
C          SORTIES:XM     VALEUR DU MOMENT A LA FIN DU PAS
C                  XP     VALEUR DE L EFFORT NORMAL A LA FIN DU PAS
C                  DELTAM  PARTIE DE L INCREMENT DE MOMENT RESTANT
C                          A ECOULER EN PLASTICITE
C                  DELTAP  PARTIE DE L INCREMENT D EFFORT NORMAL
C                          RESTANT A ECOULER EN PLASTICITE
C                  THETA   ANGLE DE LA FISSURE A LA FIN DU PAS
C=======================================================================
-INC CCREEL
-INC CCOPTIO
      DIMENSION XX(2,2),YY(2,2)
      XM1=XM+DELTAM
      XP1=XP+DELTAP
C     
C     ON CALCULE LE NOUVEL ANGLE DE LA FISSURE PAR LA METHODE
C     DE NEWTON EN APPROXIMANT LA DERIVEE
C     
      X0=THETA
      X1=THETA+XPI/180.D0
      NBIT=0
 1    CONTINUE
      F0=CRIT2(XM1,XP1,X0,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN)
      F1=CRIT2(XM1,XP1,X1,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN)
      FP=(F1-F0)/(X1-X0)
      X2=X1-F1/FP
      IF(X2.LT.THETA0) THEN
         IF(IIMPI.EQ.999)THEN
            WRITE(IOIMP,111)
 111        FORMAT(1X,'L ANGLE EST INFERIEUR A L ANGLE INITIAL')
         ENDIF
         DELTAM=0.D0
         DELTAP=0.D0
         RETURN
      ENDIF
      IF(X2.GT.360.D0) THEN
         IF(IIMPI.EQ.999)THEN
            WRITE(IOIMP,222)
 222        FORMAT(1X,'PAS DE SOLUTION REALISTE')
         ENDIF
         DELTAM=0.D0
         DELTAP=0.D0
         RETURN
      ENDIF
      NBIT=NBIT+1
      CRI=ABS((X2-X1)/X1)
      F2=CRIT2(XM1,XP1,X2,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN)
      CRU=ABS(F2)
      IF(CRI.GE.PRECIS.OR.CRU.GE.PRECIS) THEN
         X0=X1
         X1=X2
         GOTO 1
      ELSE
         THETA=X2
         IF(IIMPI.EQ.999)THEN
            WRITE(IOIMP,*)THETA
            WRITE(IOIMP,*)'NBRE ITERATIONS',NBIT
            VER=CRIT2(XM1,XP1,THETA,XJP,XJ1C,T,RAYOM,THETA0,EPAI,YOUN)
            WRITE(IOIMP,900)VER
 900        FORMAT(1X,'VERIF RESOLUTION',1X,F12.5)
         ENDIF
C     
C     THETA EST LE NOUVEL ANGLE DE LA FISSURE
C     
      ENDIF
C     
C     ON TESTE POUR SAVOIR SI ON EST TOUJOURS ELASTIQUE
C     
      CALL TUFIC1(XM1,XP1,THETA,IRETOU,PRECIS,XM0,XP0)
C     
C     SI OUI ON A TOUT ECOULE ET ON MET DELTAM ET DELTAP
C     A ZERO
C     
      IF(IRETOU.EQ.0) THEN
         XM=XM+DELTAM
         XP=XP+DELTAP
         DELTAM=0.D0
         DELTAP=0.D0
         RETURN
C     
C      SINON ON CHERCHE L ANGLE QUI CORRESPOND A LA FRONTIERE
C      DU DOMAINE D ELASTICITE
C      ON RESOUD UN SYSTEME DU TYPE F(A,B)=G(A,B)=0 D UNE
C      MANIERE ANALOGUE A LA METHODE DE NEWTON
C     
      ELSE
         IF(IIMPI.EQ.999)THEN
            WRITE(IOIMP,901)
 901        FORMAT(1X,'MAIS ON PLASTIFIE')
         ENDIF
         NBIT=0
         A1=0.D0
         B1=1.D0
         XM1=XM+A1*DELTAM
         XP1=XP+A1*DELTAP
         THETA1=B1*THETA
         A2=1.D0
         B2=0.9D0
         XM2=XM+A2*DELTAM
         XP2=XP+A2*DELTAP
         THETA2=B2*THETA
 2       CONTINUE
         F11=CRIT1(XM1,XP1,THETA1,XM0,XP0)
         F21=CRIT1(XM2,XP2,THETA1,XM0,XP0)
         XX(1,1)=(F21-F11)/(A2-A1)
         F12=CRIT1(XM1,XP1,THETA2,XM0,XP0)
         XX(1,2)=(F12-F11)/(B2-B1)
         F22=CRIT1(XM2,XP2,THETA2,XM0,XP0)
         G11=CRIT2(XM1,XP1,THETA1,XJP,XJ1C,T,RAYOM,THETA0,EPAI,
     &        YOUN)
         G21=CRIT2(XM2,XP2,THETA1,XJP,XJ1C,T,RAYOM,THETA0,EPAI,
     &        YOUN)
         XX(2,1)=(G21-G11)/(A2-A1)
         G12=CRIT2(XM1,XP1,THETA2,XJP,XJ1C,T,RAYOM,THETA0,EPAI,
     &        YOUN)
         XX(2,2)=(G12-G11)/(B2-B1)
         G22=CRIT2(XM2,XP2,THETA2,XJP,XJ1C,T,RAYOM,THETA0,EPAI,
     &        YOUN)
         DET=XX(1,1)*XX(2,2)-XX(1,2)*XX(2,1)
         YY(1,1)=XX(2,2)/DET
         YY(1,2)=-XX(1,2)/DET
         YY(2,1)=-XX(2,1)/DET
         YY(2,2)=XX(1,1)/DET
         A3=A2-YY(1,1)*F22-YY(1,2)*G22
         B3=B2-YY(2,1)*F22-YY(2,2)*G22
         XM3=XM+A3*DELTAM
         XP3=XP+A3*DELTAP
         THETA3=B3*THETA
         NBIT=NBIT+1
         CRA=ABS((A3-A2)/A2)
         CRO=ABS((B3-B2)/B2)
         CRU=CRA+CRO
         F3=CRIT1(XM3,XP3,THETA3,XM0,XP0)
         G3=CRIT2(XM3,XP3,THETA3,XJP,XJ1C,T,RAYOM,THETA0,EPAI,
     &        YOUN)
         CRI=(F3**2+G3**2)**0.5D0
         IF(CRI.GE.PRECIS.OR.CRU.GE.PRECIS) THEN
            A1=A2
            B1=B2
            A2=A3
            B2=B3
            XM1=XM2
            XP1=XP2
            THETA1=THETA2
            XM2=XM3
            XP2=XP3
            THETA2=THETA3
            GOTO 2
         ELSE
C     
C    ON A ENFIN TROUVE A ET B QUI CONVIENNENT
C     
            A=A3
            B=B3
         ENDIF
         XM=XM+A*DELTAM
         XP=XP+A*DELTAP
         DELTAM=(1.D0-A)*DELTAM
         DELTAP=(1.D0-A)*DELTAP
         THETA=B*THETA
         VER1=CRIT1(XM,XP,THETA,XM0,XP0)
         VER2=CRIT2(XM,XP,THETA,XJP,XJ1C,T,RAYOM,THETA0,EPAI,
     &        YOUN)
         IF(IIMPI.EQ.999)THEN
            WRITE(IOIMP,903)VER1,VER2
 903        FORMAT(1X,'REVERIF RESO',1X,2F12.5)
            WRITE(IOIMP,*)'NBRE ITERATIONS',NBIT
         ENDIF
         RETURN
      ENDIF
      END