// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

/************************************************************************
 *     DATE:        06/03/2023                                          *
 *                                                                $     *
 *     AUTEUR:      G RIO   (mailto:gerardrio56@free.fr)                *
 *                                                                $     *
 *     PROJET:      Herezh++                                            *
 *                                                                $     *
 ************************************************************************
 *     BUT:   Defini l'element generique de thermique.                  *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *                                                                      *
 *     VERIFICATION:                                                    *
 *                                                                      *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *     !        !            !                                    !     *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *
 *     MODIFICATIONS:                                                   *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *                                                                $     *
 ************************************************************************/
#ifndef ELEMTHERMI_H
#define  ELEMTHERMI_H

#include "Element.h"
#include "Tenseur.h"
#include "NevezTenseur.h"
#include "Deformation.h"
#include "Loi_comp_abstraite.h"
#include "Enum_calcul_masse.h"
#include "Basiques.h"
#include "Enum_dure.h"
#include "CompThermoPhysiqueAbstraite.h"
#include "CompFrotAbstraite.h"
#include "LesPtIntegThermiInterne.h"
#include "Enum_StabHourglass.h"
#include "EnergieThermi.h"


/// @addtogroup groupe_des_elements_finis
///  @{
///

class   ElemThermi : public Element
{
  public :
    // VARIABLES PUBLIQUES :
    
    // CONSTRUCTEURS :
       ElemThermi ();
       // Constructeur utile quand le numero de maillage et d'identification de l'element est connu
       ElemThermi (int num_maill,int num_id) ;
       // Constructeur utile quand le numero de maillage et d'identification et le tableau des noeuds 
       // de l'element sont connu
       ElemThermi (int num_maill,int num_id,const Tableau<Noeud *>& tab);
       // Constructeur utile quand le numero de maillage et d'identification est connu,
       // ainsi que la geometrie et le type d'interpolation de l'element	
       ElemThermi (int num_maill,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,string info="");
       // Constructeur utile quand le numero de maillage et d'identification est connu,
       // ainsi que la geometrie et le type d'interpolation de l'element
       ElemThermi (int num_maill,int num_id,char* nom_interpol,char* nom_geom,string info="");
       // Constructeur utile quand toutes les donnees de la classe Element sont connues
       ElemThermi (int num_maill,int num_id,const Tableau<Noeud *>& tab,Enum_interpol id_interp_elt,
										     Enum_geom id_geom_elt,string info="");
       // Constructeur utile quand toutes les donnees de la classe Element sont connues
       ElemThermi (int num_maill,int num_id,const Tableau<Noeud *>& tab,char* nom_interpol,
												char* nom_geom,string info="");
       // Constructeur de copie
       ElemThermi (const ElemThermi& elt);
    
    // DESTRUCTEUR :
       ~ElemThermi ();
    
    // METHODES PUBLIQUES :
        // test si l'element est complet        
        // = 1 tout est ok, =0 element incomplet
		int TestComplet();

	    // calcul si un point est a l'interieur de l'element ou non
        // il faut que M est la dimension globale
        // les trois fonctions sont pour l'etude a t=0, t et tdt
        // retour : =0 le point est externe, =1 le point est interne , 
        //          = 2 le point est sur la frontière à la précision près 
        // coor_locales : s'il est différent de NULL, est affecté des coordonnées locales calculées,
        //                 uniquement précises si le point est interne
        int Interne_0(const Coordonnee& M,Coordonnee* coor_locales=NULL);
        int Interne_t(const Coordonnee& M,Coordonnee* coor_locales=NULL); 
        int Interne_tdt(const Coordonnee& M,Coordonnee* coor_locales=NULL);
        
        // récupération des énergies intégrées sur l'éléments, résultants d'un précédent calcul
        // explicite, ou implicite
        const EnergieThermi& EnergieTotaleElement() const {return energie_totale;};

        // test pour savoir si le calcul de Flux en absolu est possible
        bool FluxAbsoluePossible();
 
    // METHODES VIRTUELLES:
        
        // retourne la liste de tous les  types de ddl interne actuellement utilisés
        // par l'élément (actif ou non), sont exclu de cette liste les ddl des noeuds
        // reliés à l'élément (ddl implique grandeur uniquement scalaire !)
        virtual List_io <Ddl_enum_etendu> Les_type_de_ddl_internes(bool absolue) const ;
        // idem pour les grandeurs évoluées c'est-à-dire directement sous forme de vecteur, tenseurs ....		 
        virtual List_io <TypeQuelconque> Les_type_evolues_internes(bool absolue) const ;
        // idem pour les données particulières
        virtual List_io <TypeQuelconque> Les_types_particuliers_internes(bool absolue) const;
        
        // retourne la liste de toutes les  grandeurs quelconques relatives aux faces de
        // l'élément (actif ou non),
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        virtual List_io <TypeQuelconque> Les_type_quelconque_de_face(bool absolue) const ;

        // retourne la liste de toutes les  grandeurs quelconques relatives aux arêtes de
        // l'élément (actif ou non),
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        virtual List_io <TypeQuelconque> Les_type_quelconque_de_arete(bool absolue) const;

