// 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/>.

//#include "Debug.h"
#include "ElemThermi.h"
#include <iomanip>
#include "ConstMath.h"
#include "Util.h"
#include "Coordonnee2.h"
#include "Coordonnee3.h"
#include "CharUtil.h"


// actualisation des ddl et des grandeurs actives de t+dt vers t
// appelé par les classes dérivées      
void ElemThermi::TdtversT_()
	{// on met à jour l'indicateur de premier calcul
	 // s'il y a des sauvegardes de grandeur aux déformations
	 // on ne regarde que le premier élément de tableau, a priori
	 // il y a toujours un pt d'integ et l'organisation est identique pour tous les pt d'integ
	 if (tabSaveDefDon(1) != NULL)
        premier_calcul_thermi_impli_expli=false;
	 // cas des énergies
	 int nbi= tab_energ.Taille();
     for (int ni=1;ni<= nbi; ni++)
        tab_energ_t(ni) = tab_energ(ni);
     E_elem_bulk_t = E_elem_bulk_tdt; // énergie due au bulk viscosity
	};
// actualisation des ddl et des grandeurs actives de t vers tdt      
// appelé par les classes dérivées      
void ElemThermi::TversTdt_()
	{// on met à jour l'indicateur de premier calcul
	 // on considère que si l'on revient en arrière, il vaut mieux re-initialiser les
	 // grandeurs correspondantes au premier_calcul
	 // s'il y a des sauvegardes de grandeur aux déformations
	 if (tabSaveDefDon(1) != NULL)
         premier_calcul_thermi_impli_expli=true;
	 // cas des énergies
	 int nbi= tab_energ.Taille();
     for (int ni=1;ni<= nbi; ni++)
        tab_energ(ni) = tab_energ_t(ni);
     E_elem_bulk_tdt = E_elem_bulk_t; // énergie due au bulk viscosity
	};
	    
// 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 ElemThermi::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)
   {  int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
//      // dimension des tenseurs
//      int dim = lesPtIntegThermiInterne->DimTens();
//      // inialisation du second membre et de la raideur
//      int nbSM = resErr.Taille(); // nombre de second membre
//      for (int j =1; j<= nbSM; j++)  //  boucle sur les seconds membres
//           (*resErr(j)).Zero();
//      raidErr.Zero();
//      // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit
//      // être compatible avec l'ordre des ddl
//      Tableau2 <int> ordre = OrdreContrainte(nbSM);
//      // création d'un tenseur au dimension absolu pour le calcul des contraintes
//      // dans la base absolue
//      TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ;
//       // ====== calcul du second membre ======= 
//      int ni; // compteur globale de point d'integration
//      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
//                                  // donc il faut calculer tous les éléments de la métrique
//      for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
//       {PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
//        // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//        // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//        const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
//        // passage dans le repère absolu du tenseur contrainte final
//        sigHH = (*(ptIntegThermi.SigHH_t())).BaseAbsolue(sigHH,*ex.giB_t);
//        for (int ne =1; ne<= nbne; ne++)  // 1ere boucle sur les noeuds
//          {  // les résidus : mais il faut suivre l'ordre de l'enregistrement des ddl
//             for (int itot = 1; itot<= nbSM; itot++)
//               { int ix = (int) (ordre(itot,1)); int iy = (int) (ordre(itot,2));
//                (*resErr(itot))(ne) += taphi(ni)(ne)*sigHH(ix,iy) * (poids(ni) * (*ex.jacobien_t));
//                }
//            }
//          }  
//       // ====== calcul de la raideur c'est à dire du hessien ========  
//       // boucle sur les pt d'integ spécifiques à l'erreur
//      for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
//       {
//       // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//       // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//        const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
//        for (int ne =1; ne<= nbne; ne++)  // 1ere boucle sur les noeuds
//           for (int me =1; me<= nbne; me++)  // 2ere boucle sur les noeuds
//                   raidErr(ne,me) += taphiEr(ni)(ne) * taphiEr(ni)(me) * poidsEr(ni) * (*ex.jacobien_t);
//          }  
//      // liberation des tenseurs intermediaires
//      TenseurHH * ptsigHH = &sigHH; delete ptsigHH;
//      LibereTenseur();
//      // calcul de l'erreur relative
    
   };

// calcul de l'erreur sur l'élément. Ce calcul n'est disponible
// qu'une fois la remontée aux contraintes effectuées sinon aucune
// action. En retour la valeur de l'erreur sur l'élément
// = 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 , les autres variables sont spécifiques
 // a l'element.
