// 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 <cmath>
#include "ElemThermi.h"
#include <iomanip>
#include "ConstMath.h"
#include "Util.h"
#include "Coordonnee2.h"
#include "Coordonnee3.h" 
#include "ParaAlgoControle.h"
#include "ExceptionsElemThermi.h"
#include "Loi_Umat.h"

// =========================== variables communes ===============================
// a priori les pointeurs de fonction pointent sur rien au debut
const Coordonnee & (Met_abstraite::*(ElemThermi::PointM))
           (const  Tableau<Noeud *>& tab_noeud,const  Vecteur& Phi) ; 
void  (Met_abstraite::*(ElemThermi::BaseND))
             (const  Tableau<Noeud *>& tab_noeud,const  Mat_pleine& dphi,
                    const Vecteur& phi,BaseB& bB,BaseH& bH); 
int ElemThermi::bulk_viscosity = 0; // indique si oui ou non on utilise le bulk viscosity
double ElemThermi::c_traceBulk = 0.; // coeff de la trace de D pour le bulk
double ElemThermi::c_carreBulk = 0.; // coeff du carré de la trace de D pour le bulk
TenseurHH* ElemThermi::sig_bulkHH = NULL;  // variable de travail pour le bulk
// ---------- pour le calcul de la masse pour la méthode de relaxation dynamique
Ddl_enum_etendu ElemThermi::masse_relax_dyn; // au début rien, def à la première utilisation
                    
// =========================== fin variables communes ===========================

ElemThermi::ElemThermi () :
Element(),tabSaveDon(),tabSaveTP(),tabSaveDefDon(),defSurf(),defArete()
,masse_volumique(masse_volumique_defaut),dilatation(false),tab_energ(),tab_energ_t(),energie_totale()
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,premier_calcul_thermi_impli_expli(true),lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
// Constructeur par defaut
{	loiComp = NULL;loiTP=NULL;loiFrot=NULL;
  met = NULL;
  def = NULL;
  defEr = NULL;
  defMas = NULL;
  fluxErreur = NULL;
  id_problem = THERMIQUE;
};

ElemThermi::ElemThermi (int num_mail,int num_id) :
Element(num_mail,num_id),tabSaveDon(),tabSaveTP(),tabSaveDefDon(),defSurf(),defArete()
,masse_volumique(masse_volumique_defaut),dilatation(false),tab_energ(),tab_energ_t(),energie_totale()
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,premier_calcul_thermi_impli_expli(true)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
// Constructeur utile quand le numero d'identification de l'element est connu
{	loiComp = NULL;loiTP=NULL;loiFrot=NULL;
  met = NULL;
  def = NULL;
  defEr = NULL;
  defMas = NULL;
  fluxErreur = NULL;
  id_problem = THERMIQUE;
};

ElemThermi::ElemThermi (int num_mail,int num_id,const Tableau<Noeud *>& tab):
// Constructeur utile quand le numero et le tableau des noeuds 
// de l'element sont connu
Element(num_mail,num_id,tab)
,tabSaveDon(),tabSaveTP(),tabSaveDefDon(),defSurf(),defArete()
,masse_volumique(masse_volumique_defaut),dilatation(false),tab_energ(),tab_energ_t(),energie_totale()
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,premier_calcul_thermi_impli_expli(true)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
  met = NULL;
  def = NULL;
  defEr = NULL;
  defMas = NULL;
  fluxErreur = NULL;
  id_problem = THERMIQUE;
};
	

ElemThermi::ElemThermi (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,string info):
// Constructeur utile quand le numero d'identification est connu,
// ainsi que la geometrie et le type d'interpolation de l'element
Element (num_mail,num_id,id_interp_elt,id_geom_elt,THERMIQUE,info)
,tabSaveDon(),tabSaveTP(),tabSaveDefDon()
,defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(false)
,tab_energ(),tab_energ_t(),energie_totale()
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,premier_calcul_thermi_impli_expli(true)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
   met = NULL;
   def = NULL;
   defEr = NULL;
   defMas = NULL;
   fluxErreur = NULL;
  id_problem = THERMIQUE;
};

ElemThermi::ElemThermi (int num_mail,int num_id,char* nom_interpol,char* nom_geom,string info):
// Constructeur utile quand le numero d'identification est connu,
// ainsi que la geometrie et le type d'interpolation de l'element
Element(num_mail,num_id,nom_interpol,nom_geom,NomElemTypeProblem(THERMIQUE),info)
,tabSaveDon()
,tabSaveTP(),tabSaveDefDon(),defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(false)
,tab_energ(),tab_energ_t(),energie_totale(),premier_calcul_thermi_impli_expli(true)
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
{  loiComp = NULL;loiTP=NULL;loiFrot=NULL;
    met = NULL;
    def = NULL;
    defEr = NULL;
    defMas = NULL;
   fluxErreur = NULL;
  id_problem = THERMIQUE;
};

