C T1T2DP    SOURCE    CHAT      05/01/13    03:29:37     5004
      SUBROUTINE T1T2DP(SIG,DSIG,YOUN,ANU,RT1,RT2,H1,H2,RDP,ADP,HDP,
     1 XLAM1,XLAM2,XLAM3,IDAM,ANG,EPPLDP,ALPHA,KERRE)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
      DIMENSION SIG(3),DSIG(3),SGG(3),DSGG(3),SFG(3),DSFG(3),
     1 A(3),B(3),IDAM(3),EPC(3),SSI(3),SI(3),EPC1(3),SG1(3),
     1 SG2(3),EC1(3),EC2(3),AT(4),BT(2),EPDP(3),EPPLDP(3)
      ITER=0
      Y=YOUN/(1.D0-ANU*ANU)
      DO 5 ITYP=1,3
      EPPLDP(ITYP)=0.D0
    5 EPDP(ITYP)=0.D0
      XT=0.D0
      IF(IIMPI.EQ.9) WRITE(IOIMP,9999)
9999  FORMAT(1X,'T1T2DP TRIPLAGE TRACTION 1 TRACTION 2 DRUCKER',/)
 9998 FORMAT(1X,'T1T2DP ',I4,'ITERATIONS INTERNES',/)
C
C------------------------------------------------
C     TRIPLAGE TRACTION TRACTION DRUCKER
C------------------------------------------------
C----------------------------------------------
C     DEFINITION D UNE CONTRAINTE DE REFERENCE
C     SI RMAX DU DRUCKER INFERIEURE A CETTE VALEUR
C     ALORS SIGMA=0
C----------------------------------------------
C
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      SIREF=1.E-6*YOUN
C
C--------------------------------------------------
C     CAS OU LE RAYON DU DRUCKER EST INFERIEUR
C     A UNE VALEUR DE REFERENCE
C--------------------------------------------------
C
      IF(SIREF.GT.RMAX) THEN
      DO 10 ITYP=1,3
  10  SIG(ITYP)=SIG(ITYP)+DSIG(ITYP)
      RDP=0.D0
      HDP=0.D0
      ADP=0.D0
      CALL NORME(SIG,DL3)
      DL3=DL3/YOUN*10.D0
      XLAM3=XLAM3+DL3
      CALL CPHOIN(SIG,EPPLDP,YOUN,ANU)
      DO 20 ITYP=1,3
      IDAM(ITYP)=0
      DSIG(ITYP)=0.D0
  20  SIG(ITYP)=0.D0
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      ENDIF
C-----------------------------------------------
C     ON SE PLACE DANS LE REPERE DE FISSURATION
C-----------------------------------------------
C
      CALL CHREP(SIG,SFG,ANG)
      CALL CHREP(DSIG,DSFG,ANG)
C
C-----------------------------------------------------
C     ESTIMATION DU PAS D'INCREMENT DE CONTRAINTES
C-----------------------------------------------------
C
      CALL SIJ(SFG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC)
      CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA)
      CALL SCAL2(EPC,SSI,VAL)
      VAL=VAL-HDP
      IF(VAL.LT.0.D0) THEN
      KERRE=459
      RETURN
      ENDIF
      CALL SCAL(SSI,DSFG,VAL)
      CALL NORME(DSFG,VA1)
      CALL NORME(SSI,VA2)
      IF(VA1.EQ.0.D0) THEN
      IDAM(1)=0
      IDAM(2)=0
      IDAM(3)=0
      RETURN
      ENDIF
      RMIN=MIN((RDP/1.73),(RDP/(1.D0+2.D0*ADP)))
      X=VAL/VA1/VA2
      IF(X.GE.1.D0) THEN
      X=1.D0
      GOTO 25
      ENDIF
      X=1.D0/SQRT(1.0001D0-X*X)*RMIN/VA1/8.D0
      IF(X.GT.1.D0) X=1.D0
C
C----------------------------------------------------------
C     ON ECOULE ET ON REGARDE LA VARIATION DE NORMALE
C----------------------------------------------------------
C
C
C---------------------------------------
C     RESOLUTION DU SYSTEME EN DL3
C---------------------------------------
   25 X=X/2.D0
      EC1(1)=1.D0
      EC1(2)=0.D0
      EC1(3)=0.D0
      EC2(1)=0.D0
      EC2(2)=1.D0
      EC2(3)=0.D0
      CALL CPHOOB(EC1,SG1,YOUN,ANU)
      CALL CPHOOB(EC2,SG2,YOUN,ANU)
   51 X=X*2.D0
   52 RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      ITER=ITER+1
  653 IF(ITER.GT.200) THEN
      KERRE=460
      RETURN
      ENDIF
      IF(SIREF.GT.RMAX) THEN
      RDP=0.D0
      ADP=0.D0
      HDP=0.D0
      DO 60 ITYP=1,3
   60 SFG(ITYP)=SFG(ITYP)+(1.D0-XT)*DSFG(ITYP)
      CALL NORME(SFG,DL3)
      DL3=DL3/YOUN*10.D0
      XLAM3=XLAM3+DL3
      CALL CPHOIN(SFG,EPPLDP,YOUN,ANU)
      DO 70 ITYP=1,3
   70 EPDP(ITYP)=EPDP(ITYP)+EPPLDP(ITYP)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      DO 80 ITYP=1,3
      IDAM(ITYP)=0
      DSIG(ITYP)=0.D0
   80 SIG(ITYP)=0.D0
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      ENDIF
      IF((XT+X).GT.1.D0) X=1.D0-XT
      CALL SIJ(SFG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC)
      CALL CPHOMO(EPC,SSI,YOUN,ANU,ALPHA)
      CALL SCAL2(EPC,SSI,VAL)
      VAL=VAL-HDP
      IF(VAL.LT.0.D0) THEN
      KERRE=459
      RETURN
      ENDIF