      // --------- calculs utils dans le cadre de la recherche du flambement linéaire		
      // dans un premier temps uniquement virtuelles, ensuite se sera virtuelle pure pour éviter
      // les oubli de définition ----> AMODIFIER !!!!!!!!    
        // Calcul de la matrice géométrique et initiale
        class MatGeomInit // pour le retour des pointeurs sur des matrices stockées
                           // une paire par classe d'éléments
         { public :  MatGeomInit(Mat_pleine * matG,Mat_pleine * matI) :
                       matGeom(matG),matInit(matI) {};
           Mat_pleine * matGeom ;Mat_pleine * matInit ; 
          };
        virtual MatGeomInit MatricesGeometrique_Et_Initiale (const ParaAlgoControle & pa) ;
 
// --------- calcul d'erreur, calculs du champs de Flux continu ---------
        // ajout des ddl de Flux pour les noeuds de l'élément
        virtual void Plus_ddl_Flux() = 0;
         // inactive les ddl du problème de recherche d'erreur : les Flux
        virtual void Inactive_ddl_Flux() = 0;
         // active les ddl du problème de recherche d'erreur : les Flux
        virtual void Active_ddl_Flux() = 0 ;
         // active le premier ddl du problème de recherche d'erreur : FLUX
        virtual void Active_premier_ddl_Flux() = 0 ;
        // retourne un tableau de ddl element, correspondant à la
        // composante de flux -> FLUX, pour chaque noeud qui contiend
        // des ddl de flux
        // -> utilisé pour l'assemblage de la raideur d'erreur
        //!!!!!!!!!!!!! pour l'instant en virtuelle il faudra après en 
        // virtuelle pure !!!!!!!!!!!!!!!!!!!!!!!!!!!
        virtual DdlElement&  Tableau_de_Flux1() const ;
        // retourne un tableau de ddl element, correspondant à la 
        // composante d'erreur -> ERREUR, pour chaque noeud 
        // -> utilisé pour l'assemblage de la raideur d'erreur
        //dans cette version tous les noeuds sont supposés avoi un ddl erreur
        // dans le cas contraire il faut redéfinir la fonction dans l'élément terminal
        virtual DdlElement  Tableau_de_ERREUR() const ;
        // calcul de l'erreur sur l'élément. Ce calcul n'est disponible
        // qu'une fois la remontée aux Flux effectuées sinon aucune
        // action. En retour la valeur de l'erreur sur l'élément
        // type indique le type de calcul d'erreur :
        // = 1 : erreur = (int (delta sigma):(delta sigma) dv)/(int sigma:sigma dv)
        // le numerateur et le denominateur sont tel que :
        // errElemRelative = numerateur / denominateur , si denominateur different de 0
        // sinon denominateur = numerateur si numerateur est different de 0, sinon
        // tous sont nuls mais on n'effectue pas la division 
		//!!!!!!!!!!!!! pour l'instant en virtuelle il faudra après en 
		// virtuelle pure !!!!!!!!!!!!!!!!!!!!!!!!!!!
        virtual void ErreurElement(int type,double& errElemRelative
                            ,double& numerateur, double& denominateur);
        // les 3 routines qui suivent sont virtuelles, car la définition
        //qui est  faite dans ElemThermi.cp considère qu'il y a un ddl erreur
        // par noeud de l'élément de manière systématique, 
        // avec le status virtuel on peut définir dans la classe dérivée
        // un cas particulier
        // ajout des ddl d'erreur pour les noeuds de l'élément
        virtual void Plus_ddl_Erreur() ;
        // inactive les ddl d'erreur 
        virtual void Inactive_ddl_Erreur() ;
        // active les ddl d'erreur
        virtual void Active_ddl_Erreur()  ;
        // test pour savoir si l'erreur a été calculée
        bool ErreurDejaCalculee()
          { if (fluxErreur == NULL)
              return false;
            else return true;};  
        // sortie de l'erreur à l'élément
        double Erreur( )
          { return (*fluxErreur);};
	    
        // lecture de données diverses sur le flot d'entrée
        // l'implantation est faite dans les classe dérivées
        virtual void LectureFlux(UtilLecture * entreePrinc) =0 ;
        
        // retour des Flux en absolu retour true si elle existe sinon false
        virtual bool FluxAbsolues(Tableau <Vecteur>& tabFlux) = 0;
        
// --------- calcul dynamique ---------

        // calcul  de la longueur d'arrête de l'élément minimal
        // divisé par la célérité dans le matériau
        virtual double Long_arrete_mini_sur_c(Enum_dure temps) = 0;
        // cas du bulk viscosity
        double E_elem_bulk_t,E_elem_bulk_tdt,P_elem_bulk;
        static void ActiveBulkViscosity(int choix) {bulk_viscosity=choix;};
        static void InactiveBulkViscosity() {bulk_viscosity=false;};
        static void ChangeCoefsBulkViscosity(const DeuxDoubles & coef) 
            {  c_traceBulk=coef.un;c_carreBulk=coef.deux;}; 
            
        // initialisation pour le calcul de la matrice masse dans le cas de l'algorithme 
        // de relaxation dynamique  avec optimisation en continu de la matrice masse
        // casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse
        void InitCalculMatriceMassePourRelaxationDynamique(int casMass_relax); 
        // phase de calcul de la matrice masse dans le cas de l'algo de relaxation dynamique
        // mi=fonction de (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises
        // ep: epaisseur, K module de compressibilite, mu: module de cisaillement, Isig trace de sigma,  
        // Sig_mises la contrainte de mises 
        // casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse
        void  CalculMatriceMassePourRelaxationDynamique 
                      (const double& alph, const double& beta, const double & lambda
                      ,const double & gamma,const double & theta, int casMass_relax);  
            
// ---------- informations annexes  ----------------

        // récupération de la base locales au noeud noe, pour le temps: temps 
        const BaseB & Gib_elemeca(Enum_dure temps, const Noeud * noe);
        
        // récupération de grandeurs particulières au numéro d'ordre  = iteg  
        // celles-ci peuvent être quelconques
        // en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        void Grandeur_particuliere (bool absolue,List_io<TypeQuelconque>& liTQ,int iteg);
                                 
        // récupération de grandeurs particulières pour une face au numéro d'ordre  = iteg
        // celles-ci peuvent être quelconques
        // en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        void Grandeur_particuliere_face (bool absolue,List_io<TypeQuelconque>& liTQ,int face, int iteg);
        
        // récupération de grandeurs particulières pour une arête au numéro d'ordre  = iteg
        // celles-ci peuvent être quelconques
        // en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        void Grandeur_particuliere_arete (bool absolue,List_io<TypeQuelconque>& liTQ,int arete, int iteg);