ElemThermi::ElemThermi (int num_mail,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 sont connues
Element (num_mail,num_id,tab,id_interp_elt,id_geom_elt,THERMIQUE,info)
,tabSaveDon(),tabSaveTP()
,tabSaveDefDon(),defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(false)
,tab_energ(),tab_energ_t(),energie_totale(),premier_calcul_thermi_impli_expli(true)
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
{	loiComp = NULL;loiTP=NULL;loiFrot=NULL;
    met = NULL;
    def = NULL;	
    defEr = NULL;
    defMas = NULL;
    fluxErreur = NULL;
  id_problem = THERMIQUE;
};

ElemThermi::ElemThermi (int num_mail,int num_id,const Tableau<Noeud *>& tab,char* nom_interpol,
												   char* nom_geom,string info):
// Constructeur utile quand toutes les donnees sont connues
Element (num_mail,num_id,tab,nom_interpol,nom_geom,NomElemTypeProblem(THERMIQUE),info)
,tabSaveDon()
,tabSaveTP(),tabSaveDefDon(),defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(false)
,tab_energ(),tab_energ_t(),energie_totale(),premier_calcul_thermi_impli_expli(true)
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
{	loiComp = NULL;loiTP=NULL;loiFrot=NULL;
    met = NULL;
    def = NULL;
    defEr = NULL;
    defMas = NULL;
    fluxErreur = NULL;
  id_problem = THERMIQUE;
};

ElemThermi::ElemThermi (const   ElemThermi& elt):
// Constructeur de copie
Element (elt),defSurf(elt.defSurf),defArete(elt.defArete)
,masse_volumique(elt.masse_volumique),dilatation(elt.dilatation)
,tab_energ(elt.tab_energ),tab_energ_t(elt.tab_energ_t),energie_totale(elt.energie_totale)
,E_elem_bulk_t(elt.E_elem_bulk_t),E_elem_bulk_tdt(elt.E_elem_bulk_tdt),P_elem_bulk(elt.P_elem_bulk)
,premier_calcul_thermi_impli_expli(true)
,lesPtIntegThermiInterne(NULL)
,type_stabHourglass(elt.type_stabHourglass),tab_elHourglass()
,coefStabHourglass(elt.coefStabHourglass),raid_hourglass_transitoire(NULL)
{ // pour la loi on pointe sur la même que l'élément de référence, car de toute
  // la loi ne contient pas les données spécifique
    loiComp = (elt.loiComp);loiTP = (elt.loiTP);loiFrot = (elt.loiFrot);
  // idem pour la métrique qui est a priori une métrique générique  
    met = (elt.met);
  id_problem = THERMIQUE;
  // la déformation est spécifique à l'élément
    if (elt.def!=NULL) def = elt.def->Nevez_deformation(tab_noeud);  // *def = *(elt.def);
    if (elt.defEr!=NULL) defEr = elt.defEr->Nevez_deformation(tab_noeud); //*defEr = *(elt.defEr);
    if (elt.defMas!=NULL) defMas = elt.defMas->Nevez_deformation(tab_noeud); //*defMas = *(elt.defMas);
  // les données spécifiques à la déformation mécanique
    int idefdtr = (elt.tabSaveDefDon).Taille();
    tabSaveDefDon.Change_taille(idefdtr);
    for (int i=1;i<= idefdtr; i++)
     { tabSaveDefDon(i) = def->New_et_Initialise(); 
       if (tabSaveDefDon(i)!= NULL) *tabSaveDefDon(i) = *(elt.tabSaveDefDon)(i);
       }
  // def: cas des arrêtes et des faces si besoin
    int tailledefSurf = defSurf.Taille();
    int itab = 1; // indice pour parcourir tabb
    for (int idsu =1;idsu<=tailledefSurf;idsu++,itab++)
      {if (defSurf(idsu) != NULL) 
          defSurf(idsu) = elt.defSurf(idsu)->Nevez_deformation(tabb(itab)->TabNoeud());}
    int tailledefArete = defArete.Taille();
    for (int idA =1;idA<=tailledefArete;idA++)
      {if (defArete(idA) != NULL) 
          defArete(idA) = elt.defArete(idA)->Nevez_deformation(tabb(itab)->TabNoeud());}
  // les donnees de la loi de comportement mécanique 
    int dtr = (elt.tabSaveDon).Taille();
    tabSaveDon.Change_taille(dtr);
    for (int i=1;i<= dtr; i++)
     { tabSaveDon(i) = loiComp->New_et_Initialise(); 
       if (tabSaveDon(i)!= NULL) *tabSaveDon(i) = *(elt.tabSaveDon)(i);
       }
  // les donnees de la loi de comportement thermo physique 
    int dtrTP = (elt.tabSaveTP).Taille();
    tabSaveTP.Change_taille(dtrTP);
    for (int i=1;i<= dtrTP; i++)
     { tabSaveTP(i) = NULL; // init
       if (loiTP != NULL ) tabSaveTP(i) = loiTP->New_et_Initialise(); 
       if (tabSaveTP(i)!= NULL) *tabSaveTP(i) = *(elt.tabSaveTP)(i);
       }
  // et autre     
    if (elt.fluxErreur != NULL)
      { fluxErreur = new double;
        fluxErreur = elt.fluxErreur;
       }
    else
      fluxErreur = NULL;
		
};

ElemThermi::~ElemThermi ()
// Destructeur
{ // les donnees pour la loi de comportement mécanique
  int imaxi = tabSaveDon.Taille();
  for (int i=1;i<=imaxi; i++)
       if (tabSaveDon(i) != NULL) {delete tabSaveDon(i) ; tabSaveDon(i)=NULL;}
  // les donnees pour la loi de comportement thermo physique
  int imaxTP = tabSaveTP.Taille();
  for (int i=1;i<=imaxTP; i++)
       if (tabSaveTP(i) != NULL) {delete tabSaveTP(i) ; tabSaveTP(i)=NULL;}
  // les donnees à la déformation mécanique
  int idefmaxi = tabSaveDefDon.Taille();
  for (int i=1;i<=idefmaxi; i++)
       if (tabSaveDefDon(i) != NULL) {delete tabSaveDefDon(i) ; tabSaveDefDon(i)=NULL;}
  // les déformations sont spécifiquent aux elements, on les detruits si necessaire
  if (def != NULL) { delete def; def = NULL;}
  if (defEr != NULL) {delete defEr; defEr = NULL;}
  if (defMas != NULL) {delete defMas; defMas = NULL;}
  int nbdefSurf = defSurf.Taille();
  for (int i = 1; i<= nbdefSurf; i++) 
     if (defSurf(i) != NULL) {delete defSurf(i); defSurf(i) = NULL;} 
  int nbdefArete = defArete.Taille();
  for (int i = 1; i<= nbdefArete; i++) 
     if (defArete(i) != NULL) {delete defArete(i);defArete(i) = NULL;}  
  // par contre la metrique et la loi de comportement sont communes a une classe d'element
  // donc ici met et loicomp ne sont que des pointeurs dont il ne faut pas supprimer les objets
  // pointes 
  if (fluxErreur!= NULL) {delete fluxErreur; fluxErreur = NULL;}
  // cas des éléments de stabilisation d'hourglass
  if (tab_elHourglass.Taille() != 0)
   { int taillhour = tab_elHourglass.Taille();
	  for (int i=1;i<= taillhour;i++)
	    { if (tab_elHourglass(i) != NULL) delete tab_elHourglass(i);};
	 };
  if (raid_hourglass_transitoire != NULL)	delete raid_hourglass_transitoire;
};

// =========================== methodes publiques ===============================

// 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
// dans le cas où l'on veut un fonctionnement différent, c'est surchargé dans les 
// éléments fils
int ElemThermi::Interne_0(const Coordonnee& M,Coordonnee* coor_locales)
  { PointM = &Met_abstraite::PointM_0;
    BaseND = &Met_abstraite::BaseND_0;
    return Interne(TEMPS_0,M,coor_locales);    
   }; 
int ElemThermi::Interne_t(const Coordonnee& M,Coordonnee* coor_locales)
  { PointM = &Met_abstraite::PointM_t;
    BaseND = &Met_abstraite::BaseND_t;
    return Interne(TEMPS_t,M,coor_locales);    
   }; 
int ElemThermi::Interne_tdt(const Coordonnee& M,Coordonnee* coor_locales)
  { PointM = &Met_abstraite::PointM_tdt;
    BaseND = &Met_abstraite::BaseND_tdt;
    return Interne(TEMPS_tdt,M,coor_locales);    
   }; 

// test si l'element est complet
int ElemThermi::TestComplet()
    { int ret = Element::TestComplet();  // teste de la class mere
      if (loiTP == NULL)
       { cout << " \n la loi thermo physique n'est pas defini dans l element " << num_elt << '\n';
         ret = 0;
        }
      else if (loiTP->TestComplet() != 1)
       { cout << " \n la loi thermo physique n'est pas complete dans l element " << num_elt << '\n';
         ret = 0;
        }
      // dans le cas d'une loi  thermo physique avec une partie mécanique on vérifie qu'elle est complète
      if (loiComp != NULL)
       { if (loiComp->TestComplet() != 1)
          { cout << " \n la loi thermo physique n'est pas complete dans l element " << num_elt << '\n';
           ret = 0;
          }
         };
      // dans le cas d'une loi  de frottement on vérifie qu'elle est complète
      if (loiFrot != NULL)
       { if (loiFrot->TestComplet() != 1)
          { cout << " \n la loi de frottement n'est pas complete dans l element " << num_elt << '\n';
           ret = 0;
          }
         };
      // suite des datas    
      if ( masse_volumique == masse_volumique_defaut)
        { cout << "\n la masse volumique n'est pas defini dans l element " << num_elt << '\n';
          ret = 0; 
        }
      if (met == NULL) 
       { cout << " \n la metrique n'est pas defini dans l element " << num_elt << '\n';
         ret = 0;
       }  
      if (def == NULL) 
       { cout << " \n le type de deformation n'est pas defini dans l element "
                   << num_elt << '\n';
         ret = 0;
       }
      return ret; 
    };
	    
// test pour savoir si le calcul de contrainte en absolu est possible
bool ElemThermi::FluxAbsoluePossible()
  { // il faut qu'il soit complet et il faut que les coordonnées a t =0 et t existe
    bool retour = true;
    if (TestComplet() !=0)
     {// on regarde la définition des coordonnées
      int nbmaxnoeud=tab_noeud.Taille();
      for (int i=1;i<=nbmaxnoeud;i++)
        if((!tab_noeud(i)->ExisteCoord0())||(!tab_noeud(i)->ExisteCoord1()))
          { retour = false; break;}
      }    
    else
      retour = false;
    return retour;       
   };  

    
      // --------- 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
ElemThermi::MatGeomInit ElemThermi::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & )
  { cout << "\n attention le calcul de la matrice g\'eom\'etrique  et de la matrice " 
         << " initiale   ne sont pas implant\'ees \n" ;
    this->Affiche(1);
    cout << "\nvoid ElemThermi::MatricesGeometrique_Et_Initiale (Mat_pleine &,Mat_pleine &)"
         << endl;
    Sortie(1);
    return MatGeomInit(NULL,NULL); // pour supprimer le warning
   };

// 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 flux).(delta flux) dv)/(int flux*flux 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 !!!!!!!!!!!!!!!!!!!!!!!!!!! donc le code est à virer
void ElemThermi::ErreurElement(int ,double& ,double& , double& )
  { cout << "\n attention le calcul de l'erreur sur l'element n'est pas implant\'e  " 
         << " pour cet element \n" ;
    this->Affiche(1);
    cout << "\nElemThermi::ErreurElement(int....."
         << endl;
    Sortie(1);
   };

// ajout des ddl d'erreur pour les noeuds de l'élément
void ElemThermi::Plus_ddl_Erreur()
  {int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) // pour chaque noeud
     { Ddl ddl(ERREUR,0.,LIBRE);
      tab_noeud(ne)->PlusDdl(ddl);
      }
   };
// inactive les ddl d'erreur 
void ElemThermi::Inactive_ddl_Erreur()
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_hors_service(ERREUR);
   };
// active les ddl d'erreur
void ElemThermi::Active_ddl_Erreur()
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(ERREUR);
   };
// 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
DdlElement  ElemThermi::Tableau_de_ERREUR() const
  { return DdlElement(tab_noeud.Taille(), DdlNoeudElement(ERREUR) );
    };
    
// =========================== methodes protegees ===============================


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

