Sources Esope | Cast3M

ecoinc
C ECOINC    SOURCE    OF166741  25/11/04    21:15:45     12349          
      SUBROUTINE ECOINC(SIG0,DEPST,DDAUX,CMATE,VALMAT,VALCAR,N2EL,
     . N2PTEL,VARIN0,VAREX0,VAREXF,XMAT,PRECIS,WORK,CARAC,TRAC,
     . SIGF,VARINF,DEFP,KERRE,    IB,IGAU,NSTRS,EPAIST,MELE,NPINT,
     . NBPGAU,              SECT,LHOOK,TXR,XLOC,XGLOB,D1HOOK,
     .                           ROTHOO,DDHOMU,CRIGI,DSIGT,INPLAS)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C---------------------------------------------------------------------
C               ECOULEMENT PLASTIQUE POUR UN POINT
C                       D 'APRES INCA
C---------------------------------------------------------------------
C
C      EN ENTREE :
C
C                  SIG0       CONTRAINTES AU DEBUT DU PAS
C                  DEPST      INCREMENT DE DEFORMATIONS TOTALES
C                               ( THERMIQUE ENLEVEE )
C                  VARIN0     VARIABLES INTERNES DEDUT DU PAS
C                  VAREX0     VARIABLES EXTERNES DEBUT DU PAS
C                  VAREXF     VARIABLES EXTERNES FIN DU PAS
C                  XMAT       COEFFICIENTS DU MATERIAU
C                  PRECIS     PRECISION POUR ITERATIONS INTERNES
C                  CARAC      DES CARACTERISTIQUES
C                  TRAC       COURBE DE TRACTION
C                  MFR        INDICE DE FORMULATION
C                  NSTRS      NOMBRE DE CONTRAINTES CA2000
C                  INPLAS     NUMERO DU MODELE DE PLASTICITE
*                  DDAUX = MATRICE DE HOOKE ELASTIQUE
*                  CMATE = NOM DU MATERIAU
*                  VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
*                  VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
*                  N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
*                  N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
*                  IFOU  = OPTION DE CALCUL
*                  IB    = NUMERO DE L ELEMENT COURANT
*                  IGAU  = NUMERO DU POINT COURANT
*                  EPAIST= EPAISSEUR
*                  NBPGAU= NBRE DE POINTS DE GAUSS
*                  MELE  = NUMERO DE L ELEMENT FINI
*                  NPINT = NBRE DE POINTS D INTEGRATION
*                  NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
*                  NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
*                  SECT  = SECTION
*                  LHOOK = TAILLE DE LA MATRICE DE HOOKE
*                  TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
*                  UTILISES POUR LE CALCUL DE LA MATRICE DE HOOKE
C
C     EN SORTIE :
C
C                  SIGF        CONTRAINTES A LA FIN DU PAS
C                  VARINF      VARIABLES INTERNES FIN DU PAS
C                  DEFP        INCR. DE DEFORMATIONS PLASTIQUES
C                  KERRE       CODE D'ERREUR
C                   = 0  SI TOUT OK
C                   = 99 SI FORMULATION NON DISPONIBLE
C    EN PLASTICITE
C                   = 1  SI DEPS EST NEGATIF
C                   = 2  SI NOMBRE MAX D'ITERATIONS INTERNES DEPASSE
C    EN FLUAGE
C                   = 3  DIMINUER LE PAS ( FLUAGE RUNGE )
C                   = 4  DEPS EST NEGATIF
C                   = 5  PAS DE CONVERGENCE INTERNE ( ANCIEN NSTOP=0 )
C                   = 6  MATRICE SINGULIERE DANS CHAFLU
C    EN PLUS
C                   = 7  PROBLEME DANS LES CARACT DU TUYAU
C                   = 49 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 51 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 52 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 53 PB DANS LES ITER. INTERNES DRUCKER-PRAGER
C                   = 75 PB AVEC LA COURBE DE TRACTION
C
C-----------------------------------------------------------------------
C      VARIABLES PASSEES PAR LES COMMON ECOU ET NECOU
C
C  IFOUR     INDICE DU TYPE DE PROBLEME
C            -2  CONTRAINTES PLANES
C            -1  DEFORMATIONS PLANES
C             0
C             N  SERIE DE FOURIER
C  IMAPLA    INDICE DE MATERIAU PLASTIQUE
C            0  MATERIAU ELASTIQUE
C            4  MODELE DE CHABOCHE
C            5  MODELE DE DRUCKER-PRAGER
C            7  MODELE DE ROUSSELIER
C  ITYP      TYPE DE FORMULATION MECANIQUE
C       ITYP=1   CAS DES ELEMENTS MASSIFS
C       ITYP=2   CAS DES COQUES
C       ITYP=3   CAS DES MEMBRANES
C       ITYP=4   CAS DES CABLES ET DES BARRES
C       ITYP=5   CAS QUELCONQUE
C       ITYP=6   CAS DES CONTRAINTES PLANES
C       ITYP=7   CAS DES COQUES A NU=0. OU CONTRAINTES PLANES
C       ITYP=8   CAS DES MEMBRANES A NU=0. OU CONTRAINTES PLANES                 C       ITYP=9   CAS DES COQUES EPAISSES
C       ITYP=10  CAS DES JOINTS
C       ITYP=11  CAS DES POUTRES
C       ITYP=12  CAS DES TUYAUX
C       ITYP=13  CAS DES COQUES AVEC CISAILLEMENT TRANSVERSE
C  IPLAST    1 SI ON A PLASTIFIE , 0 SINON
C  ICINE  = 1 ECROUISSAGE CINEMATIQUE  ( CHABOCHE COMPRIS )
C  JFLUAG = 1 ON FLUE AVEC SIGMA
C  JFLUAG = 2 ON FLUE AVEC (SIG-X)
C  LFLUAG = 0 ON ECROUIT EN CAS DE FLUAGE
C  LFLUAG = 1 ON N'ECROUIT PAS EN CAS DE FLUAGE
C  JNTRIN = 0 LE TEMPS INTRINSEQUE EST LE TEMPS REEL
C  JNTRIN = 1 ON EVOLUE LE TEMPS INTRINSEQUE EN PLASTICITE ET FLUAGE
C  JNTRIN = 2 ON N'EVOLUE LE TEMPS INTRINSEQUE QU'EN FLUAGE
C  JNTRIN = 3 ON EST EN PLASTICITE ISOCHRONE
C
C  TETI = TEMPERATURE AU DEBUT DU PAS
C  TET  = TEMPERATURE A LA FIN DU PAS
C-----------------------------------------------------------------------
 
-INC PPARAM
-INC CCOPTIO
-INC CECOU
      COMMON/ECOU/TEST,ALFAH,
     1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
     2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6),
     1 DALPHA(6),EPSPLA(6),E(12),XINV(3),
     2 SIPLAD(6),DSIGP0(6),TET,TETI