        // récupération de la loi de frottement, dans le cas où elle n'existe pas
        // retour d'un pointeur nul
        CompFrotAbstraite* LoiDeFrottement() const {return loiFrot;};

		      // --- transfert des grandeurs des points d'intégration aux noeuds 
		      // transfert de ddl des points d'intégrations (de tous) d'un éléments (on ajoute aux noeuds, on ne remplace pas)
        // les ddl doivent déjà exister aux noeuds sinon erreur
		      // il doit s'agir du même type de répartition de pt d'integ pour toutes les grandeurs
        // tab_val(i)(j) : valeur associée au i ième pt d'integ et au j ième ddl_enum_etendu
		      void TransfertAjoutAuNoeuds(const List_io < Ddl_enum_etendu >& lietendu
                                     ,const Tableau <Tableau <double> > & tab_val,int cas); 
        // transfert de type quelconque des points d'intégrations (de tous) aux noeuds d'un éléments (on ajoute aux noeuds, 
        // on ne remplace pas). Les types quelconques doivent déjà exister
        // un tableau dans tab_liQ  correspondent aux grandeurs quelconque pour tous les pt integ, 
        // tab_liQ(i) pour le pt d'integ i
	    // liQ_travail: est une liste de travail qui sera utilisée dans le transfert                        
        void TransfertAjoutAuNoeuds(const Tableau <List_io < TypeQuelconque > >& tab_liQ
                                            ,List_io < TypeQuelconque > & liQ_travail,int cas);
 
        // accumulation aux noeuds de grandeurs venant de l'éléments vers ses noeuds (exemple la pression appliquée)
	       // autres que celles aux pti classiques, mais directements disponibles
        // le contenu du conteneur stockées dans liQ est utilisé en variable intermédiaire
        void Accumul_aux_noeuds(const List_io < Ddl_enum_etendu >& lietendu
                                ,List_io < TypeQuelconque > & liQ,int cas);
        
        
        // modification de l'orientation de l'élément en fonction de cas_orientation
        // =0: inversion simple (sans condition) de l'orientation
        // si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration
        // 1.  si tous les jacobiens sont négatifs on change d'orientation
        // 2.  si tous les jacobiens sont positifs on ne fait rien
        // 3.  si certains jacobiens sont positifs et d'autres négatifs message
        //     d'erreur et on ne fait rien
        // ramène true: s'il y a eu changement effectif, sinon false
		      bool Modif_orient_elem(int cas_orientation);
  
        // calcul éventuel de la normale à un noeud
        // ce calcul existe pour les éléments 2D, 1D axi, et aussi pour les éléments 1D
        // qui possède un repère d'orientation
        // en retour coor = la normale si coor.Dimension() est = à la dimension de l'espace
        // si le calcul n'existe pas -->  coor.Dimension() = 0
        // ramène un entier :
        //    == 1 : calcul normal
        //    == 0 : problème de calcul -> coor.Dimension() = 0
        //    == 2 : indique que le calcul n'est pas licite pour le noeud passé en paramètre
        //           c'est le cas par exemple des noeuds exterieurs pour les éléments SFE
        //           mais il n'y a pas d'erreur, c'est seulement que l'élément n'est pas ad hoc pour
        //           calculer la normale à ce noeud là
        // temps: indique  à quel moment on veut le calcul
        virtual int CalculNormale_noeud(Enum_dure temps, const Noeud& noe,Coordonnee& coor);
 
	       // calcul si un point est a l'interieur de l'element ou non
        // il faut que M est la dimension globale
        // retour : =0 le point est externe, =1 le point est interne , 
        //          = 2 le point est sur la frontière à la précision près 
        // coor_locales : s'il est différent de NULL, est affecté des coordonnées locales calculées,
        //                 uniquement précises si le point est interne
        int Interne(Enum_dure temps,const Coordonnee& M,Coordonnee* coor_locales=NULL);
        // -- connaissances particulières sur l'élément
        // ramène  l'épaisseur de l'élément 
        // =0. si la notion d'épaisseurs ne veut rien dire pour l'élément
        virtual double Epaisseurs(Enum_dure , const Coordonnee& ) {return 0.;};
        // ramène  l'épaisseur moyenne de l'élément (indépendante du point)
        // =0. si la notion d'épaisseurs ne veut rien dire pour l'élément
        virtual double EpaisseurMoyenne(Enum_dure ) {return 0.;};
        // ramène  la section de l'élément 
        // =0. si la notion de section ne veut rien dire pour l'élément
        virtual double Section(Enum_dure , const Coordonnee& ) {return 0.;};
        // ramène  la section moyenne de l'élément (indépendante du point)
        // =0. si la notion de section ne veut rien dire pour l'élément
        virtual double SectionMoyenne(Enum_dure ) {return 0.;};
    
        // fonction a renseigner par les classes dérivées, concernant les répercutions
        // éventuelles due à la suppression de tous les  frontières
        // nums_i : donnent les listes de frontières supprimées
        virtual void Prise_en_compte_des_consequences_suppression_tous_frontieres();
        // idem pour une frontière (avant qu'elle soit supprimée)        
        virtual void Prise_en_compte_des_consequences_suppression_une_frontiere(ElFrontiere* elemFront);

 // -------------------- calcul de frontières -------------------

        // ramène la frontière point
        // éventuellement création des frontieres points de l'element et stockage dans l'element 
        // si c'est la première fois  sinon il y a seulement retour de l'elements
        // a moins que le paramètre force est mis a true
        // dans ce dernier cas la frontière effacéee est recréée
        // num indique le numéro du point à créer (numérotation EF)
        virtual ElFrontiere* const  Frontiere_points(int num,bool force);
        