int ElemThermi::Interne(Enum_dure temps,const Coordonnee& M,Coordonnee* coor_locales)
  { int dim_glob = M.Dimension();  // dim de M = dim de l'espace géométrique
    ElemGeomC0& elgeom = this->ElementGeometrique(); // la geometrie 
    int dim_loc = elgeom.Dimension(); // dim de l'élément de référence
    #ifdef MISE_AU_POINT	 	 
      if (ParaGlob::Dimension() != dim_glob) 
        { cout << "\n *** erreur, la dimension du point = (" << dim_glob
               << ")  et de l'espace geometrique (" << ParaGlob::Dimension() << ") sont differentes! ";
          cout << "\nElemThermi::Interne(Coordonnee& M) " << endl;
          this->Affiche();
          Sortie(1);     
         }
    #endif
    // la methode consiste tout d'abord a calculer les coordonnees locales correspondants
    // a M. Pour cela on utilise un processus iteratif
    // ---- 1 :init de la recherche
    const Coordonnee refL(dim_loc) ; // coordonnees locales du point de reference, initialisees a zero
    // - construction du point local en fonction de l'existance ou non de coor_locales
    Coordonnee* pt_ret=NULL; Coordonnee theta_local(refL);
    if (coor_locales != NULL) {pt_ret = coor_locales;} else {pt_ret = &theta_local;};
    Coordonnee& theta = *pt_ret; // init de AL qui representera les coordonnes locales de M 
    // - fin construction du point local 
    Vecteur ph_i = elgeom.Phi_point(refL); // fonctions d'interpolation en refL
    Mat_pleine dph_i = elgeom.Dphi_point(refL); // derivees des fonctions d'interpolation en refL
    // l'image de theta par l'interpolation, de dimension : espace géométrique
    Coordonnee A = (met->*PointM) (tab_noeud,ph_i); //
	// calcul de la base naturelle et duale en refL, en fonction des coord a t 
	BaseB bB(dim_glob,dim_loc);BaseH bH(dim_glob,dim_loc);
	(met->*BaseND)(tab_noeud,dph_i,ph_i,bB,bH);
	// calcul des coordonnees locales du point M dans le repere en refl
	Coordonnee AM = M - A; // en global
	for (int i=1;i<=dim_loc;i++)
	  theta(i) = AM * bH(i).Coor(); // ok même si dim_loc < dim_glob -> projection sur le plan tangent
	Coordonnee theta_repere(theta); // le point où l'on calcul le repère 
	Coordonnee theta_initial(theta); Coordonnee delta_theta(theta); 
	// dans le cas où le point est externe à l'élément, on limite le repère de calcul au point externe
	// de l'élément dans la direction de theta
    if (!(elgeom.Interieur(theta)))
     { // calcul d'un point  extreme de l'élément dans le sens de M
       theta_repere = elgeom.Maxi_Coor_dans_directionGM(theta);
      } 
	// calcul du point correspondant à theta dans l'element, qui remplace A
    ph_i = elgeom.Phi_point(theta_repere); // fonctions d'interpolation en A
    Coordonnee A_sauve= A;
	A = (met->*PointM) (tab_noeud,ph_i);
    AM = M-A; // nouveau gap
    
    // ---- 2 :iteration
    // on boucle avec un test de stabilité des coordonnées locales
    int compteur = 1; // contrôle du maxi d'itération
    int nb_exterieur=0;  // contrôle du maxi de fois que l'on trouve consécutivement le point à l'extérieur
    // - récup des paramètres de controle
    double delta_thetai_maxi = ParaGlob::param->ParaAlgoControleActifs().PointInterneDeltaThetaiMaxi();
	 int nb_boucle_sur_delta_thetai= ParaGlob::param->ParaAlgoControleActifs().PointInterneNbBoucleSurDeltaThetai();
    int nb_max_a_l_exterieur= ParaGlob::param->ParaAlgoControleActifs().PointInterneNbExterne();
    double prec_tehetai_interne= ParaGlob::param->ParaAlgoControleActifs().PointInternePrecThetaiInterne();
    // - boucle de recherche des thetai
	 while (delta_theta.Norme() >= delta_thetai_maxi)
	        
     {dph_i = elgeom.Dphi_point(theta_repere); // derivees des fonctions d'interpolation en theta_repere
      ph_i = elgeom.Phi_point(theta_repere); // fonctions d'interpolation en theta_repere	(déjà calculé la première fois)
	  // calcul de la base naturelle et duale  en fonction des coord a t 
	  (met->*BaseND)(tab_noeud,dph_i,ph_i,bB,bH);
	  // amelioration  des coordonnees locales du point M dans le repere en theta
	  for (int i=1;i<=dim_loc;i++)
	     delta_theta(i)=AM * bH(i).Coor();
//	  theta += delta_theta;  erreur je pense: le 11 juin 2007
      // le nouveau theta = la position du repère + le deltat theta
	  theta = theta_repere + delta_theta;
	   
      theta_repere=theta;
      if (!(elgeom.Interieur(theta))) // si le point est à l'extérieur
       { // calcul d'un point  extreme de l'élément dans le sens de M
         theta_repere = elgeom.Maxi_Coor_dans_directionGM(theta);
         nb_exterieur++;
        }
      else // si nb_exterieur est diff de zero, on le remet à zéro, car on est de nouveau à l'intérieur
       { if (nb_exterieur != 0) nb_exterieur=0; };
      
      // si cela fait nb_max_a_l_exterieur fois que l'on est à l'extérieur on arrête
      if ( nb_exterieur >= nb_max_a_l_exterieur)   break;  
	  // calcul du point correspondant dans l'element, qui remplace theta
      ph_i = elgeom.Phi_point(theta_repere); // fonctions d'interpolation en theta	
	  A = (met->*PointM) (tab_noeud,ph_i);
      AM = M-A; // nouveau gap
      compteur++;
      if (compteur > nb_boucle_sur_delta_thetai)
       { // cela signifie que l'on n'arrive pas à converger, petit message si besoin 
         // et on choisit le point initial pour le teste
         if (ParaGlob::NiveauImpression() >= 4)
           cout << " \n ** warning, l'algorithme de recherche d'un point a l'interieur d'un "
                << "element ne converge pas, utilisation du repere a l'origine"
                << "\nElemThermi::Interne(Coordonnee& M)" << endl;
         if (elgeom.Interieur(theta_initial))  {theta=theta_initial; }
         else                                  {theta=theta_initial; }; 
         break; 
        };          
     };
//   // à la sortie de la boucle on calcul les thétas correspondants à M final
//   dph_i = elgeom.Dphi(theta_repere); // derivees des fonctions d'interpolation en theta
//   ph_i = elgeom.Phi(theta_repere); // fonctions d'interpolation en theta	
//   // calcul de la base naturelle et duale  en fonction des coord a t 
//   (met->*BaseND)(tab_noeud,dph_i,ph_i,bB,bH);
//   // amelioration  des coordonnees locales du point M dans le repere en theta
//   for (int i=1;i<=dim_glob;i++)
//	     theta(i) += AM * bH(i).Vect();
//   
   // -- maintenant on regarde si le point est interne ou pas a l'element
   int retour=0;
   if (elgeom.Interieur(theta))
      { // on regarde si le point est intérieur d'une toute petite précision 
           // calcul d'un point  extreme de l'élément dans le sens de M
        theta_repere = elgeom.Maxi_Coor_dans_directionGM(theta);
        double diff =(theta_repere-theta).Norme();
        if (diff <= prec_tehetai_interne)
              {retour = 2;}  // !-! à l'erreur près, le point est sur la frontière
        else  {retour = 1;}; // !-! point à l'intérieur
       }
   else  {return 0;};      // !-! point à l'extérieur -> retour directe
   //      --- cas particuliers ---
   // -- dans le cas où "retour" est différent de 0, on regarde les cas particuliers
   switch (dim_glob)
   {case 3:
   	{switch (dim_loc)
        { case 3: break; // cas pas de pb
			 case 2: // cas des "coques ou plaques" on regarde si le point est réellement à l'intérieur (+- 1/2 épaisseur)
   		  { double epaisseur = Epaisseurs(temps,theta);
	          A = (met->*PointM) (tab_noeud,ph_i); // le point projeté sur la surface de l'élément
             AM = M-A; double delta_hauteur = 0.5 * epaisseur - AM.Norme();
             if (Dabs(delta_hauteur) <= prec_tehetai_interne)
                 {retour = 2;}  // !-! à l'erreur près, le point est sur la frontière
             else if (delta_hauteur < 0.)
                 {retour = 0;}  // !-! point externe
             else   {retour = 1;}; // !-! point à l'intérieur
             break;
   		   }
   	     default:
   	       cout << " \n *** erreur : cas non implemente pour l'instant "
   	            << " dim_glob= " << dim_glob << " dim_loc= " << dim_loc 
   	            << "\n ElemThermi::Interne(... ";
   	       Sortie(1);
   	   };// -- fin du switch dim_loc pour dim_glob=3
   	 break;          		
   	 } // fin dim_glob=3
	 case 2: 
   	{switch (dim_loc)
		  {
	      case 2: break; // pas de pb
   	   case 1: case 3: // cas des poutres et volume (pb non consistance)   pour l'instant non traité
					 // donc on laisse aller à la valeur par défaut
			default:
   	     cout << " \n *** erreur : cas non implemente pour l'instant "
   	          << " dim_glob= " << dim_glob << " dim_loc= " << dim_loc 
   	          << "\n ElemThermi::Interne(... ";
   	     Sortie(1);
   	   };// -- fin du switch dim_loc pour dim_glob=3
   	 break;          		
   	 } // fin dim_glob=2
   	default: // default pour le switch sur di_glob
   	     cout << " \n *** erreur : cas non implemente pour l'instant "
   	          << " dim_glob= " << dim_glob << " dim_loc= " << dim_loc 
   	          << "\n ElemThermi::Interne(... ";
   	     Sortie(1);
   	break;
   }; // -- fin du switch sur dim_glob
   // retour 
   return retour;
  };
   