void ElemThermi::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) 
  { int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
//    // dimension des tenseurs
//    int dim = lesPtIntegThermiInterne->DimTens();
//    bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
//                                  // donc il faut calculer tous les éléments de la métrique
//    switch (type)
//	{case 1 : // cas du calcul aux moindres carrés
//
//     {// création d'un tenseur au dimension absolu pour le calcul des contraintes
//      // dans la base absolue, on le choisit HB  pour le double produit contracté
//      // mais en absolu la variance n'a pas d'importance
//      TenseurHB& sigHB = *(NevezTenseurHB(dimAbsolue)) ;
//      // idem au point d'intégration et un tenseur nul pour l'initialisation
//      TenseurHB& signiHB = *(NevezTenseurHB(dimAbsolue)) ;
//      TenseurHB& sigiHB = *(NevezTenseurHB(dimAbsolue)) ; // tenseur de travail
//      TenseurHB& nulHB = *(NevezTenseurHB(dimAbsolue)) ;
//       // ====== calcul des termes de l'erreur =======
//       // ---- tout d'abord on parcourt les points d'intégration de la mécanique
//       numerateur = 0.; denominateur = 0.; // initialisation  
//      int ni; // compteur globale de point d'integration
//      for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
//       {PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
//        // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//        // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//        const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
//        // passage dans le repère absolu du tenseur contrainte initiale
//        sigHB = (*(ptIntegThermi.SigHH_t())).BaseAbsolue(sigHB,*ex.giB_t);
//        // calcul du denominateur
//        denominateur += (sigHB && sigHB) * (poids(ni) * (*ex.jacobien_t));
//        // calcul de la premiere partie du numerateur, celle qui dépend des points d'intégration
//        // mécanique. 
//           // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, 
//        signiHB = nulHB; // initialisation   
//        for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
//           signiHB += taphi(ni)(ne) * ((tab_noeud(ne))->Contrainte(sigiHB));
//           // 2) intégrale de la partie dépendant de ni
//        numerateur += denominateur - 2 * (signiHB && sigHB) * (poids(ni) * (*ex.jacobien_t))  ; 
//        } 
//       // ---- on parcourt maintenant les points d'intégration pour le calcul d'erreur
//       // ---- ce qui permet de calculer la deuxième partie du numérateur
//       // boucle sur les pt d'integ spécifiques à l'erreur
//      for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
//       {
//       // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//       // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//        const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
//           // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, 
//        signiHB = nulHB; // initialisation   
//        for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
//           signiHB += taphiEr(ni)(ne)*(tab_noeud(ne))->Contrainte(sigiHB);
//           // 2) intégrale de la partie dépendant de ni
//        numerateur += (signiHB && signiHB) * poidsEr(ni) * (*ex.jacobien_t);
//        }  
//      // liberation des tenseurs intermediaires
//      TenseurHB * ptsigHB = &sigHB; delete ptsigHB;
//      ptsigHB = &signiHB; delete ptsigHB;ptsigHB = &sigiHB; delete ptsigHB;
//      ptsigHB = &nulHB; delete ptsigHB;
//      LibereTenseur();
//      // enregistrement de l'erreur pour l'élément
//      if (sigErreur == NULL) sigErreur = new double;        
//	  *sigErreur = numerateur; 
//	  break;
//	  }
//	 default :
//		{   cout << "\nErreur : valeur incorrecte du type de calcul d'erreur !\n";
//			cout << "ElemThermi::ErreurElement(int type  , type = \n" << type;
//			Sortie(1);
//	     }		
//    };
//    // --- calcul de l'erreur relative
//    if (denominateur <= ConstMath::trespetit) 
//       // cas d'un champ de contraintes nulles initialement
//        if (numerateur <= ConstMath::trespetit) 
//           // cas également du numérateur nul
//            errElemRelative = 0.;
//        else
//           // on fait denominateur = numérateur -> erreur relative = 1
//            errElemRelative = 1.;
//    else
//      // cas du dénominateur non nul
//      errElemRelative = numerateur / denominateur;                       
   };
	    

// calcul de l'erreur aux noeuds. Contrairement au cas des contraintes
// seul le résidu est calculé. Cas d'une intégration suivant une liste
void ElemThermi::Cal_ErrAuxNoeuds(const int nbne, const Tableau <Vecteur>& taphi,
                const Vecteur& poids,Tableau <Vecteur *>&  resErr )
   {  int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
//      // inialisation du second membre, on utilise le premier vecteur uniquement
//      (*resErr(1)).Zero();
//       // ====== calcul du second membre ======= 
//      int ni; // compteur globale de point d'integration
//      double volume = 0.; // le volume de l'élément
//      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
//                                  // donc il faut calculer tous les éléments de la métrique
//      for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
//       {
//       // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//       // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//        const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
//        for (int ne =1; ne<= nbne; ne++)  // 1ere boucle sur les noeuds
//             (*resErr(1))(ne) += taphi(ni)(ne)*(*sigErreur) * (poids(ni) * (*ex.jacobien_t)); 
//       // calcul du volume
//       volume += (poids(ni) * (*ex.jacobien_t));              
//       } 
//      // on relativise par rapport au volume, car initialement sigErreur représente l'erreur
//      // totale sur tout l'élément. Donc on divise par le volume pour retrouver après résolution
//      // une valeur au noeud qui représente une valeur ponctuelle et non une valeur qui
//      // qui est relative à un volume
//      *resErr(1) /= volume;  
////      *resErr(1) = -sigErreur;  
//      LibereTenseur();
//      // calcul de l'erreur relative
    
   };
        
    