        // ramène la frontière linéique
        // éventuellement création des frontieres linéique de l'element et stockage dans l'element
        // si c'est la première fois et en 3D sinon il y a seulement retour de l'elements
        // a moins que le paramètre force est mis a true
        // dans ce dernier cas la frontière effacéee est recréée
        // num indique le numéro de l'arête à créer (numérotation EF)
        virtual ElFrontiere* const  Frontiere_lineique(int num,bool force);
		
        // ramène la frontière surfacique
        // éventuellement création des frontieres surfacique de l'element et stockage dans l'element 
        // si c'est la première fois sinon il y a seulement retour de l'elements
        // a moins que le paramètre force est mis a true
        // dans ce dernier cas la frontière effacéee est recréée
        // num indique le numéro de la surface à créer (numérotation EF)
        virtual ElFrontiere* const  Frontiere_surfacique(int num,bool force);
        
              
        // Calcul spécifiques des frontieres de l'element et retour d'un tableau tabb des frontières
        //  creation des elements frontieres et stockage dans l'element
        // la création n'a lieu qu'au premier appel
        // ou lorsque l'on force le paramètre force a true
        // dans ce dernier cas seul les frontière effacées sont recréée
        // cas :
        //  = 0 -> on veut toutes les frontières
        //  = 1 -> on veut uniquement les surfaces
        //  = 2 -> on veut uniquement les lignes
        //  = 3 -> on veut uniquement les points
        //  = 4 -> on veut les surfaces + les lignes
        //  = 5 -> on veut les surfaces + les points
        //  = 6 -> on veut les lignes + les points
//        virtual Tableau <ElFrontiere*> const & Frontiere_specifique(int cas, bool force = false)
//          {return Frontiere_elethermi(cas,force);};

        // -------------------- init éventuelle avant le chargement -------------------
        // initialisation éventuelle, nécessaire avant d'appliquer l'ensemble des charges
        // par exemple des stockages intermédiaires
        virtual void Initialisation_avant_chargement() {};

//=====================================================================================
  protected : 			
//=====================================================================================
	// METHODES  PROTEGEES utilisables par les classes derivees :
    
        // Calcul des frontieres de l'element
        //  creation des elements frontieres et retour du tableau de ces elements
        // la création n'a lieu qu'au premier appel
        // ou lorsque l'on force le paramètre force a true
        // dans ce dernier cas seul les frontière effacées sont recréée
        // cas :
        //  = 0 -> on veut toutes les frontières
        //  = 1 -> on veut uniquement les surfaces
        //  = 2 -> on veut uniquement les lignes
        //  = 3 -> on veut uniquement les points
        //  = 4 -> on veut les surfaces + les lignes
        //  = 5 -> on veut les surfaces + les points
        //  = 6 -> on veut les lignes + les points
		      Tableau <ElFrontiere*> const & Frontiere_elethermi(int cas, bool force = false);

    // ---------------------------------------------------------------------				
    // cas où l'on intègre que selon une liste (un axe, un plan, un volume)
    // ---------------------------------------------------------------------
        // Calcul du residu local et de la raideur locale,
        //  pour le schema implicite d'ou a l'instant t + dt
        // ddl represente les degres de liberte specifiques a l'element
        // tabDepsBB = vitesse de déformation, tabDeltaEpsBB = incrément de def entre t et t+dt
        // cald_Dvirtuelle = indique si l'on doit calculer la dérivée de la vitesse de déformation virtuelle
        void Cal_implicit (DdlElement & tab_ddl,Tableau <CoordonneeB >& d_gradTB
                 ,Tableau < Tableau2 <CoordonneeB> * > d2_gradTB,Tableau <CoordonneeH >& d_fluxH,int nbint
                 ,const Vecteur& poids,const ParaAlgoControle & pa,bool cald_DGradTvirtuelle);
		
        // Calcul du residu local a l'instant t ou tdt
        // atdt = true : calcul à tdt, valeur par défaut
        //      = false: calcul à t
        // ddl represente les degres de liberte specifiques a l'element
        // nbint = nb de pt d'integration , poids = poids d'integration 
        void Cal_explicit (DdlElement & ddl,Tableau <CoordonneeB >& d_gradTB,int nbint
                           ,const Vecteur& poids,const ParaAlgoControle & pa,bool atdt=true);
		
        // Calcul de la matrice géométrique et de la matrice initiale
        // cette fonction est éventuellement appelée par les classes dérivées 
        // ddl represente les degres de liberte specifiques a l'element
        // nbint = nb de pt d'integration , poids = poids d'integration 
        // cald_Dvirtuelle = indique si l'on doit calculer la dérivée de la vitesse de déformation virtuelle
        void Cal_matGeom_Init (Mat_pleine & matGeom, Mat_pleine & matInit
                 ,DdlElement & ddl,Tableau <CoordonneeB >& d_gradTB,Tableau < Tableau2 <CoordonneeB> * > d2_gradTB
                 ,Tableau <CoordonneeH>& d_fluxH,int nbint,const Vecteur& poids
                 ,const ParaAlgoControle & pa,bool cald_Dvirtuelle);
                
        // Calcul de la matrice masse selon différent choix donné par type_matrice_masse,
        // a l'instant initial. 
        void Cal_Mat_masse (DdlElement & tab_ddl,Enum_calcul_masse type_matrice_masse,
                            int nbint,const Tableau <Vecteur>& taphi,int nbne
                            ,const Vecteur& poids);
                
    // ---------------------------------------------------------------------				
    // cas où l'on intègre  selon deux listes (ex un axe et un plan, etc.)
    // ---------------------------------------------------------------------
        // Calcul du residu local et de la raideur locale,
        //  pour le schema implicite d'ou a l'instant t + dt
        // ddl represente les degres de liberte specifiques a l'element
        void Cal_implicitap (DdlElement & tab_ddl,Tableau <TenseurBB *> & d_epsBB
                ,Tableau < Tableau2 <CoordonneeB> * > d2_gradTB,Tableau <CoordonneeH>& d_fluxH
                ,int nbint1,Vecteur& poids1,int nbint2,const Vecteur& poids2
                ,const ParaAlgoControle & pa);
		
