C BEHAV1    SOURCE    CB215821  16/04/21    21:15:17     8920
      SUBROUTINE BEHAV1(SIGR,DSTRN,DEPS1,DEPS2,IPLA,
     &   SIG1,SIG2,IFIS,SIGF,DSIGT,NSTRS,IFOUB,DEP,SIGRV,SIGP,
     &   BETJEF,VISCO,NECH0,NECH1)
C
C ==================================================================
C
C             MODELE DE PLASTICITE EN TRAITEMENT POST PIC
C                Un seul critere de traction: Rankine
C
C ==================================================================
C CE SOUS-PROGRAMME EST APPELE DANS "BONE".
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DFSIG(4),SIGF(4),DSIGT(4),VEC1(4),DEPSI(4),SIGRV(4)
      DIMENSION V1(4),AC(4),SIGE(4),SIGP(4),SIGR(4),DSTRN(4)
      DIMENSION D2FSIG(4,4),DEP(4,4),P1(4,4),D(4,4),DP(4,4),A(4,4)
      DIMENSION AI(4,4),AH(4,4)
*
      SEGMENT BETJEF
       REAL*8 AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
     &          TCON,DPSTF1,DPSTF2,TETA,PDT,TP0
       INTEGER ICT,ICC,IMOD,IVIS,ITR,
     &                ISIM,IBB,IGAU,IZON
      ENDSEGMENT
*
      SEGMENT VISCO
       REAL*8 DPSTV1,DPSTV2,SIGV1,SIGV2,ENDV
      ENDSEGMENT
      SEGMENT NECH0
       REAL*8 DT,DC,ALFG,S0,ENDO
      ENDSEGMENT
      SEGMENT NECH1
       REAL*8 ENDL
      ENDSEGMENT
C
*     COMMON /DBETJEF/AA,BETA,COLI,PALF,YOUN,XNU,GFC,GFT,CAR,ETA,TDEF,
*    &              TCON,DPSTF1,DPSTF2,TETA,PDT,ICT,ICC,IMOD,IVIS,ITR,
*    &                ISIM,IBB,IGAU,IZON
C
C
      IAPEX=0
      PRB=1.D-10
      PRB2=1.D-6
      ITER=1
      ITANG=0
      IBROY=0
      Ft=PALF*COLI
      CRIMAX=0.D0
      SEQ = 0.D0
      SEQ1 = 0.D0
      CALL ZERO(SIGE,4,1)
      CALL ZERO(V1,4,1)
      CALL ZERO(P1,4,4)
      CALL ZERO(D,4,4)
C
      DO 10 I=1,NSTRS
      SIGE(I)=SIGF(I)
   10 CONTINUE
C
C      DO 11 I=1,NSTRS
C      DO 11 J=1,NSTRS
C      D(I,J)=0.D0
C   11 CONTINUE
C
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
      AD=YOUN/(1.D0-XNU*XNU)
      D(1,1)=AD
      D(2,2)=D(1,1)
      D(3,3)=AD*(1.D0-XNU)/2.D0
      D(1,2)=AD*XNU
      D(2,1)=D(1,2)
      ENDIF
C
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
      ADD=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
      D(1,1)=ADD*(1.D0-XNU)
      D(2,2)=D(1,1)
      D(3,3)=D(1,1)
      D(1,2)=ADD*XNU
      D(2,1)=D(1,2)
      D(1,3)=D(1,2)
      D(2,3)=D(1,2)
      D(3,1)=D(1,2)
      D(3,2)=D(1,2)
      D(4,4)=0.5*ADD*(1.D0-2.D0*XNU)
      ENDIF
C
C ************ Le point est fissure pour 1ere fois ***********
C
      IF (IFIS.EQ.0) THEN
        CALL PRINC(SIGF,V1,NSTRS)
        DEPS1=0.D0
        SIG1=Ft
        IF (V1(1).GT.Ft) THEN
         IFIS=1
        ENDIF
      ENDIF
      IF (IPLA.EQ.0) THEN
          DEPS2=0.D0
          SIG2=COLI*AA
      ENDIF