// Calcul du residu local et de la raideur locale,
//  pour le schema implicite
// cald_Dvirtuelle = indique si l'on doit calculer la dérivée de la vitesse de déformation virtuelle 
void ElemThermi::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)
   {  
      int nbddl = tab_ddl.NbDdl();
      volume = 0. ; // init
      Vecteur d_jacobien_tdt;
      energie_totale.Inita(0.); // init de l'énergie totale sur l'élément
      E_elem_bulk_tdt = E_elem_bulk_t; P_elem_bulk = 0.; // init pour l'énergie et la puissance associées au bulk
      // init par défaut pour un cas où la pression ne bouge pas 
      const double  P_t = 1.; // pour l'instant on considère une pression unitaire
      const double  P = 1.; // pour l'instant on considère une pression unitaire
    
      for (int ni =1; ni<= nbint; ni++)  // boucle sur les pt d'integ
       {def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
        PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
        EnergieThermi& energ =  tab_energ(ni);
        EnergieThermi& energ_t =  tab_energ_t(ni);
        CompThermoPhysiqueAbstraite::SaveResul* saveTP=tabSaveTP(ni);  // les données spécifique thermo physiques
        
        if ((loiTP->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
            || (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),ni);
        const Met_abstraite::Impli& ex = loiTP->Cal_implicit
                 (P_t,saveTP,P,*def,tab_ddl,ptIntegThermi,d_gradTB,d_fluxH,pa,loiTP
		  	             ,dilatation,energ,energ_t,premier_calcul_thermi_impli_expli);

        // on calcul éventuellement la dérivée de la vitesse de gradient virtuelle
        Tableau2 <CoordonneeB>* d2_gradTB =  d2_gradTB_(ni); // pour simplifier
        if (cald_DGradTvirtuelle)
             // si en retour le pointeur est null cela signifie qu'il ne faut pas considérer d2_gradTB
             d2_gradTB =def->Der2GradDonneeInterpoleeScalaire(d2_gradTB,TEMP);
                          
//        // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
//        DeuxDoubles delta_ener_bulk_vol; // init
//        if (bulk_viscosity) 
//                delta_ener_bulk_vol=ModifContrainteBulk(*ex.gijHH_tdt,sigHH,DepsBB_tdt,(*sig_bulkHH) );
        double rap=1.; // pour le calcul éventuel du rapport de jacobien actuel
        if (pa.MaxiVarJacobien() > 1.) // cas de la demande du calcul du rapport de jacobien 
           {if (Dabs(*(ex.jacobien_tdt)) > (*(ex.jacobien_0))) {rap = Dabs(*(ex.jacobien_tdt)) / (*(ex.jacobien_0));}
            else                                     {rap = (*(ex.jacobien_0)) / Dabs(*(ex.jacobien_tdt)) ;}
            };  
                      
        if ((*(ex.jacobien_tdt) <= 0.)|| (std::isinf(*(ex.jacobien_tdt)))||( std::isnan(*(ex.jacobien_tdt)))) // vérif du jacobien
           { if ((std::isinf(*(ex.jacobien_tdt)))||( std::isnan(*(ex.jacobien_tdt))))
               // on met un message quelque soit le niveau d'impression
                { cout << "\n **********  attention on a trouve un jacobien infini ou nan !! = ("<<*(ex.jacobien_tdt)<<")********* " << endl; };
             if (ParaGlob::NiveauImpression() >= 1)
               { cout << "\n **********  attention on a trouve un jacobien  negatif!! ("<<*(ex.jacobien_tdt)<<")********* " << endl; };
           };
        if ((*(ex.jacobien_tdt) <= 0.) && (pa.JabobienNegatif()!= 0)) // vérif du jacobien  et traitement adéquate si besoin
            { switch (pa.JabobienNegatif())
              { case 1:   
                  { cout << "\n on annulle sa contribution. Element nb: " <<  Num_elt() << " nb d'integ : " << ni;
                    break;
                   }
                case 2:
                  { // on génère une exception
                    throw ErrJacobienNegatif_ElemThermi();break;
                   }
                // dans les autres cas on ne fait rien   
                };       
              }
        else if ((pa.MaxiVarJacobien() > 1.) && ( rap > pa.MaxiVarJacobien()))
              { if (ParaGlob::NiveauImpression() >= 3)
                 { cout << "\n ***  attention la variation maximal du jacobien est atteinte  !! *** " 
                        << "\n ele= "<<  Num_elt() << " nbi=  " << ni << " jacobien= " << *(ex.jacobien_tdt) << " jacobien_0= "
                        << (*(ex.jacobien_0)) << endl; };
                // on génère une exception
                throw ErrVarJacobienMini_ElemThermi();break;
               }       
        else 
         {  // on prend en compte toutes les  variations
            double poid_jacobien= poids(ni) * *(ex.jacobien_tdt);
            energie_totale += energ * poid_jacobien;
//            E_elem_bulk_tdt += delta_ener_bulk_vol.un * poid_jacobien;
//            P_elem_bulk += delta_ener_bulk_vol.deux * poid_jacobien;
            volume += *(ex.jacobien_tdt) * poids(ni);
            if (d2_gradTB == NULL)
             {for (int j =1; j<= nbddl; j++)  // 1ere boucle sur les ddl
                {(*residu)(j) -= d_gradTB(j) * ptIntegThermi.FluxH() * poid_jacobien;
                 for (int i =1; i<= nbddl; i++)  // 2ere boucle sur les ddl
                   (*raideur)(j,i) += d_gradTB(j) * d_fluxH(i) * poid_jacobien;
  //						raideur->Affiche();
                }
             }
            else
             {for (int j =1; j<= nbddl; j++)  // 1ere boucle sur les ddl
                {(*residu)(j) -= d_gradTB(j) * ptIntegThermi.FluxH() * poid_jacobien;
                 for (int i =1; i<= nbddl; i++)  // 2ere boucle sur les ddl
                   (*raideur)(j,i) += d_gradTB(j) * d_fluxH(i) * poid_jacobien
                                      + (*d2_gradTB)(j,i) * ptIntegThermi.FluxH() * poid_jacobien;
  //						raideur->Affiche();
                }
             };
             
         };// fin du if sur la valeur non négative du jacobien 
//        if (bulk_viscosity) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
//           { sigHH -= (*sig_bulkHH); }  // ceci pour la sauvegarde  ou autre utilisation       
 
        } // fin boucle sur les points d'intégration
	  // --- intervention dans le cas d'une stabilisation d'hourglass	  
		// stabilisation pour un calcul implicit
		if (type_stabHourglass)
          Cal_implicit_hourglass();
		  
     // liberation des tenseurs intermediaires
     LibereTenseur(); 
   };
		
// Calcul uniquement du residu local a l'instant t ou pas suivant la variable atdt
// pour le schema explicit par exemple 
void ElemThermi::Cal_explicit (DdlElement & tab_ddl,Tableau <CoordonneeB >& d_gradTB
                ,int nbint,const Vecteur& poids,const ParaAlgoControle & pa,bool atdt)
            
   {  double jacobien; // def du jacobien
//      int nbddl = tab_ddl.NbDdl();
//      energie_totale.Inita(0.); // init de l'énergie totale sur l'élément
//      volume = 0. ; // init
//      E_elem_bulk_tdt = E_elem_bulk_t; P_elem_bulk = 0.; // init pour l'énergie et la puissance associées au bulk
//      // init éventuel de la contrainte de bulk viscosity
//      TenseurHH* sigHH_t_1 = (*lesPtIntegThermiInterne)(1).SigHH_t(); // simplification d'écriture
//      if (bulk_viscosity)  
//          { if (sig_bulkHH == NULL)  {sig_bulkHH = NevezTenseurHH(*sigHH_t_1);}
//            else if (sig_bulkHH->Dimension() != sigHH_t_1->Dimension())
//               {delete sig_bulkHH; sig_bulkHH = NevezTenseurHH(*sigHH_t_1);}; 
//           };               
//      for (int ni =1; ni<= nbint; ni++)  // boucle sur les pt d'integ
//        { def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
//          PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
//          TenseurHH & sigHH_t = *(ptIntegThermi.SigHH_t());
//          TenseurHH & sigHH = *(ptIntegThermi.SigHH());
//          TenseurBB & epsBB = *(ptIntegThermi.EpsBB());
//          TenseurBB & DepsBB_ = *(ptIntegThermi.DepsBB());
//          EnergieThermi& energ =  tab_energ(ni);
//          EnergieThermi& energ_t =  tab_energ_t(ni);
//          CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
//          if (loiTP != NULL) {sTP = tabSaveTP(ni);};         // au pt d'integ si elles existes 
//          if (loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
//            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),ni);
//          DeuxDoubles delta_ener_bulk_vol; // init
//          if (atdt)
//             { //const Met_abstraite::Expli_t_tdt& ex = loiComp->Cal_explicit_tdt
//               //              (tabSaveDon(ni), *def,tab_ddl,ptIntegThermi,d_epsBB,jacobien
//               //              ,sTP,loiTP,dilatation,energ,energ_t,ppremier_calcul_Thermi_impli_expli;
//               const Met_abstraite::Expli_t_tdt ex;
//               // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
//               if (bulk_viscosity) 
//                    delta_ener_bulk_vol=ModifContrainteBulk(*ex.gijHH_tdt,sigHH,DepsBB_,(*sig_bulkHH) );
//              }               
//          else   
//             { //const Met_abstraite::Expli& ex = loiComp->Cal_explicit_t
//               //              (tabSaveDon(ni), *def,tab_ddl,ptIntegThermi,d_epsBB,jacobien
//                //             ,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_Thermi_impli_expli);
//               const Met_abstraite::Expli_t_tdt ex;
//               // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
//               if (bulk_viscosity) 
//                    delta_ener_bulk_vol=ModifContrainteBulk(*ex.gijHH_t,sigHH,DepsBB_,(*sig_bulkHH) );
//              } 
//         if ((jacobien <= 0.)|| (std::isinf(jacobien))||( std::isnan(jacobien))) // vérif du jacobien
//            { if ((std::isinf(jacobien))||( std::isnan(jacobien)))
//                // on met un message quelque soit le niveau d'impression
//                { cout << "\n **********  attention on a trouve un jacobien infini ou nan !! = ("<<jacobien<<")********* " << endl; };
//              if (ParaGlob::NiveauImpression() >= 1)
//                { cout << "\n **********  attention on a trouve un jacobien  negatif!! ("<<jacobien<<")********* " << endl; };
//            };
//          if ((jacobien <= 0.) && (pa.JabobienNegatif() != 0)) // vérif du jacobien  et traitement adéquate si besoin
//            { switch (pa.JabobienNegatif())
//              { case 1:   
//                  { cout << "\n on annulle sa contribution. Element nb: " <<  Num_elt() << " nb d'integ : " << ni;
//                    break;
//                   }
//                // dans les autres cas on ne fait rien   
//                };       
//              }
//          else
//          {  double poid_jacobien= poids(ni) * jacobien;
//             energie_totale += energ * poid_jacobien;
//             E_elem_bulk_tdt += delta_ener_bulk_vol.un * poid_jacobien;
//             P_elem_bulk += delta_ener_bulk_vol.deux * poid_jacobien;
//             volume += poid_jacobien;
//             for (int j =1; j<= nbddl; j++)
//              {(*residu)(j) -=  (sigHH && (*(d_epsBB(j)))) * poid_jacobien;
//              }
//            }  
//          if (bulk_viscosity) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
//           { sigHH -= (*sig_bulkHH); }  // ceci pour la sauvegarde  ou autre utilisation       
//        }         
//	  // --- intervention dans le cas d'une stabilisation d'hourglass	  
//		// stabilisation pour un calcul implicit
//		if (type_stabHourglass)
//        Cal_explicit_hourglass(atdt);
//     // liberation des tenseurs intermediaires
     LibereTenseur();
   };
		
// 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
// epsBB = deformation, sigHH = contrainte, d_epsbb = variation des def
// 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 ElemThermi::Cal_matGeom_Init (Mat_pleine & matGeom, Mat_pleine & matInit,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_Dvirtuelle)
   {  // le début du programme est semblable au calcul implicit
      // seule la constitution des matrices finales est différente
//      int nbddl = tab_ddl.NbDdl();
//      double jacobien; // def du jacobien
//      Vecteur d_jacobien_tdt;
//      volume = 0. ; // init
//
//      for (int ni =1; ni<= nbint; ni++)  // boucle sur les pt d'integ
//       {def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
//        PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
//        TenseurHH & sigHH = *(ptIntegThermi.SigHH());
//        TenseurBB & DepsBB_tdt = *(ptIntegThermi.DepsBB());
//        EnergieThermi& energ =  tab_energ(ni);
//        EnergieThermi& energ_t =  tab_energ_t(ni);
//        CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
//        if (loiTP != NULL) {sTP = tabSaveTP(ni);};         // au pt d'integ si elles existes 
//        if (loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
//            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),ni);
////        loiComp->Cal_flamb_lin(tabSaveDon(ni), *def,tab_ddl,ptIntegThermi,d_epsBB,jacobien,d_jacobien_tdt,d_sigHH
////                              ,pa,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_Thermi_impli_expli);
//        // on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
//        Tableau2 <TenseurBB *>& d2_gradTB =  d2_gradTB_(ni); // pour simplifier
//        if (cald_Dvirtuelle) 
//             def->Cal_var_def_virtuelle(false,d2_gradTB);
//        // calcul du volume
//        volume += jacobien * poids(ni);
//        for (int j =1; j<= nbddl; j++)  // 1ere boucle sur les ddl
//           for (int i =1; i<= nbddl; i++)  // 2ere boucle sur les ddl
//              {matInit(j,i) += ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien * poids(ni);
//               matGeom(j,i) += (sigHH && (*(d2_gradTB(j,i)))) * jacobien * poids(ni);
//              }
//        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
   };
                
// Calcul de la matrice masse selon différent choix donné par type_matrice_masse,
void ElemThermi::Cal_Mat_masse (DdlElement & tab_ddl,Enum_calcul_masse type_matrice_masse,
                            int nbint,const Tableau <Vecteur>& taphi,int nbne
                            ,const Vecteur& poids)
   {  
      int nbddl = tab_ddl.NbDdl();
      double jacobien_0; // def du jacobien initial
      // initialisation de la matrice masse
      mat_masse->Initialise ();
      
      // le calcul est fonction du type de matrice masse demandé 
      switch (type_matrice_masse)
	   {case MASSE_DIAG_COEF_EGAUX  : 
	     // même masse sur chaque noeud
	     { // calcul du volume totale
	       double vol_totale =0.;
           for (int ni =1; ni<= nbint; ni++)  // boucle sur les pt d'integ
            { defMas->ChangeNumInteg(ni);
              const Met_abstraite::Dynamiq& ex = defMas->Cal_matMass();
              jacobien_0 = *ex.jacobien_0;
              vol_totale += jacobien_0 * poids(ni);          
	         }
	       // la masse élémentaire pour chaque noeud  
	       double masse_elementaire = vol_totale * masse_volumique / tab_noeud.Taille();
	       // on initialise toutes les valeurs de la matrice masse à la valeur moyenne
	       mat_masse->Initialise (masse_elementaire );  
		   break;
		  }	
		case MASSE_DIAG_COEF_VAR :
		  // masse diagonale répartie au prorata des fonctions d'interpolation
		  // (cf page 304, tome 2 Batoz)
	     { // définition d'un vecteur qui contiendra l'intégral de chaque
	       // fonction d'interpolation (ici la masse est considérée constante)
	       Vecteur fonc2(nbne);
	       // calcul de fonc2 et du volume totale
	       double vol_totale =0.;
           for (int ni =1; ni<= nbint; ni++)  // boucle sur les pt d'integ
            { defMas->ChangeNumInteg(ni);
              const Met_abstraite::Dynamiq& ex = defMas->Cal_matMass();
              jacobien_0 = *ex.jacobien_0;
              for (int ne =1;ne<=nbne;ne++)
                 fonc2(ne) += taphi(ni)(ne) * taphi(ni)(ne) * jacobien_0 * poids(ni);
              vol_totale += jacobien_0 * poids(ni);          
	         }
	       // calcul de la somme des composantes de fonc2
	       double somme_fonc2 = 0.;
	       for (int ne = 1;ne<=nbne;ne++) somme_fonc2 += fonc2(ne);
	       // calcul des termes diagonaux de la matrice de masse
           int dimension = ParaGlob::Dimension();       
           if(ParaGlob::AxiSymetrie())
           // cas d'élément axisymétrique, dans ce cas on ne prend en compte que les 
           // dimension-1 coordonnées donc on décrémente
              dimension--;
           int ix=1;       
           for (int ne =1; ne<= nbne; ne++) 
             { double inter_toto = vol_totale * masse_volumique * fonc2(ne) / somme_fonc2;
               for (int i=1;i<= dimension;i++,ix++)
                 (*mat_masse)(1,ix) = inter_toto;
              }   
		   break;
		  }	
		case MASSE_CONSISTANTE :
		  // matrice masse complète avec les mêmes fonctions d'interpolation 
		  // que pour la raideur.
	     { 
           int dimension = ParaGlob::Dimension();       
           if(ParaGlob::AxiSymetrie())
           // cas d'élément axisymétrique, dans ce cas on ne prend en compte que les 
           // dimension-1 coordonnées donc on décrémente
              dimension--;
           for (int ni =1; ni<= nbint; ni++)  // boucle sur les pt d'integ
            {defMas->ChangeNumInteg(ni); // def du pt d'intégration
             // récup des fonctions d'interpolation et du jacobien
             const Met_abstraite::Dynamiq& ex = defMas->Cal_matMass();
             double masse_i = (*ex.jacobien_0) * masse_volumique;
             // calcul de la matrice
            int ix=1;       
            for (int ne =1; ne<= nbne; ne++) 
             for (int a=1;a<= dimension;a++,ix++)
              { int jy=1; 
                for (int me =1; me<= nbne; me++) 
                  for (int b=1;b<= dimension;b++,jy++)
                    if (a==b) // seule des termes de même indice sont non nulles
                       (*mat_masse)(ix,jy) +=   taphi(ni)(ne)* taphi(ni)(me) * masse_i * poids(ni);
               }
            }   
		   break;
		  }	
		default :
			cout << "\nErreur : valeur incorrecte du type Enum_calcul_masse !\n";
			cout << "ElemThermi::Cal_Mat_masse (... \n";
			Sortie(1);
	    };

     // liberation des tenseurs intermediaires
     LibereTenseur();   
   };

  // ------------ calcul de second membre  ou de l'action des efforts extérieures -----------
  
// calcul des seconds membres suivant les chargements 
// cas d'un chargement surfacique de direction constante, 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 à tdt ou t suivant la variable atdt
Vecteur& ElemThermi::SM_charge_surf_E (DdlElement & ,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & ,bool atdt)
   {  // définition du vecteur de retour 
      Vecteur& SM = *((*res_extS)(nSurf));
      int ni; // compteur globale de point d'integration
      Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.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
        double jacobien;
        bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
        if (atdt)
          { const Met_abstraite::Expli_t_tdt& ex = defS.Cal_explicit_tdt(premier_calcul);
            jacobien = (*ex.jacobien_tdt);
            }
        else    
          { const Met_abstraite::Expli& ex = defS.Cal_explicit_t(premier_calcul);
            jacobien = (*ex.jacobien_t);
            }
        int ix=1; int dimf = force.Dimension();       
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* force(i) * (poids(ni) * jacobien);   
        }
      // retour du second membre  
      return SM;  
    };
    
