Sources Esope | Cast3M

sictens
C SICTENS   SOURCE    OF166741  25/11/04    21:16:07     12349          
 
C     Modèle ONERA pseudo-tensoriel pour le SICf/SiC tissé
C     variables en entree
C
C     WRK0,KRK1,WRK5  pointeurs sur des segments de travail
C
C
C     WRK0:
C            XMAT(NMATT): tableau composantes materiaux
C
C     WRK1:
C            SIG0(NSTRS) : contraintes au debut du pas
C            VAR0(NVARI) : variables internes au debut du pas
C            DEPST(NSTRS): increment de deformation totale
C
C     WRK22:
C            XXE: coordonnees de l'element en double precison
C
C     WTRAV:
C            TXR: cosinus directeurs des axes locaux de l'element massif
C
C
C     WRK5:
C
C
C
C     NSTRS      nombre de composantes dans les vecteurs des contraintes
C                                        et les vecteurs des deformations
C
C     NVARI      nombre de variables internes (doit etre egal a 3)
C
C     NMATT      nombre de constantes du materiau
C
C
C
C     variables en sortie
C
C     VARF      variables internes finales dans WRK0
C
C     SIGF      contraintes finales dans WRK0
C
      SUBROUTINE SICTENS(WRK52,WRK53,WRK54,WRK22,IB,IGAU,NVARI,NBPGAU
     & ,necou,iecou)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC DECHE
-INC TECOU
 
      SEGMENT WRK22
        REAL*8 XXE(3,NBNNbi)
      ENDSEGMENT
 
      DIMENSION ZIDAUX(6,6)
      DIMENSION H1 (6,6), H2(6,6), H3(6,6), H4(6,6), H5(6,6)
      DIMENSION H01 (6,6), H02(6,6), H03(6,6), H04(6,6), H05(6,6)
      DIMENSION SEFF(6,6),CEFF(6,6)
*
*      DORTH :  DEFO DAN S LE REPÈRE OrTHO
*      sigfv : CONTRAINTE RÉORGANISÉES ORDRE CASTEM
*      DEFI : DÉFORMATIONS À LA FIN DU PAS PRÉCEDENT
*      depst : INCREMENT DE DÉFO DANS LE REPÈRE GÉNÉ
*      depstV : INCREMENT DE DÉFO DANS LE REPÈRE GÉNÉ
*      epsel : DEFO ELASTIQUE CALCULÉE POUR L'OBTENTION DES CONTRAINTES
      DIMENSION SORTH(6),  SIGFV(6),EPSIN(6),EPSINV(6),
     &       DEFI(6), DEPSTV(6), DORTH(6), EPSEL(6),EPSTFV(6),
     &       SORTH1(6),SORTH2(6)
*     defoth : DEFO DANS LE REPERE ORTHO SOus FORME (3X3)
*       pOUR LE CALCUL DES DEFO PROPRES PUIS PARTIE >0 DES DÉFO PROPRES
*      xprop : DIRECTIONS PROPRES/REPERE ortho
*       xpinv : MATRICE INVERSE DE CHANGEMENT DE REPERE XPROP
*      xPINVT  : MATRICE TRANSPOSÉE DE XPINV
*      defotht : MATRICE TEMPORAIRE POUR LE CHANGEMENT DE REPÈRE
*      epsclo : DEFO DE FERMETURE (CLOSURE)
*      EPSBAR : DEFO - DEFO CLOSURE
      DIMENSION DEFOTH(3,3),PROP(3),XPROP(3,3),XPINV(3,3),
     &       DEFOTHT(3,3),XPINVT(3,3),EPSCLO(6),
     &       EPSBAR(6),EPSIM(6),SIGPLU(6)
      DIMENSION ZK1(6,6),ZK2(6,6),ZK3(6,6),ZK4(6,6),ZK5(6,6)
     &       ,TERME(6,6)
      DIMENSION EPSS1(6),EPSR1(6)
*
*     DRIVING FORCES NORMALES ET TRANSVERSALES
*     dom : VARIABLES D'ENDOMAGEMENT
      DIMENSION Y(5),DOM(5)