// ajout des ddl relatif aux contraintes pour les noeuds de l'élément
void ElemThermi::Ad_ddl_Flux(const DdlElement&  tab_ddlErr)
  { int nbne = tab_noeud.Taille();
//    for (int ne=1; ne<= nbne; ne++) // pour chaque noeud
//      { // création du tableau de ddl
//        int tab_ddlErr_Taille = tab_ddlErr(ne).tb.Taille(); // nb de ddl du noeud Thermi
//        Tableau <Ddl> ta(tab_ddlErr_Taille);
//        for (int i =1; i<= tab_ddlErr_Taille; i++)
//	        ta(i) = Ddl (tab_ddlErr(ne).tb(i),0.,LIBRE);
//	    // ajout du tableau dans le noeud    
//        tab_noeud(ne)->PlusTabDdl(ta);
//       } 
   };
       
// inactive les ddl du problème primaire de mécanique
void ElemThermi::Inact_ddl_primaire(DdlElement&  tab_ddl)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_hors_service(tab_ddl(ne).tb);
   };
// active les ddl du problème primaire de mécanique
void ElemThermi::Act_ddl_primaire(DdlElement&  tab_ddl)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(tab_ddl(ne).tb);
   };
// inactive les ddl du problème de recherche d'erreur : les contraintes
void ElemThermi::Inact_ddl_Flux(DdlElement&  tab_ddlErr)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_hors_service(tab_ddlErr(ne).tb);
   };
// active les ddl du problème de recherche d'erreur : les contraintes
void ElemThermi::Act_ddl_Flux(DdlElement&  tab_ddlErr)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(tab_ddlErr(ne).tb);
   };
// active le premier ddl du problème de recherche d'erreur : FLUX1
void ElemThermi::Act_premier_ddl_Flux()
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(FLUXD1);
   };
    
// retourne un tableau de ddl element, correspondant à la 
// composante de densité de flux -> FLUXD1, pour chaque noeud qui contiend
// des ddl de contrainte
// -> utilisé pour l'assemblage de la raideur d'erreur
DdlElement&  ElemThermi::Tableau_de_Flux1() const
  {       cout << "\n erreur, fonction non defini pour cette element "
               << "\n ElemThermi::Tableau_de_Flux1()" << endl;
          Sortie(1); 
          DdlElement * toto = new DdlElement();
          return *toto;    
   };
   
// lecture des flux sur le flot d'entrée
void ElemThermi::LectureDesFlux(bool cas,UtilLecture * entreePrinc,Tableau <CoordonneeH *>& tabfluxH)
   { // dimensionnement de la metrique identique au cas d'un calcul explicite
       if( cas)
        { Tableau<Enum_variable_metrique> tab(7);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         met->PlusInitVariables(tab) ;
        };  
     int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
     // dimension des flux
     int dim = (*(tabfluxH(1))).Dimension();
     // création d'un Coordonnee au dimension absolu pour récupérer les flux
     // dans la base absolue
     Coordonnee flux(dimAbsolue) ;
     int ni; // compteur globale de point d'integration
     bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                              // donc il faut calculer tous les éléments de la métrique
     for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
       { // calcul des éléments de la métrique,  on utilise le même calcul
         // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
         const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
         for (int i=1;i<= dimAbsolue;i++)
           // récupération des coordonnées du flux en absolu
            *(entreePrinc->entree)  >> flux(i);
         // passage dans le repère local du flux final
         CoordonneeH & fluxH = *tabfluxH(ni); // pour simplifier
         fluxH.Change_dim(dim);
         for (int i = 1;i<=dim;i++)
          fluxH(i) = ex.giH_t->Coordo(i) * flux;
        }
   };
   

// retour des flux en absolu retour true si ils existes sinon false
void  ElemThermi::FluxEnAbsolues
           (bool cas,Tableau <CoordonneeH *>& tabfluxH,Tableau <Vecteur>& tabflux)
   { // dimensionnement de la metrique identique au cas d'un calcul explicite
       if( cas)
        { Tableau<Enum_variable_metrique> tab(7);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         met->PlusInitVariables(tab) ;
        };  
     int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
     // dimension des tenseurs
     int dim = (*(tabfluxH(1))).Dimension();
	    // redimensionnement éventuel du tableau de sortie
	    int nbi = tabfluxH.Taille();
	    if (tabflux.Taille() != nbi)
	       tabflux.Change_taille(nbi);
	    for (int ni=1;ni<= nbi;ni++)
	       if ( tabflux(ni).Taille() != dimAbsolue)
	         tabflux(ni).Change_taille(dimAbsolue);
     int ni; // compteur globale de point d'integration
     bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                               // donc il faut calculer tous les éléments de la métrique
     for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
        { // calcul des éléments de la métrique,  on utilise le même calcul
         // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
         const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
         // passage dans le repère global du flux local
         Coordonnee flux; // init à 0
         for (int i=1;i<= dim;i++)
           flux += (*tabfluxH(ni))(i) * (ex.giB_t->Coordo(i));
        };
   };
   