// calcul des seconds membres suivant les chargements 
// cas d'un chargement surfacique  de direction constante, 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
// -> 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 ElemThermi::SMR_charge_surf_I (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa)
   {  bool avec_raid = pa.Var_charge_externe(); // récup indicateur
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      int ni; // compteur globale de point d'integration
      Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
      int nbddl = ddls.NbDdl();
      Vecteur& SM = (*((*res_extS)(nSurf)));  //  "      "          "
      Mat_pleine & KM = (*((*raid_extS)(nSurf))); //  "      "       "
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.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 implicit 
        const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul); 
        int ix=1; int dimf = force.Dimension();       
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
           {  SM(ix) +=  taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.jacobien_tdt));   
              // dans le cas avec_raideur on calcul la contribution à la raideur
              if (avec_raid)
                for (int j =1; j<= nbddl; j++)  
                  KM(ix,j) += 
                          taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j));
             }     
        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
    Element::ResRaid el;
    el.res = (*res_extS)(nSurf); 
    el.raid = (*raid_extS)(nSurf);       
    return el;
    };
                                
// 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& ElemThermi::SM_charge_pres_E (DdlElement & ,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,double  pression
                                ,const ParaAlgoControle & ,bool atdt)
   {  
      Deformation* definter=NULL; // variable de travail transitoire
      Vecteur* SM_P = NULL; 
      if (!(ParaGlob::AxiSymetrie())) 
           { definter = defSurf(nSurf);    // le cas normal est non axisymétrique
             SM_P = ((*res_extS)(nSurf));
           }
      else { definter = defArete(nSurf);  // en axisymétrie, c'est une def d'arête
             SM_P = ((*res_extA)(nSurf));
           };
      Deformation & defS = *definter;  // pour simplifier l'écriture
      Vecteur& SM = *SM_P;          //  "      "          "

      // définition du vecteur de retour
      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 (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.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
        Met_abstraite::Expli ex;
        if (atdt)
           ex = ((defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t());
        else    
           ex.operator=( defS.Cal_explicit_t(premier_calcul));
        // détermination de la normale à la surface
        #ifdef MISE_AU_POINT	 	 
         // test pour savoir si l'on a réellement affaire à une surface
         if ((ex.giB_t->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3)) 
          { cout << "\n erreur, il doit y avoir deux vecteurs pour la surface  de dimension 3"
               << " ElemThermi::SM_charge_pres (DdlElement & ...nbvec=  " << ex.giB_t->NbVecteur()
               << " dim = " << ParaGlob::Dimension();
            Sortie(1);     
           }
        #endif
        // la normale vaut le produit vectoriel des 2 premiers vecteurs
        Coordonnee normale = Util::ProdVec_coorBN( (*ex.giB_t)(1), (*ex.giB_t)(2));
        // calcul de la normale normée et pondérée de la pression
        // la pression est positive qu'en elle appuis (diminution de volume)
        normale.Normer(); normale *= -pression;
        int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3    
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_t));   
        }
      // retour du second membre  
      return SM;  
    };
    