        // Calcul du residu local a l'instant t ou tdt
        // atdt = true : calcul à tdt, valeur par défaut
        //      = false: calcul à t
        // ddl represente les degres de liberte specifiques a l'element
        // d_epsbb = variation des def
        // nbint = nb de pt d'integration , poids = poids d'integration 
        void Cal_explicitap (DdlElement & ddl,Tableau <TenseurBB *>& d_epsBB
                             ,int nbint,const Vecteur& poids,bool atdt=true);

        // Calcul de la matrice géométrique et de la matrice initiale
        // cette fonction est éventuellement appelée par les classes dérivées 
        // ddl represente les degres de liberte specifiques a l'element
        // d_epsbb = variation des def
        // nbint = nb de pt d'integration , poids = poids d'integration 
        void Cal_matGeom_Initap (Mat_pleine & matGeom, Mat_pleine & matInit
                ,DdlElement & ddl,Tableau <TenseurBB *>& d_epsBB,Tableau < Tableau2 <CoordonneeB> * > d2_gradTB
                ,Tableau <CoordonneeH>& d_fluxH,int nbint,const Vecteur& poids);
                
  // -------------------------- calcul de second membre ---------------------------
        // calcul des seconds membres suivant les chargements 
        // cas d'un chargement surfacique, sur les frontières des éléments
        // force indique la force surfacique appliquée
        // retourne  le second membre résultant
        // nSurf : le numéro de la surface externe
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_surf_E (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa,bool atdt=true);
	       // idem SM_charge_surf_E mais  -> implicite,
	       // pa : permet de déterminer si l'on fait ou non le calcul de la contribution à la raideur
	       // retourne le second membre et la matrice de raideur correspondant
        Element::ResRaid SMR_charge_surf_I (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa);
        // calcul des seconds membres suivant les chargements 
        // cas d'un chargement pression, sur les frontières des éléments
        // pression indique la pression appliquée
        // retourne  le second membre résultant
        // nSurf : le numéro de la surface externe
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_pres_E (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,double  pression
                                ,const ParaAlgoControle & pa,bool atdt=true);
	       // idem SM_charge_pres_E mais -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
	       // retourne le second membre et la matrice de raideur correspondant
        Element::ResRaid SMR_charge_pres_I (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,double  pression
                                ,const ParaAlgoControle & pa);
        // cas d'un chargement lineique, sur les arêtes frontières des éléments
        // force indique la force lineique appliquée
        // retourne  le second membre résultant
        // nArete : le numéro de l'arête externe
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_line_E (DdlElement & ddls,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa,bool atdt=true);
	       // idem SM_charge_line_E mais  -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
	       // retourne le second membre et la matrice de raideur correspondant
	       Element::ResRaid SMR_charge_line_I (DdlElement & ddlA,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa);

        // cas d'un chargement lineique suiveur, sur les arêtes frontières des éléments
        // pas valable pour des éléments 3D !
        // force indique la force lineique appliquée
        // retourne  le second membre résultant
        // nArete : le numéro de l'arête externe
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_line_Suiv_E (DdlElement & ddls,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa,bool atdt=true);
	       // idem SM_charge_line_E mais  -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
	       // retourne le second membre et la matrice de raideur correspondant
	       Element::ResRaid SMR_charge_line_Suiv_I (DdlElement & ddlA,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa);

 
        // cas d'un chargement surfacique suiveur, sur les surfaces de l'élément
        // la direction varie selon le système suivant: on définit les coordonnées matérielles 
        // de la direction, ce qui sert ensuite à calculer les nouvelles directions. L'intensité
        // elle  est constante. 
        // force indique la force surfacique appliquée
        // retourne  le second membre résultant
        // nSurf : le numéro de la surface externe
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_surf_Suiv_E (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa,bool atdt=true);
        // idem SM_charge_surf_Suiv_E mais  -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur  
        // retourne le second membre et la matrice de raideur correspondant
	       Element::ResRaid SMR_charge_surf_Suiv_I (DdlElement & ddlA,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa);

        // cas d'un chargement volumique, sur l'élément
        // force indique la force volumique appliquée
        // retourne  le second membre résultant
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_vol_E (DdlElement & ddls
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa,bool atdt=true);

        // idem SM_charge_vol_E mais  -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
        // retourne le second membre et la matrice de raideur correspondant
        void SMR_charge_vol_I (DdlElement & ddls
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa);

        // cas d'un chargement hydrostatique, sur les surfaces de l'élément
        // la charge dépend de la hauteur à la surface libre du liquide déterminée par un point
        // et une  direction  normale à la surface libre: 
        // nSurf : le numéro de la surface externe
        // poidvol: indique le poids volumique du liquide
        // M_liquide : un point de la surface libre
        // dir_normal_liquide : direction normale à la surface libre
        // retourne  le second membre résultant
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        Vecteur& SM_charge_hydro_E (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids
                                ,const Coordonnee& dir_normal_liquide,const double& poidvol
                                ,const Coordonnee& M_liquide
                                ,const ParaAlgoControle & pa,bool atdt=true);
        // idem SM_charge_hydro_E mais  -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur  
        // retourne le second membre et la matrice de raideur correspondant
        Element::ResRaid SMR_charge_hydro_I (DdlElement & ddlA,int nSurf
                                    ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids
                                ,const Coordonnee& dir_normal_liquide,const double& poidvol
                                ,const Coordonnee& M_liquide
                                ,const ParaAlgoControle & pa);

