C BETOKT    SOURCE    PV        11/03/07    21:15:10     6885           
      SUBROUTINE BETOKT(XMAT,TRAC,VARF,SIGF,NCURVT,NCURVC,XKTAT,KERRE)
C
C   Matrice tangente pour le modele endommageable en cyclique
C   D. Combescure 18/09/95 ELSA-Ispra
C
C ENTREES:
C -------
C     XMAT    = Caracteristiques materiaux
C
C     VARF    = Variables internes
C
C     SIGF    = Contraintes
C
C SORTIES:
C -------
C     XKTAT    = Matrice tangente
C
C================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC CCOPTIO
         DIMENSION XMAT(*)
         DIMENSION TRAC(*)
C
      DIMENSION VARF(14),SIGF(4),EPSI(2),FG(2)
      DIMENSION XNOR1(4),XNOR2(4),XCOU1(3,3),XCOU2(3,3)
      DIMENSION GAMBDA(2),EEE(2,2),SIGP(2),EEE2(2,2)
      DIMENSION UUU(3,2),VVV(2,3),UUU2(3,2),VVV2(2,3),UUU3(3,2)
      DIMENSION HHH1(3,3),HHH2(3,3),HHH3(3,3)
      DIMENSION XKTAT(4,4),RESIST(6),XKISO(2)
      PARAMETER (UN=1.D0,UNDEMI=0.5D0,EPSILO=1.D-20,XZER=0.D0)
      PARAMETER (PRECIS=1.D-4)
      KERRE=0
C================================================================
C   LECTURE DES CARACTERISTIQUES MAT. ET DES VAR. INT.
C===============================================================
C Pour la lecture des courbes
      IPT    = 1
      IPC    = IPT+2*NCURVT
C
      YOUN = XMAT(1)
      XNU  = XMAT(2)
      XGGG = (UNDEMI*YOUN/(UN+XNU))
      XHH1 = XMAT(3)
      DEGRAD    = XMAT(8)/XMAT(9)
      RESIST(1) = XMAT(4)
      RESIST(2) = XMAT(4)
      RESIST(3) = XMAT(5)
      RESIST(4) = XMAT(5)
      RESIST(5) = XMAT(6)
      RESIST(6) = XMAT(7)
      XLCAT     = XMAT(10)
      XLCAC     = XMAT(11) 
      CPHI    = VARF(5) + UN 
      SPHI    = VARF(6)
      EPSI(1) = VARF(3)
      EPSI(2) = VARF(4)
      ICAS1 = nint(VARF(11))
      ICAS2 = nint(VARF(12))
C==============================================================
C    RAIDEUR ELASTIQUE
C===============================================================
      CALL ZERO(XKTAT,4,4)
      XKTAT(1,1) = YOUN/(UN - (XNU**2))
      XKTAT(1,2) = XNU*XKTAT(1,1)
      XKTAT(2,1) = XKTAT(1,2) 
      XKTAT(2,2) = XKTAT(1,1)
      XKTAT(4,4) = YOUN/(2.D0*(UN+XNU))
C==============================================================
      IF (ICAS1.EQ.0) THEN
         GOTO 9999
      ENDIF 
C
      CALL VALEUR(EPSI(1),TRAC(IPT),NCURVT,XLCAT,
     &XKISO(1),XPEN1,KERRE)
      CALL VALEUR(EPSI(2),TRAC(IPC),NCURVC,XLCAC,
     &XKISO(2),XPEN2,KERRE)
      GAMBDA(1) = VARF(13)
      GAMBDA(2) = VARF(14)
      IF (ICAS1.LE.4) THEN
         FG(1) = RESIST(ICAS1)
      ELSE
         FG(1) = XKISO(ICAS1-4)*RESIST(ICAS1)
      ENDIF
      IF (ICAS2.LE.4) THEN
         FG(2) = RESIST(ICAS2)
      ELSE
         FG(2) = XKISO(ICAS2-4)*RESIST(ICAS2)
      ENDIF