C
      DIMENSION SIG0(*),DEPST(*),SIGF(*),VARIN0(*),VARINF(*),XMAT(*),
     .          VAREX0(*),VAREXF(*),DEFP(*),CARAC(*),TRAC(*),WORK(*)
      DIMENSION SIGMA(6),DSIGMA(6),SPHER(6),AUXIL(6),COEF(40)
      DIMENSION DIV(7),SIGFIN(6),DEFPLA(6),DSIGT(*)
      DIMENSION TABID(6),TABID2(6)
      DIMENSION DDAUX(LHOOK,*),TXR(IDIM,*),CRIGI(*)
      DIMENSION VALMAT(*),VALCAR(*),DDHOMU(LHOOK,*)
      DIMENSION XLOC(3,3),XGLOB(3,3)
      DIMENSION D1HOOK(LHOOK,*),ROTHOO(LHOOK,*)
      CHARACTER*8 CMATE
 
      REAL*8    DDOT
      EXTERNAL  DDOT
C ZZZZZZZZZZZZZZZZZZZZZZZ
C    DIMENSIONS A REVOIR
C ZZZZZZZZZZZZZZZZZZZZZ
      DIMENSION DHOOK(1),SIG(130),EPS(130),EI(1)
      DIMENSION ORMAT(1),ANORM(1)
      DIMENSION DSIGZE(6),DETFNO(1)
      DATA ITMAX/15/
      DATA PI4,R33,R22/0.785398164D0,1.732050808D0,1.414213562D0/
      DATA A,B,C,D/.577350269D0,.7071067814D0,.4082482904D0,
     . -0.8164965808D0/
      DATA A1/1.D0/
      DATA A2/.5D0/
      DATA A3/3.D0/
      DATA COEF/0.D0,1.5D0,1.5D0,3.D0,3.D0,3.D0,0.5D0,1.5D0,3.D0,0.5D0,
     . 1.5D0,3.D0,0.D0,0.D0,1.D0,0.D0,0.D0,0.D0,1.D0,
     . 1.D0,1.D0,2.D0,2.D0,2.D0,0.5D0,1.5D0,0.D0,3.D0,
     . 1.D0,0.D0,0.D0,1.D0,1.D0,1.D0,0.5D0,1.5D0,0.D0,3.D0,3.D0,3.D0/
      DATA XRF/0.0001D0/
