C BEHAV2    SOURCE    CB215821  26/08/24    21:15:14     12622          
      SUBROUTINE BEHAV2(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 PRE ET POST PIC
C      Un seul critere de compression: Von Mises ou Drucker-Prager
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),DSTRN(4)
      DIMENSION SIGR(4),AC(4),SIGE(4),SIGRV(4),SIGP(4)
      DIMENSION D(4,4),DEP(4,4),D2FSIG(4,4),P2(4,4),DP(4,4)
      DIMENSION A(4,4),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
      IAPEX=0
      PRB=1.D-10
      PRB2=1.D-6
      ITER=1
      IBROY=0
      ITANG=0
      Rb=COLI
      Ft=PALF*COLI
      IEC2=0
      CRIMAX=0.D0
      SEQ = 0.D0
      SEQ2 = 0.D0
      CALL ZERO(SIGE,4,1)
***** CALL ZERO(V1,4,1)
      CALL ZERO(P2,4,4)
      CALL ZERO(D,4,4)
C
      DO 10 I=1,NSTRS
      SIGE(I)=SIGF(I)
   10 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 comprime pour 1ere fois ***********
C
      IF (IPLA.EQ.0) THEN
       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
        CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
       ENDIF
       IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
        CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
       ENDIF
        DEPS2=0.D0
        SIG2=Rb*AA
        IF (SEQ.GT.SIG2) THEN
         IPLA=1
        ENDIF
      ENDIF
      IF (IFIS.EQ.0) THEN
        DEPS1=0.D0
        SIG1=Ft
      ENDIF
C
    8 CONTINUE
C
       IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
        CALL CRIVON(SIGE,SEQ,NSTRS,BETJEF)
       ENDIF
       IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
        CALL CRIDRU(SIGE,SEQ,NSTRS,BETJEF)
       ENDIF
      FCRIC0=SEQ-SIG2
      CRIMAX=ABS(100.D0*FCRIC0)
      IF (FCRIC0.LE.0.D0) THEN
        IEC2=1
      ENDIF
      DK=DEPS2
      DLAM0=0.D0
C
   12 CONTINUE
C
      DO 15 I=1,NSTRS
      SIGF(I)=SIGE(I)
   15 CONTINUE
C
      IF (IEC2.EQ.1) THEN
      DO 101 I=1,NSTRS
      DO 16 J=1,NSTRS
      DEP(I,J)=D(I,J)
   16 CONTINUE
  101 CONTINUE
       GOTO 100
      ENDIF
C
    9 CONTINUE
C
      IF(IVIS.LE.2) THEN
       CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
      ELSE
       CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
      ENDIF
C
C ************ Preparation des criteres ********************
C
C ------------------ Von Mises -----------------------------
      IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
        A1=0.D0
        A2=1.D0
      ENDIF
C ----------------- Drucker Prager -------------------------
      IF (IMOD.EQ.3.OR.IMOD.EQ.4) THEN
        A1=(BETA-1.D0)/(2*BETA-1.D0)
        A2=BETA/(2*BETA-1.D0)
      ENDIF
C ----------------------------------------------------------
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
      SX=SIGF(1)
      SY=SIGF(2)
      SXY=SIGF(3)
        R1=(SX-SY)*(SX-SY)+SX*SX+SY*SY+6*SXY*SXY
        R=SQRT(0.5*R1)
      IF (ABS(R).LE.1.D-10) THEN
      WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
      STOP
      ENDIF
       DFSIG(1)=(2.D0*SX-SY)/2.D0/R/A2+A1/A2
       DFSIG(2)=(2.D0*SY-SX)/2.D0/R/A2+A1/A2
       DFSIG(3)=3.D0*SXY/R/A2
      P2(1,1)=2.D0/A2/A2
      P2(1,2)=-1.D0/A2/A2
      P2(1,3)=0.D0/A2/A2
      P2(2,1)=-1.D0/A2/A2
      P2(2,2)=2.D0/A2/A2
      P2(2,3)=0.D0/A2/A2
      P2(3,1)=0.D0/A2/A2
      P2(3,2)=0.D0/A2/A2
      P2(3,3)=6.D0/A2/A2
      ENDIF