C==============================================================
C   CALCUL DE L'INVERSE
C==============================================================
        CALL ZERO(HHH1,3,3)
        HHH1(1,1) = UN/YOUN
        HHH1(1,2) = -XNU*HHH1(1,1)
        HHH1(2,1) = HHH1(1,2)
        HHH1(2,2) = HHH1(1,1)
        HHH1(3,3) = (2.D0*(UN+XNU)/YOUN)
C
        CALL CONPRI (SIGF,CPHI,SPHI,SIGP,CPHI,SPHI)
C
        IF (ICAS1.EQ.1) THEN
           ISGNTC = 1
        ELSE IF (ICAS1.EQ.2) THEN
           ISGNTC = 1
        ELSE IF (ICAS1.EQ.3) THEN
           ISGNTC = -1
        ELSE IF (ICAS1.EQ.4) THEN
           ISGNTC = -1
        ELSE IF (ICAS1.EQ.5) THEN
           ISGNTC = 1
        ELSE IF (ICAS1.EQ.6) THEN
           ISGNTC = -1
        ENDIF
        IF (ICAS1.LE.4) THEN
          IF (SIGP(2).GE.SIGP(1)) THEN
            IF ((ICAS1.EQ.1).OR.(ICAS1.EQ.3)) THEN
              ISGN2 = -1
            ELSE
              ISGN2 = 1
            ENDIF
          ELSE
            IF ((ICAS1.EQ.1).OR.(ICAS1.EQ.3)) THEN
              ISGN2 = 1
            ELSE
              ISGN2 = -1
            ENDIF         
          ENDIF
        ELSE
          ISGN2 = 1
        ENDIF 
        IF (ICAS1.LE.4) THEN
          CALL COURBU(SIGF,FG(1),ISGNTC,ISGN2,1,XCOU1)
          CALL NORMBL(SIGF,FG(1),ISGNTC,ISGN2,1,XNOR1)
        ELSE
          PHI = SQRT ((UNDEMI*(SIGF(1)-SIGF(2)))**2+SIGF(4)**2)
          IF ((PHI.LT.ABS(PRECIS*FG(1))).AND.(ICAS1.EQ.ICAS2)) THEN
C            Cas de l'apex ou on calcule directement H-1
            HHH1(1,1) = XGGG*((4.D0*PHI)+(YOUN*GAMBDA(1)))/
     & (2.D0*(UN-XNU)*((2.D0*PHI)+(XGGG*GAMBDA(1))))
            HHH1(2,2)=HHH1(1,1)
            HHH1(1,2) = XGGG*((4.D0*XNU*PHI)+(YOUN*GAMBDA(1)))/
     & (2.D0*(UN-XNU)*((2.D0*PHI)+(XGGG*GAMBDA(1))))
            HHH1(2,1)=HHH1(1,2)
            HHH1(3,3) = XGGG*PHI/(PHI+XGGG*GAMBDA(1))
            CALL NORMBL(SIGF,FG(1),ISGNTC,1,2,XNOR1)
            GOTO 1010
          ELSE
            CALL COURBU(SIGF,FG(1),ISGNTC,1,2,XCOU1)
            CALL NORMBL(SIGF,FG(1),ISGNTC,1,2,XNOR1)
          ENDIF
        ENDIF
C
C
        IF (ICAS1.EQ.ICAS2) THEN
           GOTO 1000
        ELSE
          IF (ICAS2.LE.4) THEN
             FG(2) = RESIST(ICAS2)
          ELSE
             FG(2) = XKISO(ICAS2-4)*RESIST(ICAS2)
          ENDIF
          IF (ICAS2.EQ.1) THEN
             ISGNTCB = 1
          ELSE IF (ICAS2.EQ.2) THEN
             ISGNTCB = 1
          ELSE IF (ICAS2.EQ.3) THEN
             ISGNTCB = -1
          ELSE IF (ICAS2.EQ.4) THEN
             ISGNTCB = -1
          ELSE IF (ICAS2.EQ.5) THEN
             ISGNTCB = 1
          ELSE IF (ICAS2.EQ.6) THEN
             ISGNTCB = -1
          ENDIF
          IF (ICAS2.LE.4) THEN
            IF (SIGP(2).GE.SIGP(1)) THEN
              IF ((ICAS2.EQ.1).OR.(ICAS2.EQ.3)) THEN
                ISGN2B = -1
              ELSE
                ISGN2B = 1
              ENDIF
            ELSE
              IF ((ICAS2.EQ.1).OR.(ICAS2.EQ.3)) THEN
                ISGN2B = 1
              ELSE
                ISGN2B = -1
              ENDIF
            ENDIF
          ELSE
            ISGN2B = 1
          ENDIF