C
    8 CONTINUE
C
      DO 4 I=1,NSTRS
       SIGF(I)=SIGE(I)
    4 CONTINUE
C
      CALL PRINC(SIGF,V1,NSTRS)
      TETA=V1(4)
      DK=DEPS1
      DLAM0=0.D0
      FCRI0=V1(1)-SIG1
      CRIMAX=ABS(100.D0*FCRI0)
C
      IF (FCRI0.LT.0.D0) THEN
      DO 45 I=1,NSTRS
      DO 45 J=1,NSTRS
      DEP(I,J)=D(I,J)
   45 CONTINUE
        GOTO 100
      ENDIF
C
C ************ Traitement du point deja fissure  ********************
C
    9 CONTINUE
      PI=4.D0*ATAN(1.D0)
      PHIC=V1(4)*(PI/180.D0)
      COSA=COS(PHIC)
      SINA=SIN(PHIC)
C-------------------------------------------------------------------
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
          DFSIG(1)=COSA*COSA
          DFSIG(2)=SINA*SINA
          DFSIG(3)=2.D0*SINA*COSA
      DO 12 I=1,NSTRS
      DO 12 J=1,NSTRS
       P1(I,J)=0.D0
   12 CONTINUE
      P1(1,1)=1.D0/2.D0
      P1(1,2)=-1.D0/2.D0
      P1(2,1)=-1.D0/2.D0
      P1(2,2)=1.D0/2.D0
      P1(3,3)=2.D0
      ENDIF
C-------------------------------------------------------------------
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
          DFSIG(1)=COSA*COSA
          DFSIG(2)=SINA*SINA
          DFSIG(3)=0.D0
          DFSIG(4)=2.D0*SINA*COSA
      DO 13 I=1,NSTRS
      DO 13 J=1,NSTRS
       P1(I,J)=0.D0
   13 CONTINUE
      P1(1,1)=1.D0/2.D0
      P1(1,2)=-1.D0/2.D0
      P1(2,1)=-1.D0/2.D0
      P1(2,2)=1.D0/2.D0
      P1(4,4)=2.D0
      ENDIF
C-------------------------------------------------------------------
      DO 20 I=1,NSTRS
      AC(I)=0.D0
      DO 20 J=1,NSTRS
      AC(I)=AC(I)+D(I,J)*DFSIG(I)
   20 CONTINUE
C
      FDF=0.D0
      DO 30 J=1,NSTRS
      FDF=FDF+DFSIG(J)*AC(J)
   30 CONTINUE
C--------------- Determination du parametre d'ecrouissage ----------
C
      IF(IVIS.LE.2) THEN
       CALL UNICOU(DK,PAEC,1,SEQ,BETJEF)
      ELSE
       CALL UNICO1(DK,PAEC,1,SEQ,BETJEF,NECH0,NECH1)
      ENDIF
C
C-------------------------------------------------------------------
      DJAC0=-(PAEC+FDF)
C
C ************ Debut des iterations internes ***************
C
   40 CONTINUE
C
C *************** Determination de DK  *********************
C
      DLAM1=-FCRI0/DJAC0+DLAM0
      DK=DEPS1+DLAM1
C      IF (DLAM1.LT.0.D0) THEN
C       WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1
C       WRITE(*,*)'A l iteration :',ITER
C      ENDIF
C
C--------------- Estimation contrainte quivalente ----------
C
      IF(IVIS.LE.2) THEN
       CALL UNICOU(DK,PAEC,1,SEQ1,BETJEF)
      ELSE
       CALL UNICO1(DK,PAEC,1,SEQ1,BETJEF,NECH0,NECH1)
      ENDIF
C
C ************** Determination de DPHI *********************
C
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
      AD1=YOUN/(1.D0-XNU)
      TO=SIGF(3)
      ENDIF
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
      AD1=YOUN/((1.D0+XNU)*(1.D0-2.D0*XNU))
      TO=SIGF(4)
      ENDIF
C
      DPHI=SEQ1-0.5*(SIGE(1)+SIGE(2))+0.5*AD1*DLAM1