// 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 -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
Element::ResRaid ElemThermi::SMR_charge_pres_I(DdlElement & ddlS,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,double  pression
                                ,const ParaAlgoControle & pa)
   {  
      int ni; // compteur globale de point d'integration
      
      Deformation* definter=NULL; // variable de travail transitoire
      Vecteur* SM_P = NULL; Mat_pleine* KM_P = NULL;
      if (!(ParaGlob::AxiSymetrie())) 
           { definter = defSurf(nSurf);    // le cas normal est non axisymétrique
             SM_P = ((*res_extS)(nSurf));
             KM_P = ((*raid_extS)(nSurf));
           }
      else { definter = defArete(nSurf);  // en axisymétrie, c'est une def d'arête
             SM_P = ((*res_extA)(nSurf));
             KM_P = ((*raid_extA)(nSurf));
           };
      Deformation & defS = *definter;  // pour simplifier l'écriture
      Vecteur& SM = *SM_P;          //  "      "          "
      Mat_pleine & KM = *KM_P;      //  "      "       "

      int nbddl = ddlS.NbDdl();
      // controle
      bool avec_raid = pa.Var_charge_externe();
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.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::Impli& ex = defS.Cal_implicit(premier_calcul); 
        // détermination de la normale à la surface
        #ifdef MISE_AU_POINT	 	 
         // test pour savoir si l'on a réellement affaire à une surface
         if ((ex.giB_tdt->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3)) 
          { cout << "\n erreur, il doit y avoir deux vecteurs pour la surface  de dimension 3"
               << " ElemThermi::SM_charge_pres (DdlElement & ...nbvec=  " << ex.giB_tdt->NbVecteur()
               << " dim = " << ParaGlob::Dimension();
            Sortie(1);     
           }
        #endif
        // la normale vaut le produit vectoriel des 2 premiers vecteurs
        CoordonneeB normale = Util::ProdVec_coorB( (*ex.giB_tdt)(1), (*ex.giB_tdt)(2));
        // calcul de la variation de la normale
        // 1) variation du produit vectoriel
        Tableau <CoordonneeB > D_pasnormale = 
                 Util::VarProdVect_coorB( (*ex.giB_tdt)(1),(*ex.giB_tdt)(2),(*ex.d_giB_tdt));
        // 2) de la normale        
        Tableau <CoordonneeB> D_normale = Util::VarUnVect_coorB(normale,D_pasnormale,normale.Norme());
        // calcul de la normale normée et pondérée de la pression
        // la pression est positive qu'en elle appuis (diminution de volume)
        normale.Normer(); normale *= -pression;
        // également pondération de la variation de la 
        for (int ihi =1;ihi<= nbddl;ihi++)
           D_normale(ihi) *= -pression;
        
        int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3    
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
           {  SM(ix) +=  taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_tdt)); 
              // dans le cas avec_raideur on calcul la contribution à la raideur
              if (avec_raid)
                for (int j =1; j<= nbddl; j++)  
                  KM(ix,j) -= 
                          taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j))
                         + taphi(ni)(ne)* D_normale(j)(i) * (poids(ni) * (*ex.jacobien_tdt)) ;
             }     
        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
    Element::ResRaid el;
    el.res = SM_P;
    el.raid = KM_P;
    return el;
    };
	    
// calcul des seconds membres suivant les chargements 
// cas d'un chargement lineique, sur les aretes frontières des éléments
// force indique la force lineique appliquée
// retourne  le second membre résultant
// nArete : le numéro de la surface externe
// calcul à l'instant tdt ou t en fonction de la variable atdt
Vecteur& ElemThermi::SM_charge_line_E (DdlElement & ,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & ,bool atdt)
   {  // définition du vecteur de retour 
      Vecteur& SM = *((*res_extA)(nArete));
      int ni; // compteur globale de point d'integration
      Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
      for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.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
        Met_abstraite::Expli ex;
        bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
        if (atdt)
           ex = (defA.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
        else    
           ex = defA.Cal_explicit_t(premier_calcul);
        int ix=1; int dimf = force.Dimension();       
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.jacobien_t));   
        }
      // retour du second membre  
      return SM;  
    };
    
// 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 -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
Element::ResRaid ElemThermi::SMR_charge_line_I(DdlElement & ddlA,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa)
   {  int ni; // compteur globale de point d'integration
      Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
      int nbddl = ddlA.NbDdl();
      Vecteur& SM = (*((*res_extA)(nArete)));  //  "      "          "
      Mat_pleine & KM = (*((*raid_extA)(nArete))); //  "      "       "
      // controle
      bool avec_raid = pa.Var_charge_externe();
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.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 implicit 
        const Met_abstraite::Impli& ex = defA.Cal_implicit(premier_calcul); 
        int ix=1; int dimf = force.Dimension();
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
           {  SM(ix) += 
                         taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.jacobien_tdt));   
              // dans le cas avec_raideur on calcul la contribution à la raideur
              if (avec_raid)
                for (int j =1; j<= nbddl; j++)  
                  KM(ix,j) -=
                          taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j));
             }     
        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
    Element::ResRaid el;
    el.res = (*res_extA)(nArete); 
    el.raid = (*raid_extA)(nArete);       
    return el;
    };	    
	    
// cas d'un chargement lineique suiveur, sur les aretes frontières des éléments
// force indique la force lineique appliquée
// retourne  le second membre résultant
// nArete : le numéro de la surface externe
// calcul à l'instant tdt ou t en fonction de la variable atdt
Vecteur& ElemThermi::SM_charge_line_Suiv_E (DdlElement & ,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & ,bool atdt)
   {  // définition du vecteur de retour 
      Vecteur& SM = *((*res_extA)(nArete));
      int ni; // compteur globale de point d'integration
      Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.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
        Met_abstraite::Expli ex;
        if (atdt)
           ex = (defA.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
        else    
           ex = defA.Cal_explicit_t(premier_calcul);
        #ifdef MISE_AU_POINT	 	 
         // test pour savoir si l'on est en dimension 2 ou en axy sinon pas possible
         if ((ParaGlob::Dimension() != 2) & (!(ParaGlob::AxiSymetrie()))) 
          { cout << "\n erreur, pour l'instant seul la dimension 2 est permise avec les forces suiveuses"
               << "\n si pb contacter gerard Rio ! "
               << " ElemThermi::SM_charge_line_Suiv_E ( ..."
               << " \n dim = " << ParaGlob::Dimension();
            Sortie(1);     
           }
        #endif
        // calcul du repère locale initiale orthonormée
//        Coordonnee2 v1((*ex.giB_0)(1).Vect()); v1.Normer();
        // on récupère que les deux premières coordonnées mais on travaille en 3D pour intégrer le cas axi
		  Coordonnee v1((*ex.giB_0).Coordo(1)); v1.Normer();
		  int dim = v1.Dimension();
        Coordonnee v2(dim); v2(1)=-v1(2); v2(2)=v1(1);
        // composante du vecteur force linéique dans le repère initiale
        Coordonnee force_0(v1*force,v2*force);
        //repère locale actuel orthonormée
        Coordonnee giB_t = (*ex.giB_t)(1).Coor();
        Coordonnee vf1(giB_t); 
        Coordonnee vf2(dim);vf2(1)=-vf1(2);vf2(2)=vf1(1);
        // composantes actuelles de la force
        Coordonnee force_t = force_0(1)*vf1 + force_0(2)*vf2;
        // calcul du résidu
        int ix=1; int dimf = 2 ; // ici uniquement en dim 2 sinon pb 
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_t));   
        }
      // retour du second membre  
      return SM;  
    };
    
// cas d'un chargement lineique suiveur, 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 -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
Element::ResRaid ElemThermi::SMR_charge_line_Suiv_I(DdlElement & ddlA,int nArete
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa)
   {  int ni; // compteur globale de point d'integration
      Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
      int nbddl = ddlA.NbDdl();
      Vecteur& SM = (*((*res_extA)(nArete)));  //  "      "          "
      Mat_pleine & KM = (*((*raid_extA)(nArete))); //  "      "       "
      // controle
      bool avec_raid = pa.Var_charge_externe();
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.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 implicit 
        const Met_abstraite::Impli& ex = defA.Cal_implicit(premier_calcul); 
        #ifdef MISE_AU_POINT	 	 
         // test pour savoir si l'on est en dimension 2 ou en axy sinon pas possible
         if ((ParaGlob::Dimension() != 2) & (!(ParaGlob::AxiSymetrie()))) 
          { cout << "\n erreur, pour l'instant seul la dimension 2 est permise avec les forces suiveuses"
               << "\n si pb contacter gerard Rio ! "
               << " ElemThermi::SM_charge_line_Suiv_E ( ..."
               << " \n dim = " << ParaGlob::Dimension();
            Sortie(1);     
           }
        #endif
        // calcul du repère locale initiale orthonormée
//        Coordonnee2 v1((*ex.giB_0)(1).Vect()); v1.Normer();
        // on récupère que les deux premières coordonnées mais on travaille en 3D pour intégrer le cas axi
        Coordonnee v1((*ex.giB_0).Coordo(1)); v1.Normer();
        int dim = v1.Dimension();
        Coordonnee v2(dim); v2(1)=-v1(2); v2(2)=v1(1);
        // composante du vecteur force linéique dans le repère initiale
        Coordonnee force_0(v1*force,v2*force);
        //repère locale actuel orthonormée
        Coordonnee giB_tdt = (*ex.giB_tdt)(1).Coor();
        Coordonnee vf1(giB_tdt); 
        double vf1Norme = vf1.Norme();vf1 /= vf1Norme;
        Coordonnee vf2(dim);vf2(1)=-vf1(2);vf2(2)=vf1(1);
        // composantes actuelles de la force
        Coordonnee force_t = force_0(1)*vf1 + force_0(2)*vf2;

        // calcul du résidu
        int ix=1; int dimf = 2 ; // ici uniquement en dim 2 sinon pb 
        Coordonnee dvf1,dvf2,d_giB_tdt; // def de vecteurs intermediaires     
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
           {  SM(ix) +=  taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_tdt));   
              // dans le cas avec_raideur on calcul la contribution à la raideur
              if (avec_raid)
                for (int j =1; j<= nbddl; j++)
                  { // calcul des variations de la force
                    d_giB_tdt = ((((*ex.d_giB_tdt)(j))(1)).Coor());
                    dvf1 = Util::VarUnVect_coor // variation de vf1 par rapport au ddl j
                          (giB_tdt,d_giB_tdt,vf1Norme);
                    dvf2 = -(dvf1 * vf2) * vf1;
                    KM(ix,j) -= taphi(ni)(ne) *  poids(ni)
                                * (force_t(i) *(*ex.d_jacobien_tdt)(j)
                                   + (*ex.jacobien_tdt) * ((force_0(1) * dvf1 + force_0(2) * dvf2))(i)
                                  );
                   }                
             }     
        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
    Element::ResRaid el;
    el.res = (*res_extA)(nArete); 
    el.raid = (*raid_extA)(nArete);       
    return el;
    };	    
	    