// ---------------- 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 ElemThermi::Lecture_bas_inf
  (ifstream& ent,const Tableau<Noeud  *> * tabMaillageNoeud,const int cas) 
 { // appel de la routine d'élément
   Element::Lect_bas_inf_element(ent,tabMaillageNoeud,cas);switch (cas)
  { case 1 : // ------- on récupère tout -------------------------
     { string toto,nom;
       // récup de la masse volumique
       ent >> toto >> masse_volumique ; 
       // données thermique
       ent >> toto >> dilatation;
		 // blocage éventuelle d'hourglass
		 ent >> toto >>  nom;
		 type_stabHourglass=Id_Nom_StabHourglass(nom.c_str());
	   break;
	 }
    case  2 : // ----------- lecture uniquement de se qui varie --------------------
     { 	   break;
	 }
    default :
	 { cout << "\nErreur : valeur incorrecte du type de lecture !\n";
	   cout << "ElemThermi::Lecture_bas_inf (ifstream& ent,const int cas)"
			<< " cas= " << cas << endl;
	   Sortie(1);
	  }
   }
 // ------ lecture dans tous les cas -------  	   
 // résultat d'erreur
 string toto;
 ent >> toto;
 if (toto == "erreur_de_densite_flux")
        { if (fluxErreur == NULL)
           fluxErreur = new double;
          ent  >> (*fluxErreur) ;
         } 
 // données particulière pour les lois de comportement mécanique
 int tabSaveDonTaille = tabSaveDon.Taille();
 if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) ent >> toto; 
 int num ; // numéro du pt d'integ, n'est pas vraiment utilisé mais c'est mieux que de lire un string
 for (int i=1; i<= tabSaveDonTaille; i++)
   if (tabSaveDon(i) != NULL) 
    { ent >> toto >> num; 
      tabSaveDon(i)->Lecture_base_info(ent,cas);}
 // données particulière pour les lois de comportement thermo physique
 int tabSaveTPTaille = tabSaveTP.Taille();
 if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) ent >> toto; 
 for (int i=1; i<= tabSaveTPTaille; i++)
   if (tabSaveTP(i) != NULL) 
    { ent >> toto >> num; 
      tabSaveTP(i)->Lecture_base_info(ent,cas);}
 // données particulière pour la déformation mécanique
 int tabSaveDefDonTaille = tabSaveDefDon.Taille();
 if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL))  ent >> toto;
 for (int i=1; i<= tabSaveDefDonTaille; i++)
    if (tabSaveDefDon(i) != NULL) 
    { ent >> toto >> num; 
      tabSaveDefDon(i)->Lecture_base_info(ent,cas);}
  // lecture des énergies
  int energ_taille = tab_energ_t.Taille();
  for (int i=1;i<= energ_taille;i++) ent >> tab_energ_t(i) ;
//  // énergie et puissance éventuelle de la partie bulk viscosity
//  ent >> toto >> E_elem_bulk_t >> toto >> P_elem_bulk;
//  E_elem_bulk_tdt = E_elem_bulk_t;
  };
       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void ElemThermi::Ecriture_bas_inf(ofstream& sort,const int cas)
 { // appel de la routine d'élément
   Element::Ecri_bas_inf_element(sort,cas);
   // en fait ici on sauvegarde la même chose dans tous les cas, par contre la sortie
   // totale est documentée.
   switch (cas)
  { case 1 : // ------- on sauvegarde tout -------------------------
     { // écriture de la masse volumique, 
       sort << "masse_volumique " << masse_volumique <<" " << "\n";
       // données thermique
       sort << "dilatation_thermique " << dilatation << " ";
       // blocage éventuelle d'hourglass
       sort << "\n hourglass: " << Nom_StabHourglass(type_stabHourglass) << " ";
	   break;
	 }
    case  2 : // ----------- sauvegarde uniquement de se qui varie --------------------
     { 	   break;
	 }
    default :
	 { cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
	   cout << "ElemThermi::Ecriture_bas_inf(ofstream& sort,const int cas)"
			<< " cas= " << cas << endl;
	   Sortie(1);
	  }
   }
  // --- informations à sortir dans tous les cas ----  	   
  // résultat d'erreur
  if (fluxErreur != NULL)
          sort << "erreur_de_densite_de_flux " << (*fluxErreur) <<" " << "\n";
  else
          sort << "pas_d'erreur_de_densite_de_flux \n";
  // données particulière pour les lois de comportement mécaniques
  int tabSaveDonTaille = tabSaveDon.Taille();
  if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) sort << " data_spec_loi_comp_meca ";
  for (int i=1; i<= tabSaveDonTaille; i++)
     if (tabSaveDon(i) != NULL) 
        { if (i==1) {sort << "pt_int= " <<i << " ";} else {sort << "\n                         pt_int= " <<i << " ";};
         tabSaveDon(i)->Ecriture_base_info(sort,cas);}
  // données particulière pour les lois de comportement thermo physiques
  int tabSaveTPTaille = tabSaveTP.Taille();
  if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) sort << " data_spec_loi_comp_ThemoPhysique ";
  for (int i=1; i<= tabSaveTPTaille; i++)
     if (tabSaveTP(i) != NULL) 
        { if (i==1) {sort << "pt_int= " <<i << " ";} else {sort << "\n                       pt_int= " <<i << " ";};
         tabSaveTP(i)->Ecriture_base_info(sort,cas);}
  // données particulière pour la déformation mécanique
  int tabSaveDefDonTaille = tabSaveDefDon.Taille();
  if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL)) sort << " data_spec_def ";
  for (int i=1; i<= tabSaveDefDonTaille; i++)
     if (tabSaveDefDon(i) != NULL) 
        { if (i==1) {sort << "pt_int= " <<i << " ";} else {sort << "\n               pt_int= " <<i << " ";};
          tabSaveDefDon(i)->Ecriture_base_info(sort,cas);}
  // sortie des énergies
  int energ_taille = tab_energ_t.Taille();
  sort << "\n ";
  for (int i=1;i<= energ_taille;i++) sort << tab_energ_t(i) << " " ;