C
          IF (ICAS2.LE.4) THEN
            CALL COURBU(SIGF,FG(2),ISGNTCB,ISGN2B,1,XCOU2)
            CALL NORMBL(SIGF,FG(2),ISGNTCB,ISGN2B,1,XNOR2)
          ELSE
            CALL COURBU(SIGF,FG(2),ISGNTCB,1,2,XCOU2)
            CALL NORMBL(SIGF,FG(2),ISGNTCB,1,2,XNOR2)
          ENDIF
C
           GOTO 2000
        ENDIF
C
C
C
 1000 CONTINUE 
C
C  Cas 1 mecanisme
C
C  Calcul de l'inverse de H et de H 
      CALL PLMATR(HHH1,UN,XCOU1,GAMBDA(1),3,3,HHH1)
      CALL INVALM(HHH1,3,3,IRD,1.D-8)
      IF (IRD.NE.0) THEN
        WRITE (IOIMP,*)'BETOKT: Determinant nul'
        GOTO 9999
      ENDIF
C  Point d'entree dans le cas de l'apex
 1010 CONTINUE
C  Calcul de U,V et E-1
      CALL ZERO(UUU,3,2)
      CALL ZERO(VVV,2,3)
      CALL ZERO(UUU2,3,2)
      CALL ZERO(VVV2,2,3)
      UUU(1,1) = XNOR1(1)
      UUU(2,1) = XNOR1(2)
      UUU(3,1) = XNOR1(4)
      VVV(1,1) = XNOR1(1)
      VVV(1,2) = XNOR1(2)
      VVV(1,3) = XNOR1(4)
C
      IF (ICAS1.LE.4) THEN
        EEE(1,1) = XHH1
      ELSE IF (ICAS1.EQ.5) THEN 
        EEE(1,1) = ISGNTC*XPEN1*RESIST(ICAS1)
      ELSE IF (ICAS1.EQ.6) THEN
        EEE(1,1) = ISGNTC*XPEN2*RESIST(ICAS1)
      ENDIF
C
      CALL MATMAT(HHH1,UUU,3,3,2,UUU2)
      CALL MATMAT(VVV,HHH1,2,3,3,VVV2)
C
      CALL MATMAT(VVV,UUU2,2,3,2,HHH2)
C
      EEE(1,1) = EEE(1,1) + HHH2(1,1)
      IF (ABS(EEE(1,1)).GT.EPSILO) THEN
        EEE(1,1) = -UN/EEE(1,1) 
      ELSE
        WRITE (IOIMP,*)'BETOKT: Determinant nul'
        GOTO 9999
      ENDIF
      CALL MATMAT(UUU2,VVV2,3,2,3,HHH3)
      CALL PLMATR(HHH1,UN,HHH3,EEE(1,1),3,3,HHH1)
C
      GOTO 3000 
C
 2000 CONTINUE
C
C  Cas 2 mecanismes
C
C Calcul de H-1 et H
      CALL PLMATR(HHH1,UN,XCOU1,GAMBDA(1),3,3,HHH1)
      CALL PLMATR(HHH1,UN,XCOU2,GAMBDA(2),3,3,HHH1)
      CALL INVALM(HHH1,3,3,IRD,1.D-8)
      IF (IRD.NE.0) THEN
        WRITE (IOIMP,*)'BETOKT: Determinant nul'
        GOTO 9999
      ENDIF