// 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& ElemThermi::SM_charge_surf_Suiv_E (DdlElement & ,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& forc
                                ,const ParaAlgoControle & ,bool atdt)
   {  // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
      // *** sinon c'est intractable et long
      // définition du vecteur de retour 
      Vecteur& SM = *((*res_extS)(nSurf));
      double module_f = forc.Norme();
      // dans le cas ou la force est de module nul, on retourne sans calcul
      if (module_f < ConstMath::trespetit)
           return SM;
      #ifdef MISE_AU_POINT	 	 
      // test : a priori non valide pour les éléments axi et pour un pb autre que 3D
      if (( ParaGlob::Dimension() != 3)  || ((ParaGlob::Dimension() == 3)&&(ParaGlob::AxiSymetrie())))         
          { cout << "\n erreur, ce type de chargement n'est valide sur pour la dimension 3 !!  et non axi"
                 << " ElemThermi::SM_charge_line_Suiv_E ( ..."
                 << " \n dim = " << ParaGlob::Dimension();
            Sortie(1);     
           };
      #endif
      int ni; // compteur globale de point d'integration
      Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.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
        Met_abstraite::Expli ex;
        if (atdt)
           ex = (defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
        else    
           ex = defS.Cal_explicit_t(premier_calcul);
        // détermination de la normale à la surface
        // la normale vaut le produit vectoriel des 2 premiers vecteurs
        Coordonnee normale_0 = Util::ProdVec_coorBN( (*ex.giB_0)(1), (*ex.giB_0)(2));
        // on passe en BN pour dire que c'est du B en normal
        const Coordonnee& giBN_t1 = (*ex.giB_t).Coordo(1); const Coordonnee& giBN_t2 = (*ex.giB_t).Coordo(2); 
        Coordonnee normale_t = Util::ProdVec_coor( giBN_t1, giBN_t2);

        // calcul de la normale normée 
        normale_0.Normer(); normale_t.Normer();
        // calcul des composantes de direction de la force dans le repère local initial
        Coordonnee dir_forceH_0(3);
        // ici le fait d'utiliser la méthode .Vect() supprime la variance, mais justement 
        // on ne veut pas de vérification de variance qui est ici incohérente
        dir_forceH_0(1) =  giBN_t1 * forc;  dir_forceH_0(2) =  giBN_t2 * forc;
        dir_forceH_0(3) = normale_0 * forc;
        // calcul de la direction finale
        Coordonnee dir_force_t = dir_forceH_0(1) * giBN_t1 + dir_forceH_0(2) * giBN_t2
                       + dir_forceH_0(3) * normale_t;
        dir_force_t.Normer();
        // calcul de la force finale
        Coordonnee force_t = (forc.Norme()) * dir_force_t;
        // calcul du second membre
        int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3 
        // sauf pour les éléments axisymétriques pour lesquelles on ne considère que les deux premières composantes   
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_t)); 
        }
      // retour du second membre  
      return SM;  
    };
	    
// 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 ElemThermi::SMR_charge_surf_Suiv_I (DdlElement & ddls,int nSurf
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& forc
                                ,const ParaAlgoControle & pa)
   {  // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
      // *** sinon c'est intractable et long
      double module_f = forc.Norme();
      // dans le cas ou la force est de module nul, on retourne sans calcul
      if (module_f < ConstMath::trespetit)
        {  Element::ResRaid el;    
           el.res = (*res_extS)(nSurf); 
           el.raid = (*raid_extS)(nSurf);       
           return el;
         }
      // controle
      bool avec_raid = pa.Var_charge_externe();
      int ni; // compteur globale de point d'integration
      Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
      int nbddl = ddls.NbDdl();
      Vecteur& SM = (*((*res_extS)(nSurf)));  //  "      "          "
      Mat_pleine & KM = (*((*raid_extS)(nSurf))); //  "      "       "
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.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 implicit 
        const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul); 
        #ifdef MISE_AU_POINT	 	 
         // test : a priori non valide pour les éléments axi et pour un pb autre que 3D
         if (( ParaGlob::Dimension() != 3)  || ((ParaGlob::Dimension() == 3)&&(ParaGlob::AxiSymetrie())))         
          { cout << "\n erreur, ce type de chargement n'est valide pour pour la dimension 3 !!  et non axi"
                 << " ElemThermi::SM_charge_line_Suiv_E ( ..."
                 << " \n dim = " << ParaGlob::Dimension();
            Sortie(1);     
           };
         // test pour savoir si l'on a réellement affaire à une surface
         if (ex.giB_t->NbVecteur() != 2)
          { cout << "\n erreur, pb de nombre de vecteurs pour la surface !!"
               << " ElemThermi::SM_charge_surf_Suiv_I  (....  " << ex.giB_t->NbVecteur()
               << " dim = " << ParaGlob::Dimension();
            Sortie(1);     
           }
        #endif
        // détermination de la normale à la surface
        // la normale vaut le produit vectoriel des 2 premiers vecteurs
        Coordonnee normale_0 = Util::ProdVec_coorBN( (*ex.giB_0)(1), (*ex.giB_0)(2));
        const Coordonnee& giB_tdt1 = (*ex.giB_tdt).Coordo(1); const Coordonnee& giB_tdt2 = (*ex.giB_tdt).Coordo(2); 
        Coordonnee normale_t = Util::ProdVec_coor( giB_tdt1, giB_tdt2);
        // calcul des normales normées 
        normale_0.Normer(); 
        double norme_normarle_t = normale_t.Norme();normale_t /= norme_normarle_t;
        // calcul des composantes de direction de la force dans le repère local initial
        Coordonnee dir_forceH_0(3);
        // on ne veut pas de vérification de variance
        dir_forceH_0(1) =  giB_tdt1 * forc;  dir_forceH_0(2) =  giB_tdt2 * forc;
        dir_forceH_0(3) = normale_0 * forc;
        // calcul de la direction finale
        Coordonnee dir_f_t = dir_forceH_0(1) * giB_tdt1 + dir_forceH_0(2) * giB_tdt2
                           + dir_forceH_0(3) * normale_t;
        Coordonnee dir_f_t_normer =  dir_f_t/(dir_f_t.Norme());             
        // calcul de la force finale
        Coordonnee force_t = module_f * dir_f_t_normer;
        // calcul du résidu et éventuellement de la raideur
        int ix=1; int dimf = forc.Dimension();  // normalement ne fonctionne qu'en 3D     
        // sauf pour les éléments axisymétriques pour lesquelles on ne considère que les deux premières composantes   
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_tdt));   
        // dans le cas avec_raideur on calcul la contribution à la raideur
        Coordonnee d_giB_tdt1,d_giB_tdt2,d_dir_f_t; // def de vecteurs intermediaires 
        Coordonnee d_dir_f_t_normer;            //       "   " 
        Coordonnee D_normale_pasnormer,D_normale_t;
        if (avec_raid)
         {for (int ne =1; ne<= nbne; ne++)
            for (int j =1; j<= nbddl; j++)  
               { // calcul de la variation de la normale
                 // 1) variation du produit vectoriel pour la normale finale
                 D_normale_pasnormer = Util::ProdVec_coor(giB_tdt1,(*ex.d_giB_tdt)(j).Coordo(2)) 
                                       + Util::ProdVec_coor((*ex.d_giB_tdt)(j).Coordo(1),giB_tdt2); 
                 // 2) de la normale  finale      
                 D_normale_t = Util::VarUnVect_coor(normale_t,D_normale_pasnormer,norme_normarle_t);
                 // calcul des variations de la force
                 d_giB_tdt1 = ((((*ex.d_giB_tdt)(j))(1)).Coor());
                 d_giB_tdt2 = ((((*ex.d_giB_tdt)(j))(2)).Coor());
                 d_dir_f_t = dir_forceH_0(1) * d_giB_tdt1 + dir_forceH_0(2) * d_giB_tdt2
                       + dir_forceH_0(3) * D_normale_t;
                 d_dir_f_t_normer = Util::VarUnVect_coor( dir_f_t,d_dir_f_t,module_f);
                 for (int i=1;i<= dimf;i++)
                    { ix = i+ (ne-1)*dimf;
                      KM(ix,j) -= taphi(ni)(ne) *  poids(ni)
                                 * (force_t(i) *(*ex.d_jacobien_tdt)(j)
                                     + (*ex.jacobien_tdt) * d_dir_f_t(i)  );
                     }                
                }                
           }
       }
    // liberation des tenseurs intermediaires
    LibereTenseur();   
    Element::ResRaid el;
    el.res = (*res_extS)(nSurf); 
    el.raid = (*raid_extS)(nSurf);       
    return el;
    };