//  // énergie et puissance éventuelle de la partie bulk viscosity
//  sort << "\n E_el_bulk= " << E_elem_bulk_t << " P_el_bulk= " << P_elem_bulk;
  };
 
// 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 ElemThermi::Interne_Long_arrete_mini_sur_c(Enum_dure temps,int nb_noeud)
  { int nbne = tab_noeud.Taille(); // récup du nombre de noeud
    if (nb_noeud != 0)
        nbne = nb_noeud;
    // tout d'abord l'objectif est de déterminer la distance minimum
    // entre les différents noeuds
    // initialisation de la distance entre les deux noeuds
    double dist = ConstMath::tresgrand;
    for (int i=1;i<= nbne;i++)
      // on itère sur les noeuds restants
      switch (temps)
      { case TEMPS_0:
          for (int j=i+1;j<= nbne;j++)
            { double dist_new = (tab_noeud(i)->Coord0() - tab_noeud(j)->Coord0()).Norme();
              if (dist_new < dist) dist = dist_new;
             }
          break;   
        case TEMPS_t:
          for (int j=i+1;j<= nbne;j++)
            { double dist_new = (tab_noeud(i)->Coord1() - tab_noeud(j)->Coord1()).Norme();
              if (dist_new < dist) dist = dist_new;
             }
          break;   
        case TEMPS_tdt:
          for (int j=i+1;j<= nbne;j++)
            { double dist_new = (tab_noeud(i)->Coord2() - tab_noeud(j)->Coord2()).Norme();
              if (dist_new < dist) dist = dist_new;
             }
          break;   
        default :
          cout << "\n cas du temps non implante temps= " << Nom_dure(temps)
               << "\n ElemThermi::Interne_Long_arrete_mini_sur_c(Enum_dure temps)";
          Sortie(1);                       
       }      
    // traitement d'une erreur éventuelle
    if (dist <= ConstMath::petit)
     { cout << "\n **** ERREUR une longueur d'arrete de l'element est nulle"
            << "\n ElemThermi::Interne_Long_arrete_mini_sur_c(..."
            << "\n element: "; this->Affiche(1);
       #ifdef MISE_AU_POINT
        cout << "\n *** version mise au point: on continue neanmoins avec une longueur "
             << " arbitrairement tres petite (" <<ConstMath::petit <<")  ";
      #else
        Sortie(1);
      #endif
     };
    // on calcul la célérité
    double cc = sqrt(ElemThermi::Calcul_maxi_E_qui(temps)/masse_volumique);
    // calcul du résultat
    double l_sur_c = dist/cc;
    return l_sur_c;
   };
   

            
// 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 ElemThermi::InitCalculMatriceMassePourRelaxationDynamique(int )
{ // on active le calcul des invariants de contraintes
  int taille = lesPtIntegThermiInterne->NbPti();
 
  // on définit  le ddl étendu correspondant, s'il n'existe pas 
  // on définit le  Ddl_enum_etendu correspondant à la masse au noeud pour la méthode
  if (masse_relax_dyn.Nom_vide() )
  	  masse_relax_dyn = (Ddl_enum_etendu("masse_relax_dyn"));        		
};

