C ISSGRA    SOURCE    PV        17/12/08    21:17:29     9660           
C
      SUBROUTINE ISSGRA(WRK52,WRK53,WRK54,WRK27,IB,IGAU,NBPGAU)
C
C______________________________________________________________________
C
C    MODELE  ISSGRA
C______________________________________________________________________
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
c inutile? 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC DECHE
*
C     GLOBAL MODEL
        DIMENSION XMKel(5,5)
C        DIMENSION XMCel(5,5)
        DIMENSION uo(5), duo(5), Fo(5)
        DIMENSION XMK(5,5), Ftr(5), force(5)
*
        DIMENSION upl(5), uel(5)
        DIMENSION uoi(5), duoi(5), dueli(5), force_preci(5)
C        DIMENSION ueli(5)

C     Modèle de Plasticité
        DIMENSION ddlambda(2), yield(2)
        DIMENSION Ho(2,2), H(2,2), Hinv(2,2)
        DIMENSION duP(5)
*
C     Modèle de plasticité 1
        DIMENSION q1(5), upl1(5)
C        DIMENSION q1_prec(5)
        DIMENSION Gradgs(5),Gradfs(5),Gradfq(5)
        DIMENSION hh1(5), XKGradgs(5)
*
C     Modèle de plasticité vertical
        DIMENSION upl2(5), Gradf2s(5), XKGradf2s(5)
*
C     Modèle de Décollement
        DIMENSION dforce(5), forcepred(5), force_prec(5), duel(5)
*
C     Amortissment
        DIMENSION upl_preci(5),q1_preci(5)
*
        DIMENSION da(20)
C        DIMENSION TIME(2)
*
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C     Affectation des données du modèle
        da(1) =XMAT(13)
        da(2) =XMAT(14)
        da(3) =XMAT(15)
        da(4) =XMAT(7)
        da(5) =XMAT(8)
        da(6) =XMAT(12)
        da(7) =XMAT(11)
        da(20)=XMAT(10)
        da(8) =XMAT(16)
        da(9) =XMAT(17)
        da(10)=XMAT(18)
        da(11)=XMAT(19)
        da(12)=XMAT(20)
        da(13)=XMAT(21)
        da(14)=XMAT(22)
        da(15)=XMAT(23)
        da(16)=XMAT(24)
        da(17)=XMAT(25)
        da(18)=XMAT(26)
        da(19)=XMAT(27)
*
        Diam=da(1)
        XLx=da(2)
        XLy=da(3)
        XKelz=da(4)
        XKelh=da(5)
        XKelry=da(6)
        XKelrx=da(7)
C        a=da(8)
C        b=da(9)
C        c=da(10)
C        d=da(11)
C        e=da(12)
C        f=da(13)
        qmax=da(14)
        a9=da(15)
C        a6=da(16)
        eta3=da(17)
        dtime=da(18)
        a8=da(19)
        XKelt=da(20)
*
*
        upl1(1)=VAR0(2)
        upl1(2)=VAR0(3)
        upl1(3)=VAR0(4)
        upl1(4)=VAR0(5)
        upl1(5)=VAR0(6)
        q1(1)=VAR0(7)
        q1(2)=VAR0(8)
        q1(3)=VAR0(9)
        q1(4)=VAR0(10)
        q1(5)=VAR0(11)
        force_prec(1)=-SIG0(1)
        force_prec(2)=-SIG0(2)
        force_prec(3)=SIG0(6)
        force_prec(4)=SIG0(3)
        force_prec(5)=-SIG0(5)
        delta_prec=VAR0(12)
        Xkapa_prec=VAR0(13)
        Xksi_prec=VAR0(14)
        distmin_prec=VAR0(15)
        upl2(1)=VAR0(17)
        upl2(2)=VAR0(18)
        upl2(3)=VAR0(19)
        upl2(4)=VAR0(20)
        upl2(5)=VAR0(21)
        q2=VAR0(22)
       DO I=1,5
        upl(I)=upl1(I)+upl2(I)
       ENDDO
*
*deplx joint --> -deplz ISS
        uo(1)=-epst0(1)