C
C
   90 ITER=ITER+1
      IF(ITER.GT.201) GOTO 653
      AT(1)=Y-H1
      AT(4)=Y-H2
      AT(2)=Y*ANU
      AT(3)=AT(2)
      DO 95 ITYP=1,3
   95 DSGG(ITYP)=X*DSFG(ITYP)
      BT(1)=SFG(1)+DSGG(1)-RT1
      BT(2)=SFG(2)+DSGG(2)-RT2
      CALL SYLIN2(AT,BT,DL11,DL21)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,DL12,DL22)
      DO 100 ITYP=1,3
      A(ITYP)=SFG(ITYP)+DSGG(ITYP)-DL11*SG1(ITYP)-DL21*SG2(ITYP)
  100 B(ITYP)=-DL12*SG1(ITYP)-DL22*SG2(ITYP)-SSI(ITYP)
      CALL DLAMDP(A,B,DL3,RDP,ADP,HDP,ITEST)
      IF(ITEST.EQ.1) THEN
      X=X/2.D0
      GOTO 90
      ENDIF
      DL1=DL11+DL3*DL12
      DL2=DL21+DL3*DL22
      DO 105 ITYP=1,3
      SGG(ITYP)=SFG(ITYP)-DL1*SG1(ITYP)-DL2*SG2(ITYP)
  105 SGG(ITYP)=SGG(ITYP)-DL3*SSI(ITYP)+X*DSFG(ITYP)
      IF(DL1.LT.-1.D-10) IDAM(1)=0
      IF(DL2.LT.-1.D-10) IDAM(2)=0
      IF(DL3.LT.-1.D-10) IDAM(3)=0
      IF(IDAM(1).EQ.0.OR.IDAM(2).EQ.0.OR.IDAM(3).EQ.0) THEN
      DL1=0.D0
      DL2=0.D0
      DL3=0.D0
      X=0.D0
      GOTO 2000
      ENDIF
C
C-------------------------------------------------------
C     ON  VERIFIE SI LA NORMALE DU DRUCKER VARIE PEU
C-------------------------------------------------------
C
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO=VAL/VA1/VA2
      IF(CO.LT.0.99) THEN
      X=X/2.D0
      GOTO 90
      ENDIF
C
C-------------------------------------------------
C     CAS OU LA TRACTION PASSE EN DESSOUS DE 0
C-------------------------------------------------
C
      IF((RT1-H1*DL1).GE.-1.E-10) GOTO 200