// phase de calcul de la matrice masse dans le cas de l'algo de relaxation dynamique
// mi= lambda * delta_t**2 / 2 * (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises)
// avec delta_t=1 par défaut a priori, donc n'intervient pas ici
// 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  ElemThermi::CalculMatriceMassePourRelaxationDynamique
                      (const double& alph, const double& beta, const double & lambda
                      ,const double & gamma,const double & theta, int casMass_relax)
 {  // 1)
    // on calcule la moyenne des grandeurs qui servent pour le calcul de la pseudo-masse que l'on distribue 
    // également sur les noeuds
    double trace_moyenne= 0.; double  mises_moyen= 0.; 
    double compress_moy=0.; double cisaille_moy = 0.;
    int taille = lesPtIntegThermiInterne->NbPti();
  
cout << "\n stop *** pour l'instant la methode : ElemThermi::CalculMatriceMassePourRelaxationDynamique"
     << " n'est pas encore operationnelle " << endl;
Sortie(1);

//	 // on boucle sur les points d'intégrations
//    for (int i= 1; i<= taille; i++)
//        {PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(i);
//    	 const Vecteur & invariant = ptIntegThermi.SigInvar();
//    	 trace_moyenne += invariant(1) ;  // le premier invariant c'est la trace
//    	 if (invariant.Taille() > 1)
//    	    mises_moyen += sqrt((3.*invariant(2)-invariant(1)*invariant(1))/2.);
//    	 compress_moy +=  ptIntegThermi.ModuleCompressibilite();
//    	 cisaille_moy +=  ptIntegThermi.ModuleCisaillement();
//    	};
  
    // --- tout d'abord la partie contrainte
//    double ensemble = (val_propre_sup + val_propre_inf + val_cisaillement) / taille;
    // on divise par taille, ce qui conduit à la moyenne sur les points d'intégration	
    double ensemble = (gamma*trace_moyenne /3.+ 0.5*theta*mises_moyen ) / taille;	
    // --- maintenant la partie raideur
    ensemble += (alph * compress_moy + beta * cisaille_moy)  / taille;
    // --- pour avoir l'intégration totale on utilise le volume
//    ensemble *= volume;
    // --- et maintenant on distribue sur tous les noeuds de manière identique
    // ne fonctionne que si tous les noeuds sont actifs !!! (ce qui est un cas particulier courant)
    // la masse élémentaire pour chaque noeud  
//	double masse_elementaire = ensemble * lambda * 0.5 / tab_noeud.Taille();
   // pour des membrannes on a: 
	// mi= lambda * delta_t**2 / 2 * epaisseur/4 * (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises)
	// avec delta_t=1 par défaut a priori, donc n'intervient pas ici
	// d'où le facteur 0.125, que l'on garde également pour le volume
	double masse_elementaire = ensemble * lambda * 0.125 ;
	
	// prise en compte du type d'élément: 
	double ParNoeud = 1.; // n'est utilisé que par la méthode de Barnes classique
	int nbn=tab_noeud.Taille();
	switch (Type_geom_generique(ElementGeometrique_const().TypeGeometrie()))
	{  case VOLUME : // par rapport à la surface, on utilise la racine cubique du volume
	     	{  double long_caracteristique = pow(Volume(),1./3.); 
	     	   masse_elementaire *= long_caracteristique; 
				// ParNoeud = volume / nombre de noeuds, donc c'est la partie du volume 
				// attribuée à chaque noeud
	     	   ParNoeud = Volume() / nbn;
	     	   break;
	     	  }
	   case SURFACE : // l'algorithme historique
	     	{  double epaisseur_moyenne = EpaisseurMoyenne(TEMPS_tdt);
	     	   masse_elementaire *= epaisseur_moyenne; 
				// ParNoeud = section / nombre de noeuds, donc c'est la partie de la section 
				// attribuée à chaque noeud
	     	   ParNoeud = Volume() / epaisseur_moyenne/ nbn;
	     	   break;
	     	  }
	   case LIGNE : 
	     	{  double section_moyenne = SectionMoyenne(TEMPS_tdt);
			   double long_caracteristique = Volume()/section_moyenne;
	     	   masse_elementaire *= long_caracteristique; 
				// ParNoeud = longueur / nombre de noeuds, donc c'est la partie de la longueur 
				// attribuée à chaque noeud
	     	   ParNoeud = Volume() / section_moyenne/ nbn;
	     	   break;
	     	  }
		default :
			cout << "\nErreur : pour l'instant les types autres que volume, surface, ligne ne sont pas pris en compte dans la relaxation "
			     << " dynamique selon barnes !\n";
			cout << "\n ElemThermi::CalculMatriceMassePourRelaxationDynamique(.. \n";
			Sortie(1);
	};	
	// on parcours les noeuds de l'élément et on remplace le cas échéent, la valeur de la masse au noeud.
    switch (casMass_relax)
	 { case 1: case 3: // cas 1 : on cumule aux noeuds, cas 3 on fait la moyenne (pas effectué ici, on ne fait 
	                   //         que préparer le travail pour le cas 3)
	    {for (int inoe=1;inoe<=nbn;inoe++)
		  {Noeud& noe = *tab_noeud(inoe);
		   // on cumule
		   noe.ModifDdl_etendu(masse_relax_dyn).Valeur() += masse_elementaire;
		   };
		 break;
	    }
	   case 2: // cas 2 : on prend la masse maxi
	    {for (int inoe=1;inoe<=nbn;inoe++)
		  {Noeud& noe = *tab_noeud(inoe);
		   const Ddl_etendu& masse_actuel =noe.DdlEtendue(masse_relax_dyn);
		   // dans le cas où la masse actuelle est plus petite on la remplace
		   if (masse_actuel.ConstValeur() <= masse_elementaire)
		      noe.ModifDdl_etendu(masse_relax_dyn).Valeur() = masse_elementaire;
		   };
		 break;
	    }
	   case 4: case 5: //  on cumule (puis on moyennera pour le cas 4, autre part) et on divise par la section comme dans la formule de barnes
	    {for (int inoe=1;inoe<=nbn;inoe++)
		  {Noeud& noe = *tab_noeud(inoe);
		   const Ddl_etendu& masse_actuel =noe.DdlEtendue(masse_relax_dyn);
		   // on cumule
		   if (masse_actuel.ConstValeur() <= masse_elementaire)
			   // !!!! je pense que ce cas est très bizarre il faudra revoir ce que cela signifie....  !!
		      noe.ModifDdl_etendu(masse_relax_dyn).Valeur() += (masse_elementaire/ParNoeud);
		   };
		 break;
	    }
	 }; //-- fin du switch
	};  
            
    // METHODES VIRTUELLES: 
        