        // cas d'un chargement aero-hydrodynamique, sur les frontières de l'élément 
        // Il y a trois forces: une suivant la direction de la vitesse: de type traînée aerodynamique
        // Fn = poids_volu * fn(V) * S * (normale*u) * u, u étant le vecteur directeur de V (donc unitaire)
        // une suivant la direction normale à la vitesse de type portance
        // Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire, normal à V, et dans le plan n et V
        // une suivant la vitesse tangente de type frottement visqueux
        // T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
        // retourne  le second membre résultant
        // calcul à l'instant tdt ou t en fonction de la variable atdt
        // coef_mul: est un coefficient multiplicateur global (de tout)
        Vecteur& SM_charge_hydrodyn_E (const double& poidvol,const Tableau <Vecteur>& taphi,int nbne
                                         ,Courbe1D* frot_fluid,const Vecteur& poids
                                         ,Courbe1D* coef_aero_n,int numfront,const double& coef_mul
                                         ,Courbe1D* coef_aero_t,const ParaAlgoControle & ,bool atdt=true);
        // idem SM_charge_hydrodyn_E mais  -> implicite,
        // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
        // retourne le second membre et la matrice de raideur correspondant
        Element::ResRaid SM_charge_hydrodyn_I (const double& poidvol,const Tableau <Vecteur>& taphi,int nbne
                                         ,Courbe1D* frot_fluid,const Vecteur& poids,DdlElement & ddls
                                         ,Courbe1D* coef_aero_n,int numfront,const double& coef_mul
                                         ,Courbe1D* coef_aero_t,const ParaAlgoControle & pa);

 // -------------------- calcul de frontières en protected -------------------
    
        //  --- fonction nécessaire pour la construction des Frontières linéiques ou surfaciques particulière à l'élément
        // adressage des frontières linéiques et surfacique
        // définit dans les classes dérivées, et utilisées pour la construction des frontières
        virtual ElFrontiere* new_frontiere_lin(int num,Tableau <Noeud *> & tab, DdlElement& ddelem) = 0;
        virtual ElFrontiere* new_frontiere_surf(int num,Tableau <Noeud *> & tab, DdlElement& ddelem) = 0;
                 
 // -------------------- stabilisation d'hourglass -------------------
        // calcul d'élément de contrôle d'hourglass associée à un comportement
        void Init_hourglass_comp(const ElemGeomC0& elgeHour, const string & str_precision
                                 ,LoiAbstraiteGeneral * loiHourglass,const BlocGen & bloc);
         
        // stabilisation pour un calcul implicit
        void Cal_implicit_hourglass();
       
        // stabilisation pour un calcul explicit
        void Cal_explicit_hourglass(bool atdt);
         
      public :
          // récupération de l'energie d'hourglass éventuelle
        double Energie_Hourglass();
           // idem pour l'énergie et la puissance de bulk viscosity
              double Energie_Bulk() const {return E_elem_bulk_tdt;};
              double Puissance_Bulk() const {return P_elem_bulk;};
      protected :
 // -------------------- calcul d'erreur, remontée des Flux -------------------
 
        // calcul du résidu et de la matrice de raideur pour le calcul d'erreur
        // cas d'une intégration suivant une seule liste
        void FluxAuNoeud_ResRaid(const int nbne,const Tableau <Vecteur>& taphi
                ,const Vecteur& poids,Tableau <Vecteur *>&  resErr,Mat_pleine&  raidErr
                ,const Tableau <Vecteur>& taphiEr,const Vecteur& poidsEr );

        // calcul de l'erreur sur l'élément. Ce calcul n'est disponible
        // qu'une fois la remontée aux Flux effectuées sinon aucune
        // action. pour les autres arguments de retour, idem ErreurElement 
        // qui est la fonction generique, les autres variables sont spécifiques
        // a l'element.
        void Cal_ErrElem(int type,double& errElemRelative
                ,double& numerateur, double& denominateur,const int nbne,const Tableau <Vecteur>& taphi
                ,const Vecteur& poids,const Tableau <Vecteur>& taphiEr,const Vecteur& poidsEr);

        // calcul de l'erreur aux noeuds. Contrairement au cas des Flux
        // seul le résidu est calculé. Cas d'une intégration suivant une liste
        void Cal_ErrAuxNoeuds(const int nbne, const Tableau <Vecteur>& taphi, 
                const Vecteur& poids,Tableau <Vecteur *>&  resErr );
                
                
 // ---------------------------------------------------------------------
        // affichage dans la sortie transmise, des variables duales "nom"
        // aux differents points d'integration
        // dans le cas ou nom est vide, affichage de "toute" les variables
        // cas =1 -> premier passage pour de l'implicit
        // cas = 2 -> premier passage pour de l'explicit
        // cas = 11 -> passage autre que le premier pour de l'implicit
        // cas = 12 -> passage autre que le premier pour de l'explicit
        void VarDualSort(ostream& sort, Tableau<string>& nom,int nbint,int cas);
        // utilitaires de VarDualSort 
        // affiche en fonction d'indic les differentes variables et appel
        // AffDefContiD  en fonction de la dimension i  
        //                  
        void AffDefCont( ostream& sort,CompThermoPhysiqueAbstraite::SaveResul * saveDon,
                         CoordonneeB& gradTB,Coordonnee& gradT,double& norme_gradT,
                         CoordonneeB& DgradTB,Coordonnee& DgradT,double& norme_dGradT,
                         CoordonneeH& fluxDH,Coordonnee & fluxD,double& norme_flux,
                         double& temperature, int indic);
        // cas 1D                               
        void AffDefCont1D( ostream& sort,CompThermoPhysiqueAbstraite::SaveResul * saveDon,
                         CoordonneeB& gradTB,Coordonnee& gradT,double& norme_gradT,
                         CoordonneeB& DgradTB,Coordonnee& DgradT,double& norme_dGradT,
                         CoordonneeH& fluxDH,Coordonnee & fluxD,double& norme_flux,
                         double& temperature, int indic);
        // cas 2D                               
        void AffDefCont2D( ostream& sort,CompThermoPhysiqueAbstraite::SaveResul * saveDon,
                         CoordonneeB& gradTB,Coordonnee& gradT,double& norme_gradT,
                         CoordonneeB& DgradTB,Coordonnee& DgradT,double& norme_dGradT,
                         CoordonneeH& fluxDH,Coordonnee & fluxD,double& norme_flux,
                         double& temperature, int indic);
        // cas 3D                               
        void AffDefCont3D( ostream& sort,CompThermoPhysiqueAbstraite::SaveResul * saveDon,
                         CoordonneeB& gradTB,Coordonnee& gradT,double& norme_gradT,
                         CoordonneeB& DgradTB,Coordonnee& DgradT,double& norme_dGradT,
                         CoordonneeH& fluxDH,Coordonnee & fluxD,double& norme_flux,
                        double& temperature, int indic);
       
       
  // ---------------------------------------------------------------------
       