C
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
      SX=SIGF(1)
      SY=SIGF(2)
      SZ=SIGF(3)
      SXY=SIGF(4)
      TX=(SX-SY)*(SX-SY)
      TY=(SX-SZ)*(SX-SZ)
      TZ=(SY-SZ)*(SY-SZ)
         R1=TX+TY+TZ+6*SXY*SXY
         R=SQRT(0.5*R1)
      IF (ABS(R).LE.1.D-10) THEN
      WRITE(*,*)'INDETERMINATION de DFSIG ds BEHAV2'
      ENDIF
        DFSIG(1)=(2.D0*SX-SY-SZ)/2.D0/R/A2+A1/A2
        DFSIG(2)=(2.D0*SY-SX-SZ)/2.D0/R/A2+A1/A2
        DFSIG(3)=(2.D0*SZ-SX-SY)/2.D0/R/A2+A1/A2
        DFSIG(4)=3.D0*SXY/R/A2
      P2(1,1)=2.D0/A2/A2
      P2(1,2)=-1.D0/A2/A2
      P2(1,3)=-1.D0/A2/A2
      P2(1,4)=0.D0/A2/A2
      P2(2,1)=-1.D0/A2/A2
      P2(2,2)=2.D0/A2/A2
      P2(2,3)=-1.D0/A2/A2
      P2(2,4)=0.D0/A2/A2
      P2(3,1)=-1.D0/A2/A2
      P2(3,2)=-1.D0/A2/A2
      P2(3,3)=2.D0/A2/A2
      P2(3,4)=0.D0/A2/A2
      P2(4,1)=0.D0/A2/A2
      P2(4,2)=0.D0/A2/A2
      P2(4,3)=0.D0/A2/A2
      P2(4,4)=6.D0/A2/A2
      ENDIF
C
      DO 102 I=1,NSTRS
      AC(I)=0.D0
      DO 25 J=1,NSTRS
      AC(I)=AC(I)+D(I,J)*DFSIG(J)
   25 CONTINUE
  102 CONTINUE
C
C ************ Preparation du jacobien  ***********************
C
      FDF=0.D0
      DO 30 J=1,NSTRS
      FDF=FDF+DFSIG(J)*AC(J)
   30 CONTINUE
C
      DJAC0=-(FDF+PAEC)
C
C ************ Debut des iterations internes *******************
C
   40 CONTINUE
C
C *************** Determination de DK et DLAM ******************
C
C
      DLAM1=DLAM0-FCRIC0/DJAC0
C
      IF (DLAM1.LT.0.D0) IEC2=1
      IF (IEC2.EQ.1) THEN
C       WRITE(*,*)'Dans BEHAV2, DLAMDA est negatif:',DLAM1
C       WRITE(*,*)'A l iteration :',ITER
        GOTO 12
      ENDIF
C
      DK=DEPS2+DLAM1
C
      IF(IVIS.LE.2) THEN
       CALL UNICOU(DK,PAEC,2,SEQ2,BETJEF)
      ELSE
       CALL UNICO1(DK,PAEC,2,SEQ2,BETJEF,NECH0,NECH1)
      ENDIF
C
C ************** Determination de DPHI2 *********************
C
      IF (IMOD.EQ.1.OR.IMOD.EQ.2) THEN
       DPHI2=SEQ2
      ENDIF
      IF (IMOD.EQ.3) THEN
       SOM=SIGF(1)+SIGF(2)+SIGF(3)
       DPHI2=SEQ2-A1*SOM/A2
      ENDIF
      IF (IMOD.EQ.4) THEN
      SOM=SIGE(1)+SIGE(2)+SIGE(3)
      DPHI2=A1*SOM/A2-A1*A1*DLAM1*3.D0*ADD*(1.D0+XNU)/A2/A2
      DPHI2=SEQ2-DPHI2
      ENDIF
C      IF (DPHI2.LT.0.D0) THEN
C       WRITE(*,*)'ATTENTION DPHI2 ds BEHAV2 NEGATIF'
C       WRITE(*,*)'DPHI2=',DPHI2
C      ENDIF
C
C---------- Cas du critere reduit a un point ----------------
C
      IF (ABS(DPHI2).LE.10E-10) THEN
      IAPEX=1
C      WRITE(*,*)'IAPEX ds BEHAV2 =',IAPEX
      DO 103 I=1,NSTRS
      DO 50 J=1,NSTRS
      AI(I,J)=0.D0
   50 CONTINUE
  103 CONTINUE
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
      DO 104 I=1,3
      DO 53 J=1,3
      AI(I,J)=1.D0/3.D0
   53 CONTINUE
  104 CONTINUE
      ENDIF
        GOTO 75
      ENDIF
C
C ************** Mise a jour des contraintes ***************
C
C ---------------- calcul de la matrice A ------------------
C
      DO 105 I=1,NSTRS
      DO 60 J=1,NSTRS
      A(I,J)=0.D0
   60 CONTINUE
  105 CONTINUE
C
      DG=YOUN/2.D0/(1.D0+XNU)
C
      IF (IMOD.EQ.1.OR.IMOD.EQ.3) THEN
      DW1=AD*(2.D0-XNU)
      DW2=AD*(2.D0*XNU-1.D0)
      A(1,1)=1.D0+(DLAM1*DW1)/2.D0/DPHI2/A2/A2
      A(2,2)=A(1,1)
      A(1,2)=(DLAM1*DW2)/2.D0/DPHI2/A2/A2
      A(2,1)=A(1,2)
      A(3,3)=1.D0+(DLAM1*3.D0*DG)/DPHI2/A2/A2
      ENDIF
      IF (IMOD.EQ.2.OR.IMOD.EQ.4) THEN
      A(1,1)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
      A(2,2)=A(1,1)
      A(3,3)=1.D0+(DLAM1*2.D0*DG)/DPHI2/A2/A2
      A(1,2)=-(DLAM1*DG)/DPHI2/A2/A2
      A(2,1)=A(1,2)
      A(1,3)=-(DLAM1*DG)/DPHI2/A2/A2
      A(3,1)=A(1,3)
      A(3,2)=A(1,3)
      A(2,3)=A(1,3)
      A(4,4)=1.D0+(DLAM1*DG*3.D0)/DPHI2/A2/A2
      ENDIF