*deply joint --> -deplx ISS
        uo(2)=-epst0(2)
*rotz  joint --> -roty  ISS
        uo(3)=-epst0(6)
*deplz joint -->  deply ISS
        uo(4)=epst0(3)
*roty  joint --> -rotx  ISS
        uo(5)=-epst0(5)
*
        duo(1)=-DEPST(1)
        duo(2)=-DEPST(2)
        duo(3)=-DEPST(6)
        duo(4)=DEPST(3)
        duo(5)=-DEPST(5)

C Correction données d'entrées Castem
       DO I=1,5
        uo(I)=uo(I)+duo(I)
*        epstf(I) = epst0(I)+DEPST(I)
       ENDDO
*
C   on adimensionne les déplacements
        if (a9.eq.2) then
       DO I=1,5
        uo(I)    = uo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy)
       ENDDO
        uo(3) = uo(3)*XLx
        uo(5) = uo(5)*XLy
       DO I=1,5
        duo(I)   = duo(I)*sqrt(XLx**2+XLy**2)/(XLx*XLy)
       ENDDO
        duo(3)= duo(3)*XLx
        duo(5)= duo(5)*XLy
        else
       DO I=1,5
        uo(I)    = uo(I)/Diam
       ENDDO
        uo(3) = uo(3)*Diam
        uo(5) = uo(5)*Diam
       DO I=1,5
        duo(I)   = duo(I)/Diam
       ENDDO
        duo(3)= duo(3)*Diam
        duo(5)= duo(5)*Diam
        endif
*
C       Adimensionnement des efforts
        if (a9.eq.1) then
C semelle filante
       DO I=1,5
        force_prec(I)=force_prec(I)/(Diam*qmax)
       ENDDO
        force_prec(3)=force_prec(3)/Diam
        force_prec(5)=force_prec(5)/Diam
        else
        if (a9.eq.2) then
C semelle rectangulaire
       DO I=1,5
        force_prec(I)=force_prec(I)/(XLx*XLy*qmax)
       ENDDO
        force_prec(3)=force_prec(3)/XLx
        force_prec(5)=force_prec(5)/XLy
        else
C semelle circulaire
       DO I=1,5
        force_prec(I)=force_prec(I)/(XPI*Diam**2/4*qmax)
       ENDDO
        force_prec(3)=force_prec(3)/Diam
        force_prec(5)=force_prec(5)/Diam
        endif
        endif
C     ***********************************************************************
C     Macro-élément d'ISS
C adimensionnement de la matrice de rigidité élastique suivant la forme de la fondation (paramètre a9)
      call MKADIM(da,XMKel)
         do i=1,5
              do j=1,5
                 XMK(i,j)=XMKel(i,j)
              enddo
         enddo

C     Calcul de l'effort de prédiction élastique en considérant le décollement élastique non-linéaire.
        duP(1)=0
        duP(2)=0
        duP(3)=0
        duP(4)=0
        duP(5)=0

C       Réduction du pas de calcul
        dtimei=dtime/10.
       DO I=1,5
        duoi(I)=1./10.*duo(I)
        force_preci(I)=force_prec(I)
       ENDDO
        delta_preci=delta_prec
        Xkapa_preci=Xkapa_prec
        Xksi_preci=Xksi_prec
        distmin_preci=distmin_prec
      do k=1,10
        q2_preci=q2
*
        Xkapa_k=Xkapa_preci
        Xksi_k=Xksi_preci
       DO I=1,5
        q1_preci(I)=q1(I)
        upl_preci(I)=upl(I)
        uoi(I)=uo(I)-(10.-k)*duoi(I)
        uel(I)=uoi(I)-upl(I)
        duel(I)=duoi(I)-duP(I)
       ENDDO
*si a8=1 on désactive le décollement
       if (a8.eq.1) then
        call MATVE1 (XMK, uel,5,5,Ftr,2)
        delta=delta_prec
        deltamax=1
       else
        call MATVE1 (XMK, duel,5,5,dforce,2)
       DO I=1,5
        force(I)=force_preci(I)+dforce(I)
       ENDDO
        call ELNLIN (da, uel, force, XMKel,
     & delta_preci, Ftr, delta, deltamax)
       endif