        // --- méthodes relatives aux calculs d'erreurs -----------
        // ---- utilisées par les classes dérivées ---------
         // ajout des ddl relatif aux Flux pour les noeuds de l'élément
        void Ad_ddl_Flux(const DdlElement&  tab_ddlErr);
         // inactive les ddl du problème primaire de mécanique
        void Inact_ddl_primaire(DdlElement&  tab_ddl);
          // active les ddl du problème primaire de mécanique
        void Act_ddl_primaire(DdlElement&  tab_ddl);
        // inactive les ddl du problème de recherche d'erreur : les Flux
        void Inact_ddl_Flux(DdlElement&  tab_ddlErr);
        // active les ddl du problème de recherche d'erreur : les Flux
        void Act_ddl_Flux(DdlElement&  tab_ddlErr);
        // active le premier ddl du problème de recherche d'erreur : FLUX
        void Act_premier_ddl_Flux();
       
 // ---------------------------------------------------------------------
        // lecture des Flux sur le flot d'entrée
        void LectureDesFlux
              (bool cas,UtilLecture * entreePrinc,Tableau <CoordonneeH *>& tabfluxH);

        // retour des Flux en absolu retour true si ils existent sinon false
        void  FluxEnAbsolues
           (bool cas,Tableau <CoordonneeH *>& tabfluxH,Tableau <Vecteur>& tabflux);
       
	 
	// ---------------- lecture écriture dans base info  ----------------
	//  programmes utilisés par les classes dérivées
	
        // cas donne le niveau de la récupération
           // = 1 : on récupère tout
           // = 2 : on récupère uniquement les données variables (supposées comme telles)
        void Lecture_bas_inf
          (istream& ent,const Tableau<Noeud  *> * tabMaillageNoeud,const int cas) ;
           // cas donne le niveau de sauvegarde
           // = 1 : on sauvegarde tout
           // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
        void Ecriture_bas_inf(ostream& sort,const int cas) ;
 
        // --------- utilitaires pour la dynamique
        // calcul  de la longueur d'arrête de l'élément minimal
        // divisé par la célérité  la plus rapide dans le matériau
        // appelé par les classes dérivées
        // nb_noeud : =0 indique que l'on utilise tous les noeuds du tableau de noeuds
        //        = un nombre > 0, indique le nombre de noeuds à utiliser au début du tableau
        double Interne_Long_arrete_mini_sur_c(Enum_dure temps,int nb_noeud=0);
       
        // recuperation des coordonnées du point d'intégration numéro = iteg pour 
        // la grandeur enu
		      // temps: dit si c'est à 0 ou t ou tdt
        // si erreur retourne erreur à true
        Coordonnee CoordPtInt(Enum_dure temps,Enum_ddl enu,int iteg,bool& erreur); 
        
        // recuperation des coordonnées du point d'intégration numéro = iteg pour
        // la face : face
        // temps: dit si c'est à 0 ou t ou tdt
        // si erreur retourne erreur à true
        Coordonnee CoordPtIntFace(int face, Enum_dure temps,int iteg,bool& erreur);
        // recuperation des coordonnées du point d'intégration numéro = iteg pour
        // la face : face
        // temps: dit si c'est à 0 ou t ou tdt
        // si erreur retourne erreur à true
        Coordonnee CoordPtIntArete(int arete, Enum_dure temps,int iteg,bool& erreur);


        // procesure permettant de completer l'element, en ce qui concerne les variables gérés
        // par ElemThermi, apres sa creation avec les donnees du bloc transmis
        // peut etre appeler plusieurs fois
        Element* Complete_ElemThermi(BlocGen & bloc,LesFonctions_nD*  lesFonctionsnD);
        
        // retourne le numero du pt d'ing le plus près ou est exprimé la grandeur enum
        // temps: dit si c'est à 0 ou t ou tdt
        int PtLePlusPres(Enum_dure temps,Enum_ddl enu, const Coordonnee& M);
                 
        // récupération des  valeurs au numéro d'ordre  = iteg pour
        // les grandeur enu
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        Tableau <double> Valeur_multi(bool absolue,Enum_dure enu_t,const List_io<Ddl_enum_etendu>& enu
                                      ,int iteg,int cas ) ; 
                                
        // récupération des  valeurs Tensorielles (et non scalaire comme avec Valeur_multi)
        // au numéro d'ordre  = iteg pour les grandeur enu
        // enu contiend les grandeurs de retour
        // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
        void Valeurs_Tensorielles(bool absolue,Enum_dure enu_t,List_io<TypeQuelconque>& enu
                                        ,int iteg,int cas ) ;               	

        // ==== >>>> methodes virtuelles definis dans les classes mères ============  
        // ramene l'element geometrique correspondant au ddl passé en paramètre
        virtual ElemGeomC0& ElementGeometrie(Enum_ddl  ddl) const ; 
        // ramène le nombre de grandeurs génératrices pour un pt d'integ, correspondant à un type enuméré
        // peut-être surchargé pour des éléments particuliers 
        virtual int NbGrandeurGene(Enum_ddl enu) const;
        
