Sources Esope | Cast3M

sigelp
C SIGELP    SOURCE    PV        18/11/20    21:15:05     1001           
      SUBROUTINE SIGELP(SIGMA,DSIGMA,SIGEL,DSIGP,STOT,STEST,ITYPE,
     1   SEL,SPLA,IBOU,S0,SELAS,XMAT,COVNMS,ALFAH,IMAPLA,SSTAR)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SIGMA(*),DSIGMA(*),STOT(*),SIGEL(*),DSIGP(*)
     1,XMAT(*),COVNMS(*)
      DATA UN/1.0D0/,ZERO/0.0D0/
-INC CCREEL
C
CCCCCCCCCCCCCCC CALCUL DE L INTERSECTION AVEC LE CONVEXE
C
      CALL ZDANUL(SIGEL,6)
      CALL ZDANUL(DSIGP,6)
      AA=VONMIS(DSIGMA,ITYPE,ALFAH,COVNMS)
      BB=VNMS12(DSIGMA,SIGMA,ITYPE,ALFAH,COVNMS)
      AA=AA*AA
      AA1=AA
      BB1=BB
      TST=STEST*STEST
      IF(IMAPLA.EQ.11) GO TO 20
      C=S0*S0
      CC=C-SELAS*SELAS
      IF(IMAPLA.NE.5) GO TO 11
C
CCCCCCCCCCCCCC   CAS DRUCKER PRAGER
C
      EE=TRACE(DSIGMA)
      FF=TRACE(SIGMA)
      XMAT22=XMAT(2)*XMAT(2)
      XMAT11=XMAT(1)*XMAT(1)
      XMAT13=XMAT(1)*SELAS
      AA=AA*XMAT22-EE*EE*XMAT11
      BB=BB*XMAT22-EE*FF*XMAT11+XMAT13*EE
      GO TO 11
  20  CONTINUE
C
CCCCCCCCCCCCAS CHEN ET CHEN
C
      IN=0
  90  CONTINUE
      IF (IN.EQ.0) GO TO 120
      AAA=AA
      BBB=BB
      CCC=CC
      DDD=DD
      BB=BB1
      AA=AA1
 120  CONTINUE
      DINV1=TRACE(DSIGMA)
      AA=AA/3+IN*DINV1*DINV1/6
      CC=VONMIS(SIGMA,ITYPE,ALFAH,COVNMS)
      DINVO=TRACE(SIGMA)
      CC=CC*CC/3+IN*DINVO*DINVO/6+XMAT(2)*DINVO/3
      CC=CC+(-UN+XMAT(1)*DINVO/3)*SELAS*SELAS
      BB=BB*2/3+IN*DINV1*DINVO/3+DINV1*XMAT(2)/3
      BB=BB+XMAT(1)*SELAS*SELAS*DINV1/3
      BB=BB/2
CCCCC
   11 CONTINUE
      IF(ABS(AA).GT.TST) GO TO 1
C
CCCCCCCCCC  AA EST NUL UNE SEULE SOLUTION
C
      IF(ABS(BB).GT.TST) GO TO 10
C
CCCCCCCCCC   BB EST NUL &lt;=> DSIGMA=0 SIGMA0=STOT
C
  110 CONTINUE
             if (s0.lt.-xpetit/xzprec.or.s0.gt.xpetit/xzprec)then
               XX=SELAS/S0
             else
               xx=xgrand*xzprec
             endif
      CALL AEQBX(SIGEL,SIGMA,XX,IBOU)
      CALL AEQBPC(DSIGP,SIGMA,SIGEL,UN,-UN,IBOU)
      SEL=0.
      SPLA=1.
      GO TO 555
  10  CONTINUE
C
CCCCCCCCCCCCCCAS BB NON NUL
C
      SEL=CC/(2.*BB)
      GO TO 333
  333 CALL AEQBPC(SIGEL,SIGMA,DSIGMA,UN,SEL,IBOU)
      SPLA=1.-SEL
      CALL AEQBX(DSIGP,DSIGMA,SPLA,IBOU)
      GO TO 555
    1 CONTINUE
C
CCCCCCCCCCCCCC  AA NON NUL ON A 2 SOLUTIONS
C
      IF(ABS(CC).GT.TST) GO TO 2
C
CCCCCCCCCCCCCC 1 SEULE SOLUTION CC NUL
C
      SEL=0.
      GO TO 333
   2  CONTINUE
C
CCCCCCCCCCCC   2 SOLUTIONS
C
      DD=4.*BB*BB-4.*AA*CC
      B=2.*BB
      IF(AA*CC.GT.0.) GO TO 3
C
CCCC          CC EST NEGATIF ON A 2 SOLUTIONS 1 NEGATIVE 1 POSITIVE
C      LA BONNE EST LA POSITIVE
C
      SEL=(-B+SQRT(DD))/(2.*AA)
      GO TO 333
    3 CONTINUE
C
CCCCCCCC    ON A 2 SOLUTIONS POSITIVES OU NEGATIVES
C
      IF(DD.GT.0.) GO TO 5
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCC LA DROITE S0 STOT NE COUPE PAS LE CONVEXE   C
C              ON CHERCHE LE POINT DU CONVEXE LE          C
CCC             PLUS PRES DE STOT ON L APPELLE SIGEL      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  130 CONTINUE
        if (sSTAR.lt.-xpetit/xzprec.or.SSTAR.gt.xpetit/xzprec)then
          XX=SELAS/SSTAR
        else
          xx=xgrand*xzprec
        endif
      CALL AEQBX(SIGEL,STOT,XX,IBOU)
      SEL=0.D0
      SPLA=1.D0
      CALL AEQBPC(DSIGP,STOT,SIGEL,UN,-UN,IBOU)
      GO TO 555
    5 CONTINUE
C
CCCCCCCCCCCCCCCCCCC   ON A 2 SOLUTIONS POSITIVES OU NEGATIVES
C
      SEL=(-B+SQRT(DD))/(2.*AA)
      GO TO 333
 555  CONTINUE
      IF(IMAPLA.NE.11.OR.IN.EQ.1) GO TO 666
      IF(IN.EQ.0) GO TO 80
      IF(ABS(SELMIN).GE.ABS(SEL)) GO TO 666
      SEL=SELMIN
      IF(ABS(AAA).LT.TST.AND.ABS(BBB).LT.TST) GO TO 110
      IF(ABS(AAA).GT.TST.AND.ABS(CCC).GT.TST.AND.(AAA*CCC).GT.ZERO.
     .   AND.DDD.LT.ZERO) GO TO 130
      IN=1
      GO TO 333
  80  SELMIN=SEL
      IN=-1
      GO TO 90
 666  CONTINUE
      RETURN
      END