C     Calcul des fonctions de charge au point de charge élastique
        call LOADSF (da, Ftr,q1,yield1)
        call LOASFV (Ftr,q2,yield2)
       DO I=1,5
        force(I)=Ftr(I)
       ENDDO
C     Boucle de plasticité (Return Mapping)
        epsilon=1E-9
       if ((yield1.le.epsilon) .and. (yield2.le.epsilon)) then
               DO I=1,5
                force(I)=Ftr(I)
               ENDDO
                Xkapa=Xkapa_k
                Xksi=Xksi_k
                distmin=distmin_preci
       else
            m=0
            do while ((yield1.gt.epsilon) .or. (yield2.gt.epsilon))
               m=m+1
               IF (m.gt.1000) THEN
                 KERRE = 22
                 RETURN
               ENDIF
               IF ((force(1).LT.0).or.(force(1).GE.1)) THEN
                 KERRE = 23
                 RETURN
               ENDIF
                call GRADIF(da,force,q1,Gradfs,Gradfq)
                call GRADIG(da,force, q1, force_preci, Xkapa_k,
     & Xksi_k, Xkapa_preci, Xksi_preci, Gradgs, Xkapa, Xksi)
                call FVEH(da,force,q1, XMKel,
     & Gradgs,distmin_preci,hh1,distmin)
                call GRADFV(force,q2,Gradf2s,Gradf2q)
                call FVEHV(force,q2, XMKel,Gradf2s,hh2)

                call MATVE1(XMK,Gradgs,5,5,XKGradgs,2)
                call MATVE1(XMK,Gradf2s,5,5,XKGradf2s,2)

                call SCALPR(5,Gradfs,XKGradgs,Ho11)
                call SCALPR(5,Gradfs,XKGradf2s,Ho12)
                call SCALPR(5,Gradf2s,XKGradgs,Ho21)
                call SCALPR(5,Gradf2s,XKGradf2s,Ho22)

                call SCALPR(5,Gradfq,hh1,H1)
                H2=Gradf2q*hh2
                H(1,1)=Ho11+H1
                H(1,2)=Ho12
                H(2,1)=Ho21
                H(2,2)=Ho22+H2
                yield(1)=yield1
                yield(2)=yield2
*initialisation de ces variables necessaire pour les tests sur les mécanismes
                ddlambda(1)=1
                ddlambda(2)=1
               if ((yield(1).gt.epsilon).
     & and.(yield(2).gt.epsilon)) then
                call HINV22 (H,Hinv)
                call MATVE1(Hinv,yield,2,2,ddlambda,2)
                ddlambda1=ddlambda(1)
                ddlambda2=ddlambda(2)
               endif
C       Tests sur les mécanismes
               if ((yield(1).le.epsilon).or.(ddlambda(1).le.0)) then
                ddlambda1=0
                ddlambda2=yield2/H(2,2)
                       if (ddlambda2.le.0) then
                        ddlambda2=0
                       endif
               endif
               if ((yield(2).le.epsilon).or.(ddlambda(2).le.0)) then
                ddlambda2=0
                ddlambda1=yield1/H(1,1)
                       if (ddlambda1.le.0) then
                        ddlambda1=0
                       endif
               endif

                ddlambda(1)=ddlambda1
                ddlambda(2)=ddlambda2


                q2=q2-ddlambda(2)*hh2
             DO I=1,5
                q1(I)=q1(I)-ddlambda(1)*hh1(I)
                upl1(I)=upl1(I)+ddlambda(1)*Gradgs(I)
                upl2(I)=upl2(I)+ddlambda(2)*Gradf2s(I)
                upl(I) =upl1(I)+upl2(I)
*
                uel(I)=uoi(I)-upl(I)
                duP(I)=upl(I)-upl_preci(I)
                duel(I)=duoi(I)-duP(I)
             ENDDO
