Sources Esope | Cast3M

chafon
C CHAFON    SOURCE    CHAT      05/01/12    21:54:13     5004
      SUBROUTINE CHAFON(EPST,F1,F2,F3, DYK,SI,
     . C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP,R0,RM,B,ICAS,WEP,
     . SINT,X1INT,X2INT,PSI,OME,ecou)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION F1(*),F2(*),F3(*),DYK(*),SINT(*),X1INT(*),X2INT(*)
C
*  Commun ECOU: sert de fourre-tout pour les tableaux
*
      SEGMENT ECOU
***     COMMON/ECOU/TEST,ALFAH,
       REAL*8 TEST, ALFAH,
     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),X1(6),X2(6),
     1 DALPHA(6),DSIGO(6),S(12),XINV(3),
     2 SIPLAD(6),DSIGP0(6),TET,TETI
      ENDSEGMENT
C      COMMON/ECOU/TEST,ALFAH,
C     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
C     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),X1(6),X2(6),
C     1 DALPHA(6),DSIGO(6),S(12),XINV(3),
C     2 SIPLAD(6),DSIGP0(6),TET,TETI
C
C  LES CONTRAINTES SONT DANS S ET LE VON-MISES CORRESPONDANT EST
C  DANS SI
C  ICAS = 1  ON CHERCHE UNE PREMIERE ESTIMATION DU HDEP
C  ICAS = 2  ON CALCULE LE VECTEUR INCREMENT DYK ET ON MET A JOUR
C            S ,SI,X1,X2
C  ICAS = 3  ON CALCULE LA PENTE ( STOCKEE DANS WEP ) POUR LA
C            DETERMINATION DU DEPSI
C  DANS DYK ON STOCKE DANS L ORDRE : SIG-X1-X2,X1,X2
C  R EST CALCULE EN ENTREE POUR UN EPST DONNE
C
      IF(IDIAM.NE.0) R=RM-(RM-R0)*EXP(-B*EPST)
      PHI=1.D00
      IF(PSI.NE.1.D00) PHI=1.D00+(PSI-1.D00)*EXP(-OME*EPST)
      UNSUR=1./SI
      H=B*(RM-R)
      GO TO (101,102,103,104,105,106,107,108,109,999,
     .       999,999,113,114),ITYP
101   F1(1)=(S(1)-0.5*(S(2)+S(3)))*UNSUR
      F1(2)=(S(2)-0.5*(S(1)+S(3)))*UNSUR
      F1(3)=(S(3)-0.5*(S(1)+S(2)))*UNSUR
      DO 10 I=4,6
      F1(I)=3.*S(I)*UNSUR
      F2(I)=F1(I)*0.5
10    F3(I)=G*F1(I)
      DO 11 I=1,3
11    F2(I)=F1(I)
      F3(1)=S(7)*F1(1)+S(8)*(F1(2)+F1(3))
      F3(2)=S(7)*F1(2)+S(8)*(F1(1)+F1(3))
      F3(3)=S(7)*F1(3)+S(8)*(F1(1)+F1(2))
      DUM=H+1.5*ELT
      GO TO (111,300,300),ICAS
111   DUM=DUM+1.5*S(9)*PHI
      GO TO 200
102   IBO=2
      I3=3
      AL1=SQRT(ALFAH)
      GO TO 121
103   I3=3
      GO TO 110
106   I3=4
110   IBO=1
121   DUM=H
      DO 125 IPP=1,IBO
      IP=3*(IPP-1)
      IP1=IP+1
      IP2=IP+2
      IP3=IP+I3
      F1(IP1)=(S(IP1)-0.5*S(IP2))*UNSUR
      F1(IP2)=(S(IP2)-0.5*S(IP1))*UNSUR
      IF(ITYP.NE.6) GO TO 127
      F1(3)=0.
      F2(3)=0.
      F3(3)=0.
127   CONTINUE
      F1(IP3)=3.*S(IP3)*UNSUR
      F2(IP1)=S(IP1)*UNSUR
      F2(IP2)=S(IP2)*UNSUR
      F2(IP3)=S(IP3)*UNSUR
      F3(IP1)=S(7)*F1(IP1)+S(8)*F1(IP2)
      F3(IP2)=S(7)*F1(IP2)+S(8)*F1(IP1)
      F3(IP3)=G*F1(IP3)
      DUM2=F1(IP1)*F3(IP1)+F1(IP2)*F3(IP2)+F1(IP3)*F3(IP3)
      IF(ITYP.NE.2.OR.IPP.NE.1) GO TO 128
      DUM2=DUM2*ALFAH
      F1(IP1)=F1(IP1)*AL1
      F1(IP2)=F1(IP2)*AL1
      F1(IP3)=F1(IP3)*AL1
