C RUPTHE    SOURCE    PV        17/12/08    21:17:45     9660           
C--------------------------------------------------------------------
        subroutine RUPTHE(WRK52,WRK53,WRK54,WRK27,IB,IGAU,NBPGAU)

C     ***************************************************
C       Elément macro de rupteur (derniÚre version)
C          crée par Huyen Nguyen - décembre 2010
C          *************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

-INC PPARAM
-INC CCOPTIO
-INC DECHE
C-----------------------------------------------------------------------
C     GLOBAL MODEL
        DIMENSION XMKBel(6,6)
        DIMENSION XMKAel(6,6)
*       real*8 XKBel
*       real*8 XKAel

C     Endommagement du béton
*       real*8 Zplus, Zmoin
*       real*8 Dplus, Dmoin, Dh
*       real*8 Ym, Ymp, Fonc1

C     Frottement du béton
         logical conv
*        integer iter

*       real*8 up, Ecrb, Fp
*       real*8 Fonc2, Xmul1, Fonc3, Xlamda1
*       real*8 dfdfb, dpidfb, dfdxb, dpidxb

C     Frottement du béton

*       real*8 ua, Ecra, Fa
*       real*8 Fonc4, Xmul2, Fonc5, Xlamda2
*       real*8 dfdfa, dpidfa, dfdxa, dpidxa

C
        DIMENSION uo(6)
        DIMENSION duo(6)
        DIMENSION force(6)
        DIMENSION da(10)


C----------------------------------------------------------------------
C     Affectation des données du modÚle


        da(1) =XMAT(13)
        da(2) =XMAT(14)
        da(3) =XMAT(15)
        da(4) =XMAT(16)
        da(5) =XMAT(17)
        da(6) =XMAT(18)
        da(7) =XMAT(19)
        da(8) =XMAT(20)
        da(9) =XMAT(21)
        da(10) =XMAT(22)

        XKBel=da(1)
        Ybo=da(2)
        C1=da(3)
        D1=da(4)
        alphab=da(5)
        betab=da(6)
        XKAel=da(7)
        Yao=da(8)
        alphaa=da(9)
        betaa=da(10)
C---------------------------------------------------------------------
C     Affectation des données historiques courantes

        force(1)=SIG0(1)
        force(2)=SIG0(2)
        force(3)=SIG0(3)
        force(4)=SIG0(4)
        force(5)=SIG0(5)
        force(6)=SIG0(6)

        Dplus=VAR0(2)
        Dmoin=VAR0(3)
        Dh=VAR0(4)
        Zplus=VAR0(5)
        Zmoin=VAR0(6)
        up=VAR0(7)
        Ecrb=VAR0(8)
        Fp=VAR0(9)
        ua=VAR0(10)
        Ecra=VAR0(11)
        Fa=VAR0(12)
C---------------------------------------------------------------------

        uo(1)=epst0(1)
        uo(2)=epst0(2)
        uo(3)=epst0(3)
        uo(4)=epst0(4)
        uo(5)=epst0(5)
        uo(6)=epst0(6)


        duo(1)=DEPST(1)
        duo(2)=DEPST(2)
        duo(3)=DEPST(3)
        duo(4)=DEPST(4)
        duo(5)=DEPST(5)
        duo(6)=DEPST(6)

C Correction données d'entrées Castem

        DO I=1,6
          uo(I)=uo(I)+duo(I)
        ENDDO

C       PRINT*, uo
C       PRINT*, duo



C     *****************************************************************


        do i=1,6
           do j=1,6
              XMKBel(i,j)=0.
           enddo
        enddo
        XMKBel(2,2)=XKBel
C       XMKBel(4,4)=XKBel

        do i=1,6
           do j=1,6
              XMKAel(i,j)=0.
           enddo
        enddo
        do i=1,6
              XMKAel(i,i)=XMAT(i+6)
        enddo
        XMKAel(2,2)=XKAel
C       XMKAel(4,4)=XKAel

C     ************************Endomagement du béton

        Ym=5.d-1*XMKBel(2,2)*(uo(2)**2.d0)
        Ymp=Ym

            IF (uo(2) .ge. 0.d0) THEN

        Fonc1=(Ymp-Zplus-Ybo)

                IF (Fonc1 .gt. 0.d0) THEN

        Dplus=1.d0-(1.d0/(1.d0+C1*(Ym**D1)))
        Dmoin=Dmoin
        Dh=max(Dplus,Dmoin)
        Zplus=Zplus+(5.d-1*XMKBel(2,2)*duo(2)*(2.d0*uo(2)-duo(2)))
        Zmoin=Zmoin

                ELSE
        Dh=Dh
        Dplus=Dplus
        Dmoin=Dmoin
        Zplus=Zplus
        Zmoin=Zmoin

                ENDIF

                ELSE

         Fonc1=(Ymp-Zmoin-Ybo)

