C ACTI1     SOURCE    PV        22/04/19    16:17:58     11344          
C ACTI1    SOURCE
      SUBROUTINE ACTI1(SIG0,SIGF,D,NSTRSS,BETINSA)
C
C ===============================================================
C           Un seul critere de traction: DRUCKER PRAGER 3D
C ================================================================
C CE SOUS-PROGRAMME EST APPELE DANS "BONE3D".
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION SIG0(NSTRSS),SIGF(NSTRSS)
      DIMENSION DFSIG(6),VEC1(6),DEPSI(6),SIGE(6)
      DIMENSION D(6,6),AC(6),A(6,6),AI(6,6),AIM(4,4)
C
      SEGMENT BETINSA
       REAL*8 RT,RC,YOUN,XNU,GFT,GFC,CAR
       REAL*8 DKT,DKC,SEQT,SEQC,ENDT,ENDC
       INTEGER IFIS,IPLA,IBB,IGAU
      ENDSEGMENT
C
      CALL ZERO(SIGE,6,1)
      DO 2 I=1,NSTRSS
       SIGE(I)=SIGF(I)
    2 CONTINUE
      IAPEX=0
      PRB=1.D-5
      PRB2=1.D-2
      ITER=1
      ITR = 1500
      ITANG=0
      IBROY=0
      CRIMAX=0.D0
      SEQ = 0.D0
      SEQ1 = 0.D0
C
    8 CONTINUE
C
      DO 4 I=1,NSTRSS
       SIGF(I)=SIGE(I)
    4 CONTINUE
      DK=DKT
      DLAM0=0.D0
      CALL DRUTRA(SIGF,SEQTE,BETINSA)
      FCRI0 = SEQTE - SEQT
      CRIMAX=ABS(100.D0*FCRI0)
    9 CONTINUE
C
C ************ CALCUL DU JACOBIEN INITIAL  ********************
C
      CALL DRUTR1(SIGF,DFSIG,BETINSA)
C
      DO    I=1,NSTRSS
      AC(I)=0.D0
      DO    J=1,NSTRSS
      AC(I)=AC(I)+D(I,J)*DFSIG(J)
      enddo
      enddo
      FDF=0.D0
      DO 30 J=1,NSTRSS
      FDF=FDF+DFSIG(J)*AC(J)
   30 CONTINUE
C
      CALL ENDAME(1,BETINSA)
      CALL FORECR(DK,PAEC,1,SEQ,BETINSA)
C
      DJAC0=-(PAEC+FDF)
C
C ************ DEBUT ITERATION INTERNES ************************
C
   40 CONTINUE
C
C *************** Determination de DK  *************************
C
      DLAM1=-FCRI0/DJAC0+DLAM0
      DK=DKT+DLAM1
C
C *********** Estimation contrainte quivalente *****************
C
      CALL ENDAME(1,BETINSA)
      CALL FORECR(DK,PAEC,1,SEQ,BETINSA)
C
C ************** Determination de DPHI *************************
C
      CALL DRUTR2(SIGE,SEQ,DPHI,DLAM1,VEC1,BETINSA)
C
C **************** Cas de l'apex *******************************
C
      IF (ABS(DPHI).LE.10E-10) THEN
      IAPEX=1
      DO    I=1,NSTRSS
      DO    J=1,NSTRSS
      AI(I,J)=0.D0
      enddo
      enddo
        AI(1,1)=1./3.
        AI(1,2)=AI(1,1)
        AI(1,2)=AI(1,1)
        AI(1,3)=AI(1,1)
        AI(2,1)=AI(1,1)
        AI(2,2)=AI(1,1)
        AI(2,3)=AI(1,1)
        AI(3,1)=AI(1,1)
        AI(3,2)=AI(1,1)
        AI(3,3)=AI(1,1)
        GOTO 75
      ENDIF
C
C ************** Mise a jour des contraintes ***************
C
C ---------------- calcul de la matrice A ------------------
C
      DO    I=1,NSTRSS
      DO    J=1,NSTRSS
      A(I,J)=0.D0
      enddo
      enddo
C
      DG=YOUN/(1.D0+XNU)
C
      A(1,1)=1.D0+2.*(DLAM1*DG)/2.D0/DPHI
      A(2,2)=A(1,1)
      A(3,3)=A(1,1)
      A(1,2)=-(DLAM1*DG)/2.D0/DPHI
      A(1,3)=A(1,2)
      A(2,1)=A(1,2)
      A(2,3)=A(1,2)
      A(3,1)=A(1,2)
      A(3,2)=A(1,2)
      A(4,4)=1.D0+3.*(DLAM1*DG)/2.D0/DPHI
      A(5,5)=A(4,4)
      A(6,6)=A(4,4)
C
C -------------- invertion de la matrice A -----------------
C
      CALL ZERO(AI,6,6)
      CALL ZERO(AIM,4,4)
C
      DO    I=1,3
      DO    J=1,3
      AIM(I,J)=A(I,J)
      enddo
      enddo
      CALL INVMA2(AIM,3,ISING)
      IF (ISING.EQ.1) THEN
       WRITE(*,*)'MATRICE AIM singuliere ds ACTI1'
      ENDIF
      DO    I=1,3
      DO    J=1,3
      AI(I,J)=AIM(I,J)
      enddo
      enddo
      AI(4,4) = 1./A(4,4)
      AI(5,5) = 1./A(5,5)
      AI(6,6) = 1./A(6,6)
C
C -------------- mise a jour des contraintes  ------------
C
   75 CONTINUE
C
      DO 80 I=1,NSTRSS
      DEPSI(I)=SIGE(I)-DLAM1*VEC1(I)
   80 CONTINUE
C
      DO    I=1,NSTRSS
      SIGF(I)=0.
      DO    J=1,NSTRSS
      SIGF(I)=SIGF(I)+AI(I,J)*DEPSI(J)
      enddo
      enddo
C
C ******** VERIFICATION DU CRITERE ****************
C
      CALL DRUTRA(SIGF,SEQI,BETINSA)
      FCRI1 = SEQI - SEQ
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
       GOTO 40
       ENDIF
       IF (IBROY.EQ.1.AND.ITANG.EQ.1) THEN
       DLAM0=DLAM1
       FCRI0=FCRI1
       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 ACTI1'
        WRITE(*,*)'FCRI=',FCRI1
        WRITE(*,*)'Dans l element',IBB
        WRITE(*,*)'et au point d intégration',IGAU
        STOP
      ENDIF
C
C *******  convergence **************************
C
      SEQT = SEQ
      DKT = DK
C
C ***********************************************
C
  100 CONTINUE
C
      RETURN
      END