C
C
      IF(H1.EQ.0.D0) H1=1.D0
      DL1=RT1/H1
      AT(1)=-DSFG(1)
      AT(3)=-DSFG(2)
      AT(4)=Y-H2
      AT(2)=Y*ANU
      BT(1)=SFG(1)-DL1*SG1(1)
      BT(2)=SFG(2)-DL1*SG1(2)-RT2
      CALL SYLIN2(AT,BT,X1,DL21)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,X2,DL22)
      DO 110 ITYP=1,3
      A(ITYP)=SFG(ITYP)+X1*DSFG(ITYP)-DL1*SG1(ITYP)-DL21*SG2(ITYP)
  110 B(ITYP)=X2*DSFG(ITYP)-DL22*SG2(ITYP)-SSI(ITYP)
      CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST)
      DO 111 ITYP=1,3
      SGG(ITYP)=A(ITYP)+DL31*B(ITYP)
  111 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP)
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO1=VAL/VA1/VA2
      CALL SIJ(DSGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO2=VAL/VA1/VA2
      XIN1=X1+X2*DL31
      XIN2=X1+X2*DL32
      DLIN1=DL21+DL31*DL22
      DLIN2=DL21+DL32*DL22
      IF(XIN1.GT.-1.E-10.AND.DLIN1.GT.-1.E-10.AND.XIN1.LE.X
     1 .AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      DL2=DLIN1
      GOTO 200
      ENDIF
      IF(XIN2.GT.-1.E-10.AND.DLIN2.GT.-1.E-10.AND.XIN2.LE.X
     1 .AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      DL2=DLIN2
      GOTO 200
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,10101)
10101 FORMAT(1X,'T1T2DP  ERREUR TRAC 1')
      DL3=DL31
      X=XIN1
      DL2=DLIN1
  200 IF((RT2-H2*DL2).GE.-1.E-10) GOTO 300
C
C-------------------------------------------------
C     CAS OU LA TRACTION 2 PASSE EN DESSOUS DE 0
C-------------------------------------------------
C
C
C
      IF(H2.EQ.0.D0) H2=1.D0
      DL2=RT2/H2
      AT(1)=-DSFG(1)
      AT(3)=-DSFG(2)
      AT(2)=Y-H1
      AT(4)=Y*ANU
      BT(1)=SFG(1)-DL2*SG2(1)-RT1
      BT(2)=SFG(2)-DL2*SG2(2)
      CALL SYLIN2(AT,BT,X1,DL11)
      BT(1)=-SSI(1)
      BT(2)=-SSI(2)
      CALL SYLIN2(AT,BT,X2,DL12)
      DO 210 ITYP=1,3
      A(ITYP)=SFG(ITYP)+X1*DSFG(ITYP)-DL11*SG1(ITYP)-DL2*SG2(ITYP)
  210 B(ITYP)=X2*DSFG(ITYP)-DL12*SG1(ITYP)-SSI(ITYP)
      CALL DLAMD(A,B,DL31,DL32,RDP,ADP,HDP,ITEST)
      DO 211 ITYP=1,3
      SGG(ITYP)=A(ITYP)+DL31*B(ITYP)
  211 DSGG(ITYP)=A(ITYP)+DL32*B(ITYP)
      CALL SIJ(SGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO1=VAL/VA1/VA2
      CALL SIJ(DSGG,SI,SIETOI)
      CALL EPSDP(SI,SIETOI,ADP,EPC1)
      CALL SCAL(EPC,EPC1,VAL)
      CALL NORME(EPC,VA1)
      CALL NORME(EPC1,VA2)
      CO2=VAL/VA1/VA2
      XIN1=X1+X2*DL31
      XIN2=X1+X2*DL32
      DLIN1=DL11+DL31*DL12
      DLIN2=DL11+DL32*DL12
      IF(XIN1.GT.-1.E-10.AND.DLIN1.GT.-1.E-10.AND.XIN1.LE.X
     1 .AND.CO1.GT.0.9) THEN
      DL3=DL31
      X=XIN1
      DL1=DLIN1
      GOTO 300
      ENDIF
      IF(XIN2.GT.-1.E-10.AND.DLIN2.GT.-1.E-10.AND.XIN2.LE.X
     1 .AND.CO2.GT.0.9) THEN
      DL3=DL32
      X=XIN2
      DL1=DLIN2
      GOTO 300
      ENDIF
      IF(IIMPI.EQ.9) WRITE(IOIMP,20202)
20202 FORMAT(1X,'T1T2DP  ERREUR TRAC 2')
      DL3=DL31
      X=XIN1
      DL1=DLIN1
  300 DO 310 ITYP=1,3
      SGG(ITYP)=SFG(ITYP)-DL1*SG1(ITYP)-DL2*SG2(ITYP)
  310 SGG(ITYP)=SGG(ITYP)-DL3*SSI(ITYP)+X*DSFG(ITYP)
      RDP=RDP-DL3*HDP
      RMAX=MAX((RDP/1.73),(RDP/(1.D0-2.D0*ADP)))
      IF(RMAX.LT.SIREF) THEN
      RDP=0.D0
      GOTO 52
      ENDIF
      XT=XT+X
      DO 330 ITYP=1,3
      EPDP(ITYP)=EPDP(ITYP)+DL3*EPC(ITYP)
  330 SFG(ITYP)=SGG(ITYP)
      XLAM1=XLAM1+DL1
      XLAM2=XLAM2+DL2
      XLAM3=XLAM3+DL3
      RT1=RT1-H1*DL1
      IF(RT1.LT.1.E-10) THEN
      H1=0.D0
      RT1=0.D0
      ENDIF
      RT2=RT2-H2*DL2
      IF(RT2.LT.1.E-10) THEN
      H2=0.D0
      RT2=0.D0
      ENDIF
      C=1.D0-1.D-10
      IF(XT.GE.C) THEN
      DO 350 ITYP=1,3
      IDAM(ITYP)=0
  350 DSIG(ITYP)=0.D0
      CALL CHREP(SFG,SIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      ENDIF
      GOTO 51
 2000 CONTINUE
      CALL CHREP(SFG,SIG,-ANG)
      DO 2010 ITYP=1,3
 2010 DSFG(ITYP)=(1.D0-XT)*DSFG(ITYP)
      CALL CHREP(DSFG,DSIG,-ANG)
      CALL CHREP(EPDP,EPPLDP,-ANG)
      IF(IIMPI.EQ.9) WRITE(IOIMP,9998) ITER
      RETURN
      END