// recuperation des coordonnées du point de numéro d'ordre = iteg pour 
// la grandeur enu
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
// utilisé par les classes dérivées
Coordonnee ElemThermi::CoordPtInt(Enum_dure temps,Enum_ddl enu,int iteg,bool& err)
  { Coordonnee ptret;err = false;
    // récupération de l'élément géométrique correspondant à Enu
    ElemGeomC0& ele = this->ElementGeometrie(enu);
    // récupération du tableau des fonctions d'interpolations
    const Tableau <Vecteur> & tabphi =   ele.TaPhi();
    if ((iteg < 0) || (iteg > tabphi.Taille()))
     { err = true;}
    else
     { switch (temps)
       { case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
         case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
         case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
        }
     };
    return ptret; // retour
   }; 
                
                
 // 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 ElemThermi::CoordPtIntFace(int face, Enum_dure temps,int iteg,bool& err)
 { Coordonnee ptret;err = false;
  if (SurfExiste(face))
   {// récupération de l'élément géométrique correspondant à la face
    const ElFrontiere* elf = this->Frontiere_surfacique(face,false);
    if (elf != NULL) // cas où la frontière existe déjà
     { const ElemGeomC0 & elegeom = elf->ElementGeometrique();
       // récupération du tableau des fonctions d'interpolations
       const Tableau <Vecteur> & tabphi =   elegeom.TaPhi();
       if ((iteg < 0) || (iteg > tabphi.Taille()))
        { err = true;}
       else
        { switch (temps)
          { case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
            case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
            case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
          };
        };
     }
   }
  else // cas où la frontière n'existe pas, on renvoie une erreur
       // et on laisse la valeur par défaut pour ptret
   {err = true;
   };
  return ptret; // retour
 };

 // 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 ElemThermi::CoordPtIntArete(int arete, Enum_dure temps,int iteg,bool& err)
 { Coordonnee ptret;err = false;
   if (AreteExiste(arete))
    {// récupération de l'élément géométrique correspondant à l'arête
     const ElFrontiere* elf = this->Frontiere_lineique(arete,false);
     if (elf != NULL) // cas où la frontière existe déjà
      { const ElemGeomC0 & elegeom = elf->ElementGeometrique();
        // récupération du tableau des fonctions d'interpolations
        const Tableau <Vecteur> & tabphi =   elegeom.TaPhi();
        if ((iteg < 0) || (iteg > tabphi.Taille()))
         { err = true;}
        else
         { switch (temps)
           { case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
             case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
             case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
           };
         };
      }
   }
  else // cas où la frontière n'existe pas, on renvoie une erreur
       // et on laisse la valeur par défaut pour ptret
   {err = true;
   };
  return ptret; // retour
 };

// 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
// utilisé par les classes dérivées
int ElemThermi::PtLePlusPres(Enum_dure temps,Enum_ddl enu, const Coordonnee& M)
  { int iret;
    // récupération de l'élément géométrique correspondant à Enu
	ElemGeomC0& ele = this->ElementGeometrie(enu);
	// récupération du tableau des fonctions d'interpolations
	const Tableau <Vecteur> & tabphi =   ele.TaPhi();
	// on boucle sur les pt d'integ pour avoir le point le plus près
	int tabphitaille = tabphi.Taille();
	Coordonnee P; iret=1; double dist= ConstMath::tresgrand;
	for (int ipt = 1;ipt<=tabphitaille;ipt++)
	 { switch (temps)
	   { case TEMPS_0 : P = met->PointM_0(tab_noeud,tabphi(ipt)); break;
	     case TEMPS_t : P = met->PointM_t(tab_noeud,tabphi(ipt)); break;
	     case TEMPS_tdt : P = met->PointM_tdt(tab_noeud,tabphi(ipt)); break;
	    }
	   double di=(M-P).Norme();
	   if (di < dist) { dist = di; iret = ipt;}; 
     }
    return iret; // retour
   }; 
	     
// ==== >>>> methodes virtuelles redéfini éventuellement dans les classes dérivées  ============  