C Calcul de U ,V et E
      UUU(1,1) = XNOR1(1)
      UUU(2,1) = XNOR1(2)
      UUU(3,1) = XNOR1(4)
      UUU(1,2) = XNOR2(1)
      UUU(2,2) = XNOR2(2)
      UUU(3,2) = XNOR2(4)
      VVV(1,1) = XNOR1(1)
      VVV(1,2) = XNOR1(2)
      VVV(1,3) = XNOR1(4)
      VVV(2,1) = XNOR2(1)
      VVV(2,2) = XNOR2(2)
      VVV(2,3) = XNOR2(4)
C
      IF (ICAS1.LE.4) THEN
        EEE(1,1) = XHH1
        IF (ICAS2.EQ.6) THEN
          EEE(2,1) = DEGRAD*ISGNTCB*RESIST(6)
     &*ISGNTC*RESIST(ICAS1)*XPEN2
        ELSE
          EEE(2,1) = XZER
        ENDIF
      ELSE IF (ICAS1.EQ.5) THEN
        EEE(1,1) = ISGNTC*XPEN1*RESIST(ICAS1)
        EEE(2,1) = XZER
      ELSE IF (ICAS1.EQ.6) THEN
        EEE(1,1) = ISGNTC*XPEN2*RESIST(ICAS1)
        EEE(2,1) = XZER
      ENDIF
      IF (ICAS2.LE.4) THEN
        EEE(2,2) = XHH1
        IF (ICAS1.EQ.6) THEN
          EEE(1,2) = DEGRAD*ISGNTC*RESIST(6)
     &*ISGNTCB*RESIST(ICAS2)*XPEN2
        ELSE
          EEE(1,2) = XZER
        ENDIF
      ELSE IF (ICAS2.EQ.5) THEN
        EEE(2,2) = ISGNTCB*XPEN1*RESIST(ICAS2)
        EEE(1,2) = XZER
      ELSE IF (ICAS2.EQ.6) THEN
        EEE(2,2) = ISGNTCB*XPEN2*RESIST(ICAS2)
        EEE(1,2) = XZER
      ENDIF
C
      CALL MATMAT(HHH1,UUU,3,3,2,UUU2)
      CALL MATMAT(VVV,HHH1,2,3,3,VVV2)
C
      CALL MATMAT(VVV,UUU2,2,3,2,HHH2)
C
      CALL PLMATR(EEE,UN,HHH2,UN,2,2,EEE)
      XDET = EEE(1,1)*EEE(2,2) - EEE(1,2)*EEE(2,1)
      IF (ABS(XDET).GT.EPSILO) THEN
        XDET = UN/XDET
        EEE2(1,1) = XDET*EEE(2,2)
        EEE2(2,2) = XDET*EEE(1,1)
        EEE2(2,1) = -XDET*EEE(1,2)
        EEE2(1,2) = -XDET*EEE(2,1)
      ELSE
        WRITE (IOIMP,*)'BETOKT: Determinant nul'
        GOTO 9999
      ENDIF
      CALL MATMAT(UUU2,EEE2,3,2,2,UUU3)
      CALL MATMAT(UUU3,VVV2,3,2,3,HHH3)
      CALL PLMATR(HHH1,UN,HHH3,-UN,3,3,HHH1)
C
      GOTO 3000
C
C==============================================================
 3000   CONTINUE
        CALL ZERO(XKTAT,4,4)
        XKTAT(1,1) = HHH1(1,1)
        XKTAT(1,2) = HHH1(1,2)
        XKTAT(2,1) = HHH1(2,1)
        XKTAT(2,2) = HHH1(2,2)
        XKTAT(4,4) = HHH1(3,3)
        XKTAT(4,1) = HHH1(3,1)
        XKTAT(4,2) = HHH1(3,2)
        XKTAT(1,4) = HHH1(1,3)
        XKTAT(2,4) = HHH1(2,3)
 9999 CONTINUE
C===============================================================
C                 Fin de la routine
C===============================================================
      RETURN
      END