C
C -------------- invertion de la matrice A -----------------
C
      DO 106 I=1,NSTRS
      DO 70 J=1,NSTRS
      AI(I,J)=A(I,J)
   70 CONTINUE
  106 CONTINUE
      CALL INVMA2(AI,NSTRS,ISING)
      IF (ISING.EQ.1) THEN
       WRITE(*,*)'MATRICE AI singuliere ds BEHAV2'
       WRITE(*,*)'NSTRS=',NSTRS
      ENDIF
C
C -------------- mise a jour des contraintes  --------------
C
   75 CONTINUE
      IF(IMOD.EQ.3) THEN
       VEC1(1)=AD*(1.D0+XNU)
       VEC1(2)=AD*(1.D0+XNU)
       VEC1(3)=0.D0
       VEC1(4)=0.D0
      ENDIF
      IF(IMOD.EQ.4) THEN
       VEC1(1)=ADD*(1.D0+XNU)
       VEC1(2)=ADD*(1.D0+XNU)
       VEC1(3)=ADD*(1.D0+XNU)*2.D0*XNU
       VEC1(4)=0.D0
      ENDIF
      IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
       VEC1(1)=0.D0
       VEC1(2)=0.D0
       VEC1(3)=0.D0
       VEC1(4)=0.D0
      ENDIF
C
      DO 80 I=1,NSTRS
      DEPSI(I)=SIGE(I)-A1*DLAM1*VEC1(I)/A2
   80 CONTINUE
C
      DO 107 I=1,NSTRS
      SIGF(I)=0.0D+00
      DO 90 J=1,NSTRS
      SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
   90 CONTINUE
  107 CONTINUE
C
C ************ Verification du critere ********************
C
      IF(IMOD.EQ.1.OR.IMOD.EQ.2) THEN
       CALL CRIVON(SIGF,SEQ,NSTRS,BETJEF)
      ENDIF
      IF(IMOD.EQ.3.OR.IMOD.EQ.4) THEN
       CALL CRIDRU(SIGF,SEQ,NSTRS,BETJEF)
      ENDIF
      FCRIC1=SEQ-SEQ2
C
      IF (IBROY.EQ.0.AND.ABS(FCRIC1).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
        IAPEX=0
        IEC2=0
        ITER=1
        GOTO 8
      ENDIF
C
C ******* non convergence **************************
C
      IF (ABS(FCRIC1).GT.PRB.AND.ITER.LT.ITR) THEN
       IF (IBROY.EQ.0) THEN
       DJAC0=(FCRIC0-FCRIC1)/(DLAM0-DLAM1)
       DLAM0=DLAM1
       FCRIC0=FCRIC1
       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=',FCRIC1
C        WRITE(*,*)'***********************'
       ENDIF
       GOTO 40
      ENDIF
      IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
       DLAM0=DLAM1
       FCRIC0=FCRIC1
       ITER=ITER+1
       IF (ITER.GE.(ITR-5)) THEN
C        WRITE(*,*)'ITER=',ITER
C        WRITE(*,*)'FCRI=',FCRIC1
       ENDIF
       GOTO 9
      ENDIF
      ENDIF
      IF (ITER.GE.ITR.AND.ABS(FCRIC1).GT.PRB2) THEN
        WRITE(*,*)'NON CONVERGENCE INTERNE dans BEHAV2'
        WRITE(*,*)'FCRI=',FCRIC1
        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,DPHI2,P2,SIGF,NSTRS,BETJEF)
      CALL LADEP(SIGF,DEP,PAEC,DPHI,NSTRS,IFOU,
     & DFSIG,D2FSIG,DLAM1,D,DP,BETJEF)
      ELSE IF (IAPEX.EQ.1) THEN
C      WRITE(*,*)'ELAS2',IAPEX
      DO 109 I=1,NSTRS
      DO 108 J=1,NSTRS
         DEP(I,J)=0.D0
      DO 96 K=1,NSTRS
         DEP(I,J)=DEP(I,J)+AI(I,K)*D(K,J)
   96 CONTINUE
  108 CONTINUE
  109 CONTINUE
      ENDIF
      IAPEX=0
      CRIMAX=0.D0
*      DO 97 I=1,NSTRS
*      DO 97 J=1,NSTRS
*         DEP(I,J)=D(I,J)
*   97 CONTINUE
C
C ***********************************************
C
      DEPS2=DK
      SIG2=SEQ2
      IF(IVIS.LE.2) THEN
       CALL INDICA(0.D0,DK,IFIS,IPLA,2,BETJEF)
      ELSE
       CALL INDIC1(0.D0,DK,IFIS,IPLA,2,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
      IEC2=0
C
      RETURN
      END