// ramene l'element geometrique correspondant au ddl passé en paramètre
ElemGeomC0& ElemThermi::ElementGeometrie(Enum_ddl  ddl) const
   { Enum_ddl enu = PremierDdlFamille(ddl);
     switch (enu)
     { case X1 : return ElementGeometrique(); break;
       case TEMP : return ElementGeometrique(); break;
       case FLUXD1 : return ElementGeometrique(); break;
       case GRADT1 : return ElementGeometrique(); break;
       case DGRADT1 : return ElementGeometrique(); break;
       case ERREUR : return ElementGeometrique(); break;
       default :
        {cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl) 
              << "\n ElemThermi::ElementGeometrie(Enum_ddl  ddl) " ;
         Sortie(1);}
      }        
     return ElementGeometrique(); // pour taire le compilo 
    }; 

// 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 
int ElemThermi::NbGrandeurGene(Enum_ddl ddl) const
   { Enum_ddl enu = PremierDdlFamille(ddl);
     int nbGG = 0.; // par défaut
     switch (enu)
     { case TEMP :  case FLUXD1 :case GRADT1 :case DGRADT1 :
         { // cas d'un calcul de thermique classique
           switch (Type_geom_generique(id_geom))
	            {   case LIGNE : case POINT_G :   nbGG = 1; break; // FLUXD1
		               case SURFACE : nbGG = 2; break; // FLUXD1, FLUXD2
		               case VOLUME :  nbGG = 3; break; // FLUXD1 FLUXD2 FLUXD3
		               default :
                  cout << "\nErreur : cas non traite, id_geom= :" << Nom_type_geom(Type_geom_generique(id_geom))
                       << "ElemThermi::NbGrandeurGene(.. \n";
                  Sortie(1);
	         };
        	 break;
        	} 
       default :
        {cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl) 
              << "\n ElemThermi::NbGrandeurGene(Enum_ddl  ddl) " ;
         Sortie(1);}
      };        
     return nbGG;  
    }; 
        
// 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 ElemThermi::Modif_orient_elem(int cas_orientation)
 { // retour: 
   bool retour=false; // par défaut pas de changement
   if (cas_orientation == 0) // cas où on inverse l'orientation sans condition particulière
    { // on change l'orientation de l'élément
      retour = true;
      int nbnoe = tab_noeud.Taille();
      Tableau<Noeud *> tab_inter(tab_noeud); // on crée un tableau intermédiaire
      // on récupère la numérotation locale inverse
      const Tableau<int>  tabi =  ElementGeometrique().InvConnec();
		// on met à jour le tableau actuel
      for (int n=1;n<=nbnoe;n++)
      	 tab_noeud(n)=tab_inter(tabi(n));
	 }
   else 
	 { // 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
      int cas=1; // a priori tout est ok
      // boucle sur les pt d'integ
      for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg())
        { double jacobien_0 = def->JacobienInitial();
          
          if (jacobien_0 < 0)
          	{if (cas ==1) // on a trouvé un jacobien négatif  
          	  {cas =0; }
          	}// si c'était positif --> négatif
          else // cas positif
           	{if (cas == 0) // on a déjà changé
           	  { cas = 2; break;} // on sort de la boucle
          	};
        };
      // gestion de pb
      if (cas == 2)
      	{ cout << "\n **** Attention **** element nb= "<< this->Num_elt() << " peut-etre trop distordu ?" 
      	       << " pt d'integ Thermi positif et  negatif !! ";
          if (ParaGlob::NiveauImpression() > 7)
          {// on va essayer de sortir plus d'info
           int pti=1;cout << "\n les jacobiens initiaux \n";
           // on sort les valeurs des jacobiens
           for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg(),pti++)
            { cout << " pti= "<<pti << ", J= "<< def->JacobienInitial();};
           // les coordonnées des noeuds en suivant l'ordre de la table de connection locale
           int nbnoe = tab_noeud.Taille();
           cout << "\n les coordonnees de l'element \n";
           for (int ino=1;ino<=nbnoe;ino++)
            { cout << "\n "<< ino << " " << tab_noeud(ino)->Coord0();
            };
          
          };
      	}
      else if (cas == 0)
      	{ switch (cas_orientation)
           { case 1: // on change l'orientation de l'élément
      	     {retour = true;
      	      int nbnoe = tab_noeud.Taille();
      	      Tableau<Noeud *> tab_inter(tab_noeud); // on crée un tableau intermédiaire
      	      // on récupère la numérotation locale inverse
              const Tableau<int>  tabi =  ElementGeometrique().InvConnec();
      	      for (int n=1;n<=nbnoe;n++)
      	       tab_noeud(n)=tab_inter(tabi(n));
      	      break;
      	     }
           	 case -1: // on sort une information à l'écran
      	      { cout << "\n element nb= "<< this->Num_elt() << "jacobien negatif " ;
      	        break;;
      	      }
             default:
              cout << "\n erreur le cas : " << cas_orientation
                   << " n'est pas actuellement pris en compte"
                   << "\n ElemThermi::Modif_orient_elem(...";
              Sortie(1);     
           };
      	};
	 };	
	// retour
	return retour;	
 };