Sources Esope | Cast3M

gdp
C GDP       SOURCE    CHAT      05/01/13    00:16:43     5004
      SUBROUTINE GDP(TENS,DTENS,RDP,ADP,GAMDP)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION TENS(3),DTENS(3)
C
C----------------------------------------------
C     CETTE ROUTINE PERMET DE CALCULER LE
C     GAMMA POUR ARRIVER SUR LA LIMITE DU
C     CRITERE DRUCKER PRAGER
C----------------------------------------------
C
      GAMDP=0.
      A=-ADP*(DTENS(1)+DTENS(2))*ADP*(DTENS(1)+DTENS(2))
     1 +DTENS(1)*DTENS(1)+DTENS(2)*DTENS(2)+3.D0*DTENS(3)*DTENS(3)
     1 -DTENS(1)*DTENS(2)
      B=2.D0*ADP*(DTENS(1)+DTENS(2))*(RDP-ADP*(TENS(1)+TENS(2)))
     1 +2.D0*DTENS(1)*TENS(1)+2.D0*DTENS(2)*TENS(2)+
     1 6.D0*DTENS(3)*TENS(3)-DTENS(1)*TENS(2)-DTENS(2)*TENS(1)
      C=RDP-ADP*(TENS(1)+TENS(2))
      C=C*C
      C=-C+TENS(1)*TENS(1)+TENS(2)*TENS(2)+3.D0*TENS(3)*TENS(3)
     1 -TENS(2)*TENS(1)
      X1=0.D0
      X2=0.D0
      IF(A.EQ.0.D0) GOTO 1000
      B=B/A
      C=C/A
      A=1.D0
      DIS=B*B-4.D0*A*C
      ADIS=ABS(DIS)
      IF(DIS.GT.0.D0) THEN
      DIS=SQRT(ADIS)
      X2=(-B+DIS)/2.D0
      X1=(-B-DIS)/2.D0
      GOTO 2000
      ENDIF
      IF(ADIS.LE.1.E-10) THEN
      X1=-B/2.D0
      X2=X1
      GOTO 2000
      ENDIF
      GAMDP=0.D0
      RETURN
 1000 IF(B.EQ.0.D0) THEN
      X1=0.D0
      X2=0.D0
      GOTO 2000
      ENDIF
      X1=-C/B
      X2=X1
 2000 CONTINUE
      S1=RDP-ADP*(TENS(1)+TENS(2)+X1*(DTENS(1)+DTENS(2)))
      S2=RDP-ADP*(TENS(1)+TENS(2)+X2*(DTENS(1)+DTENS(2)))
      IF(RDP.NE.0.D0) THEN
      S1=S1/RDP
      S2=S2/RDP
      ENDIF
      IF(S1.GE.-1.E-10.AND.S2.GE.-1.E-10) THEN
      GAMDP=MAX(X1,X2)
      GOTO 3000
      ENDIF
      IF(S1.GE.-1.E-10.AND.S2.LT.-1.E-10) THEN
      GAMDP=X1
      GOTO 3000
      ENDIF
      IF(S1.LT.-1.E-10.AND.S2.GE.-1.E-10) THEN
      GAMDP=X2
      GOTO 3000
      ENDIF
 3000 IF(GAMDP.LT.0.D0) GAMDP=0.D0
      IF(GAMDP.GT.1.D0) GAMDP=1.D0
      RETURN
      END