         // ==== >>>> methodes virtuelles definis dans les classes dérivées ============  
        // ramene la dimension des vecteurs flux et gradient de température de l'élément
        virtual int Dim_flux_gradT() const = 0;
 
    // VARIABLES PROTEGEES :
    
        //---- variables: Flux, déformations etc.. aux points d'intégrations -------
        // lesPtIntegThermiInterne pointe vers une entitée entièrement gérés par les classes dérivées
        // contenant les grandeurs mécaniques (Flux, déformations etc.) stockées aux points d'intégration
        LesPtIntegThermiInterne * lesPtIntegThermiInterne;
    
        Loi_comp_abstraite * loiComp; // loi de comportement mécanique defini dans les classes derivees
        CompThermoPhysiqueAbstraite * loiTP; // éventuellement une loi de comportement thermo physique 
                                             // defini dans les classes derivees
        CompFrotAbstraite * loiFrot;  // éventuellement une loi de comportement au frottement pour ses frontières
                                             
        bool dilatation; // indique si oui ou non on tiend compte de la dilatation thermique                                     
        Met_abstraite * met; // definition specifique dans les classes derivees
        Deformation * def; // definition specifique dans les classes derivees 
                           // mais relative à la déformation mécanique
        Deformation * defEr; // idem que def mais pour la remonte aux Flux
        Tableau <Deformation *> defSurf; // idem mais pour les déformations des surfaces
              // externes (frontières) si cela est nécessaire
        Tableau <Deformation *> defArete; // idem mais pour les déformations des arretes
              // externes (frontières) si cela est nécessaire
        Deformation * defMas; // idem que def mais pour le calcul de la matrice masse
        double* fluxErreur; // erreur sur l'élément dans le calcul des densités de flux;
        Tableau <Loi_comp_abstraite::SaveResul *> tabSaveDon; // donnée particulière à la loi mécanique
        Tableau <CompThermoPhysiqueAbstraite::SaveResul *> tabSaveTP; // donnée parti à la loi thermo physique
        Tableau <Deformation::SaveDefResul *> tabSaveDefDon; // donnée particulière pour la déf
        Tableau <EnergieThermi > tab_energ,tab_energ_t; //les différentes énergies mécaniques mises en jeux
                                          // dimensionné dans les classes dérivées, a t+dt, et a t
        EnergieThermi  energie_totale; // le bilan sur l'élément de l'énergie
        bool premier_calcul_thermi_impli_expli; // au premier calcul soit en implicite, explicite, mat geom
           // toutes les grandeurs à t=0 de la métrique sont calculées et stockées dans la déformation
           // ensuite elles ne sont plus calculées -> premier_calcul_thermi_impli_expli sert pour l'indiquer
        double masse_volumique; // masse volumique de l'élément
		  //--- stabilisation d'hourglass éventuelle
        Enum_StabHourglass type_stabHourglass; // méthode de stabilisation	
        double E_Hourglass;

      private:
		     Tableau <ElemThermi*> tab_elHourglass; // tableau de pointeurs éventuels sur des éléments servant à la stabilisation
									// les éléments sont définies au travers de l'appel à la méthode Cal_hourglass_comp, par une classe dérivée,		                   								 
        double coefStabHourglass ; // coef de stabilisation
        Mat_pleine*  raid_hourglass_transitoire; // raideur transitoire d'hourglass, éventuellement définie
		  		  
        // --------  pour faciliter les routines Interne_t  _0 et _tdt + paramètres de contrôle ----
        static const Coordonnee & (Met_abstraite::*PointM) 
              (const Tableau<Noeud *>& tab_noeud,const Vecteur& phi); 
        static void  (Met_abstraite::*BaseND)
             (const Tableau<Noeud *>& tab_noeud,
                      const Mat_pleine& dphi,const Vecteur& phi,BaseB& bB,BaseH& bH); 
        // ---------- cas du bulk viscosity -------------
      private:  // pour éviter les modifications par les classes dérivées     
        static int bulk_viscosity; // indique si oui ou non on utilise le bulk viscosity
                                   // en plus = 1: fonctionnement normal, =2 fonctionnement quelque soit 
                                   // le signe de IDeps
        static double c_traceBulk; // coeff de la trace de D pour le bulk
        static double c_carreBulk; // coeff du carré de la trace de D pour le bulk
        static TenseurHH * sig_bulkHH; // variable de travail pour le bulk
        // ---------- pour le calcul de la masse pour la méthode de relaxation dynamique
        static Ddl_enum_etendu masse_relax_dyn; // simplement pour le définir une fois pour toute

      protected:
        // pour une utilisation par les classes dérivées
        inline bool Bulk_visco_actif() {return bulk_viscosity;};
        inline double& C_traceBulk() {return c_traceBulk;};
        inline double& C_carreBulk() {return c_carreBulk;};
//        // calcul du Flux modifiée en fonction du bulk viscosity
//        // ramène l'énergie et la puissance par unité de volume, dépensée par le bulk
//        DeuxDoubles  ModifContrainteBulk(const TenseurHH& gijHH_tdt,TenseurHH & sigHH,const TenseurBB & DepsBB_tdt
//                                ,TenseurHH & sigbulkHH); 
        // actualisation des ddl et des grandeurs actives de t+dt vers t
        void TdtversT_();
        // actualisation des ddl et des grandeurs actives de t vers tdt      
        void TversTdt_();

       // ---- utilitaires interne 
      private:
        // calcul du mini  du module d'young équivalent pour tous les points d'intégrations
        double  Calcul_mini_E_qui(Enum_dure temps);
        // calcul  du maxi du module d'young équivalent pour tous les points d'intégrations
        double  Calcul_maxi_E_qui(Enum_dure temps);    
 };
/// @}  // end of group


 
#endif  