C      IF (ITER.EQ.1) THEN
C      DPHI=0.25*(SIGF(1)-SIGF(2))*(SIGF(1)-SIGF(2))
C      DPHI=DPHI+TO*TO
C      DPHI=SQRT(DPHI)
C      ENDIF
      IF (DPHI.LT.0.D0) THEN
C       WRITE(*,*)'ATTENTION DPHI NEGATIF'
C       WRITE(*,*)'DPHI=',DPHI
      ENDIF
C
C--------------- Cas de l'apex -----------------------------
C
      IF (ABS(DPHI).LE.10E-10) THEN
      IAPEX=1
C      WRITE(*,*)'IAPEX ds BEHAV1 =',IAPEX
C      WRITE(*,*)'Dans l element',IBB
C      WRITE(*,*)'et au point d intégration',IGAU
      DO 50 I=1,NSTRS
      DO 50 J=1,NSTRS
      AI(I,J)=0.D0
   50 CONTINUE
        AI(1,1)=0.5
        AI(1,2)=0.5
        AI(2,1)=AI(1,2)
        AI(2,2)=AI(1,1)
        IF (IMOD.EQ.1.OR.IMOD.EQ.3) AI(3,3)=0.D0
        IF (IMOD.EQ.2.OR.IMOD.EQ.4) AI(3,3)=1.D0
        GOTO 75
      ENDIF
C
C ************** Mise a jour des contraintes ***************
C
C ---------------- calcul de la matrice A ------------------
C
      DO 60 I=1,NSTRS
      DO 60 J=1,NSTRS
      A(I,J)=0.D0
   60 CONTINUE
C
      DG=YOUN/2.D0/(1.D0+XNU)
C
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
      A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI
      A(2,2)=A(1,1)
      A(1,2)=-(A(1,1)-1.D0)
      A(2,1)=A(1,2)
      A(3,3)=1.D0+(DLAM1*DG)/DPHI
      ELSE
      A(1,1)=1.D0+(DLAM1*DG)/2.D0/DPHI
      A(2,2)=A(1,1)
      A(1,2)=-(A(1,1)-1.D0)
      A(2,1)=A(1,2)
      A(3,3)=1.D0
      A(4,4)=1.D0+(DLAM1*DG)/DPHI
      ENDIF
C
C -------------- invertion de la matrice A -----------------
C
      DO 70 I=1,NSTRS
      DO 70 J=1,NSTRS
      AI(I,J)=A(I,J)
   70 CONTINUE
      CALL INVMA2(AI,NSTRS,ISING)
      IF (ISING.EQ.1) THEN
       WRITE(*,*)'MATRICE AI singuliere ds BEHAV1'
      ENDIF
C
C -------------- mise a jour des contraintes  ------------
C
   75 CONTINUE
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
      VEC1(1)=AD1
      VEC1(2)=AD1
      VEC1(3)=0.D0
      VEC1(4)=0.D0
      ENDIF
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
      VEC1(1)=AD1
      VEC1(2)=AD1
      VEC1(3)=AD1*2.D0*XNU
      VEC1(4)=0.D0
      ENDIF
C
      DO 80 I=1,NSTRS
      DEPSI(I)=SIGE(I)-0.5*DLAM1*VEC1(I)
   80 CONTINUE
C
      DO 90 I=1,NSTRS
      SIGF(I)=0.0D+00
      DO 90 J=1,NSTRS
      SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
   90 CONTINUE
C
C ******** Verification du critere ****************
C
      CALL PRINC(SIGF,V1,NSTRS)
      FCRI1=V1(1)-SEQ1
C
      IF (IBROY.EQ.0.AND.ABS(FCRI1).GE.CRIMAX) THEN
C      WRITE(*,*)'****************************************'
C      WRITE(*,*)'LE RESIDU DIVERGE AVEC BROYDEN'
C      WRITE(*,*)'on passe donc a la secante'
C      WRITE(*,*)'Dans l element',IBB
C      WRITE(*,*)'et au point d intégration',IGAU
C      WRITE(*,*)'CRIMAX=',CRIMAX
C      WRITE(*,*)'****************************************'
      ITER=ITR
      ENDIF