128   CONTINUE
      DUM=DUM+DUM2
125   CONTINUE
      GO TO (124,300,300),ICAS
124   DUM=DUM+S(9)*PHI
      GO TO 200
104   I1=3
      GO TO 140
107   AL1=SQRT(ALFAH)
108   I1=1
140   DUM=H
      DO 141 I=1,IBOU
      F1(I)=0.
      F2(I)=0.
141   F3(I)=0.
      F1(I1)=S(I1)*UNSUR
      F2(I1)=F1(I1)
      F3(I1)=S(7)*F1(I1)
      DUM2=F1(I1)*F3(I1)
      IF(ITYP.NE.7) GO TO 142
      DUM2=DUM2*ALFAH
      F1(I1)=F1(I1)*AL1
      F1(4)=S(4)*UNSUR
      F2(4)=F1(4)
      F3(4)=S(7)*F1(4)
      DUM=DUM+F1(4)*F3(4)
142   DUM=DUM+DUM2
      GO TO (124,300,300),ICAS
113   DUM=H
      F1(1)=(S(1)-0.5*S(2))*UNSUR
      F1(2)=(S(2)-0.5*S(1))*UNSUR
      F1(3)=0.
      F2(3)=0.
      F3(3)=0.
      DO 1113 IB=4,6
      F1(IB)=3.*S(IB)*UNSUR
      F2(IB)=S(IB)*UNSUR
      F3(IB)=G*F1(IB)
1113  CONTINUE
      F2(1)=S(1)*UNSUR
      F2(2)=S(2)*UNSUR
      F3(1)=S(7)*F1(1)+S(8)*F1(2)
      F3(2)=S(7)*F1(2)+S(8)*F1(1)
      DUM2=F1(1)*F3(1)+F1(2)*F3(2)+F1(4)*F3(4)+F1(5)*F3(5)+F1(6)*F3(6)
      DUM=DUM+DUM2
      GO TO (124,300,300),ICAS
C= Modes de calcul UNIDIMENSIONNELs (1D)
114   F1(1)=(S(1)-0.5*(S(2)+S(3)))*UNSUR
      F1(2)=(S(2)-0.5*(S(1)+S(3)))*UNSUR
      F1(3)=(S(3)-0.5*(S(1)+S(2)))*UNSUR
      F2(1)=F1(1)
      F2(2)=F1(2)
      F2(3)=F1(3)
      F3(1)=S(7)*F1(1)+S(8)*(F1(2)+F1(3))
      F3(2)=S(7)*F1(2)+S(8)*(F1(1)+F1(3))
      F3(3)=S(7)*F1(3)+S(8)*(F1(1)+F1(2))
      DUM=H+1.5*ELT
      GO TO (1140,300,300),ICAS
1140  DUM=DUM+1.5*S(9)*PHI
      GO TO 200
105   CONTINUE
109   CONTINUE
C-----------------------------------------------------------------------
C                       ICAS = 1
C-----------------------------------------------------------------------
 200  CONTINUE
      DO 1 I=1,IBOU
      DUM=DUM-C1*F1(I)*X1(I)
      IF(ICENT2.EQ.0) GO TO 1
      DUM=DUM-C2*F1(I)*X2(I)
   1  CONTINUE
      HDEP=(SI-R)/DUM
      RETURN
C-----------------------------------------------------------------------
C             ICAS = 2 OU 3
C-----------------------------------------------------------------------
 300  CONTINUE
      DO 2 I=1,IBOU
      DYK(6+I)=(S(10)*PHI*F2(I)-C1*X1(I))*HDEP
      X1(I)=X1INT(I)+DYK(6+I)*WEP
      IF(ICENT2.EQ.0) GO TO 2
      DYK(12+I)=(S(11)*PHI*F2(I)-C2*X2(I))*HDEP
      X2(I)=X2INT(I)+DYK(12+I)*WEP
2     CONTINUE
      DO 3 I=1,IBOU
      DYK(I)=-F3(I)*HDEP -DYK(6+I)
      IF(ICENT2.EQ.0) GO TO 33
      DYK(I)=DYK(I) -DYK(12+I)
  33  S(I)=SINT(I)+DYK(I)*WEP
   3  CONTINUE
      SI=VONMIS(S,ITYP,ALFAH,COVNMS)
      IF(ICAS.EQ.2) RETURN
C
C  CALCUL DE LA PENTE POUR LA DETERMINATION DU DEPSI
C
      WEP=H
      DO 631 IB=1,IBOU
631   WEP=WEP-F1(IB)*DYK(IB)/HDEP
      RETURN
999   WRITE(6,7999)
7999  FORMAT('0 CHAFON - CAS NON IMPLEMENTE '/)
      RETURN
      END