// 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& ElemThermi::SM_charge_vol_E (DdlElement &
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & ,bool atdt )
   {  
      #ifdef MISE_AU_POINT	 	 
      // axisymétrie: pour l'instant non testé a priori avec les forces hydro
      if(ParaGlob::AxiSymetrie())
          { cout << "\n erreur, les charges volumique ne sont utilisable avec les elements axisymetriques"
              << "\n utilisez à la place les charges surfaciques !! "
               << " ElemThermi::SM_charge_vol_E ( ...";
            Sortie(1);     
           };
      #endif
   
      // définition du vecteur de retour 
      Vecteur& SM = * residu;
      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, 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
        Met_abstraite::Expli ex;
        if (atdt)
           ex = (def->Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
        else    
           ex = def->Cal_explicit_t(premier_calcul);
        int ix=1; int dimf = force.Dimension();       
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
             SM(ix) +=  taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.jacobien_t));   
        }
      // retour du second membre  
      return SM;  
    }; 
    
// cas d'un chargement volumique, sur l'élément
// force indique la force volumique appliquée
// retourne  le second membre résultant -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
void ElemThermi::SMR_charge_vol_I(DdlElement & tab_ddl
                                ,const Tableau <Vecteur>& taphi,int nbne
                                ,const Vecteur& poids,const Coordonnee& force
                                ,const ParaAlgoControle & pa)
   {  
      #ifdef MISE_AU_POINT	 	 
      // axisymétrie: pour l'instant non testé a priori avec les forces hydro
      if(ParaGlob::AxiSymetrie())
          { cout << "\n erreur, les charges volumique ne sont utilisable avec les elements axisymetriques"
              << "\n utilisez à la place les charges surfaciques !! "
               << " ElemThermi::SMR_charge_vol_I ( ...";
            Sortie(1);     
           };
      #endif
      int ni; // compteur globale de point d'integration
      int nbddl = tab_ddl.NbDdl();
      Vecteur& SM = *residu;  
      Mat_pleine & KM = *raideur; 
      // controle
      bool avec_raid = pa.Var_charge_externe();
      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 implicit 
        const Met_abstraite::Impli& ex = def->Cal_implicit(premier_calcul); 
        int ix=1; int dimf = force.Dimension();       
        if(ParaGlob::AxiSymetrie())
          // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les 
          // dimf-1 coordonnées donc on décrémente
          dimf--;
        for (int ne =1; ne<= nbne; ne++) 
          for (int i=1;i<= dimf;i++,ix++)
           {  SM(ix) += 
                         taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.jacobien_tdt));   
              // dans le cas avec_raideur on calcul la contribution à la raideur
              if (avec_raid)
                for (int j =1; j<= nbddl; j++)  
                  KM(ix,j) -=
                          taphi(ni)(ne)* force(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j));
             }     
        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
    };
	    
// 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& ElemThermi::SM_charge_hydro_E (DdlElement & ,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 & ,bool atdt)
   {  
      #ifdef MISE_AU_POINT	 	 
      // axisymétrie: pour l'instant non testé a priori avec les forces hydro
      if(ParaGlob::AxiSymetrie())
          { cout << "\n erreur, les charges hydro ne sont pas testees avec les elements axisymetriques"
               << " ElemThermi::SM_charge_hydro_E ( ...";
            Sortie(1);     
           };
      #endif

      // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
      // *** sinon c'est intractable et long
      // définition du vecteur de retour 
      Vecteur& SM = *((*res_extS)(nSurf));
      // dans le cas ou le poid volumique est nul, on retourne sans calcul
      if (poidvol < ConstMath::trespetit)
           return SM;
      int ni; // compteur globale de point d'integration
      Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
       {
        // calcul de la pression au niveau de la position du point d'intégration
        //  1) on récupère la position dans I_a du point d'intégration
        const Coordonnee & M = defS.Position_tdt();
        //  2) calcul de la distance à la surface libre
        Coordonnee MMp= M_liquide - M; double dis = dir_normal_liquide * MMp;
        //  3) calcul de la pression effective uniquement si le point est dans l'eau
        // en fait si le point est hors de l'eau il n'y a pas de pression hydrostatique
        if (dis > 0.)
         { double pression = -dis * poidvol;
           // -- 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
           Met_abstraite::Expli ex;
           if (atdt)
              ex = (defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
           else    
              ex = defS.Cal_explicit_t(premier_calcul);
           // détermination de la normale à la surface
           #ifdef MISE_AU_POINT	 	 
             // test pour savoir si l'on a réellement affaire à une surface
             if ((ex.giB_t->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3)) 
             { cout << "\n erreur, il doit y avoir deux vecteurs pour la surface  de dimension 3"
                    << " ElemThermi::SM_charge_hydro_E (DdlElement & ...nbvec=  " << ex.giB_t->NbVecteur()
                    << " dim = " << ParaGlob::Dimension();
               Sortie(1);     
              }
           #endif
           // la normale vaut le produit vectoriel des 2 premiers vecteurs
           CoordonneeB normale = Util::ProdVec_coorB( (*ex.giB_t)(1), (*ex.giB_t)(2));
           // calcul de la normale normée et pondérée de la pression
           normale.Normer();  normale *= pression;
           int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3    
           for (int ne =1; ne<= nbne; ne++) 
             for (int i=1;i<= dimf;i++,ix++)
                SM(ix) +=  taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_t));   
           }// fin test : point dans le liquide ou pas
         } // fin boucle sur les points d'intégrations  
      // retour du second membre  
      return SM;  
    };
    
    	    
// idem SMR_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 ElemThermi::SMR_charge_hydro_I (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)
   {  
      #ifdef MISE_AU_POINT	 	 
      // axisymétrie: pour l'instant non testé a priori avec les forces hydro
      if(ParaGlob::AxiSymetrie())
          { cout << "\n erreur, les charges hydro ne sont pas testees avec les elements axisymetriques"
               << " ElemThermi::SM_charge_hydro_I ( ...";
            Sortie(1);     
           };
      #endif
   
      // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
      // *** sinon c'est intractable et long
      // dans le cas ou le poid volumique est nul, on retourne sans calcul
      if (poidvol < ConstMath::trespetit)
        {  Element::ResRaid el;    
           el.res = (*res_extS)(nSurf); 
           el.raid = (*raid_extS)(nSurf);       
           return el;
         }
      // controle
      bool avec_raid = pa.Var_charge_externe();
      int ni; // compteur globale de point d'integration
      Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
      int nbddl = ddls.NbDdl();
      Vecteur& SM = (*((*res_extS)(nSurf)));  //  "      "          "
      Mat_pleine & KM = (*((*raid_extS)(nSurf))); //  "      "       "
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
       {
        // calcul de la pression au niveau de la position du point d'intégration
        //  1) on récupère la position dans I_a du point d'intégration
        const Coordonnee & M = defS.Position_tdt();
        //  2) calcul de la distance à la surface libre
        Coordonnee MMp= M_liquide - M; double dis = dir_normal_liquide * MMp;
        //  3) calcul de la pression effective uniquement si le point est dans l'eau
        // en fait si le point est hors de l'eau il n'y a pas de pression hydrostatique
        if (dis > 0.)
         { double pression = -dis * poidvol;
           // dans le cas avec raideur on calcul la variation de la pression
           const Tableau <Coordonnee> * d_M_pointe = NULL;
           if (avec_raid)
                 d_M_pointe=&(defS.Der_Posi_tdt());
           // --- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
           // que pour un calcul primaire en implicit 
           const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul); 
           #ifdef MISE_AU_POINT	 	 
            // test pour savoir si l'on a réellement affaire à une surface
            if ((ex.giB_t->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3)) 
              { cout << "\n erreur, il doit y avoir deux vecteurs pour la surface  de dimension 3"
                  << " ElemThermi::SMR_charge_hydro_I  (....  " << ex.giB_t->NbVecteur()
                  << " dim = " << ParaGlob::Dimension();
                Sortie(1);     
               }
           #endif
           // la normale vaut le produit vectoriel des 2 premiers vecteurs
           CoordonneeB normale = Util::ProdVec_coorB( (*ex.giB_tdt)(1), (*ex.giB_tdt)(2));
           // calcul de la variation de la normale
           // 1) variation du produit vectoriel
           Tableau <CoordonneeB > D_pasnormale = 
                 Util::VarProdVect_coorB( (*ex.giB_tdt)(1),(*ex.giB_tdt)(2),(*ex.d_giB_tdt));
           // 2) de la normale        
           Tableau <CoordonneeB> D_normale = Util::VarUnVect_coorB(normale,D_pasnormale,normale.Norme());
           // calcul de la normale normée et pondérée de la pression
           normale.Normer(); normale *= pression;
           // également pondération de la variation de la 
           if (avec_raid)
            { for (int ihi =1;ihi<= nbddl;ihi++)
              { double d_pression = poidvol * (dir_normal_liquide * (*d_M_pointe)(ihi)); 
                // le - devant dis disparait avec le - devant M
                D_normale(ihi) = pression * D_normale(ihi)  + d_pression * normale ;
               }
             }        
           int ix=1; int dimf = 3;   // ne fonctionne qu'en dim 3    
           for (int ne =1; ne<= nbne; ne++) 
             for (int i=1;i<= dimf;i++,ix++)
              {  SM(ix) +=  taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_tdt)); 
                 // dans le cas avec_raideur on calcul la contribution à la raideur
                 if (avec_raid)
                   for (int j =1; j<= nbddl; j++)  
                    KM(ix,j) -= 
                             taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j))
                            + taphi(ni)(ne)* D_normale(j)(i) * (poids(ni) * (*ex.jacobien_tdt)) ;
                } 
           }// fin test : point dans le liquide ou pas
         } // fin boucle sur les points d'intégrations  
                
     // liberation des tenseurs intermediaires
     LibereTenseur();   
    Element::ResRaid el;
    el.res = (*res_extS)(nSurf); 
    el.raid = (*raid_extS)(nSurf);       
    return el;
    };