C FLUI3D    SOURCE    FD218221  24/02/07    21:15:12     11834          
       subroutine flui3d(XMAT,NMAT,sig0,sigf,deps,nstrs,VAR0,VARF,
     # NVARI,teta1,teta2,dt,ierr1,iso,mfr,ifou,istep,epst0,epstf,
     # tetaref,NBELAS3D,NB_MAT_BASE,NB_MAT_SUPP,NB_VAR_BASE,NB_VAR_SUPP,
     # NB_DECAL,NBNMAX3D,NBNB3D,IDIMB3D,XE3D,NBVIA3D,INLVIA3D,
     # NMAT0,NMAT1,NMAT2,NMAT3,NVAR0,NVAR1,NVAR2,NVAR3,NB_PARA_FIBRE_U)

c     A Sellier Nov 2023

C     MFR           =  NUMERO DE LA FORMULATION DE L'ELEMENT FINI
C                   = 1       MASSIF
C                   = 33      POREUX
C     IFOU    INDICE DU TYPE DE PROBLEME                                     
C            -1  DEFORMATIONS PLANES                                       
C             0  AXISYMETRIQUE                                                                                     
C             2  TRIDIMENSIONNEL 

C     ISTEP=0 : CALCUL LOCAL
C     ISTEP=1 : PREPARATION DU 1ER PASSGE NON LOCAL
C     ISTEP=2 : DERNIER PASSAGE NON LOCAL EFFECTUE
C     ISTEP=3 : PASSAGE NON LOCAL INTERMEDIAIRE EN NON LOCAL ITERATIF

C     TABLES DE DIMENSION FIXE POUR RESOLUTION DES SYTEMES LINEAIRES 
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

C***********************************************************************
C     PARAMETRES OBLIGATOIRES POUR POINTER CORRECTEMENT
 
C     NMAT0 FIN DES PARAMETRES ELASTIQUES
C     NMAT1 FIN DU MODELE DE LA MATRICE 
C     NMAT2 FIN DES PARAMETRES POUR LES MODULES COMPLEMENTAIRES 
C     NMAT3 FIN DES PARAMETRES POUR LES RENFORTS REPARTIS
C     NMAT  FIN DES PARAMETRES POUR LES CHAMPS DE PHASE VIA HELMHOLTZ 
      INTEGER NMAT0,NMAT1,NMAT2,NMAT3,NMAT

C     NVAR0 0
C     NVAR1 FIN DES VARIABLES INTERNES POUR LA MATRICE
C     NVAR2 FIN DES VARIABLES INTERNES POUR LES MODULES COMPLEMENTAIRES
C     NVAR3 FIN DES VARIABLES INTERNES POUR LES RENFORTS REPARTIS
C     NVARI FIN DES VARIABLES INTERNES POUR LES CHAMPS DE PHASES VIA HELMHOLTZ
      INTEGER NVAR0,NVAR1,NVAR2,NVAR3,NVARI
C***********************************************************************

C******** DECLARATION DES VARIABLES EXTERNES ***************************
      INTEGER NSTRS,NBELAS3D,IERR1,MFR
      REAL*8 XMAT(NMAT),SIG0(NSTRS),SIGF(NSTRS),DEPS(NSTRS)
      REAL*8 EPST0(NSTRS),EPSTF(NSTRS)
      REAL*8 VAR0(NVARI),VARF(NVARI)
      REAL*8 DT,TETA1,TETA2,TETAREF  
C     VARIABLE LOGIQUE POUR MATRICE INIT ISOTROPE
      LOGICAL ISO  
      INTEGER IFOU,ISTEP
      
C     *** NOMBRE DE PARAMETRES MATERIAUX *******************************
      INTEGER NB_MAT_BASE,NB_MAT_SUPP
      
C     *** NOMBRE DE VARIABLES INTERNES ********************************* 
C     NOMBRE DE VARIABLES INTERNES
      INTEGER NB_VAR_BASE,NB_VAR_SUPP

C     *** NOMBRE DE VARIABLES FACULTATIVE ENTRE LES OBLIGATOIRES ET LES
c     HELMHOLTZ
      INTEGER NB_DECAL  

C     *********COORDONNES DES NOEUDS DE L ELEMENT FINI *****************        
C     LIGNE A INCLURE DANS LES MODELES UTILISANT LES NOEUDS     
      INTEGER NBNMAX3D,NBNB3D,IDIMB3D
      REAL*8  XE3D(3,NBNMAX3D)

c     *********nombre de sous types de modeles a coupler ***************
c     nombre de sous type de modele et de variable par tenseur (cf idvar4)
      integer NSTYPE,NDTENS
      parameter (NSTYPE=9,NDTENS=12)
c     nombre de parametres scalaires par sous type de modele
      integer VNMAT(NSTYPE) 
c     nombre de vari scalaires par sous type de modele
      integer VNVARI(NSTYPE,2)

C***********************************************************************
 
c***********************************************************************
c     declarations locales
c     indicateur 1er pas
      logical PPAS
c     teneur en eau3d
      real*8 vw0,vwf
c     module 1er passage
      real*8 E,nu,k0,mu0,k1,mu1
c     numero vari
      integer ivar,numx,numk,numm,nx
c     matrices de Hoocke debut et fin de pas
      real*8 raideur66(6,6),souplesse66(6,6)
c     contraintes initiales dans la matrice
      real*8 sigm03(3),vsigm033(3,3)
c     invariant contraintes initiales
      real*8 sigm0d3(3),devis0,tracs0
c     contraintes effectives
      real*8 sigeff6(6)
c     influence de la temperature 
      real*8 dtherm,kthermv,kthermp
c     increment total de deformation mecanique
      real*8 epst06(6),epstf6(6),deps6(6)
c     influence de l eau capillaire et nano
      real*8 keau,knanov,kmecap
c     coeff de fluage de Maxwell
      real*8 kflum
c     consolidation et temps de maxwell
      real*8 taumc3(3),taum,CC3(3)
c     increment deformation totale
      real*8 deps33(3,3),deps3(3),Vdeps33(3,3)
c     tenseur d orientations
      real*8 O333(3,3,3),DX3(3),VDX33(3,3)
c     increment de deformations elastiques et de contraintes
      real*8 depse3(3),dsig3(3)
c     contraintes
      real*8 sigm06(6),sigmf6(6)
c     increments de deformations anelastiques
      real*8 depsx3(NSTYPE,3)
c     coeff de fluage
      real*8 cflux3(NSTYPE,3)

      
C************* PARAMETRES POUR LES RENFORTS ****************************
c     lignes a inclure dans les modeles utilisants les renforts
c     nb_renf
c     nb_para_par_renf
c     nb_para_renf
c     nb_vari_par_renf
c     nb_vari_renf
C-INC HNBRREN
C     ADAPTE POUR MC3D
      PARAMETER(NB_RENF=0)
      PARAMETER(NB_PARA_PAR_RENF=0)
      PARAMETER(NB_PARA_RENF=0)
      PARAMETER(NB_VARI_PAR_RENF=0)
      PARAMETER(NB_VARI_RENF=0)
c     declaration des renforts du fluendo3d si on veut les utiliser
C-INC HDECREN
C***********************************************************************


C************** DECLARATIONS POUR LES FIBRES ROMAIN ********************
c      declaration des fibres de romain gontero si on veut les utiliser
C-INC HDECFIB
C***********************************************************************


C************** PARAMETRES POUR HELMHOLTZ ******************************
c     lignes a inclure dans les modeles utilisant helmholtz
c     recuperation des dimensions des tableaux pour helmholtz
-INC HNBRHEL
c     declaration des tableaux pour helmholtz
-INC HDECHEL
c     recuperation des parametres d helmholtz
-INC HMATHEL    
      DO NL=1,NB_HELM
c       lecture des variable internes et rangements dans helm0(nl,ii)
c       et helm1(nl,ii) avec ii=1,nb_vari_par_helm
-INC HLVIHEL
      END DO
C***********************************************************************


c******recuperation des parametres pour les fibres si modele de fibre*** 
C-INC HMATFIB
c***********************************************************************


c*************** recuperation des parametres pour les renforts *********
c-INC HMATREN
c***********************************************************************


c*************** NOMBRE DE PARAMETRES PAR SOUS TYPE ********************

c     a actualiser en fonction du contenu de idvisc
c     numero du sous type de modele  1   2   3   4   5   6  7  8   9
c     type de sous modele :          ELA,MAX,FLU,KEL,TRA,DP,CC,RAG,RSI 
                         data VNMAT /8  ,6  ,4  ,2  ,5  ,6 ,3 ,9  ,9  /

c************** NOMBRE DE VARI ET TYPE DE VARI PAR SOUS TYPE ***********

c     nombre de variables scalaire par sous type
c     numero du sous type de modele  1   2   3   4   5   6  7  8   9
c     type de sous modele :          ELA,MAX,FLU,KEL,TRA,DP,CC,RAG,RSI 
      data (VNVARI(I,1),I=1,NSTYPE) /9  ,0  ,2  ,0  ,3  ,2 ,3 ,5  ,5  /,
c     nombre de variables tensorielles par sous type
c     numero du sous type de modele  1   2   3   4   5   6  7  8   9
c     type de sous modele :          ELA,MAX,FLU,KEL,TRA,DP,CC,RAG,RSI 
     #     (VNVARI(I,2),I=1,NSTYPE) /1  ,1  ,1  ,1  ,1  ,1 ,1 ,1  ,1  /

c      print*,'debut de mc3d'
c      verification des varf
C      do i=1,nvari
c         print*,'var0',10,'=',var0(10)
C      end do
c      read*


c******* INITIALISATIONS ***********************************************

c     indicateur d erreur     
      ierr1=0 

c     copie des var0 sur les varf avec diluation des vari scalaire
      call init3d(xmat,nmat,vnmat,nstype,nbelas3d,var0,varf,
     # nvari,vnvari,ndtens,ierr1)
c      print*, 'varf apres init'
c      print*,'varf',10,'=',varf(10)
         
c     initialisation de l indicateur de premier passage dans le modele
      if (var0(1).ne.1.) then
         ppas=.true.
      else
         ppas=.false.
      end if 
c     traitements particuliers lors du premier passage
      if(ppas) then
c         controle des dimensions des listes de parametres 
          ntest=nbelas3d
          do i=1,NSTYPE
            ntest=ntest+VNMAT(i)
          end do
          if(ntest.ne.nmat1) then
             print*,'ERREUR0 dans mc3d'
             print*,'nmat1',nmat1
             print*,'ntest',ntest
             ierr1=1
             return
          end if
c         controle des dimensions des listes de variables internes
          ntest=0
          do i=1,NSTYPE
            ntest=ntest+VNvari(i,1)+VNvari(i,2)*ndtens
          end do
          if(ntest.ne.nvar1) then
             print*,'ERREUR1 dans mc3d'
             ierr1=1
             return
          end if
c         initialisation des module dans les vari pour le suivi d hydratation
          call xmat3d(E,xmat,nmat,vnmat,nstype,nbelas3d,0,1)
          call xmat3d(nu,xmat,nmat,vnmat,nstype,nbelas3d,0,2)
          k1=E/(3.d0*(1.d0-2.d0*nu))
          mu1=E/(2.d0*(1.d0+nu))
c          print*,E,nu,k1,mu1
          call xvar3d(k0,numk,var0,nvari,vnvari,nstype,ndtens,1,7)
          call xvar3d(mu0,numm,var0,nvari,vnvari,nstype,ndtens,1,8)
          var0(numk)=k1
          var0(numm)=mu1
      end if

C******* TENEUR EN EAU ************************************************* 

      if(mfr.ne.33) then
c       actualisation de la teneur en eau debut de pas pour le fluage
c       et la capillarite
        call xvar3d(vw0,numx,var0,nvari,vnvari,nstype,ndtens,1,2)
        call xmat3d(vwf,xmat,nmat,vnmat,nstype,nbelas3d,1,3)
        if (ppas) then
c           on initialise la teneur en eau pour le 1er pas
            var0(numx)=vwf
            vw0=vwf
        end if
      else
        call xvar3d(vw0,numx,var0,nvari,vnvari,nstype,ndtens,1,2)
        print*,'ERREUR2 Implementation poreux pour calcul vw ds mc3d'
        print*,'Actualiser la vari numro:', numx,' avec vwf'
        ierr1=1
        return
      end if
c     stockage pour le prochain pas
      varf(numx)=vwf

c******* CONTRAINTES INITIALES *****************************************

c     ** calcul elastique des contraintes principales debut de pas *****

c     matrice de Hoocke
      call hook3d(iso,xmat,nmat,nbelas3d,ierr1,raideur66,souplesse66)

c     contraintes principales en fonction des deformations elastiques
      call sigp3d(nstype,var0,nvari,vnvari,ndtens,iso,raideur66,
     # sigm03,vsigm033)

c      print*,'mc3d sigm03 av relaxation',sigm03
c      print*,'mc3d vsigm033',vsigm033

c     relaxation de la nouvelle contrainte principale en cas 
c     d augmentation des modules et actualisation des def visco elastiques 
      call hydr3d(xmat,nmat,vnmat,nstype,nbelas3d,var0,varf,
     # nvari,vnvari,ndtens,iso,sigm03,souplesse66)

c      print*,'mc3d sigm03 apres relaxation',sigm03

c     ** 1-1 actualisation des parametres de fluage ********************

c     prise en compte de la temperature
      call ther3d(xmat,nmat,vnmat,nstype,var0,varf,nvari,vnvari,
     # ndtens,nbelas3d,teta1,teta2,tetaref,dtherm,kthermv,kthermp)

c     prise en compte de l eau capillaire
      call eau3d(xmat,nmat,vnmat,nstype,var0,varf,nvari,vnvari,
     # ndtens,nbelas3d,teta1,teta2,tetaref,mfr,keau)

c     prise en compte de la surpression dans la nanoporosité
      call nano3d(xmat,nmat,vnmat,nstype,var0,varf,nvari,vnvari,
     # ndtens,nbelas3d,teta1,teta2,tetaref,ppas,dt,kthermv,knanov)

c     Amplification du fluage entre Rc/3 et Rc
      call kmkf3d(xmat,nmat,vnmat,nstype,nbelas3d,sigm03,kmecap)

c     ** 1-2 relaxation visco elastique a deformation constante ********

c     calcul des coefficients de fluage dans la direction des contraintes
c     initiales principales
      call cofl3d(xmat,nmat,vnmat,nstype,var0,varf,nvari,vnvari,
     # ndtens,nbelas3d,vsigm033,dt,kthermv,kthermp,kmecap,keau,
     # knanov,cflux3)
C       print*,'mc3d avant rela coeff de fluage '
C       do ityp=1,4
C         do j=1,3
C             print*,ityp,j,cflux3(ityp,j)
C         end do
C       end do


c     relaxation visco elastique
c      print*, 'varf avant relaxation'
c      print*,'varf',10,'=',varf(10)
      call rela3d(xmat,nmat,vnmat,nstype,nbelas3d,var0,varf,nvari,
     # vnvari,ndtens,ierr1,iso,cflux3,raideur66,sigm03,vsigm033)
c      print*, 'varf apres relaxation'
c      print*,'varf',10,'=',varf(10)

c     expression de sigm0 en base fixe
      call orie3d(vsigm033,O333)
      do i=1,6
        sigm06(i)=0.d0
      end do
c      print*,'aff0 ds mc3d',sigm06,sigm03,O333
      call incv3d(sigm06,sigm03,O333)
c      print*,'mc3d sigm03 ap relaxation des contraintes init',sigm03
c      print*, 'varf apres reeevaluation contrainte 1'
c      print*,'varf',10,'=',varf(10)


c******* INCREMENT TOTAL DE DEFORMATION MECANIQUE **********************

c     chargement increment de deformation imposee et deformation finale
c     remarque 4-5-6 sont des gama jusqu ici, ensuite des epsilons
      if(nstrs.lt.6) then
         do i=1,nstrs
             deps6(i)=deps(i)
             epstf6(i)=epstf(i)
             epst06(i)=epst0(i)             
         end do
         do i=nstrs+1,6
             deps6(i)=0.d0
             epstf6(i)=0.d0
             epst06(i)=0.d0
         end do
      else
         do i=1,6
             deps6(i)=deps(i)
             epstf6(i)=epstf(i)
             epst06(i)=epst0(i)
         end do      
      end if
c     passage gama en epsilon
      do i=4,6
         deps6(i)=0.5d0*deps6(i)
         epstf6(i)=0.5d0*epstf6(i)
         epst06(i)=0.5d0*epst06(i)         
      end do

c******* TIR VISCO ELASTIQUE *******************************************

c     diagonalisation du tenseur des increments de deformation
      call x6x33(deps6,deps33)
c     valeurs principales des increments
      call b3_v33(deps33,deps3,Vdeps33)
c     tenseurs d orientations des increments
      call orie3d(Vdeps33,O333)
c      print*,'aff1 cm3d',deps33,deps3,Vdeps33

c     coeff visco elastique dans la base principele de l increment
      call cofl3d(xmat,nmat,vnmat,nstype,var0,varf,nvari,vnvari,
     # ndtens,nbelas3d,Vdeps33,dt,kthermv,kthermp,kmecap,keau,
     # knanov,cflux3)
C       print*,'aff2 cm3d'
C       do ityp=1,4
C         do j=1,3
C             print*,ityp,j,cflux3(ityp,j)
C         end do
C       end do
C       read*
c     resolution visco elastique en base principales
      do i=1,3
c        deformation elastique (ityp=1)
         depsx3(1,i)=deps3(i)/(1.d0+cflux3(1,i))
c        deformations visqueuses (itype=2,3,4)
         do ityp=2,4
            depsx3(ityp,i)=depsx3(1,i)*cflux3(ityp,i)
         end do
c        deformations plastiques (itype=5,..,nstype)
         do ityp=5,nstype
            depsx3(ityp,i)=0.d0
         end do
      end do
c     mise a jour des increments de deformation visco elastiques
      do ityp=1,4
        do i=1,3
            dx3(i)=depsx3(ityp,i)
        end do
c        print*,'info envoyee a majt3d',ityp,1
c        print*,'valeurs',DX3
c        print*,'orientatio',O333
        call majt3d(DX3,O333,varf,nvari,vnvari,nstype,
     #  ndtens,ityp,1)
      end do
c      print*, 'varf apres increment visco elastique'
c      print*,'varf',10,'=',varf(10)
c     increments de contraintes
      do i=1,3
         dsig3(i)=0.d0
         do j=1,3
            dsig3(i)=dsig3(i)+raideur66(i,j)*depsx3(1,j)
         end do
      end do

c     nouvelle valeur des contraintes
      do i=1,6
        sigmf6(i)=sigm06(i)
      end do
      call incv3d(sigmf6,dsig3,O333)
c     print*,'aff3 mc3d apres tir visco elastique sigmf6',sigmf6
c     print*, 'varf apres increment contrainte visco elastique'
c     print*,'varf',10,'=',varf(10)

c******* EVALUATION DES CRITERES ***************************************

c     chargement des tenseurs de plasticité pour les ecrouissages
c      do ityp=5,nstype
c        call tvar3d(dx3,vdx33,nx,var0,nvari,vnvari,nstype,ndtens,ityp,1)        
c      end do

c     *** re-mise a zero de tous les increments de deformation *********
      do ityp=1,nstype
        do i=1,3
            depsx3(ityp,i)=0.d0
        end do
      end do

c     *** evaluation et tri des criteres actifs ************************
      na=0

c     *** gestion de l ecoulement plastique ****************************
      if(na.ne.0) then
c       *** ecoulement plastique ***************************************
        print*,'Ecoulement plastique'

c       *** mise a jour des deformations *******************************
        do ityp=1,nstype
            do i=1,3
                dx3(j)=depsx3(ityp,i)
            end do
            call majt3d(dx3,O333,varf,nvari,vnvari,nstype,
     #      ndtens,ityp,1)
        end do
      end if
c     print*,'varf apres plasticite'
c     print*,'varf',10,'=',varf(10)

c****** FIN DE L ECOULEMENT PLASTIQUE **********************************

c****** TRANSFERT DES CONTRAINTES VERS LES POINTS DE GAUSS *************
      do i=1,nstrs
             sigf(i)=sigmf6(i)
c             print*,'sigf(',i,')=',sigf(i)
      end do 
c      read*
   
     
c     indicateur de fin de 1er pas 
      if(ppas.and.((istep.eq.2).or.(istep.eq.0))) then     
             varf(1)=1.d0
      else
             if(ppas) then
                 varf(1)=0.d0
             else
                 varf(1)=var0(1)
             end if
      end if

c***********************************************************************
c     VARIABLES INTERNES A ACTUALISER POUR LES FORMULATIONS DE HELMHOLTZ
C***********************************************************************
c     sauvegarde des vari de helmholtz HELM1(NL,ii=1..NBR_VARI_HELM) 
      do NL=1,NB_HELM
-INC HEVIHEL   
      end do
c      print*,'fin de mc3d'
c      verification des varf
C       do i=1,nvari
c         print*,'varf',10,'=',varf(10)
C       end do
c       read*
        
      return
      end
c     fin         
c***********************************************************************