*si a8=1 on désactive le décollement
              if (a8.eq.1) then
                call MATVE1 (XMK, uel,5,5, force,2)
                delta=delta_prec
              else
                call MATVE1 (XMK, duel, 5,5, dforce,2)
               DO I=1,5
                forcepred(I)=force_preci(I)+dforce(I)
               ENDDO
*forcepred et delta_prec sont là pour une initialisation à la boucle de Newton
                call ELNLIN (da, uel, forcepred, XMKel,
     & delta_preci, force, delta, deltamax)
              endif
                Xkapa_k=Xkapa
                Xksi_k=Xksi
                call LOADSF (da,force,q1,yield1)
                call LOASFV (force,q2,yield2)
           enddo
       endif
*
C       Correction visqueuse
* calcul de la norme de la force dans l'hyperplan (pour voir si on est dans la phase de chargement statique poids propre)
       Xnorm_f25=sqrt(force(2)**2+force(3)**2+force(4)**2+force(5)**2)
      if (dtime.ne.0) then
       if ((eta3.ne.0) .and. (Xnorm_f25.gt.0.001)) then
       DO I=1,5
        force(I) = (Ftr(I)+dtime/eta3*force(I))/(1+dtime/eta3)
        upl(I) = (upl_preci(I)+dtime/eta3*upl(I))/(1+dtime/eta3)
        q1(I)  = (q1_preci(I)+dtime/eta3*q1(I))/(1+dtime/eta3)
       ENDDO
        q2  = (q2_preci+dtime/eta3*q2)/(1+dtime/eta3)
       endif
      endif
*
        IF (force(1).GE.1) THEN
           KERRE = 25
           RETURN
        ENDIF
        delta_preci=delta
       DO I=1,5
        force_preci(I)=force(I)
       ENDDO
        Xkapa_preci=Xkapa
        Xksi_preci=Xksi
        distmin_preci=distmin
*
      enddo
*
        call SCALPR(5,upl,upl,EPSE)
        EPSE=sqrt(EPSE)
*
        VARF(1)=EPSE
        VARF(2)=upl1(1)
        VARF(3)=upl1(2)
        VARF(4)=upl1(3)
        VARF(5)=upl1(4)
        VARF(6)=upl1(5)
        VARF(7)=q1(1)
        VARF(8)=q1(2)
        VARF(9)=q1(3)
        VARF(10)=q1(4)
        VARF(11)=q1(5)
        VARF(12)=delta
        VARF(13)=Xkapa
        VARF(14)=Xksi
        VARF(15)=distmin
*deltamax ne sert pas, on la garde juste pour la tracer (provisoire)
        VARF(16)=deltamax
        VARF(17)=upl2(1)
        VARF(18)=upl2(2)
        VARF(19)=upl2(3)
        VARF(20)=upl2(4)
        VARF(21)=upl2(5)
        VARF(22)=q2

C       Redimensionnement des variables
        if (a9.eq.1) then
C semelle filante
       DO I=1,5
        force(I)=force(I)*Diam*qmax
       ENDDO
        force(3)=force(3)*Diam
        force(5)=force(5)*Diam

        else
        if (a9.eq.2) then
C semelle rectangulaire
       DO I=1,5
        force(I)=force(I)*XLx*XLy*qmax
       ENDDO
        force(3)=force(3)*XLx
        force(5)=force(5)*XLy
        else
C semelle circulaire
       DO I=1,5
        force(I)=force(I)*XPI*Diam**2/4*qmax
       ENDDO
        force(3)=force(3)*Diam
        force(5)=force(5)*Diam
        endif
        endif
*
        call MKREDI(da,XMKel)
         do i=1,5
              do j=1,5
                 XMK(i,j)=XMKel(i,j)
              enddo
         enddo
*
*Fz ISS --> -Fx joint
       SIGF(1)=-force(1)
*Fx ISS --> -Fy joint
       SIGF(2)=-force(2)
*Fy ISS -->  Fz joint
       SIGF(3)=force(4)
*Mz ISS --> -Mx joint (torsion linéaire)
       SIGF(4)=-XKelt*epst0(4)
*Mx ISS --> -My joint
       SIGF(5)=-force(5)
*My ISS -->  Mz joint
       SIGF(6)=force(3)
*
      return
      end