C
C *******  Compteur sur la methode de resolution ****
C
      IF (IBROY.EQ.0.AND.ITER.EQ.ITR) THEN
        IBROY=1
        ITANG=1
        ITER=1
        IAPEX=0
        GOTO 8
      ENDIF
C
C ******* non convergence **************************
C
      IF (ABS(FCRI1).GT.PRB.AND.ITER.LT.ITR) THEN
       IF (IBROY.EQ.0) THEN
       DJAC1=(FCRI0-FCRI1)/(DLAM0-DLAM1)
       DLAM0=DLAM1
       DJAC0=DJAC1
       FCRI0=FCRI1
       ITER=ITER+1
       IF (ITER.GE.(ITR-1)) THEN
C        WRITE(*,*)'***********************'
C        WRITE(*,*)'BROYDEN n a pas aboutit'
C        WRITE(*,*)'Dans l element',IBB
C        WRITE(*,*)'et au point d intégration',IGAU
C        WRITE(*,*)'ITER=',ITER
C        WRITE(*,*)'FCRI=',FCRI1
C        WRITE(*,*)'***********************'
       ENDIF
       GOTO 40
      ENDIF
      IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
       DLAM0=DLAM1
       FCRI0=FCRI1
       CALL PRINC(SIGF,V1,NSTRS)
       ITER=ITER+1
       IF (ITER.GE.(ITR-5)) THEN
        WRITE(*,*)'ITER=',ITER
        WRITE(*,*)'FCRI=',FCRI1
       ENDIF
       GOTO 9
      ENDIF
      ENDIF
      IF (ITER.GE.ITR.AND.ABS(FCRI1).GT.PRB2) THEN
        WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV1'
        WRITE(*,*)'FCRI=',FCRI1
        WRITE(*,*)'Dans l element',IBB
        WRITE(*,*)'et au point d intégration',IGAU
        WRITE(*,*)'IPLA=',IPLA
        WRITE(*,*)'IFIS=',IFIS
C        STOP
      ENDIF
C
C ************* Calcul de la DEP ****************
C
      IF (IAPEX.EQ.0) THEN
      CALL DERI2(D2FSIG,DPHI,P1,SIGF,NSTRS,BETJEF)
      CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
     & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
      ELSE
C      WRITE(*,*)'ELAS',IAPEX
      DO 95 I=1,NSTRS
      DO 95 J=1,NSTRS
         AH(I,J)=0.D0
      DO 95 K=1,NSTRS
         AH(I,J)=AH(I,J)+AI(I,K)*D(K,J)
   95 CONTINUE
       DFSIG(1)=0.5*SQRT(2.D0)
       DFSIG(2)=0.5*SQRT(2.D0)
       DFSIG(3)=0.D0
       DFSIG(4)=0.D0
       CALL CREDEP(AH,DFSIG,PAEC,NSTRS,DEP,BETJEF)
      ENDIF
      IAPEX=0
      CRIMAX=0.D0
C
C ***********************************************
C
      IF (DK.LT.0.D0) THEN
       WRITE(*,*)'Dans BEHAV1, DLAMDA est negatif:',DLAM1
       WRITE(*,*)'A l iteration :',ITER
C       STOP
      ENDIF
      DEPS1=DK
      SIG1=V1(1)
      IF (V1(1).LT.0.D0) V1(1)=0.D0
      CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF)
      IF(IVIS.LE.2) THEN
       CALL INDICA(DK,0.0D0,IFIS,IPLA,1,BETJEF)
      ELSE
       CALL INDIC1(DK,0.0D0,IFIS,IPLA,1,BETJEF,NECH0,NECH1)
      ENDIF
C ********************** CALCUL VISCOPLASTIQUE ***********************
      IF (IVIS.EQ.1) THEN
        CALL VISPLA(SIGR,SIGF,DSIGT,NSTRS,DEPS1,DEPS2,
     &                   SIGP,SIGRV,DSTRN,D,BETJEF,VISCO)
      ENDIF
C ********************************************************************
  100 CONTINUE
C
      RETURN
      END