C          print*, 'Déplacement négatif'

                IF (Fonc1 .gt. 0.d0) THEN

        Dplus=Dplus
        Dmoin=1.d0-(1.d0/(1.d0+C1*(Ym**D1)))
        Dh=max(Dplus,Dmoin)
        Zmoin=Zmoin+(5.d-1*XMKBel(2,2)*duo(2)*(2.d0*uo(2)-duo(2)))
        Zplus=Zplus

                ELSE

        Dh=Dh
        Dplus=Dplus
        Dmoin=Dmoin
        Zplus=Zplus
        Zmoin=Zmoin

                ENDIF
                ENDIF

C***************Frottement du beton

        Fp=Dh*XMKBel(2,2)*(uo(2)-up)

        Fonc2=ABS(Fp-Ecrb)

        IF (Fonc2 .gt. 0.d0) THEN
        conv=.false.
        iter=0
                DO WHILE (.not. conv .and. iter .le. 500)

        iter=iter+1

                IF ( (Fp-Ecrb) .GE. 0.d0) THEN

        Xmul1=1.d0
                ELSE
        Xmul1=-1.d0
                ENDIF

        Fonc3=ABS(Fp-Ecrb)

        dfdfb=Xmul1
        dpidfb=Xmul1
        dfdxb=-1.0*Xmul1
        dpidxb=-1.0*Xmul1+alphab*Ecrb

        Xlamda1=Fonc3/(XMKBel(2,2)*Dh*dfdfb*dpidfb+betab*dfdxb*dpidxb)

        up=up+Xlamda1*dpidfb

        Ecrb=Ecrb-betab*Xlamda1*dpidxb
        Fp=Fp-XMKBel(2,2)*Dh*Xlamda1*dpidfb

        Fonc3=ABS(Fp-Ecrb)

        conv=((ABS(Fonc3/Fonc2) .LE. 1.d-7). or.
     & (Xlamda1 .LE. 1.d-7))

                ENDDO

                IF (.not. conv) THEN
                 KERRE = 22
                 RETURN
                ENDIF
                ENDIF


C**************Acier plastique

        Fa=XMKAel(2,2)*(uo(2)-ua)
        Fonc4=ABS(Fa-Ecra)-Yao

        IF (Fonc4 .gt. 0.d0) THEN
        conv=.false.
        iter=0
                DO WHILE (.not. conv .and. iter .le. 500)
        iter=iter+1

                IF ( (Fa-Ecra) .GE. 0.d0) THEN

        Xmul2=1.d0
                ELSE
        Xmul2=-1.d0
                ENDIF

        Fonc5=ABS(Fa-Ecra)-Yao

        dfdfa=Xmul2
        dpidfa=Xmul2
        dfdxa=-1.0*Xmul2
        dpidxa=-1.0*Xmul2+alphaa*Ecra

        Xlamda2=Fonc5/(XMKAel(2,2)*dfdfa*dpidfa+betaa*dfdxa*dpidxa)

        ua=ua+Xlamda2*dpidfa
        Ecra=Ecra-betaa*Xlamda2*dpidxa
        Fa=Fa-XMKAel(2,2)*Xlamda2*dpidfa

        Fonc5=ABS(Fa-Ecra)-Yao

        conv=((ABS(Fonc5/Fonc4) .LE. 1.d-7). or.
     & (Xlamda2 .LE. 1.D-7))

        ENDDO

        IF (.not. conv) THEN
                 KERRE = 22
                 RETURN
        ENDIF
        ENDIF

        Fa=XMKAel(2,2)*(uo(2)-ua)


C*****************************************************************
                IF (uo(2) .ge. 0.d0) THEN
        force(2)=XMKBel(2,2)*uo(2)*(1.d0-Dplus)
     &  +Dh*XMKBel(2,2)*(uo(2)-up)+Fa
                ELSE
        force(2)=XMKBel(2,2)*uo(2)*(1.d0-Dmoin)
     &  +Dh*XMKBel(2,2)*(uo(2)-up)+Fa
                ENDIF

        force(1)=XMAT(7)*uo(1)
        force(3)=XMAT(9)*uo(3)
        force(4)=XMAT(10)*uo(4)
        force(5)=XMAT(11)*uo(5)
        force(6)=XMAT(12)*uo(6)

C*****************************************************************
        EPSE=sqrt(EPSE)
        VARF(1)=EPSE
        VARF(2)=Dplus
        VARF(3)=Dmoin
        VARF(4)=Dh
        VARF(5)=Zplus
        VARF(6)=Zmoin
        VARF(7)=up
        VARF(8)=Ecrb
        VARF(9)=Fp
        VARF(10)=ua
        VARF(11)=Ecra
        VARF(12)=Fa

        SIGF(1)=force(1)
        SIGF(2)=force(2)
        SIGF(3)=force(3)
        SIGF(4)=force(4)
        SIGF(5)=force(5)
        SIGF(6)=force(6)

      return

C#####################################################################
C#####################################################################

      end