C-----------------------------------------------------------------------
C         CONVENTION DE REMPLISSAGE DES MEMOIRES
C-----------------------------------------------------------------------
C
C     XMAT(1)    MODULE D'YOUNG
C     XMAT(2)    COEFFICIENT DE POISSON
C
C     TRAC(1 A 2*NCOURB) COURBE DE TRACTION
C     CARAC(  " +1) ALFAH POUR COQUES PLASTICITE GLOBALE
C     CARAC(  " +..) DONNEES POUR CRITERE POUTRES GLOBALES
C
C    MODELE ISOTROPE
C            VARIN0(1)  EPS*
C
C    MODELE CINEMATIQUE LINEAIRE
C            VARIN0(1)  EPS*
C            VARIN0(2) A VARIN0(1+IBOU)  DEPLACEMENT DE LA SPHERE
C
C    MODELE CHABOCHE
C            XMAT(5) ....      COEFFICIENTS A,C,...
C            VARIN0(1)  EPS*
C            VARIN0(2) A VARIN0(1+IBOU) DEPLACEMENT DE LA SPHERE 1
C            VARIN0(2+IBOU) A VARIN0(1+2*IBOU)  "    "  "   "    2
C            ICENT2 = 1 SI ON A 2 CENTRES
C
C    MODELE DRUCKER-PRAGER
C            XMAT(5) .... XMAT(13) = LES CONSTANTES DU MODELE
C            VARIN0(1)  EPS*
C
C    MODELE ROUSSELIER
C
C    MODELE DE CYCLAGE
C
C    MODELE A TEMPS INTRINSEQUE ( FLUAGE ET PLASTICITE ISOCHRONE )
C            VARIN0(.)  TIMEXI   TEMPS INTRINSEQUE DE DEBUT DE PAS
C
C--------------------------------------------------------------------
C                     REMPLISSAGE
C-----------------------------------------------------------------------
      YUNG=XMAT(1)
      XNU =XMAT(2)
C
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C    REMPLISSAGE A ACTUALISER SELON LES MODELES
CZZZZZZZZZZZZZZZZZZZZZZZZZZZ
CZZZZZZZZZZZZZZZZZZZZZZ
C      VARIABLES A RECUPERER
CZZZZZZZZZZZZZZZZZZZZZZ
C  ORMAT  MATERIAU ORTHOTROPE      DANS XMAT
C  ANORM    A VOIR         DANS XMAT OU DANS VARIN
C  EPENT    DANS VARIN
C  VECPRO,VALPRO  DANS XMAT
C  COVNMS,CVNMSD  DANS XMAT ( SELON LES CAS )
C
C  IT    DANS VAREXT
C  EPSFLU DANS VARINT
C  SIG,EPS DANS XMAT
C   DETFNO  DANS VARIN ????
C  DSIGZE DANS VARIN  ( POUR FLUAGE RELAX )
C  SIPLAD DANS VARIN  ( POUR FLUAGE RELAX  )  ...
C  ... ( CE SONT LES DEFORMATIONS OU CONTRAINTES (A VOIR) PLASTIQUES
C                            AU DEBUT DU PAS )
C  TEMPS , HPAS DANS VAREXT
C  TET,TETI DANS VAREXT
C
      EPSM1=VARIN0(1)
      TIMEXI=VARIN0(1)
      DPSM1=VARIN0(1)
      DPSM2=VARIN0(1)
      EPENT=VARIN0(1)
      EPSFLU=VARIN0(1)
      DETFNO(1)=VARIN0(1)
      DSIGZE(1)=VARIN0(1)
      TEMPS0=VAREX0(1)
      IT=nint(VAREX0(3))
      TEMPS=VAREXF(1)
      TETI=VAREX0(2)
      TET =VAREXF(2)
      DO I=1,NCOURB
        SIG(I)=TRAC(2*I-1)
        EPS(I)=TRAC(2*I)
      ENDDO
      ORMAT(1)=XMAT(1)
      DHOOK(1)=XMAT(1)
C----------------------------------------------------------------------
C                       INITIALISATIONS
C----------------------------------------------------------------------
      KERRE=0
      JA=1
      JC=1
      IA=1
C
C  PETIT TEST SUR NU POUR CERTAINS CAS
C
      IF(MFR.EQ.3.AND.IFOURB.EQ.-2.AND.XNU.NE.0.D0) THEN
        KERRE=38
        RETURN
      ENDIF
CZZZZZZZZZZZZZ
C  PROVISOIRE
CZZZZZZZZZZZZ
      ANORM(1)=XMAT(5)
      ICENT2=0
      IF(INPLAS.EQ.12.OR.INPLAS.EQ.13) ICENT2=1
      IF(INPLAS.EQ.7) NUMCHA=1
      IF(INPLAS.EQ.11) NUMCHA=2
      IF(INPLAS.EQ.12) NUMCHA=3
      IF(INPLAS.EQ.13) NUMCHA=4
C
      TEST=0.5D00*PRECIS
      HPAS=TEMPS-TEMPS0
      MCOD=1
*
*   CALCUL DE DSIGT
*
      CALL CALSIG(DEPST,DDAUX,NSTRS,CMATE,VALMAT,VALCAR,
     1        N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST,
     2      NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR,
     3    XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,DSIGT,IRTD)
*
      IF(IRTD.NE.1) THEN
         KERRE=69
         GOTO 1000
      ENDIF
*
*    TRAITEMENT DES COQ2 ET MASSIF 2D EN CAS DE MATERIAU
*    UNIDIRECTIONNEL
*
*      LES VARIABLES INTERNES SONT DANS LE REPERE UNIDIRECTIONNEL
*      ON SE LIMITE AU CAS AXISYMETRIQUE ?
*
            IF (CMATE.EQ.'UNIDIREC'.AND.INPLAS.NE.0) THEN
              IF(MELE.EQ.44.OR.(MFR.EQ.1.AND.IFOURB.LE.0)) THEN
                 APHI=0.D0
*
*                   CAS MASSIF : ON AJOUTE L'ANGLE DU REPERE
*                   DE L'ELEMENT AVEC LE REPERE GLOBAL
*
                 IF(MFR.EQ.1) APHI=ATAN2(TXR(2,1),TXR(1,1))
                 DO 1996 I=1,NSTRS
                    TABID(I)=SIG0(I)
                    TABID2(I)=DSIGT(I)
 1996            CONTINUE
                 ANG = ATAN2(XMAT(4),XMAT(3))  + APHI
                 CALL CHREP2(TABID,SIG0,ANG,NSTRS,MFR,1,KERRE)
                 CALL CHREP2(TABID2,DSIGT,ANG,NSTRS,MFR,1,KERRE)
                 IF(KERRE.NE.0) RETURN
               ELSEIF(MFR.EQ.1.AND.IFOURB.EQ.2) THEN
                 DO 996 I=1,NSTRS
                    TABID(I)=SIG0(I)
                    TABID2(I)=DSIGT(I)
 996             CONTINUE
                 CALL CHREP3(TABID,SIG0,XMAT,TXR,NSTRS,1,1,KERRE)
                 CALL CHREP3(TABID2,DSIGT,XMAT,TXR,NSTRS,1,1,KERRE)
 
              ENDIF
            ENDIF
*
      iecou=0
      inecou=0
      CALL VISAVI(SIG0,DSIGT,VARIN0,SIGMA,DSIGMA,SPHER,AUXIL,
     .            SIGF,DEFP,VARINF,SIGFIN,DEFPLA,
     .            DSIGZE,ICENT2,MCOD,IBOU,MFR,NSTRS,CARAC,CMATE
     >     ,iecou,inecou)
 
      IF(ITYP.EQ.0) THEN
      KERRE=99
      RETURN
      ENDIF
C-----------------------------------------------------------------------
      EI(1)=YUNG
      UNSYG=1.D00/YUNG
      U5SNU=1.5D00/(1.D00+XNU)
      U5DMU=U5SNU*YUNG
      SREF=YUNG*XRF
      IPLAST=0
      LCGDF=36
      LPLUS=6
      ITER=1
      LAPOIN=0
C
C  CAS DES COQUES EN GLOBAL - ON RECUPERE LE ALFAH
C
      ALFAH=1.D0
      IF(ITYP.EQ.2) ALFAH=CARAC(2)**2
      IF(ITYP.EQ.7) ALFAH=CARAC(2)**2
      UNALF=0.D0
      IF(ALFAH.GE.1.D-12) UNALF=1.D0/ALFAH
C
C  CAS DES BARRES ET DES CABLES
C
      IF(ITYP.EQ.4) THEN
      DIV(1)=1.D0
      DIV(2)=1.D0
      DIV(3)=1.D0/CARAC(1)
      ENDIF
C
C  CAS DES POUTRES
C
      IF(ITYP.NE.11) GO TO 841
      DIV(1)=1.D0/CARAC(4)
      DIV(2)=1.D0
      DIV(3)=1.D0
      DIV(4)=CARAC(10)/CARAC(1)
      DIV(5)=CARAC(11)/CARAC(2)
      DIV(6)=CARAC(12)/CARAC(3)
      IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(CARAC(1)*CARAC(4))
      IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(CARAC(2)*CARAC(4))
      IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(CARAC(3)*CARAC(4))
      GO TO 761
 841  CONTINUE
C
C  CAS DES TUYAUX
C
      IF(ITYP.NE.12) GO TO 842
      EPAIS=CARAC(1)
      REXT=CARAC(2)
      RMOY=REXT-EPAIS*0.5D0
      RACO=CARAC(3)
      GAM=1.D0
      XK=1.D0
      IF(RACO.EQ.0.D0) GO TO 765
      XLAM=RMOY*RMOY/EPAIS/RACO
      XK=1.65D0*XLAM
      IF(XK.LT.1.D0) XK=1.D0
      GAM=0.8888888888888889D0*(XLAM)**0.6666666666666667D0
      IF(GAM.LT.1.D0) GAM=1.D0
 765  CONTINUE
C
C     NB 23/09/98
C     VALEURS PAR DEFAUT POUR LES CFFX CFMX CFMY
C     CFMZ CFPR ( COEFFICIENTS POUR CALCULER LES
C     CONTRAINTES DE MEMBRANE, TORSION, FLEXIONS
C     DANS LE PLAN, HORS PLAN ET CIRCONFERENTIELLE
C     DUE A LA PRESSION )
C     POUR L'INSTANT PAS DE CONTRAINTE CIRCONFERENTIELLE
C     DUE A LA PRESSION ON N'UTILISE DONC PAS DIV(7)
C
C
      DIV(1)=1.D0
      DIV(2)=1.D0
      DIV(3)=1.D0
      DIV(4)=R33
      DIV(5)=PI4*GAM
      DIV(6)=DIV(5)
      DIV(7)=0.D0
      IF(IDIM.EQ.2) THEN
        PRES1=CARAC(6)
        CISA1=CARAC(7)
        IXCAR1=12
        IDEB1=8
      ELSE IF(IDIM.EQ.3) THEN
        PRES1=CARAC(7)
        CISA1=CARAC(8)
        IXCAR1=13
        IDEB1=9
      ENDIF
C
      JDIV1=2
      DO 1515 IBA=IDEB1,IXCAR1
      JDIV1=JDIV1+1
      VCAR1=CARAC(IBA)
      IF (VCAR1.NE.-1.D0) DIV(JDIV1)=CARAC(IBA)
 1515 CONTINUE
C
C     NB 23/09/98
C     TRANSFERT DE CFFX DANS DIV(1) ET REMISE A
C     1.D0 DE DIV(3)
C
      DIV(1) = DIV(3)
      DIV(3)=1.D0
C
      VX=CARAC(4)
      VY=CARAC(5)
      VZ=CARAC(6)
      CALL TUYCAR(CARAC,CISA1,VX,VY,VZ,KERRE,2)
      IF(KERRE.EQ.0) GO TO 843
      KERRE=7
      RETURN
  843 CONTINUE
C
C     NB 23/09/98
C     DIVISION DE DIV(5) ET DIV(6) PAR XK
C     CAR TUYCAR SORT LES INERTIES MODIFIEES
C     POUR LE CALCUL DES DEFORMATIONS PAR EPSIG
C
      DIV(5)=DIV(5)/XK
      DIV(6)=DIV(6)/XK
C
      DIV(1)=DIV(1)/CARAC(4)
      DIV(4)=DIV(4)*RMOY/CARAC(1)
      DIV(5)=DIV(5)*RMOY/CARAC(2)
      DIV(6)=DIV(6)*RMOY/CARAC(3)
761   CONTINUE
      E(1)=YUNG
      E(2)=0.D0
      E(3)=0.D0
      E(4)=YUNG*1.5D0/(1.D0+XNU)
      E(5)=YUNG
      E(6)=YUNG
  842 CONTINUE
C
      IF(ICINE.EQ.0.OR.JFLUAG.EQ.1) GO TO 204
C
C       ON EST EN CINEMATIQUE ( PLASTIQUE OU FLUAGE )
C       DANS CE CAS ON FLUE AVEC SIGMA-ALPHA QUE L'ON CALCULE
C
      DO IBA=1,IBOU
        SIGMA(IBA)=SIGMA(IBA)-SPHER(IBA)
      ENDDO
 204  CONTINUE
C
      IF(JNTRIN.EQ.0) GO TO 830
      TIMEXF=TIMEXI
      GO TO 831
830   TIMEXI=TEMPS-HPAS
      TIMEXF=TEMPS
831   CONTINUE
      IF(ITYP.NE.11.AND.ITYP.NE.12.AND.ITYP.NE.4) GO TO 844
      DO IBA=1,IBOU
        SIGMA(IBA)=SIGMA(IBA)*DIV(IBA)
        IF(ICINE.EQ.0) GO TO 845
        IF(JFLUAG.EQ.1.AND.LFLUAG.EQ.1) GO TO 845
        SPHER(IBA)=SPHER(IBA)*DIV(IBA)
 845    CONTINUE
        DSIGMA(IBA)=DSIGMA(IBA)*DIV(IBA)
      ENDDO
844   CONTINUE
      DO IBA=1,IBOU
        STOT(IBA)=SIGMA(IBA)+DSIGMA(IBA)
      ENDDO
C----------------------------------------------------------------------
C               CALCUL DE LA LIMITE ELASTIQUE
C----------------------------------------------------------------------
      IF(LFLUAG.EQ.1) GO TO 261
      IF(IMAPLA.NE.4) GO TO 262
      BPSTAR=EPSM1
      ICOD=1
      CALL CHALIM(BPSTAR,SELAS,XMAT,TET,ICOD,
     . BID,BID,BID,BID,BID1,BID2,BI3,BI4,BI5,BI6,IBID,IBID,NUMCHA)
      GO TO 261
  262 IF(IMAPLA.NE.5) GO TO 263
*                     REMPLACEMENT  XMAT(7) PAR XMAT(12) AM  8/3/90
*                     TEST SELON EPSM1                   AM 28/1/91
      IF(EPSM1.EQ.0.D0) THEN
      SELAS = XMAT(12)
      ELSE
        SELAS = XMAT(7) + EPSM1 * XMAT(13)
*
*                     AM  24/5/93  TEST SUR SELAS
*
        IF( SELAS.LT.0.D0) THEN
          SELAS = 0.D0
        ENDIF
*
      ENDIF
      GO TO 261
  263 CONTINUE
      EPSTAR=EPSM1
      IF(ICINE.EQ.1) EPSTAR=0.D0
      CALL TRACTI(SELAS,EPSTAR,SIG,EPS,NCOURB,2,izz)
      IF(izz.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
261   CONTINUE
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C   A VOIR    CE QU'IL Y A DANS ANORM
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      IF(IMAPLA.EQ.7) COVNMS(1)=ANORM(2*JC-1)
      IF(IMAPLA.EQ.7.AND.IT.EQ.1) ANORM(2*JC)=1.D-20
*
*      SI MATERIAU DRUCKER PRAGER ON CHERCHE LE CRITERE
*               AVEC LEQUEL ON DOIT FAIRE LA PROJECTION ET LE CRITERE
      IXMAT=5
      IF(IMAPLA.EQ.5.AND.EPSM1.EQ.0.D0) IXMAT=10
      CALL CKRIT(IMAPLA,STOT,ITYP,XMAT(IXMAT),ALFAH,COVNMS,XINV,SSTAR)
      SJ1=XINV(1)
      SJ2=XINV(2)
      IF(IMAPLA.EQ.0) GO TO 3
      IF(JFLUAG.NE.0) GO TO 4
C---------------------------------------------------------------------
C          A T'ON DEPASSE LA LIMITE ELASTIQUE
C---------------------------------------------------------------------
      PETI=1.1D00*TEST*SSTAR
      CALL EPSPRE(SSTAR,SELAS,PETI,ITRY)
      IF(ITRY.EQ.0.AND.SSTAR.GT.SELAS) GO TO 4
C----------------------------------------------------------------------
C         ON N'A PAS PLASTIFIE
C----------------------------------------------------------------------
    3 CONTINUE
      IF(ITYP.EQ.11.OR.ITYP.EQ.12.OR.ITYP.EQ.4) THEN
      DO IBA=1,IBOU
        IF(ICINE.NE.0)  SPHER(IBA)=SPHER(IBA)/DIV(IBA)
        STOT(IBA)=STOT(IBA)/DIV(IBA)
      ENDDO
      ENDIF
      DO IBA=1,IBOU
        SIGFIN(IBA)=STOT(IBA)
        DEFPLA(IBA)=0.D0
        DSIGMA(IBA)=0.D0
      ENDDO
      STARF=SSTAR
      EPST=EPSM1
      GO TO 31
C----------------------------------------------------------------------
C                ON A PLASTIFIE
C----------------------------------------------------------------------
 4    CONTINUE
      IPLAST=1
      PENTE=0.D0
      SLOPE=0.D0
      DO IBA=1,IBOU
        W1(IBA)=SIGMA(IBA)
      ENDDO
      CALL CKRIT(IMAPLA,W1,ITYP,XMAT(5),ALFAH,COVNMS,XINV,S0)
      EPSTAR=EPSM1
C
C    EN FLUAGE SI ON LA VITESSE DE FLUAGE NE DEPEND QUE
C    DU EPSILON EQUIVALENT DE FLUAGE CUMULE ON VA LE LIRE DANS EPSFLU
C
      IF(JFLUAG.NE.0.AND.IFLUPL.EQ.1) EPSTAR=EPSFLU
      EPSIN=EPSTAR
      SEL=0.D0
      SPLA=1.D0
      IF(JFLUAG.NE.0) GO TO 70
      CALL SIGELP(SIGMA,DSIGMA,SIGEL,DSIGP,STOT,PETI,
     .     ITYP,SEL,SPLA,IBOU,S0,SELAS,XMAT(5),COVNMS,ALFAH,
     .     IMAPLA,SSTAR)
      CALL SHIFTD(SIGEL,DALPHA,IBOU)
      GO TO 2333
C
   70 CONTINUE
C
      DO IBA=1,IBOU
        SIGEL(IBA)=SIGMA(IBA)+SEL*DSIGMA(IBA)
        DALPHA(IBA)=SIGEL(IBA)
        DSIGP(IBA)=SPLA*DSIGMA(IBA)
      ENDDO
 2333 CONTINUE
      DTIMEX=SPLA*HPAS
      IF(JNTRIN.EQ.1.OR.JNTRIN.EQ.3) TIMEXF=TIMEXI+DTIMEX
C
C           CALCUL DES PENTES
C
      IF(ABS(SPLA).LT.0.95D00) IELPLA=0
      IF(ICYCL.NE.1) GO TO 71
C
C            CAS DES MODELES DE CYCLAGE
C
      CALL NORMAL(A1,A2,A3,SIGEL,W1,ALFAH,SELAS,IBOU,ITYP)
CZZZZZZZZZZZZZZZZZZZZZZ
C    ON RECUPERE LA NORMALE DE LA DERNIERE SORTIE
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      DO IBA=1,IBOU
        W2(IBA)=ANORM(IBA)
      ENDDO
C
C   CALCUL DU PRODUIT SCALAIRE DE LA NORMALE AU POINT OU ON SORT
C   AVEC LA DERNIERE NORMALE OU ON EST SORTI
C
      X=DDOT(IBOU,W1,1,W2,1)
C                          ON STOCKE LA NORMALE
      DO IBA=1,IBOU
        ANORM(IBA)=W1(IBA)
      ENDDO
CZZZZZZZZZZZZZZZZZZZZZZZ
C   IL FAUDRA LA STOCKER DANS VARINF AU LIEU DE ANORM
CZZZZZZZZZZZZZZZZZZZZZZZZZ
C
C           CALCUL DU EPSI AVEC QUOI CALCULER LA PENTE
C
      EPSEN=0.5D00*(X+1.D00)*EPENT
      EPSIN=EPSEN
   71 CONTINUE
      IF(JFLUAG.EQ.0) GO TO 40
      IF(IMAPLA.EQ.4) GO TO 5
C
C         INITIALISATIONS POUR LE FLUAGE
C         QUEL EST LE EPSILON CUMULE DE FLUAGE
C
      SLOPE=0.D0
      PENTE=0.D0
      IF(ICINE.EQ.0) GO TO 5
      IF(JFLUAG.EQ.2) GO TO 40
C
C      CAS CINE + FLUAGE : CALCUL DE LA PENTE A LA COURBE DE TRACTION
C
      IF(ICINE.EQ.2) EPSIN=EPS(NCOURB)
      CALL TRACTI(PENTE,EPSIN,SIG,EPS,NCOURB,1,izz)
      IF (izz.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      GO TO 5
   40 CONTINUE
      IF(IMAPLA.EQ.5.OR.IMAPLA.EQ.4) GO TO 5
C
C CALCUL DE LA PENTE DE LA COURBE DE TRACTION AU POINT X-EPSTAR
C
      SLOPE=0.D0
      CALL TRACTI(PENTE,EPSIN,SIG,EPS,NCOURB,1,izz)
      IF(izz.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      IF(ICINE.EQ.0) GO TO 5
 6    SLOPE=PENTE
      PENTE=0.
      PANTIN=SLOPE
      IF(ICINE.EQ.1) GO TO 5
      EPSIN=EPS(NCOURB)
      CALL TRACTI(PANTIN,EPSIN,SIG,EPS,NCOURB,1,izz)
      IF(izz.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      PENTE=SLOPE-PANTIN
      SLOPE=PANTIN
 5    CONTINUE
      IF(IMAPLA.NE.4) GO TO 5621
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                       MODELE CHABOCHE
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      iecou=0
      inecou=0
      IF(JFLUAG.EQ.0)
     .  CALL CHABOK(YUNG,XNU,IA,EI,
     .  XMAT,SPHER,JA,IBOU,SI,DEPS,EPST,EPSTAR,WORK,KERRE,SN,
     .  AUXIL,NUMCHA,iecou,inecou)
      IF(JFLUAG.NE.0)
     .  CALL CHAFLU(YUNG,XNU,IA,EI,SSTAR,
     1  XMAT,SPHER,IBOU,SI,DEPS,EPST,EPSTAR,WORK,
     2  AUXIL,DPSM1,DPSM2,KERRE,NUMCHA,iecou,inecou)
      IF(KERRE.NE.0) RETURN
      IF(JFLUAG.NE.0.AND.IMAPLA.EQ.4) GO TO 689
      IF(IMAPLA.EQ.4) GO TO 570
 5621 CONTINUE
      IF(IMAPLA.NE.5) GO TO 9997
 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                MODELE DRUCKER-PRAGER  -  PRELIMINAIRES
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     ON EST DANS LE CAS DU MATERIAU DRUCKER-PRAGER.ON VA FAIRE UN
C     TEST POUR SAVOIR SI L'ON DOIT PROJETER AU SOMMET DU CONE ET
C     EVENTUELLEMENT FAIRE CETTE PROJECTION.
C     DANS CETTE VERSION,LE TEST CONSISTE A SAVOIR SI L'ETAT DE
C     CONTRAINTES CALCULE ELASTIQUEMENT SE TROUVE D'UN COTE OU DE
C     L'AUTRE D'UN PLAN ORTHOGONAL A L'AXE ET PASSANT PAR LE SOMMET
C     DU CONE,PUIS A SAVOIR SI L'ON EST A L'INTERIEUR OU A
C     L'EXTERIEUR DU CONE DES NORMALES AU CRITERE DE DRUCKER.
C     ON CALCULE LE I1 QUI DELIMITE LES DEUX REGIONS ET ON LE
C     COMPARE AVEC LE PREMIER INVARIANT DU TENSEUR CALCULE ELASTIQUEMENT
C
C     SNN=XMAT(7)+XMAT(13)*EPSTAR
*                     AM  24/5/93  TEST SUR SNN
*
*     IF( SNN*XMAT(7).LT.0.D0) THEN
*       SNN=0.D0
*       XMAT(13)=0.D0
*     ENDIF
C
C     MILL  17/3/93
C     SI XMAT(5)=0. PAS DE TEST PAR RAPPORT AU SOMMET
C
C     IF(XMAT(5).EQ.0.D0) GO TO 99970
C
C     XI1LIM=SNN/XMAT(5)
C     IF(SJ1.LE.XI1LIM) GO TO 99970
C
C     ON TESTE ENSUITE PAR RAPPORT AU CONE DES NORMALES
C
C     TESTDR=-XMAT(9)*SJ1/(2.D0*XMAT(8))+SJ2+XI1LIM*XMAT(9)/
C    .                    (2.D0*XMAT(8)*XMAT(5))
C     IF(TESTDR.GT.0.D00)GO TO 99970
C
C     ON EST MAINTENANT DANS LE CAS OU L'ON DOIT PROJETER AU SOMMET
C
C     CALL SOMDRU(IBOU,SI,DEPS,EPST,EPSTAR,SNN,XMAT(5),YUNG,XNU,KERRE)
C     IF(KERRE.NE.0) RETURN
C     LAPOIN=1
C     GO TO 570
99970 CONTINUE
      CALL SHIFTD(DHOOK,WORK(109),36)
      iecou=0
      inecou=0
      CALL PRJDRU(SSTAR,SELAS,PENTE,IBOU,SI,DEPS,EPST,EPSTAR,ITER,
     .   SN,WORK(109),WORK,WORK(37),WORK(73),XMAT(5),YUNG,
     .   XNU,LAPOIN,KERRE,iecou,inecou)
      IF(KERRE.NE.0) RETURN
      GO TO 570
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C                        MODELE ROUSSELIER
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9997  IF(IMAPLA.EQ.7) THEN
      CALL FISSUR
C     CALL FISSUR(YUNG,XNU,U5DMU,XMAT(5),A,B,C,ANORM,
C    1 JC,SIG,EPS,SELAS,SSTAR,IBOU,JKI,PENTE,EPSTAR,EPST,ITMAX,COEF,
C    1 IA,S0,SI,DEPS,SN,IMPRI)
      GO TO 570
      ENDIF
C
      IF(IRELAX.EQ.0) GO TO 42
      DO IBA=1,IBOU
        SIPLAD(IBA)=0.D00
        DSIGP0(IBA)=DSIGZE(IBA)
      ENDDO
42    CONTINUE
      IF(JGRDEF.LE.1) GO TO 420
C
C       INTEGRATION EN GRANDES DEFORMATIONS D N'EST PLUS ISOTROPE
C ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C            A VOIR|||||||||||||||||||||||||
      CALL SHIFTD(DHOOK,WORK,36)
      RAI=DETFNO(1)
      CALL T6CDCT
C     CALL T6CDCT(CGDFNO((IA-1)*LCGDF+1),WORK,4,W1,RAI)
      CALL PRJGDF
C     CALL PRJGDF(SSTAR,SELAS,PENTE,IBOU,SI,DEPS,EPST,EPSTAR,
C    1 ITER,SN,TET,WORK,WORK(37),WORK(73),PANTIN,TIMEXF,SIG,EPS,IA,
C    1 WORK(109))
      CALL SHIFTD(W1,DEFPLA,6)
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      GO TO 570
  420 CONTINUE
C----------------------------------------------------------------------
C                 ECRITURE SOUS FORME DIAGONALE
C               DIAGONALISATION DE SIGEL ET DSIGP
C----------------------------------------------------------------------
      NK=1
      CALL XYZDIA(DSIGP,SIGEL,W1,W2,A,B,C,D,ITYP,VECPRO)
      IF(IRELAX.NE.0)
     1 CALL XYZDIA(SIPLAD,DSIGP0,W1,W2,A,B,C,D,ITYP,VECPRO)
C
      GO TO (60,61,61,63,64,65,66,66,68,69,78,79,80),ITYP
C
   80 CONTINUE
C                 CONTRAINTES PLANES AVEC TOUS LES CISAILLEMENTS
C
      RAI=YUNG/(1.D0-XNU*XNU)
      E(1)=RAI*(1.D0+XNU)*0.5D0
      E(2)=RAI*(1.D0-XNU)*1.5D0
      E(3)=0.D0
      E(4)=E(2)
      E(5)=E(2)
      E(6)=E(2)
      JKI=34
      GO TO 615
C
  79  CONTINUE
C                   TUYAUX
  78  CONTINUE
C                   POUTRES
C     LE REMPLISSAGE DE E EST FAIT EN AMONT
C
      JKI=28
      GO TO 615
   69 CONTINUE
      GO TO 640
   68 CONTINUE
      GO TO 640
   66 CONTINUE
C                   COQUE OU MEMBRANE A NU =0 OU CP
C
      RAI=YUNG
      DO IBA=1,6
        E(IBA)=0.D00
      ENDDO
      E(1)=RAI
      IF(ICINE.NE.0) E(1)=E(1)+SLOPE
      E(4)=E(1)*ALFAH
      W1(1)=DSIGP(1)/(E(1)+PENTE)
      IF(ITYP.EQ.7) THEN
        IF(ALFAH.NE.0.D0) THEN
        W1(4)=DSIGP(4)/(E(4)+PENTE)
        ELSE
        W1(4)=0.D0
        ENDIF
      ENDIF
      JKI=18
      DEPS=VNMISD(W1,ITYP,ALFAH,COVNMS)
      GO TO 640
   65 CONTINUE
C                      CONTRAINTES  PLANES
C
      RAI=YUNG/(1.D0-XNU*XNU)
      E(1)=RAI*(1.D0+XNU)*0.5D0
      E(2)=RAI*(1.D0-XNU)*1.5D0
      E(3)=0.D0
      E(4)=E(2)
      E(5)=E(3)
      E(6)=E(3)
      JKI=24
      GO TO 615
   64 CONTINUE
C                   MATER QUELCONQUE
C
      JKI=18
      NK=1
      DO IBA=1,6
        E(IBA)=VALPRO(IBA)
      ENDDO
      IF(ICINE.EQ.0) GO TO 6401
      DO IBA=1,6
        E(IBA)=E(IBA)+SLOPE
      ENDDO
 6401 CONTINUE
      DO 6406 IBA=1,6
      IF(E(IBA)+PENTE.EQ.0.D00) GO TO 6406
      W1(IBA)=DSIGP(IBA)/(E(IBA)+PENTE)
 6406 CONTINUE
      DEPS=VNMISD(W1,ITYP,ALFAH,CVNMSD)
      GO TO 640
C
   61 RAI=YUNG/(1.D0-XNU*XNU)
      E(1)=RAI*(1.D0+XNU)*0.5D0
      E(2)=RAI*(1.D0-XNU)*1.5D0
      E(3)=E(2)
      E(4)=E(1)*ALFAH
      E(5)=E(2)*ALFAH
      E(6)=E(3)*ALFAH
      JKI=6
  615 CONTINUE
      IF(ICINE.EQ.0)  GOTO 610
      DO IK=1,6
        E(IK)=E(IK)+SLOPE
      ENDDO
  610 CONTINUE
      IBO=2
      IF(ITYP.EQ.3.OR.ITYP.EQ.6) IBO=1
      IF(ITYP.NE.11.AND.ITYP.NE.12.AND.ITYP.NE.4) W1(3)=0.D0
      DO 614 I=1,IBOU
      IF(ITYP.EQ.6.AND.I.EQ.3) GO TO 614
      W1(I)=0.D0
      IF(PENTE+E(I).EQ.0.) GO TO 614
      W1(I)=DSIGP(I)/(E(I)+PENTE)
  614 CONTINUE
      DEPS=VNMISD(W1,ITYP,ALFAH,COVNMS)
      GO TO 640
C
   63 E(3)=U5DMU/U5SNU
      IF(ICINE.NE.0) E(3)=E(3)+SLOPE
      NK=3
      JKI=12
      E(2)=0.D0
      E(1)=0.D0
      E(4)=0.D0
      E(5)=0.D0
      E(6)=0.D0
      DEPS=DSIGP(3)/(E(3)+PENTE)
      GO TO 640
C
   60 CONTINUE
      JKI=0
      E(1)=0.D0
      DO I=2,6
        E(I)=U5DMU
      ENDDO
      IF(ICINE.EQ.0)  GO TO 603
      DO I=1,6
        E(I)=E(I)+SLOPE
      ENDDO
  603 CONTINUE
      EMU=PENTE+U5DMU
      RAI=YUNG
      IF(IMAPLA.EQ.5) EMU=U5DMU+3*XMAT(5)*XMAT(5)*RAI/(1.D0-2.D0*XNU)
      DEPS=(SSTAR-SELAS)/EMU
C
  640 CONTINUE
C
                  IF(JFLUAG.NE.0) GO TO 41
C-----------------------------------------------------------------------
C       PLASTICITE  - PROJECTION SUR LA SURFACE DE CHARGE
C-----------------------------------------------------------------------
      iecou=0
      inecou=0
      CALL PROJVM(SIG,EPS,EPST,EPSTAR,DEPS,PENTE,SN,PANTIN,
     1 SELAS,ITER,SO,SSTAR,SI,NK,IBOU,COEF,JKI,KERRE,iecou,inecou)
C
      GO TO 57
C-----------------------------------------------------------------------
C                       FLUAGE
C-----------------------------------------------------------------------
   41 CONTINUE
      iecou=0
      inecou=0
      CALL FLUAGE(ITER,EPSTAR,S0,EPST,SSTAR,SI,DPSM1,DPSM2,
     . IBOU,DEPS,ICONV,TIMEXI,TIMEXF,YUNG,KERRE,iecou,inecou)
      IF(KERRE.EQ.0) GO TO 8693
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C   TESTER  KERRE EN SORTIE DE FLUAGE
C  ON RECUPERE AUUSI ICONV EN SORTIE
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
8693  CONTINUE
689   CONTINUE
C-----------------------------------------------------------------------
C          ON STOKE LE EPS DE FLUAGE DANS EPSFLU DANS TOUS LES CAS
C-----------------------------------------------------------------------
      EPSFLU=EPST
C           ON CALCULE LE EPSILON INELASTIQUE CUMULE
C
      EPST=EPSM1
      IF(LFLUAG.EQ.1) GO TO 414
C             ON ECROUIT LE RAYON DE LA SURFACE DE CHARGE
C
      EPST=DEPS+EPSM1
  414 CONTINUE
      IF(IMAPLA.EQ.4) GO TO 570
      IF(ICINE.EQ.0) GO TO 413
      IF(LFLUAG.EQ.1) GO TO 57
C            CAS CINEMATIQUE
C
      SLIM=SELAS
      GO TO 57
  413 CONTINUE
C         CAS ISOTROPE
C
      IF(LFLUAG.EQ.1) GO TO 57
C            ON CHERCHE LE NOUVEAU RAYON DE LA SPHERE
C
      CALL TRACTI(SF,EPST,SIG,EPS,NCOURB,2,izz)
      IF (izz.EQ.1) THEN
        KERRE=75
        GOTO 1000
      ENDIF
*
      SLIM=SF
C
  57   CONTINUE
C
C----------------------------------------------------------------------
C                         RETOUR AUX AXES X Y Z
C----------------------------------------------------------------------
      CALL DIAXYZ(ITYP,SIGT,SIGEL,VECPRO,IBO,A,B,C,B1,B2)
      IF(JFLUAG.EQ.0) GO TO 10
      IF((IT.EQ.1.AND.IRELAX.EQ.0).OR.(ICONV.NE.0.AND.IRELAX.NE.0))
     . CALL DIAXYZ(ITYP,DSIGP,STOT,VECPRO,IBO,A,B,C,B1,B2)
10    CONTINUE
C-----------------------------------------------------------------------
C                     TRAITEMENT FINAL
C-----------------------------------------------------------------------
  570 CONTINUE
      EPSM1F=EPST
      IF(ICYCL.EQ.1) EPENTF=EPSEN+DEPS
      IF(IMAPLA.NE.7) GO TO 9100
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C     VOIR SI ON A BESOIN DE SIGT
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      DO IBA=1,IBOU
        SIGT(IBA)=SIGT(IBA)-W1(IBA)
      ENDDO
 9100 CONTINUE
      IF(ICINE.EQ.0) GO TO 913
      IF(LFLUAG.EQ.1) GO TO 913
      IF(JFLUAG.EQ.1) GO TO 94
      IF(IMAPLA.EQ.4) GO TO 90
C----------------------------------------------------------------------
C         CAS CINEMATIQUE - CALCUL DE DALPHA
C----------------------------------------------------------------------
      Z=PANTIN*DEPS/SELAS
      IF(Z.LT.0.01D00) GO TO 910
      EZZ=1.D0-EXP(-Z)
      Z=1.D0-EZZ/Z
      EZZ=EZZ/(1.D0-Z)
      Z=Z/(1.D0-Z)
      DO 91 IBA=1,IBOU
      DALPHA(IBA)=DALPHA(IBA)*EZZ+Z*(SIGEL(IBA)-DALPHA(IBA))
   91 CONTINUE
      GO TO 90
  910 DO 911 IBA=1,IBOU
      DALPHA(IBA)=Z*0.5D0*(DALPHA(IBA)+SIGEL(IBA))
  911 CONTINUE
      GO TO 90
   94 CONTINUE
      IF(IMAPLA.EQ.4) GO TO 913
      CALL MOVCEN(PENTE,DEPS,SELAS,JA,IBOU,DALPHA,SIGMA,SPHER,SIGEL)
  90  CONTINUE
 913  CONTINUE
C----------------------------------------------------------------------
C         CAS CINEMATIQUE    COORDONNEES DU NOUVEAU CENTRE DE LA SPHERE
C----------------------------------------------------------------------
      IF(ICINE.EQ.0) GO TO 29
C
C  MODELE CHABOC  MISE A JOUR DE SPHER DANS CHABOC
C  SPHER CONTIENT LA SOMME DES DEPLACEMENTS DES SPHERES
C
       IF(IMAPLA.EQ.4) GO TO 29
      IF(LFLUAG.EQ.1.AND.JFLUAG.NE.0) GO TO 29
      DO IBA=1,IBOU
        SPHER(IBA)=SPHER(IBA)+DALPHA(IBA)
      ENDDO
   29 CONTINUE
C-----------------------------------------------------------------------
C  CALCUL DU SIGMA FINAL ET DE LA NOUVELLE LIMITE ELASTIQUE
C-----------------------------------------------------------------------
      IF(ITYP.NE.11.AND.ITYP.NE.12.AND.ITYP.NE.4) GO TO 846
      DO IBA=1,IBOU
      STOT(IBA)=STOT(IBA)/DIV(IBA)
      IF(ICINE.EQ.0) GO TO 847
      IF(LFLUAG.EQ.1.AND.JFLUAG.NE.0) GO TO 847
      SPHER(IBA)=SPHER(IBA)/DIV(IBA)
      DALPHA(IBA)=DALPHA(IBA)/DIV(IBA)
847   CONTINUE
        SIGEL(IBA)=SIGEL(IBA)/DIV(IBA)
      ENDDO
846   CONTINUE
      DO IBA=1,IBOU
        SIGFIN(IBA)=SIGEL(IBA)
      ENDDO
      STARF=SI
      IF(JFLUAG.EQ.0) SLIMF=SN
      IF(JNTRIN.NE.0) EPENTF=TIMEXF
C
      IF(DEPS.LT.0.D00) KERRE=1
C          CALCUL DU NOUVEAU SIGMAP
      DO IBA=1,IBOU
         DSIGMA(IBA)=STOT(IBA)-SIGEL(IBA)
      ENDDO
      IF(ICINE.EQ.0) GO TO 240
      IF(JFLUAG.EQ.1) GO TO 240
      IF(LFLUAG.EQ.1) GO TO 240
      DO IBA=1,IBOU
        DSIGMA(IBA)=DSIGMA(IBA)-DALPHA(IBA)
      ENDDO
  240 CONTINUE
*
      IF(ITYP.EQ.6.AND.IMAPLA.EQ.5.AND.LAPOIN.EQ.0) THEN
      DO IBA=1,IBOU
        DEFPLA(IBA)=EPSPLA(IBA)
      ENDDO
      EPST=EPSM1 + VONEPS(DEFPLA,ITYP,ALFAH,COVNMS)
      GO TO 31
      ENDIF
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
C    ATTENTION LE EPSIG NE MARCHE QUE POUR CERTAINS CAS
CZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ
      CALL EPSIG(DSIGMA,DEFPLA,IFOURB,YUNG,XNU,ITYP,ORMAT,CARAC)
      IF(ITYP.EQ.6.AND.IMAPLA.NE.5) DEFPLA(3)=-DEFPLA(1)-DEFPLA(2)
  31  CONTINUE
C-----------------------------------------------------------------------
C                  REMPLISSAGE
C-----------------------------------------------------------------------
      VARINF(1)=EPST
CZZZZZZZZZZZ
C PROVISOIRE
CZZZZZZZZZZZ
      IF(JFLUAG.EQ.1) VARINF(1)=EPSFLU
      IF(ICINE.EQ.0) GO TO 280
C              ON CALCULE LA VRAIE CONTRAINTE EN CINEMATIQUE
      DO 281 IBA=1,IBOU
      IF(JFLUAG.EQ.1) GO TO 284
      SIGFIN(IBA)=SIGFIN(IBA)+SPHER(IBA)
284   IF(LFLUAG.EQ.1) GO TO 281
      IF(ICENT2.NE.0) SPHER(IBA)=SPHER(IBA)-AUXIL(IBA)
  281 CONTINUE
  280 CONTINUE
*
*  MISE A 0. DES COMPOSANTES NON CALCULEES SI ITYP=8
*
      IF(ITYP.EQ.8.AND.IBOU.LT.NSTRS) THEN
         DO 285 IBA=2,NSTRS
            DEFPLA(IBA)=0.D0
            SIGFIN(IBA)=0.D0
            IF(ICENT2.NE.0) SPHER(IBA)=0.D0
 285     CONTINUE
      ENDIF
*
      MCOD=2
      iecou=0
      inecou=0
      CALL VISAVI(SIG0,DSIGT,VARIN0,SIGMA,DSIGMA,SPHER,AUXIL,
     .            SIGF,DEFP,VARINF,SIGFIN,DEFPLA,
     .            DSIGZE,ICENT2,MCOD,IBOU,MFR,NSTRS,CARAC,CMATE,
     .            iecou,inecou)
C
C      RETRAITEMENT POUR LES COQ2 ET MASSIFS  EN MATERIAU
C      UNIDIRECTIONNEL
C
 
              IF (CMATE.EQ.'UNIDIREC'.AND.INPLAS.NE.0) THEN
                 IF(MELE.EQ.44.OR.(MFR.EQ.1.AND.IFOURB.LE.0)) THEN
                  DO 1997 I=1,NSTRS
                    TABID(I)=SIGF(I)
                    TABID2(I)=DEFP(I)
 1997             CONTINUE
                  CALL CHREP2(TABID,SIGF,-ANG,NSTRS,MFR,1,KERRE)
                  CALL CHREP2(TABID2,DEFP,-ANG,NSTRS,MFR,2,KERRE)
                  IF(KERRE.NE.0) RETURN
                 ELSEIF(MFR.EQ.1.AND.IFOURB.EQ.2) THEN
                  DO 997 I=1,NSTRS
                    TABID(I)=SIGF(I)
                    TABID2(I)=DEFP(I)
 997              CONTINUE
                  CALL CHREP3(TABID,SIGF,XMAT,TXR,NSTRS,1,2,KERRE)
                  CALL CHREP3(TABID2,DEFP,XMAT,TXR,NSTRS,2,2,KERRE)
                 ENDIF
              ENDIF
*
 1000 RETURN
      END