*     DEPSIF (DELTA EPSILON I f :
*     LA FERMETURE COMMENCE POUR epsbar = DELTA EPSILON
*     finit POUR  epsbar = DELTA EPSILON
*     ZNUACT  INDICE D'ACTIVATION =1 : ACTIVE
*                                =0 : PASSIVE
      DIMENSION DEPSIF(5),ZNUACT(5)
*
* Pseudo contraintes pour le calcul des
* incréments de defo residuelles EPSR1 et stockées EPS1
      DIMENSION SOMS(6,6),SOMR(6,6),
     $ SISIS(6),SISIR(6)
      DIMENSION Y0(5),AIF(5),PY(5),
     $ YC(5),DC(5),DEPSIF0(5)
      DATA ( ( H1(I,J), J=1,6), I=1,6)
     &    /1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 1.0/
      DATA ( ( H2(I,J), J=1,6), I=1,6)
     &    /0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 1., 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 1./
*
      DATA ( ( H3(I,J), J=1,6), I=1,6)
     &    /0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
*     + 45°
      DATA ( ( H4(I,J), J=1,6), I=1,6)
     &    /1.0, 0.0, 0.0, 0.0, 0.0, 0.5,
     &     0.0, 1.0, 0.0, 0.0, 0.0, 0.5,
     &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     &     0.0, 0.0, 0.0, 1.0, 1.0, 0.0,
     &     0.0, 0.0, 0.0, 1.0, 1.0, 0.0,
     &     0.5, 0.5, 0.0, 0.0, 0.0, 1.0/
*    -45°
      DATA ( ( H5(I,J), J=1,6), I=1,6)
     &    /1.0,   0.0, 0.0,  0.0,  0.0, -0.5,
     &     0.0,   1.0, 0.0,  0.0,  0.0, -0.5,
     &     0.0,   0.0, 0.0,  0.0,  0.0,  0.0,
     &     0.0,   0.0, 0.0,  1.0, -1.0,  0.0,
     &     0.0,   0.0, 0.0, -1.0,  1.0,  0.0,
     &     -0.5, -0.5, 0.0,  0.0,  0.0,  1.0/
 
*
*  Attention on rentre les racines carrée des valeurs
*  Y0 et YC telles que définies formules (2.21) et (2.24)
*  dans SEMT/LM2S/RT/05-034
*   PHENOMENOLOGICAL COEFFICIENTS
*  Coefficients pour les matrices H
*
* H1N, H1HP, H1P
* H2N, H2HP, H2P
* H3N, H3HP, H3P
* H4N, H4HP, H4P
*
* scalar damage indicators
* DOM1,DOM2,DOM3,DOM4,DOM5
* YGR1, YGR2, YGR3, YW(6)
* YEQ1, YEQ2, YEQ3
*
***** INITIALIZATION VARIABLES
*
      lhookl = WRK53.LHOOK
* Proprietes materiau
       YG1   = XMAT(1)
       YG2   = XMAT(2)
       YG3   = XMAT(3)
       XNU12 = XMAT(4)
       XNU23 = XMAT(5)
       XNU13 = XMAT(6)
       G12   = XMAT(7)
       G23   = XMAT(8)
       G13   = XMAT(9)
* Nombre de IIPARAM obligatoires pour le modèle 15 + ALP1 ALP2 ALP3 et RHO =19
       NOBL = 19
* On extrait les valeurs facultatives du modèle
      H1N = XMAT(NOBL + 1)
      IF(IVALMA(NOBL+1).EQ.0.) H1N = 1.
      H1HP = XMAT(NOBL + 2)
      IF(IVALMA(NOBL+2).EQ.0.) H1HP = 0.7
      H1P = XMAT(NOBL + 3)
      IF(IVALMA(NOBL+3).EQ.0.) H1P = 0.45
      H2N = XMAT(NOBL + 4)
      IF(IVALMA(NOBL+4).EQ.0.) H2N = 1.
      H2HP = XMAT(NOBL + 5)
      IF(IVALMA(NOBL+5).EQ.0.) H2HP = 0.7
      H2P = XMAT(NOBL + 6)
      IF(IVALMA(NOBL+6).EQ.0.) H2P = 0.45
      H3N = XMAT(NOBL+7)
      IF(IVALMA(NOBL+7).EQ.0.) H3N = 1.
      H3P = XMAT(NOBL + 8)
      IF(IVALMA(NOBL+8).EQ.0.) H3P = 0.7
      H4N = XMAT(NOBL + 9)
      IF(IVALMA(NOBL+9).EQ.0.) H4N = 1.
      H4HP = XMAT(NOBL + 10)
      IF(IVALMA(NOBL+10).EQ.0.) H4HP = 0.7
      H4P = XMAT(NOBL + 11)
      IF(IVALMA(NOBL+11).EQ.0.) H4P = 1.2
      Y0(1) = XMAT(NOBL + 12)
      IF(IVALMA(NOBL+12).EQ.0.) Y0(1) =173.205
      Y0(2) = XMAT(NOBL + 13)
      IF(IVALMA(NOBL+13).EQ.0.) Y0(2) =173.205
      Y0(3) = XMAT(NOBL + 14)
      IF(IVALMA(NOBL+14).EQ.0.) Y0(3) =173.205
      Y0(4) = XMAT(NOBL + 15)
      IF(IVALMA(NOBL+15).EQ.0.) Y0(4) =100.00
      Y0(5) = XMAT(NOBL + 16)
      IF(IVALMA(NOBL+16).EQ.0.) Y0(5) =100.00
      YC(1) = XMAT(NOBL + 17)
      IF(IVALMA(NOBL+17).EQ.0.) YC(1) =1870.83
      YC(2) = XMAT(NOBL + 18)
      IF(IVALMA(NOBL+18).EQ.0.) YC(2) =1870.83
      YC(3) = XMAT(NOBL + 19)
      IF(IVALMA(NOBL+19).EQ.0.) YC(3) =1870.83
      YC(4) = XMAT(NOBL + 20)
      IF(IVALMA(NOBL+20).EQ.0.) YC(4) =3464.1
      YC(5) = XMAT(NOBL + 21)
      IF(IVALMA(NOBL+21).EQ.0.) YC(5) =3464.1
      DC(1) = XMAT(NOBL + 22)
      IF(IVALMA(NOBL+22).EQ.0.) DC(1) =4.
      DC(2) = XMAT(NOBL + 23)
      IF(IVALMA(NOBL+23).EQ.0.) DC(2) =4.
      DC(3) = XMAT(NOBL + 24)
      IF(IVALMA(NOBL+24).EQ.0.) DC(3) =4.
      DC(4) = XMAT(NOBL + 25)
      IF(IVALMA(NOBL+25).EQ.0.) DC(4) =4.
      DC(5) = XMAT(NOBL + 26)
      IF(IVALMA(NOBL+26).EQ.0.) DC(5) =4.
      PY(1) = XMAT(NOBL + 27)
      IF(IVALMA(NOBL+27).EQ.0.) PY(1) =1.
      PY(2) = XMAT(NOBL + 28)
      IF(IVALMA(NOBL+28).EQ.0.) PY(2) =1.
      PY(3) = XMAT(NOBL + 29)
      IF(IVALMA(NOBL+29).EQ.0.) PY(3) =1.
      PY(4) = XMAT(NOBL + 30)
      IF(IVALMA(NOBL+30).EQ.0.) PY(4) =1.2
      PY(5) = XMAT(NOBL + 31)
      IF(IVALMA(NOBL+31).EQ.0.) PY(5) =1.2
      B1 = XMAT(NOBL + 32)
      IF(IVALMA(NOBL+32).EQ.0.) B1 =1.
      B2 = XMAT(NOBL + 33)
      IF(IVALMA(NOBL+33).EQ.0.) B2 =1.
      B3 = XMAT(NOBL + 34)
      IF(IVALMA(NOBL+34).EQ.0.) B3 =1.
      DELTAL = XMAT(NOBL + 35)
      IF(IVALMA(NOBL+35).EQ.0.) DELTAL =0.
      TZERO = XMAT(NOBL + 36)
      IF(IVALMA(NOBL+36).EQ.0.) TZERO =0.
      DEPSIF0(1) = XMAT(NOBL + 37)
      IF(IVALMA(NOBL+37).EQ.0.) DEPSIF0(1) =0.0003
      DEPSIF0(2) = XMAT(NOBL + 38)
      IF(IVALMA(NOBL+38).EQ.0.) DEPSIF0(2) =0.0003
      DEPSIF0(3) = XMAT(NOBL + 39)
      IF(IVALMA(NOBL+39).EQ.0.) DEPSIF0(3) =0.0003
      DEPSIF0(4) = XMAT(NOBL + 40)
      IF(IVALMA(NOBL+40).EQ.0.) DEPSIF0(4) =0.0005
      DEPSIF0(5) = XMAT(NOBL + 41)
      IF(IVALMA(NOBL+41).EQ.0.) DEPSIF0(5) =0.0005
      AIF(1) = XMAT(NOBL + 42)
      IF(IVALMA(NOBL+42).EQ.0.) AIF(1) =0.5
      AIF(2) = XMAT(NOBL + 43)
      IF(IVALMA(NOBL+43).EQ.0.) AIF(2) =0.5
      AIF(3) = XMAT(NOBL + 44)
      IF(IVALMA(NOBL+44).EQ.0.) AIF(3) =0.5
      AIF(4) = XMAT(NOBL + 45)
      IF(IVALMA(NOBL+45).EQ.0.) AIF(4) =1.
      AIF(5) = XMAT(NOBL + 46)
      IF(IVALMA(NOBL+46).EQ.0.) AIF(5) =1.
      ETA1 = XMAT(NOBL + 47)
      IF(IVALMA(NOBL+47).EQ.0.) ETA1  =0.1
      ETA2 = XMAT(NOBL + 48)
      IF(IVALMA(NOBL+48).EQ.0.) ETA2  =0.1
      ETA3 = XMAT(NOBL + 49)
      IF(IVALMA(NOBL+49).EQ.0.) ETA3 =0.0
      ETA4 = XMAT(NOBL + 50)
      IF(IVALMA(NOBL+50).EQ.0.) ETA4 =0.1
      ETA5 = XMAT(NOBL + 51)
      IF(IVALMA(NOBL+51).EQ.0.) ETA5 =0.1
 
* Pour debuggage
      IIIPARAM=0
* Variables internes du modele
      DOM1=VAR0(2)
      DOM2=VAR0(3)
      DOM3=VAR0(4)
      DOM4=VAR0(5)
      DOM5=VAR0(6)
*
****** CALCUL DE LA DEFORMATION INITIALE A PARTIR
****** DES CONTRAINTES INITIALES
 
*
      IF (IIIPARAM.EQ.1) THEN
         WRITE(IOIMP,66770) IB,IGAU
66770    format(////////2x,'element  ',i6,2x,'point  ',i3//)
 
         WRITE (IOIMP,*) 'Increment des deformations '
         WRITE (IOIMP,99999) (DEPST(ILOOP), ILOOP=1,6)
 
         WRITE (IOIMP,*) 'Contraintes au debut de l''iteration'
         WRITE (IOIMP,99999) (SIG0(ILOOP), ILOOP=1,6)
      ENDIF
 
* Controle si il faut calculer la matrice de hook
 
      IF (IB.EQ.1.AND.IGAU.EQ.1) THEN
         GOTO 100
      ELSEIF (N2PTEL.EQ.1.AND.N2EL.EQ.1) THEN
         GOTO 200
      ELSEIF (N2PTEL.EQ.1.AND.N2EL.NE.1)  THEN
         IF (IGAU.EQ.1) THEN
            GOTO 100
         ELSE
            GOTO 200
         ENDIF
      ELSE
         GOTO 100
      ENDIF
 
100   CONTINUE
      IPERR=1
* Calcul de la matrice de hook pour un materiau orthotrope: C0
      CALL ZERO (DDAUX,LHOOKL,LHOOKL)
      XNU21=(YG2/YG1)*XNU12
      XNU32=(YG3/YG2)*XNU23
      XNU31=(YG3/YG1)*XNU13
      AUX=(1.D0-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     &    -2.D0*XNU21*XNU32*XNU13)
      AUX1=AUX/YG1
      AUX2=AUX/YG2
      AUX3=AUX/YG3
      DDAUX(1,1)=(1.D0-XNU23*XNU32)/AUX1
      DDAUX(1,2)=(XNU21+XNU31*XNU23)/AUX1
      DDAUX(2,1)=DDAUX(1,2)
      DDAUX(1,3)=(XNU31+XNU21*XNU32)/AUX1
      DDAUX(3,1)=DDAUX(1,3)
      DDAUX(2,2)=(1.D0-XNU13*XNU31)/AUX2
      DDAUX(2,3)=(XNU32+XNU12*XNU31)/AUX2
      DDAUX(3,2)=DDAUX(2,3)
      DDAUX(3,3)=(1.D0-XNU12*XNU21)/AUX3
      DDAUX(4,4)=G23
      DDAUX(5,5)=G13
      DDAUX(6,6)=G12
* On recopie la matrice de Hooke pour l'inverser
200   CONTINUE
      DO JLOOP=1,LHOOKL
         DO ILOOP=1,LHOOKL
            ZIDAUX(ILOOP,JLOOP)= DDAUX(ILOOP,JLOOP)
         ENDDO
      ENDDO
* Inversion de la matrice de Hooke
      TPREC= 0.D1
      IPERR=0
      CALL INVALM(ZIDAUX,LHOOKL,LHOOKL,IPERR,TPREC)
      IF (IPERR.NE.0) THEN
        WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK'
      ENDIF
* Calcul des tenseurs d'endommagement H01, H02, H03
* selon les equations 4.2.3 à 4.2.5
      CALL ZERO (H01,LHOOKL,LHOOKL)
      CALL ZERO (H02,LHOOKL,LHOOKL)
      CALL ZERO (H03,LHOOKL,LHOOKL)
      CALL ZERO (H04,LHOOKL,LHOOKL)
      CALL ZERO (H05,LHOOKL,LHOOKL)
 
      H01(1,1)=H1(1,1)*H1N*ZIDAUX(1,1)
      H01(5,5)=H1(5,5)*H1HP*ZIDAUX(5,5)
      H01(6,6)=H1(6,6)*H1P*ZIDAUX(6,6)
      H02(2,2)=H2(2,2)*H2N*ZIDAUX(2,2)
      H02(4,4)=H2(4,4)*H2HP*ZIDAUX(4,4)
      H02(6,6)=H2(6,6)*H2P*ZIDAUX(6,6)
      H03(3,3)=H3(3,3)*H3N*ZIDAUX(3,3)
      H03(4,4)=H3(4,4)*H3P*ZIDAUX(4,4)
      H03(5,5)=H3(5,5)*H3P*ZIDAUX(5,5)
*      direction 4 : + 45   5 :-45
* 1 ere ligne
      H04(1,1)=H4(1,1)*H4N*ZIDAUX(1,1)
      H05(1,1)=H5(1,1)*H4N*ZIDAUX(1,1)
      H04(1,6)=H4(1,6)*H4P*ZIDAUX(6,6)
      H05(1,6)=H5(1,6)*H4P*ZIDAUX(6,6)
* 2 eme ligne
      H04(2,2)=H4(2,2)*H4N*ZIDAUX(1,1)
      H05(2,2)=H5(2,2)*H4N*ZIDAUX(1,1)
      H04(2,6)=H4(2,6)*H4P*ZIDAUX(6,6)
      H05(2,6)=H5(2,6)*H4P*ZIDAUX(6,6)
* 4 eme ligne
      H04(4,4)=H4(4,4)*H4HP*ZIDAUX(5,5)
      H05(4,4)=H5(4,4)*H4HP*ZIDAUX(5,5)
      H04(4,5)=H4(4,5)*H4HP*ZIDAUX(5,5)
      H05(4,5)=H5(4,5)*H4HP*ZIDAUX(5,5)
* 5 eme ligne
      H04(5,4)=H4(5,4)*H4HP*ZIDAUX(5,5)
      H05(5,4)=H5(5,4)*H4HP*ZIDAUX(5,5)
      H04(5,5)=H4(5,5)*H4HP*ZIDAUX(5,5)
      H05(5,5)=H5(5,5)*H4HP*ZIDAUX(5,5)
* 6 eme ligne
      H04(6,1)=H4(6,1)*H4P*ZIDAUX(6,6)
      H05(6,1)=H5(6,1)*H4P*ZIDAUX(6,6)
      H04(6,2)=H4(6,2)*H4P*ZIDAUX(6,6)
      H05(6,2)=H5(6,2)*H4P*ZIDAUX(6,6)
      H04(6,6)=H4(6,6)*H4P*ZIDAUX(6,6)
      H05(6,6)=H5(6,6)*H4P*ZIDAUX(6,6)
*
       IF (IIIPARAM.EQ.1) THEN
       WRITE (IOIMP,*) 'Matrice H01'
       DO 250 ILOOP= 1,LHOOKL
          WRITE (IOIMP, 99999 ) (H01(ILOOP,J), J=1,LHOOKL)
250    CONTINUE
       WRITE (IOIMP,*) 'Matrice H02'
       DO 251 ILOOP= 1,LHOOKL
          WRITE (IOIMP, 99999 ) (H02(ILOOP,J), J=1,LHOOKL)
251    CONTINUE
       WRITE (IOIMP,*) 'Matrice H03'
       DO 252 ILOOP= 1,LHOOKL
          WRITE (IOIMP, 99999 ) (H03(ILOOP,J), J=1,LHOOKL)
252    CONTINUE
       ENDIF
***** Recuperation deformation initiale dans le repere orth.
***** EN plus les VI sont rangées dans l'ordre de l'ONERA cad
****** 12 ==> 6,    13 ==> 5 et 23 ==> 4
        DEFI(1)=VAR0(7)
        DEFI(2)=VAR0(8)
        DEFI(3)=VAR0(9)
        DEFI(4)=VAR0(10)
        DEFI(5)=VAR0(11)
        DEFI(6)=VAR0(12)
**  repère global
       IF (IIIPARAM.EQ.1) THEN
           WRITE (IOIMP,*) 'DEFORMATION RECUPEREE'
           WRITE (IOIMP,99999) (DEFI(ILOOP), ILOOP=1,6)
      ENDIF
* Calcul de la nouvelle défo repère glob ordre cast
      DO 451 ILOOP=1,6
        EPSTFV(ILOOP)=EPST0(ILOOP)+DEPST(ILOOP)
*        EPSTF(ILOOP)=EPSTFV(ILOOP)
 451   CONTINUE
*           WRITE (IOIMP,*) 'DEF  IMPOSEE glob cast3m'
*           WRITE (IOIMP,99999) (EPSTFV(ILOOP), ILOOP=1,6)
* Par contre les depstv cad après les DORTH sont rangés dans l'ordre castem
*  : 12 ==> 4 13 ==> 5 23 ==>6   -
      CALL CICROT (wrk52,wrk53,wrk54,1,EPSTFV,DORTH)
*
* On les réorganise donc pour etre dans l'ordre ONERA
      VARTMP= DORTH(4)
      DORTH(4)= DORTH(6)
      DORTH(6)= VARTMP
*
*     Calcul de la nouvelle deformation totale dans le repère ortho
*     Utilisées pour le calcul des valeurs propres
      DEFOTH(1,1) =  DORTH(1)
      DEFOTH(2,2) =  DORTH(2)
      DEFOTH(3,3) =  DORTH(3)
      DEFOTH(1,2) = DORTH(6)
      DEFOTH(1,3) =  DORTH(5)
      DEFOTH(2,1) = DEFOTH(1,2)
      DEFOTH(2,3) = DORTH(4)
      DEFOTH(3,1) = DEFOTH(1,3)
      DEFOTH(3,2) = DEFOTH(2,3)
*
       CALL JACOBA(DEFOTH,3,PROP,XPROP)
C
C  partie positive des valeurs propres
       IF(PROP(1).LT.0.) PROP(1)=0.D0
       IF(PROP(2).LT.0.) PROP(2)=0.D0
       IF(PROP(3).LT.0.) PROP(3)=0.D0
C Je remet ma matrice de valeurs propres corrigée axes ortho
      DEFOTH(1,1) = PROP(1)
      DEFOTH(2,2) = PROP(2)
      DEFOTH(3,3) = PROP(3)
      DEFOTH(1,2) = 0.D0
      DEFOTH(1,3) = 0.D0
      DEFOTH(2,1) = 0.D0
      DEFOTH(2,3) = 0.D0
      DEFOTH(3,1) = 0.D0
      DEFOTH(3,2) = 0.D0
C XPROP(3,3) est la matrice donnant les coordonnées des vecteurs propre
C par rapport aux axes d'orthotropie
C 1 er indice direction 2° n° vecteur propre
C j'ai besoin de la matrice donnant les axes ortho dans le repère des vecteurs propres
C           PRINT *,'AVANT INVER3'
           CALL INVER3(XPROP,XPINV)
C           PRINT *,'APRES INVER3'
           CALL TRSPOD(XPINV,3,3,XPINVT)
           CALL MULMAT(DEFOTHT,DEFOTH,XPINV,3,3,3)
           CALL MULMAT(DEFOTH,XPINVT,DEFOTHT,3,3,3)
C
C en principe on a la matrice de defo efficasse dans les axes ortho defoth
C calcul des driving forces dans le cas scalaire
C mode normal
C calcul des driving forces dans le cas pseudo-tensoriel
 
          Y1 =0.5*(DDAUX(1,1)*DEFOTH(1,1)*DEFOTH(1,1)+b1*
     &         DDAUX(6,6)*DEFOTH(1,2)*DEFOTH(1,2)+b2*
     &         DDAUX(5,5)*DEFOTH(1,3)*DEFOTH(1,3))
          Y2 =0.5*(DDAUX(2,2)*DEFOTH(2,2)*DEFOTH(2,2)+b1*
     &         DDAUX(6,6)*DEFOTH(1,2)*DEFOTH(1,2)+b2*
     &         DDAUX(4,4)*DEFOTH(2,3)*DEFOTH(2,3))
          Y6 =0.25*(DDAUX(1,1)*DEFOTH(1,1)*DEFOTH(1,2)+
     &         DDAUX(2,2)*DEFOTH(2,2)*DEFOTH(1,2)+b1*
     &         DDAUX(6,6)*DEFOTH(1,2)*
     &         (DEFOTH(1,1)+DEFOTH(2,2)) )
C driving forces utilisées pour le calcul des dommages
          Y1 = Y1 - (ABS(Y6))
          Y2 = Y2 - (ABS(Y6))
          IF(Y6.GT.0.) THEN
          Y(4)=SQRT(Y6)
          Y(5)= 0.
          ELSE
          Y(4)=0.
          Y(5)=SQRT(-Y6)
          ENDIF
          Y3 = 0.5*(DDAUX(3,3)*DEFOTH(3,3)*DEFOTH(3,3)+b3*
     &         DDAUX(4,4)*DEFOTH(2,3)*DEFOTH(2,3)+b3*
     &         DDAUX(5,5)*DEFOTH(1,3)*DEFOTH(1,3))
C
 
          IF(Y1.LT.0.) THEN
           Y(1)=0.
          ELSE
           Y(1)=SQRT(Y1)
          ENDIF
          IF(Y2.LT.0.) THEN
           Y(2)=0.
          ELSE
           Y(2)=SQRT(Y2)
          ENDIF
          IF(Y3.LT.0.) THEN
           Y(3)=0.
          ELSE
           Y(3)=SQRT(Y3)
          ENDIF
C Calcul des variables d'endomagenement
        Do 722 IDRF=1,5
        FAFA =Y(IDRF)-Y0(IDRF)
        IF (FAFA.LT.0.) FAFA=0.
        DOM(IDRF)= DC(IDRF)*(1.D0 - EXP(-1.*
     &        ((FAFA/YC(IDRF))**PY(IDRF))))
 722    CONTINUE
      DOM1= MAX(DOM(1),VAR0(2))
      DOM2= MAX(DOM(2),VAR0(3))
      DOM3= MAX(DOM(3),VAR0(4))
      DOM4= MAX(DOM(4),VAR0(5))
      DOM5= MAX(DOM(5),VAR0(6))
C dans WRK52 (DECHE ) on a tur0(1) et turf(1)
 
      IF (IIIPARAM.EQ.1) THEN
        WRITE (IOIMP,*) 'Increment des deformations orth. '
        WRITE (IOIMP,99999) (DORTH(ILOOP), ILOOP=1,6)
      ENDIF
*
      DO 500 ILOOP=1,iecou.NSTRSS
*       On garde les déformations totales dans les var. internes
        VARF(ILOOP+6)=DORTH(ILOOP)
500   CONTINUE
*
* On peut calculer EPSILON CLosure
* EN fonction de deltal et tzero
* et de la temperature TURE0(1)ou TUREF(1)du segment WRK52 de l'include DECHE
*
      DO 701 IDOU=1,6
       EPSCLO(IDOU) = DELTAL*(TURE0(1)-TZERO)
       EPSBAR(IDOU) = DORTH(IDOU)- EPSCLO(IDOU)
       EPSIM(IDOU) = (DORTH(IDOU)+VAR0(IDOU+6))/2.
 701  CONTINUE
* On redefini les epsbar pour les direction 45 et -45 ° (cf 4.2.8)
       EPSBAR(4) = (EPSBAR(1)+EPSBAR(2)+(2.D0*EPSBAR(6)))*0.5
       EPSBAR(5) = (EPSBAR(1)+EPSBAR(2)-(2.D0*EPSBAR(6)))*0.5
 
* Calcul des valeurs délimitant les  déformations pour
* lesquelles les fissures commencent et finissent à/de se fermer
* en vue du calcul de l'index de désactivation
      DEPSIF(1)= (1.D0 + (AIF(1)*DOM1))*DEPSIF0(1)
      DEPSIF(2)= (1.D0 + (AIF(2)*DOM2))*DEPSIF0(2)
      DEPSIF(3)= (1.D0 + (AIF(3)*DOM3))*DEPSIF0(3)
      DEPSIF(4)= (1.D0 + (AIF(4)*DOM4))*DEPSIF0(4)
      DEPSIF(5)= (1.D0 + (AIF(5)*DOM5))*DEPSIF0(5)
*
* Calcul de l'index d'activation 1 si actif 0 si passif
*
      DO 703 IDO2=1,5
       IF(EPSBAR(IDO2).GE.DEPSIF(IDO2)) THEN
       ZNUACT(IDO2)= 1.D0
       ELSEIF(EPSBAR(IDO2).GT.(-1.*(DEPSIF(IDO2)))) THEN
       ZNUACT(IDO2)= 0.5D0*(1.D0-(COS(0.5D0*XPI*(EPSBAR(IDO2)
     &              +DEPSIF(IDO2))/DEPSIF(IDO2))))
       ELSE
       ZNUACT(IDO2)= 0.D0
       ENDIF
 703   CONTINUE
*
*  Calcul des déformations résiduelles et stockées
*
* Calcul des variations du dommage
      DELTDO1 = DOM1 - VAR0(2)
      DELTDO2 = DOM2 - VAR0(3)
      DELTDO3 = DOM3 - VAR0(4)
      DELTDO4 = DOM4 - VAR0(5)
      DELTDO5 = DOM5 - VAR0(6)
* Calcul des  variations d'indice d'activation
      DELTNU1 = ZNUACT(1) - VAR0(13)
      DELTNU2 = ZNUACT(2) - VAR0(14)
      DELTNU3 = ZNUACT(3) - VAR0(15)
      DELTNU4 = ZNUACT(4) - VAR0(16)
      DELTNU5 = ZNUACT(5) - VAR0(17)
*
* Calcul de Seff =somme sur i de nui di hi0
*
       DO 744 INDJ4=1,6
         Do 745 INDJ3=1,6
           SEFF(INDJ3,INDJ4) = ZIDAUX(INDJ3,INDJ4)+
     &      ZNUACT(1)*DOM1*H01(INDJ3,INDJ4)+
     &      ZNUACT(2)*DOM2*H02(INDJ3,INDJ4)+
     &      ZNUACT(3)*DOM3*H03(INDJ3,INDJ4)+
     &      ZNUACT(4)*DOM4*H04(INDJ3,INDJ4)+
     &      ZNUACT(5)*DOM5*H05(INDJ3,INDJ4)
           CEFF(INDJ3,INDJ4) = SEFF(INDJ3,INDJ4)
 745     CONTINUE
 744   CONTINUE
* inversion de la matrice de rigidité
      TPREC= 0.D1
      IPERR=0
      CALL INVALM(CEFF,LHOOKL,LHOOKL,IPERR,TPREC)
      IF (IPERR.NE.0) THEN
        WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK'
      ENDIF
      CALL MULMAT (TERME,H01,CEFF,6,6,6)
      CALL MULMAT (ZK1,CEFF,TERME,6,6,6)
      CALL MULMAT (TERME,H02,CEFF,6,6,6)
      CALL MULMAT (ZK2,CEFF,TERME,6,6,6)
      CALL MULMAT (TERME,H03,CEFF,6,6,6)
      CALL MULMAT (ZK3,CEFF,TERME,6,6,6)
      CALL MULMAT (TERME,H04,CEFF,6,6,6)
      CALL MULMAT (ZK4,CEFF,TERME,6,6,6)
      CALL MULMAT (TERME,H05,CEFF,6,6,6)
      CALL MULMAT (ZK5,CEFF,TERME,6,6,6)
*
*   Déformation stockées (S) + résiduelle (R) à mettre dans les VI
*
       DO 741 INDI=1,6
       SISIS(INDI)=0.D0
       SISIR(INDI)=0.D0
        DO 742 INDJ=1,6
        SOMS(INDI,INDJ)= DELTNU1*DOM1*ZK1(INDI,INDJ)+
     $              DELTNU2*DOM2*ZK2(INDI,INDJ)+
     $              DELTNU3*DOM3*ZK3(INDI,INDJ)+
     $              DELTNU4*DOM4*ZK4(INDI,INDJ)+
     $              DELTNU5*DOM5*ZK5(INDI,INDJ)
        SISIS(INDI)=SISIS(INDI)-SOMS(INDI,INDJ)*EPSIM(INDJ)
        SOMR(INDI,INDJ)= ETA1*ZNUACT(1)*DELTDO1*ZK1(INDI,INDJ)+
     $              ETA2*ZNUACT(2)*DELTDO2*ZK2(INDI,INDJ)+
     $              ETA3*ZNUACT(3)*DELTDO3*ZK3(INDI,INDJ)+
     $              ETA4*ZNUACT(4)*DELTDO4*ZK4(INDI,INDJ)+
     $              ETA5*ZNUACT(5)*DELTDO5*ZK5(INDI,INDJ)
        SISIR(INDI)=SISIR(INDI)+SOMR(INDI,INDJ)*EPSIM(INDJ)
 
  742   CONTINUE
  741   CONTINUE
C
      CALL MULMAT (EPSR1,ZIDAUX,SISIR,6,1,6)
      CALL MULMAT (EPSS1,ZIDAUX,SISIS,6,1,6)
       DO 743 INDI=1,6
        EPSIN(INDI)=VAR0(17+INDI)+EPSR1(INDI)+EPSS1(INDI)
        EPSEL(INDI) = DORTH(INDI)-EPSIN(INDI)
  743   CONTINUE
*
* SORTH contient les contraintes a la fin du pas
* dans le repère orthotrope
       CALL MULMAT (SORTH1,CEFF,DORTH,6,1,6)
       CALL MULMAT (SORTH2,DDAUX,EPSIN,6,1,6)
       DO 749 INDIA=1,6
       SORTH(INDIA) = SORTH1(INDIA) - SORTH2(INDIA)
 749   CONTINUE
      IF (IIIPARAM.EQ.1) THEN
         WRITE (IOIMP,*) 'Contraintes calculees (repere orth.)'
         WRITE (IOIMP,99999) (SORTH(ILOOP), ILOOP=1,6)
      ENDIF
* On réorganise les contraintes en cisaillement pour
* calculer les contraintes dans le repere global (ON REMET DANS L'ORDRE castem)
      SIGFV(1)=SORTH(1)
      SIGFV(2)=SORTH(2)
      SIGFV(3)=SORTH(3)
      SIGFV(4)=SORTH(6)
      SIGFV(5)=SORTH(5)
      SIGFV(6)=SORTH(4)
      EPSINV(1)=EPSIN(1)
      EPSINV(2)=EPSIN(2)
      EPSINV(3)=EPSIN(3)
      EPSINV(4)=EPSIN(6)
      EPSINV(5)=EPSIN(5)
      EPSINV(6)=EPSIN(4)
      VARF(1)= 0.D0
      DO 186 INDI = 1,6
        VARF(17+INDI)= EPSIN(INDI)
        VARF(1)= VARF(1)+ EPSIN(INDI)
 186  CONTINUE
* On appele SICROT pour trouver les contraintes dans le repère global
* On utilise SORTH pour garder le resultat
* (on peut pas passer SIGF à la subroutine car il est dans un segment)
      CALL CICROT (wrk52,wrk53,wrk54,0,SIGFV,SORTH)
 
* On recopie SORTH in SIGF
      DO 1000 ILOOP=1,6
         SIGF(ILOOP)=SORTH(ILOOP)
1000  CONTINUE
       IF (IIIPARAM.EQ.1) THEN
         WRITE (IOIMP,*) 'SIG  REP GLOB cast'
         WRITE (IOIMP,99999) (SIGF(ILOOP), ILOOP=1,6)
       ENDIF
      VARF(2)=  DOM1
      VARF(3)=  DOM2
      VARF(4)=  DOM3
      VARF(5)=  DOM4
      VARF(6)=  DOM5
      VARF(13)=  ZNUACT(1)
      VARF(14)=  ZNUACT(2)
      VARF(15)=  ZNUACT(3)
      VARF(16)=  ZNUACT(4)
      VARF(17)=  ZNUACT(5)
99999    format(2x,'  ',(6(1x,3pe12.5),/))
99998    format(2x,'  ',(4(1x,3pe12.5),/))
99997    format(2x,'  ',(1x,3pe12.5),/)
99996    format(2x,'  ',(3(1x,1pe12.3),/))
      RETURN
      END