
// 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 "ElemMeca.h"
#include <iomanip>
#include "ConstMath.h"
#include "Util.h"
#include "Coordonnee2.h"
#include "Coordonnee3.h"
#include "CharUtil.h"
#include "TypeQuelconqueParticulier.h"



// actualisation des ddl et des grandeurs actives de t+dt vers t
// appelé par les classes dérivées      
void ElemMeca::TdtversT_()
	{// on met à jour l'indicateur de premier calcul
	 // s'il y a des sauvegardes de grandeur aux déformations
	 // on ne regarde que le premier élément de tableau, a priori
	 // il y a toujours un pt d'integ et l'organisation est identique pour tous les pt d'integ
	 if (tabSaveDefDon(1) != NULL)
        premier_calcul_meca_impli_expli=false; 
	 // cas des énergies
	 int nbi= tab_energ.Taille();
  for (int ni=1;ni<= nbi; ni++)
     tab_energ_t(ni) = tab_energ(ni);
  E_elem_bulk_t = E_elem_bulk_tdt; // énergie due au bulk viscosity
  energie_totale_t = energie_totale; // les énergies globalisées sur l'élément
  // cas des intégrales de volumes et temps éventuelles
  if (integ_vol_typeQuel != NULL)
   {int taille = integ_vol_typeQuel->Taille();
    Tableau <TypeQuelconque> & tab_integ = *integ_vol_typeQuel; // pour simplifier
    Tableau <TypeQuelconque> & tab_integ_t = *integ_vol_typeQuel_t; // ""
    for (int il =1;il <= taille ;il++)
      {tab_integ_t(il) = tab_integ(il);
//       tab_integ(il).Grandeur_pointee()->InitParDefaut(); // on initialise pour le prochain calcul
      };
   };
  if (integ_vol_t_typeQuel != NULL)
   {int taille = integ_vol_t_typeQuel->Taille();
    Tableau <TypeQuelconque> & tab_integ = *integ_vol_t_typeQuel; // pour simplifier
    Tableau <TypeQuelconque> & tab_integ_t = *integ_vol_t_typeQuel_t; // ""
    for (int il =1;il <= taille ;il++)
      {tab_integ_t(il) = tab_integ(il);
      };
   };
  // force de stabilisation éventuelle
  if (pt_StabMembBiel != NULL)
    pt_StabMembBiel->TdtversT();
	};

// actualisation des ddl et des grandeurs actives de t vers tdt
// appelé par les classes dérivées      
void ElemMeca::TversTdt_()
	{// on met à jour l'indicateur de premier calcul
	 // on considère que si l'on revient en arrière, il vaut mieux re-initialiser les
	 // grandeurs correspondantes au premier_calcul
	 // s'il y a des sauvegardes de grandeur aux déformations
	 if (tabSaveDefDon(1) != NULL)
         premier_calcul_meca_impli_expli=true;  		
	 // cas des énergies
	 int nbi= tab_energ.Taille();
  for (int ni=1;ni<= nbi; ni++)
     tab_energ(ni) = tab_energ_t(ni);
  E_elem_bulk_tdt = E_elem_bulk_t; // énergie due au bulk viscosity
  energie_totale = energie_totale_t; // les énergies globalisées sur l'élément
  // cas des intégrales de volumes et temps éventuelles
// il n'y a rien à faire
//  if (integ_vol_typeQuel != NULL)
//   {int taille = integ_vol_typeQuel->Taille();
//    Tableau <TypeQuelconque> & tab_integ = *integ_vol_typeQuel; // pour simplifier
//    Tableau <TypeQuelconque> & tab_integ_t = *integ_vol_typeQuel_t; // ""
//    for (int il =1;il <= taille ;il++)
//      tab_integ(il) = tab_integ_t(il);
////      tab_integ(il).Grandeur_pointee()->InitParDefaut(); // on initialise pour le prochain calcul
//   };
  if (integ_vol_t_typeQuel != NULL)
   {int taille = integ_vol_t_typeQuel->Taille();
    Tableau <TypeQuelconque> & tab_integ = *integ_vol_t_typeQuel; // pour simplifier
    Tableau <TypeQuelconque> & tab_integ_t = *integ_vol_t_typeQuel_t; // ""
    for (int il =1;il <= taille ;il++)
      {tab_integ(il) = tab_integ_t(il); // on met à jour le cumul
      };
   };
  // force de stabilisation éventuelle
  if (pt_StabMembBiel != NULL)
    pt_StabMembBiel->TdtversT();

	};

// calcul du résidu et de la matrice de raideur pour le calcul d'erreur
// cas d'une intégration suivant une seule liste
void ElemMeca::SigmaAuNoeud_ResRaid(const int nbne,const Tableau <Vecteur>& taphi
                ,const Vecteur& poids,Tableau <Vecteur *>&  resErr,Mat_pleine&  raidErr
                ,const Tableau <Vecteur>& taphiEr,const Vecteur& poidsEr)
   {  int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
      // dimension des tenseurs
//      int dim = lesPtIntegMecaInterne->DimTens();
      // inialisation du second membre et de la raideur
      int nbSM = resErr.Taille(); // nombre de second membre
      for (int j =1; j<= nbSM; j++)  //  boucle sur les seconds membres
           (*resErr(j)).Zero();
      raidErr.Zero();
      // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit
      // être compatible avec l'ordre des ddl
      Tableau2 <int> ordre = OrdreContrainte(nbSM);
      // création d'un tenseur au dimension absolu pour le calcul des contraintes
      // dans la base absolue
      TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ;
       // ====== calcul du second membre ======= 
      int ni; // compteur globale de point d'integration
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
       {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
        // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
        // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
        const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
        // passage dans le repère absolu du tenseur contrainte final
        sigHH = (*(ptIntegMeca.SigHH_t())).BaseAbsolue(sigHH,*ex.giB_t);
        for (int ne =1; ne<= nbne; ne++)  // 1ere boucle sur les noeuds
          {  // les résidus : mais il faut suivre l'ordre de l'enregistrement des ddl
             for (int itot = 1; itot<= nbSM; itot++)
               { int ix = (int) (ordre(itot,1)); int iy = (int) (ordre(itot,2));
                (*resErr(itot))(ne) += taphi(ni)(ne)*sigHH(ix,iy) * (poids(ni) * (*ex.jacobien_t));
               };
          };
       };
       // ====== calcul de la raideur c'est à dire du hessien ========  
       // boucle sur les pt d'integ spécifiques à l'erreur
      for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
       {
         // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
         // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
         const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
         for (int ne =1; ne<= nbne; ne++)  // 1ere boucle sur les noeuds
            for (int me =1; me<= nbne; me++)  // 2ere boucle sur les noeuds
                    raidErr(ne,me) += taphiEr(ni)(ne) * taphiEr(ni)(me) * poidsEr(ni) * (*ex.jacobien_t);
       };
      // liberation des tenseurs intermediaires
      TenseurHH * ptsigHH = &sigHH; delete ptsigHH;
      LibereTenseur();
      // calcul de l'erreur relative
              
   };

// calcul de l'erreur sur l'élément. Ce calcul n'est disponible
// qu'une fois la remontée aux contraintes effectuées sinon aucune
// action. En retour la valeur de l'erreur sur l'élément
// = 1 : erreur = (int (delta sigma):(delta sigma) dv)/(int sigma:sigma dv)
// le numerateur et le denominateur sont tel que :
// errElemRelative = numerateur / denominateur , si denominateur different de 0
// sinon denominateur = numerateur si numerateur est different de 0, sinon
// tous sont nuls mais on n'effectue pas la division , les autres variables sont spécifiques
 // a l'element.
void ElemMeca::Cal_ErrElem(int type,double& errElemRelative,double& numerateur
                , double& denominateur,const int nbne,const Tableau <Vecteur>& taphi
                ,const Vecteur& poids,const Tableau <Vecteur>& taphiEr,const Vecteur& poidsEr) 
  { int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
    // dimension des tenseurs
//    int dim = lesPtIntegMecaInterne->DimTens();
    bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
    double domaine_integration=0.; // volume ou surface on longueur
    switch (type)
     {case 1 : // cas du calcul aux moindres carrés

         {// création d'un tenseur au dimension absolu pour le calcul des contraintes
          // dans la base absolue, on le choisit HB  pour le double produit contracté
          // mais en absolu la variance n'a pas d'importance
          TenseurHB& sigHB = *(NevezTenseurHB(dimAbsolue)) ;
          // idem au point d'intégration et un tenseur nul pour l'initialisation
          TenseurHB& signiHB = *(NevezTenseurHB(dimAbsolue)) ;
          TenseurHB& sigiHB = *(NevezTenseurHB(dimAbsolue)) ; // tenseur de travail
          TenseurHB& nulHB = *(NevezTenseurHB(dimAbsolue)) ;
          // ====== calcul des termes de l'erreur =======
          // ---- tout d'abord on parcourt les points d'intégration de la mécanique
          numerateur = 0.;  // initialisation
          double numerateur_bis = 0.;
          denominateur = 0.;// intégrale double de la contrainte élément fini
          int ni; // compteur globale de point d'integration
          for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
           {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
            // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
            // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
            const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
            // passage dans le repère absolu du tenseur contrainte initiale
            // car les contraintes aux noeuds on été calculées et stockées en absolue
            sigHB = (*(ptIntegMeca.SigHH_t())).BaseAbsolue(sigHB,*ex.giB_t);
////---debug
//cout << "\n debug ElemMeca::Cal_ErrElem: sigHB= ";
// sigHB.Ecriture(cout);
////---fin debug
            // calcul du denominateur
            double integ_sigHB_carre_du_pti = (sigHB && sigHB) * (poids(ni) * (*ex.jacobien_t));
            denominateur += integ_sigHB_carre_du_pti;
            
            // calcul de la premiere partie du numerateur, celle qui dépend des points d'intégration
            // mécanique. 
               // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, 
            signiHB = nulHB; // initialisation   
            for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
               signiHB += taphi(ni)(ne) * ((tab_noeud(ne))->Contrainte(sigiHB));
               // 2) intégrale de la partie dépendante de ni
            numerateur += integ_sigHB_carre_du_pti - ((signiHB && sigHB)+(sigHB && signiHB)) * (poids(ni) * (*ex.jacobien_t))  ;
            numerateur_bis += (signiHB-sigHB)&&(signiHB-sigHB)* (poids(ni) * (*ex.jacobien_t));
           };
          // ---- on parcourt maintenant les points d'intégration pour le calcul d'erreur
          // ---- ce qui permet de calculer la deuxième partie du numérateur
          // boucle sur les pt d'integ spécifiques à l'erreur
          double denominateur_continu = 0.;// intégrale double de la contrainte continue
          domaine_integration=0.; // init
          for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
           {
            // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
            // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
            const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
               // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds, 
            signiHB = nulHB; // initialisation   
            for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
               signiHB += taphiEr(ni)(ne)*(tab_noeud(ne))->Contrainte(sigiHB);
               // 2) intégrale de la partie dépendant de ni
////---debug
//cout << "\n debug ElemMeca::Cal_ErrElem: signiHB= ";
// signiHB.Ecriture(cout);
//cout << "\n debug ElemMeca::Cal_ErrElem: (tab_noeud(ne))->Contrainte(sigiHB)= ";
//for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
//   {cout <<"\n";(tab_noeud(ne))->Contrainte(sigiHB).Ecriture(cout);}
//
////---fin debug
            double integ_signiHB_carre_du_pti = (signiHB && signiHB) * poidsEr(ni) * (*ex.jacobien_t);
            denominateur_continu += integ_signiHB_carre_du_pti;
            numerateur += integ_signiHB_carre_du_pti;
               // 3) intégrale du domaine d'intégration
            domaine_integration += poidsEr(ni) * (*ex.jacobien_t);
           };
          // liberation des tenseurs intermediaires
          TenseurHB * ptsigHB = &sigHB; delete ptsigHB;
          ptsigHB = &signiHB; delete ptsigHB;ptsigHB = &sigiHB; delete ptsigHB;
          ptsigHB = &nulHB; delete ptsigHB;
          LibereTenseur();
          // enregistrement de l'erreur pour l'élément
          if (sigErreur == NULL) {sigErreur = new double;sigErreur_relative = new double;};
          // comme le calcul précédent correspond à l'intégrale d'un carré,
          // donc c'est analogue à "une grandeur moyenne au carré " * le volume
          //on prend la racine carrée de:
          //  (" grandeur moyenne au carré " * le volume) / le volume c-a-d |grandeur moyenne|
          // qui est donc pondèré par le domaine ce qui nous donne une grandeur analogue à
          // un delta contrainte
          
          // on sauvegarde avant
 //         double numerateur_initale = numerateur;
          double denominateur_initiale = denominateur;
// modif 28 nov 2024
//          *sigErreur = sqrt(Dabs(numerateur/domaine_integration));
          *sigErreur = sqrt(numerateur_bis/domaine_integration);
          // *sigErreur a la dimension de l'écart de contrainte, il s'agit
          // de l'écart moyenné sur tout l'élément

          // pour calculer une erreur relative on commence par calculer une grandeur moyenne
          // qui sera ici analogue à une contrainte
// modif 28 nov 2024
//          denominateur = sqrt(Dabs(denominateur_initiale)/domaine_integration);
          denominateur = sqrt(denominateur_initiale/domaine_integration);
          //denominateur a la dimention de la contrainte moyenne sur tout l'élément

// cf. pense b^te : la suite doit-êtr supprimé ... ok fait
// sigErreur_relative est calculée dans LesMaillages::ErreurSurChaqueElement
// après le calcul de toutes les erreurs, donc on ne fait rien ici
// mais on laisse en ayant en tête que c'est une erreur relative valable uniquement
// avant l'appel de LesMaillages::ErreurSurChaqueElement

          // on essaie de gérer le cas des contraintes petites
          if (denominateur < ConstMath::unpeupetit)
            // on passe en erreur absolue si le dénominateur est trop petit
            { *sigErreur_relative = *sigErreur; }
          else
            {*sigErreur_relative = *sigErreur/(denominateur);};

          numerateur = *sigErreur;
          
//////------ debug
//if (numerateur > 90)
//  {cout << "\n debug ElemMeca::Cal_ErrElem: *sigErreur= "<< *sigErreur
//       << " *sigErreur_relative = " << *sigErreur_relative
//       << "\n num élément: " <<   num_elt
//       << flush;
//   // création d'un tenseur au dimension absolu pour le calcul des contraintes
//   // dans la base absolue, on le choisit HB  pour le double produit contracté
//   // mais en absolu la variance n'a pas d'importance
//   TenseurHB& sigHB = *(NevezTenseurHB(dimAbsolue)) ;
//   TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ; // pour voir si identique à sigHB
//      // idem au point d'intégration et un tenseur nul pour l'initialisation
//   TenseurHB& signiHB = *(NevezTenseurHB(dimAbsolue)) ;
//   TenseurHB& sigiHB = *(NevezTenseurHB(dimAbsolue)) ; // tenseur de travail
//   TenseurHB& nulHB = *(NevezTenseurHB(dimAbsolue)) ;
//   // ====== calcul des termes de l'erreur =======
//   // ---- tout d'abord on parcourt les points d'intégration de la mécanique
//   double numerateur = 0.;  // initialisation
//   double numerateur_bis = 0.;
//   double denominateur = 0.;// intégrale double de la contrainte élément fini
//   int ni; // compteur globale de point d'integration
//   for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
//    {PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
//     // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//     // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//     const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
//     // passage dans le repère absolu du tenseur contrainte initiale
//     // car les contraintes aux noeuds on été calculées et stockées en absolue
//     sigHB = (*(ptIntegMeca.SigHH_t())).BaseAbsolue(sigHB,*ex.giB_t);
//     (ptIntegMeca.SigHH_t())->BaseAbsolue(sigHH,*ex.giB_t);
////---debug
//cout << "\n debug ElemMeca::Cal_ErrElem: num_pti: " << ni
//     << "\n contrainte initiale sigHB= ";
// sigHB.Ecriture(cout);
//// cout << "\n et la version en HH : sigHH= ";
////sigHH.Ecriture(cout);
////---fin debug
//      // calcul du denominateur
//      double integ_sigHB_carre_du_pti = (sigHB && sigHB) * (poids(ni) * (*ex.jacobien_t));
//      denominateur += integ_sigHB_carre_du_pti;
//
//      // calcul de la premiere partie du numerateur, celle qui dépend des points d'intégration
//      // mécanique.
//         // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds,
//      signiHB = nulHB; // initialisation
//      for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
//         signiHB += taphi(ni)(ne) * ((tab_noeud(ne))->Contrainte(sigiHB));
//         // 2) intégrale de la partie dépendante de ni
//      numerateur += integ_sigHB_carre_du_pti - ((signiHB && sigHB)+(sigHB && signiHB)) * (poids(ni) * (*ex.jacobien_t))  ;
//      numerateur_bis += (signiHB-sigHB)&&(signiHB-sigHB)* (poids(ni) * (*ex.jacobien_t));
//     };
//   // ---- on parcourt maintenant les points d'intégration pour le calcul d'erreur
//   // ---- ce qui permet de calculer la deuxième partie du numérateur
//   // boucle sur les pt d'integ spécifiques à l'erreur
//   double denominateur_continu = 0.;// intégrale double de la contrainte continue
//   domaine_integration=0.; // init
//   for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
//    {
//     // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
//     // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
//     const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
//        // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds,
//     signiHB = nulHB; // initialisation
//     for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
//        signiHB += taphiEr(ni)(ne)*(tab_noeud(ne))->Contrainte(sigiHB);
//        // 2) intégrale de la partie dépendant de ni
////---debug
//cout << "\n\n  debug ElemMeca::Cal_ErrElem: contraintes aux noeuds ramené au pti: "<<ni
//     << " signiHB=  ";
// signiHB.Ecriture(cout);
//cout << "\n debug ElemMeca::Cal_ErrElem: (tab_noeud(ne))->Contrainte(sigiHB)= ";
//for (int ne =1; ne<= nbne; ne++)  //  boucle sur les noeuds
//   {cout <<"\n";(tab_noeud(ne))->Contrainte(sigiHB).Ecriture(cout);}
//
////---fin debug
//      double integ_signiHB_carre_du_pti = (signiHB && signiHB) * poidsEr(ni) * (*ex.jacobien_t);
//      denominateur_continu += integ_signiHB_carre_du_pti;
//      numerateur += integ_signiHB_carre_du_pti;
//         // 3) intégrale du domaine d'intégration
//      domaine_integration += poidsEr(ni) * (*ex.jacobien_t);
//     };
//cout << "\n debug ElemMeca::Cal_ErrElem: domaine_integration= "
//     << domaine_integration << flush;
//   // liberation des tenseurs intermediaires
//   TenseurHB * ptsigHB = &sigHB; delete ptsigHB;
//   ptsigHB = &signiHB; delete ptsigHB;ptsigHB = &sigiHB; delete ptsigHB;
//   ptsigHB = &nulHB; delete ptsigHB;
//   TenseurHH * ptsigHH = &sigHH; delete ptsigHH;
//   LibereTenseur();
//   }
////------ fin debug

//          *sigErreur = Dabs(numerateur);
//
//          numerateur = sqrt(Dabs(numerateur/domaine_integration));
//          *sigErreur_relative = numerateur;
////          numerateur = sqrt(Dabs(numerateur)/domaine_integration)*domaine_integration;
//          *sigErreur = numerateur;
          
          
////---essai sur des grandeurs relatives, mais en fait ce n'est pas une bonne idée, car ici c'est
// la différence totale qui nous intéresse, et en relativ cela donne n'importe quoi
//          *sigErreur = numerateur;
            {// on commence par calculer le maxi des contraintes moyennes au carré, soient via le calcul élément fini initial
             // soit via la contrainte continu, et on évite d'avoir 0
////             double sig_moy_2 = DabsMaX(denominateur_continu,denominateur,ConstMath::trespetit);
//             double sig_moy_2 = DabsMaX(denominateur,ConstMath::trespetit);
//////---debug
//cout << "\n debug ElemMeca::Cal_ErrElem: "
//     << "\n Int_sig2_continu=" << denominateur_continu <<", Int_sig2_EF=" << denominateur
//     <<", rapport="<<denominateur/denominateur_continu
//     << "\n Int_(delta_sig)^2 opti="<<numerateur<<", Int_(delta_sig)^2 brut="<<numerateur_bis
//     << ", rapport="<< numerateur/numerateur_bis
//     << "\n estimateur= Int_(delta_sig)^2 opti/Int_sig2_EF= "<<numerateur/denominateur<<endl;
//     denominateur = sqrt(denominateur/domaine_integration)*domaine_integration;
//     numerateur = sqrt(numerateur/domaine_integration)*domaine_integration;
//     << "\n denominateur_continu=" << denominateur_continu <<", denominateur=" << denominateur
//     << ", numerateur="<<numerateur<<", domaine_integration="<<domaine_integration
//     << ", sig_moy="<< sqrt(Dabs(numerateur)/domaine_integration) <<endl;
//---fin debug
////---fin essai sur des grandeurs relatives
             
            }
          break;
       }
      default :
       {cout << "\nErreur : valeur incorrecte du type de calcul d'erreur !\n";
        cout << "ElemMeca::ErreurElement(int type  , type = \n" << type;
        Sortie(1);
       }
     };
    // --- calcul de l'erreur relative
    
//    if (denominateur <= ConstMath::trespetit)
//       // cas d'un champ de contraintes nulles initialement
//       {if (numerateur <= ConstMath::trespetit)
//           // cas également du numérateur nul
//            errElemRelative = 0.;
//        else
//           // on fait denominateur = numérateur -> erreur relative = 1
//            errElemRelative = 1./domaine_integration;
//       }
//    else
//      // cas du dénominateur non nul
//      { errElemRelative = numerateur / denominateur;}//(denominateur*domaine_integration);};
    errElemRelative = *sigErreur_relative;
  };
	    

// calcul de l'erreur aux noeuds. Contrairement au cas des contraintes
// seul le résidu est calculé. Cas d'une intégration suivant une liste
void ElemMeca::Cal_ErrAuxNoeuds(const int nbne, const Tableau <Vecteur>& taphi, 
                const Vecteur& poids,Tableau <Vecteur *>&  resErr )
   {  //int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
      // inialisation du second membre, on utilise le premier vecteur uniquement
      (*resErr(1)).Zero();
       // ====== calcul du second membre ======= 
      int ni; // compteur globale de point d'integration
      double volume = 0.; // le volume de l'élément
      bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
      for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
       {
       // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
       // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
        const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
        for (int ne =1; ne<= nbne; ne++)  // 1ere boucle sur les noeuds
//             (*resErr(1))(ne) += taphi(ni)(ne)*(*sigErreur) * (poids(ni) * (*ex.jacobien_t));
           (*resErr(1))(ne) += taphi(ni)(ne)*(*sigErreur_relative) * (poids(ni) * (*ex.jacobien_t));
//       // calcul du volume
//       volume += (poids(ni) * (*ex.jacobien_t));
       } 
      // on relativise par rapport au volume, car initialement sigErreur représente l'erreur
      // totale sur tout l'élément. Donc on divise par le volume pour retrouver après résolution
      // une valeur au noeud qui représente une valeur ponctuelle et non une valeur qui
      // qui est relative à un volume
// modif 28 nov 2024 : (*sigErreur_relative) est déjà une grandeur sans dimension
//      *resErr(1) /= volume;
      // ensuite on prend la racine carrée pour que ce soit homogène à une contrainte
// non car 1) il y a une erreur dans la boucle on reste sur la première valeur !! qui est la seule
// à etre changée
//         2) les grandeurs manipulées sont des rapports de delta_contrainte / contraintes maxi
// il ne faut pas prendre de racine carré

//      int tail = (*resErr(1)).Taille();
//      Vecteur& resErr_1 = *resErr(1); // pour simplifier
//      for (int i=1;i<= tail;i++)
//        {if (resErr_1(1) >= 0.)
//           resErr_1(1) = sqrt(resErr_1(1));
//         else
//           { if (ParaGlob::NiveauImpression() > 3)
//               {cout << "\n *** warning : bizarre la composante "<<i<<" du vecteur erreur "
//                    << " est negative, on retient sa valeur absolue " << endl;
//// --- debug
//cout << "\n ElemMeca::Cal_ErrAuxNoeuds ";
//cout << "\n volume = "<< volume;
//
//// fin debug
//
//               }
//             if (ParaGlob::NiveauImpression() > 5)
//               cout << "\n ElemMeca::Cal_ErrAuxNoeuds ";
//             resErr_1(1) = sqrt(-resErr_1(1));
//            };
//        };

////---debug
// cout << "\n debug ElemMeca::Cal_ErrAuxNoeuds ";
// resErr_1.Affiche();
// cout << endl;
////---fin debug
    
            
//      *resErr(1) = -sigErreur;
      LibereTenseur();
      // calcul de l'erreur relative
              
   };
        
    
// ajout des ddl relatif aux contraintes pour les noeuds de l'élément
void ElemMeca::Ad_ddl_Sigma(const DdlElement&  tab_ddlErr)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) // pour chaque noeud
      { // création du tableau de ddl
        int tab_ddlErr_Taille = tab_ddlErr(ne).tb.Taille(); // nb de ddl du noeud meca
        Tableau <Ddl> ta(tab_ddlErr_Taille);
        for (int i =1; i<= tab_ddlErr_Taille; i++)
	        ta(i) = Ddl (tab_ddlErr(ne).tb(i),0.,LIBRE);
	       // ajout du tableau dans le noeud    
        tab_noeud(ne)->PlusTabDdl(ta);
       } 
   };
       
// inactive les ddl du problème primaire de mécanique
void ElemMeca::Inact_ddl_primaire(DdlElement&  tab_ddl)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_hors_service(tab_ddl(ne).tb);
   };
// active les ddl du problème primaire de mécanique
void ElemMeca::Act_ddl_primaire(DdlElement&  tab_ddl)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(tab_ddl(ne).tb);
   };
// inactive les ddl du problème de recherche d'erreur : les contraintes
void ElemMeca::Inact_ddl_Sigma(DdlElement&  tab_ddlErr)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_hors_service(tab_ddlErr(ne).tb);
   };
// active les ddl du problème de recherche d'erreur : les contraintes
void ElemMeca::Act_ddl_Sigma(DdlElement&  tab_ddlErr)
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(tab_ddlErr(ne).tb);
   };
// active le premier ddl du problème de recherche d'erreur : SIGMA11
void ElemMeca::Act_premier_ddl_Sigma()
  { int nbne = tab_noeud.Taille();
    for (int ne=1; ne<= nbne; ne++) 
        tab_noeud(ne)->Met_en_service(SIG11);
   };
    
// retourne un tableau de ddl element, correspondant à la 
// composante de sigma -> SIG11, pour chaque noeud qui contiend
// des ddl de contrainte
// -> utilisé pour l'assemblage de la raideur d'erreur
DdlElement&  ElemMeca::Tableau_de_Sig1() const
  {       cout << "\n erreur, fonction non defini pour cette element "
               << "\n ElemMeca::Tableau_de_Sig1()" << endl;
          Sortie(1); 
          DdlElement * toto = new DdlElement();
          return *toto;    
   };
   
// lecture des contraintes sur le flot d'entrée
void ElemMeca::LectureDesContraintes(bool cas,UtilLecture * entreePrinc,Tableau <TenseurHH *>& tabSigHH)
   { // dimensionnement de la metrique identique au cas d'un calcul explicite
       if( cas)
        { Tableau<Enum_variable_metrique> tab(7);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         met->PlusInitVariables(tab) ;
        };  
     int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
     // dimension des tenseurs
     //int dim = (*(tabSigHH(1))).Dimension();
	    // récup du nombre de composantes du tenseur
	    int nbcomposante= ParaGlob::NbCompTens();
     // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit
     // être compatible avec l'ordre des ddl
     Tableau2  <int> ordre = OrdreContrainte(nbcomposante);
     // création d'un tenseur au dimension absolu pour récupérer les contraintes
     // dans la base absolue
     TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ;
     int ni; // compteur globale de point d'integration
     bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
     for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
       { // calcul des éléments de la métrique,  on utilise le même calcul
         // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
         const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
         for (int i=1;i<= nbcomposante;i++)
           // récupération des coordonnées du tenseur en absolu
            *(entreePrinc->entree)  >> sigHH.Coor(ordre(i,1),ordre(i,2));
         // passage dans le repère local du tenseur contrainte final
         (*(tabSigHH(ni))) = sigHH.Baselocale((*(tabSigHH(ni))),*ex.giH_t);
       };
   };
   

// retour des contraintes en absolu retour true si elle existe sinon false
void  ElemMeca::ContraintesEnAbsolues
          (bool cas,Tableau <TenseurHH *>& tabSigHH,Tableau <Vecteur>& tabSigAbs)
   { // dimensionnement de la metrique identique au cas d'un calcul explicite
       if( cas)
        { Tableau<Enum_variable_metrique> tab(7);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         met->PlusInitVariables(tab) ;
        };  
     int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
     // dimension des tenseurs
     //int dim = (*(tabSigHH(1))).Dimension();
     // récup du nombre de composantes du tenseur
     int nbcomposante= ParaGlob::NbCompTens(); 
     // redimensionnement éventuel du tableau de sortie 
     int nbi = tabSigHH.Taille();
     if (tabSigAbs.Taille() != nbi)
       tabSigAbs.Change_taille(nbi);
     for (int ni=1;ni<= nbi;ni++)
	     if ( tabSigAbs(ni).Taille() != nbcomposante)
	         tabSigAbs(ni).Change_taille(nbcomposante);
     // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit
     // être compatible avec l'ordre des ddl
     Tableau2 <int> ordre = OrdreContrainte(nbcomposante);
     // création d'un tenseur au dimension absolu pour récupérer les contraintes
     // dans la base absolue
     TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ;
     int ni; // compteur globale de point d'integration
     bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
                                  // donc il faut calculer tous les éléments de la métrique
     for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
        { // calcul des éléments de la métrique,  on utilise le même calcul
         // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
         const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);        
         // passage dans le repère global du tenseur contrainte local
         sigHH = (*(tabSigHH(ni))).BaseAbsolue(sigHH,*ex.giB_t);
         // remplissage du vecteur pour le retour d'info
         for (int i=1;i<= nbcomposante;i++)
          tabSigAbs(ni)(i) = sigHH(ordre(i,1),ordre(i,2));  
        };
   };
   
// ---------------- lecture écriture dans base info  ----------------
	//  programmes utilisés par les classes dérivées
	
	   // cas donne le niveau de la récupération
       // = 1 : on récupère tout
       // = 2 : on récupère uniquement les données variables (supposées comme telles)
void ElemMeca::Lecture_bas_inf 
  (istream& ent,const Tableau<Noeud  *> * tabMaillageNoeud,const int cas) 
{ // appel de la routine d'élément
  Element::Lect_bas_inf_element(ent,tabMaillageNoeud,cas);
  switch (cas)
   { case 1 : // ------- on récupère tout -------------------------
      { string toto,nom;
        // récup de la masse volumique
        ent >> toto >> masse_volumique ; 
        // données thermique
        ent >> toto >> dilatation;
        // blocage éventuelle d'hourglass
        ent >> toto >>  nom;
        type_stabHourglass=Id_Nom_StabHourglass(nom.c_str());
        // blocage éventuel transversal
        ent >> nom;
        if (nom == "stabilisation_transversale:")
         { pt_StabMembBiel = new StabMembBiel;
           pt_StabMembBiel->Lecture_bas_inf(ent,cas);
//         string nom1; string nom_fonction;
//           ent >> nom1;
//           double valeur = 0.;
//           if (nom1=="par_valeur") {ent >> valeur;}
//           else {ent >> nom_fonction;};
//           // maintenant il faut mettre à jour le conteneur
//           if ((pt_StabMembBiel == NULL)&&(nom1=="par_valeur"))
//             {pt_StabMembBiel=new StabMembBiel(valeur,NULL);}
//           else
//             {pt_StabMembBiel=new StabMembBiel(valeur,&nom_fonction);};
//           // ensuite il faudra définir la fonction de stabilisation
//           // au moment de compléter l'élément
         }
        else // cas sans stabilisation
         {if (pt_StabMembBiel != NULL)
                {delete pt_StabMembBiel; pt_StabMembBiel=NULL;}
         };
        break;
      }
     case  2 : // ----------- lecture uniquement de se qui varie --------------------
      {break;
      }
     default :
      { cout << "\nErreur : valeur incorrecte du type de lecture !\n";
        cout << "ElemMeca::Lecture_bas_inf (istream& ent,const int cas)"
             << " cas= " << cas << endl;
        Sortie(1);
      };
   };
  // ------ lecture dans tous les cas -------
  // résultat d'erreur
  string toto;
  ent >> toto;
  if (toto == "erreur_de_contrainte")
    { if (sigErreur == NULL)
       sigErreur = new double;
      ent  >> (*sigErreur) ;
    };
  // données particulière pour les lois de comportement mécanique
  int tabSaveDonTaille = tabSaveDon.Taille();
  if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) ent >> toto; 
  int num ; // numéro du pt d'integ, n'est pas vraiment utilisé mais c'est mieux que de lire un string
  for (int i=1; i<= tabSaveDonTaille; i++)
    if (tabSaveDon(i) != NULL) 
     { ent >> toto >> num; 
       tabSaveDon(i)->Lecture_base_info(ent,cas);
     };
  // données particulière pour les lois de comportement thermo physique
  int tabSaveTPTaille = tabSaveTP.Taille();
  if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) ent >> toto; 
  for (int i=1; i<= tabSaveTPTaille; i++)
    if (tabSaveTP(i) != NULL) 
     { ent >> toto >> num; 
       tabSaveTP(i)->Lecture_base_info(ent,cas);
     };
  // données particulière pour la déformation mécanique
  int tabSaveDefDonTaille = tabSaveDefDon.Taille();
  if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL))  ent >> toto;
  for (int i=1; i<= tabSaveDefDonTaille; i++)
     if (tabSaveDefDon(i) != NULL) 
     { ent >> toto >> num; 
       tabSaveDefDon(i)->Lecture_base_info(ent,cas);
     };
  // ---- lecture des énergies
  // énergie totale
  ent >> toto >>  energie_totale_t ;
  // par point d'intégration
  int energ_taille = tab_energ_t.Taille();
  ent >> toto;
  for (int i=1;i<= energ_taille;i++) ent >> tab_energ_t(i) ;
  // énergie et puissance éventuelle de la partie bulk viscosity
  ent >> toto >> E_elem_bulk_t >> toto >> P_elem_bulk;
  E_elem_bulk_tdt = E_elem_bulk_t;

  // blocage éventuel transversal
  {int taille = 0;
   ent >> toto >> taille;
   if (taille != 0)
    {// cas avec stabilisation enregistrée
     // on lit systématiquement le cas 2
     pt_StabMembBiel->Lecture_bas_inf(ent,2) ;
//     // on considère qu'il faut la récupérer
//     pt_StabMembBiel->Change_nb_pti(taille);
//     for (int i=1;i<= taille;i++)
//        {ent >> pt_StabMembBiel->FF_StabMembBiel(i);
//         pt_StabMembBiel->FF_StabMembBiel_t(i)=pt_StabMembBiel->FF_StabMembBiel(i);
//        };
    };
   // dans le cas où la taille est nulle, cela signifie qu'il n'y avait pas de stabilisation
   // auparavant, on laisse en l'état, car l'utilisateur a peut-être déjà définit une stabilisation
  };
  
//  if (pt_StabMembBiel == NULL)
//    {// pas de stabilisation active, on passe la grandeur
//     ent >> toto >> toto;
//    }
//  else
//    {// on a une stabilisation enregistrée, on lit la stabilisation
//     int taille=0; // init
//     ent >> toto >> taille;
//     pt_StabMembBiel->Change_nb_pti(taille);
//     for (int i=1;i<= taille;i++)
//       {ent >> pt_StabMembBiel->FF_StabMembBiel(i);
//        pt_StabMembBiel->FF_StabMembBiel_t(i)=pt_StabMembBiel->FF_StabMembBiel(i);
//       };
//    };

  // les forces externes éventuelles
  int indic=0;
  ent >> toto >> indic;
  if (indic == 0)
   {if (lesChargeExtSurEle != NULL)
     delete lesChargeExtSurEle;
   }
  else
   {if (lesChargeExtSurEle == NULL)
     {lesChargeExtSurEle = new LesChargeExtSurElement;};
    lesChargeExtSurEle->Lecture_base_info(ent,cas);
   };
  
};
       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void ElemMeca::Ecriture_bas_inf(ostream& sort,const int cas)
 { // appel de la routine d'élément
   Element::Ecri_bas_inf_element(sort,cas);
   // en fait ici on sauvegarde la même chose dans tous les cas, par contre la sortie
   // totale est documentée.
   switch (cas)
    { case 1 : // ------- on sauvegarde tout -------------------------
       { // écriture de la masse volumique, 
         sort << "masse_volumique " << masse_volumique <<" " << "\n";
         // données thermique
         sort << "dilatation_thermique " << dilatation << " ";
         // blocage éventuel d'hourglass
         sort << "\n hourglass: " << Nom_StabHourglass(type_stabHourglass) << " ";
         // blocage éventuel transversal
         if (pt_StabMembBiel == NULL)
           {sort << "\n pas_de_stabilisation_transversale ";}
         else
           {sort << "\n stabilisation_transversale: ";
            pt_StabMembBiel->Ecriture_bas_inf(sort,cas); 
           };
         break;
       }
      case  2 : // ----------- sauvegarde uniquement de se qui varie --------------------
       { 	   break;
       }
      default :
        { cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
          cout << "ElemMeca::Ecriture_bas_inf(ostream& sort,const int cas)"
               << " cas= " << cas << endl;
          Sortie(1);
        };
    };
   // --- informations à sortir dans tous les cas ----
   // résultat d'erreur
   if (sigErreur != NULL)
           sort << "erreur_de_contrainte " << (*sigErreur) <<" " << "\n";
   else
           sort << "pas_d'erreur_contrainte \n";   
   // données particulière pour les lois de comportement mécaniques
   int tabSaveDonTaille = tabSaveDon.Taille();
   if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) sort << "\n\n +-+-+-data_spec_loi_comp_meca-+-+-+ ";
   for (int i=1; i<= tabSaveDonTaille; i++)
      if (tabSaveDon(i) != NULL) 
         { if (i==1) {sort << "\n...pt_int= " <<i << " ";} else {sort << "\n...pt_int= " <<i << " ";};
          tabSaveDon(i)->Ecriture_base_info(sort,cas);}
   // données particulière pour les lois de comportement thermo physiques
   int tabSaveTPTaille = tabSaveTP.Taille();
   if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) sort << "\n\n +-+-+-data_spec_loi_comp_TP-+-+-+ ";
   for (int i=1; i<= tabSaveTPTaille; i++)
      if (tabSaveTP(i) != NULL) 
         { if (i==1) {sort << "\n...pt_int= " <<i << " ";} else {sort << "\n...pt_int= " <<i << " ";};
          tabSaveTP(i)->Ecriture_base_info(sort,cas);}
   // données particulière pour la déformation mécanique
   int tabSaveDefDonTaille = tabSaveDefDon.Taille();
   if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL)) sort << "\n\n +-+-+-data_spec_def-+-+-+ ";
   for (int i=1; i<= tabSaveDefDonTaille; i++)
      if (tabSaveDefDon(i) != NULL) 
         { if (i==1) {sort << "\n...pt_int= " <<i << " ";} else {sort << "\n...pt_int= " <<i << " ";};
           tabSaveDefDon(i)->Ecriture_base_info(sort,cas);
         };
   // --- sortie des énergies
   // énergie totale
   sort << "\n\n E_elem_totale: " << energie_totale << "  ";
   // par point d'intégration
   int energ_taille = tab_energ_t.Taille();
   sort << "\n   par_pti: ";
   for (int i=1;i<= energ_taille;i++) sort << tab_energ_t(i) << " " ;
   // énergie et puissance éventuelle de la partie bulk viscosity
   sort << "\n E_el_bulk= " << E_elem_bulk_t << " P_el_bulk= " << P_elem_bulk;
   // blocage éventuel transversal
   if (pt_StabMembBiel == NULL)
     {sort << "\n StabTrans 0 ";}
   else
     {sort << "\n StabTrans 1 ";
      // on sort systématiquement le cas 2
      pt_StabMembBiel->Ecriture_bas_inf(sort,2) ;
     };
  
   // les forces externes éventuelles
   if (lesChargeExtSurEle != NULL)
     {sort << "\n lesChargesExt 1 "; lesChargeExtSurEle->Ecriture_base_info(sort,cas);}
   else {sort << "\n lesChargesExt 0 ";};
 };
 
// calcul  de la longueur d'arrête de l'élément minimal
// divisé par la célérité la plus rapide dans le matériau
// appelé par les classes dérivées
// nb_noeud : =0 indique que l'on utilise tous les noeuds du tableau de noeuds
//        = un nombre > 0, indique le nombre de noeuds à utiliser au début du tableau
double ElemMeca::Interne_Long_arrete_mini_sur_c(Enum_dure temps,int nb_noeud)
  { int nbne = tab_noeud.Taille(); // récup du nombre de noeud
    if (nb_noeud != 0)
        nbne = nb_noeud;
    // tout d'abord l'objectif est de déterminer la distance minimum
    // entre les différents noeuds
    // initialisation de la distance entre les deux noeuds
    double dist = ConstMath::tresgrand;
    for (int i=1;i<= nbne;i++)
      // on itère sur les noeuds restants
      switch (temps)
      { case TEMPS_0:
          for (int j=i+1;j<= nbne;j++)
            { double dist_new = (tab_noeud(i)->Coord0() - tab_noeud(j)->Coord0()).Norme();
              if (dist_new < dist) dist = dist_new;
             }
          break;   
        case TEMPS_t:
          for (int j=i+1;j<= nbne;j++)
            { double dist_new = (tab_noeud(i)->Coord1() - tab_noeud(j)->Coord1()).Norme();
              if (dist_new < dist) dist = dist_new;
             }
          break;   
        case TEMPS_tdt:
          for (int j=i+1;j<= nbne;j++)
            { double dist_new = (tab_noeud(i)->Coord2() - tab_noeud(j)->Coord2()).Norme();
              if (dist_new < dist) dist = dist_new;
             }
          break;   
        default :
          cout << "\n cas du temps non implante temps= " << Nom_dure(temps)
               << "\n ElemMeca::Interne_Long_arrete_mini_sur_c(Enum_dure temps)";
          Sortie(1);                       
       }      
    // traitement d'une erreur éventuelle
    if (dist <= ConstMath::petit)
     { cout << "\n **** ERREUR une longueur d'arrete de l'element est nulle"
            << "\n ElemMeca::Interne_Long_arrete_mini_sur_c(..."
            << "\n element: "; this->Affiche(1);
       #ifdef MISE_AU_POINT
        cout << "\n *** version mise au point: on continue neanmoins avec une longueur "
             << " arbitrairement tres petite (" <<ConstMath::petit <<")  ";
      #else
        Sortie(1);
      #endif
    };
    // on calcul la célérité
    double cc = sqrt(ElemMeca::Calcul_maxi_E_qui(temps)/masse_volumique);
    // calcul du résultat
    double l_sur_c = dist/cc;
// debug
//if (l_sur_c ==0.00000000013879728102415373)
//cout << "\n ElemMeca::Interne_Long_arrete_mini_sur_c "<< ElemMeca::Calcul_maxi_E_qui() << " l " << dist
//     << " l_sur_c " << l_sur_c << endl;
//fin debug
    return l_sur_c;
   };
   
// activation du calcul des invariants de contraintes, qui seront calculé à chaque
// fois que l'on calcul les contraintes au travers de la loi de comportement
void ElemMeca::ActivCalculInvariantsContraintes()
 { int taille = lesPtIntegMecaInterne->NbPti();
   for (int i= 1; i<= taille; i++)
     (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_contrainte(true);
 };
// idem pour la déformation	  
void ElemMeca::ActivCalculInvariantsDeformation()
 { int taille = lesPtIntegMecaInterne->NbPti();
   for (int i= 1; i<= taille; i++)
     (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_deformation(true);
 };
// idem pour la vitesse de déformation	  
void ElemMeca::ActivCalculInvariantsVitesseDeformation()
 { int taille = lesPtIntegMecaInterne->NbPti();
   for (int i= 1; i<= taille; i++)
     (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_vitesseDeformation(true);
 };

            
// initialisation pour le calcul de la matrice masse dans le cas de l'algorithme 
// de relaxation dynamique  avec optimisation en continu de la matrice masse
// casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse
void ElemMeca::InitCalculMatriceMassePourRelaxationDynamique(int )
{ // on active le calcul des invariants de contraintes
  int taille = lesPtIntegMecaInterne->NbPti();
  for (int i= 1; i<= taille; i++)
     (*lesPtIntegMecaInterne)(i).Change_statut_Invariants_contrainte(true);
  
  // on définit  le ddl étendu correspondant, s'il n'existe pas 
  // on définit le  Ddl_enum_etendu correspondant à la masse au noeud pour la méthode
  if (masse_relax_dyn.Nom_vide() )
  	  masse_relax_dyn = (Ddl_enum_etendu("masse_relax_dyn"));
};

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

    // METHODES VIRTUELLES:
        

// recuperation des coordonnées du point de numéro d'ordre = iteg pour 
// la grandeur enu
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
// utilisé par les classes dérivées
Coordonnee ElemMeca::CoordPtInt(Enum_dure temps,Enum_ddl enu,int iteg,bool& err) 
  { Coordonnee ptret;err = false;
    // récupération de l'élément géométrique correspondant à Enu
    ElemGeomC0& ele = this->ElementGeometrie(enu);
    // récupération du tableau des fonctions d'interpolations
    const Tableau <Vecteur> & tabphi =   ele.TaPhi(); 
    if ((iteg < 0) || (iteg > tabphi.Taille()))
     { err = true;}
    else
     { switch (temps)
       { case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
         case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
         case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
       };
     };
    return ptret; // retour
   }; 
        
// recuperation des coordonnées du point d'intégration numéro = iteg pour
// la face : face
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
Coordonnee ElemMeca::CoordPtIntFace(int face, Enum_dure temps,int iteg,bool& err)
{ Coordonnee ptret(ParaGlob::Dimension());err = false;

  Deformation* definter=NULL; // variable de travail transitoire
  if (!(ParaGlob::AxiSymetrie()))
    {if (defSurf.Taille())
          {definter = defSurf(face);}    // le cas normal est non axisymétrique
     else {err = true;};
    }
  else // en axisymétrie, c'est une def d'arête
   { if (defArete.Taille())
          {definter = defArete(face);}
     else {err = true;};
   };
  if (definter == NULL)
   err = true;
  
  if (!err)
   { Deformation & defS = *definter;  // pour simplifier l'écriture
     defS.ChangeNumInteg(iteg);
     switch (temps)
        { case TEMPS_0 : ptret = defS.Position_0(); break;
          case TEMPS_t : ptret = defS.Position_t(); break;
          case TEMPS_tdt : ptret = defS.Position_tdt(); break;
        };
   };

  return ptret; // retour
 };

// recuperation des coordonnées du point d'intégration numéro = iteg pour
// la face : face
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
Coordonnee ElemMeca::CoordPtIntArete(int arete, Enum_dure temps,int iteg,bool& err)
{ Coordonnee ptret(ParaGlob::Dimension());err = false;

  Deformation* definter=NULL; // variable de travail transitoire
  if (defArete.Taille())
          {definter = defArete(arete);}
  else {err = true;};
  if (definter == NULL)
   err = true;

  if (!err)
   { Deformation & defA = *definter;  // pour simplifier l'écriture
     defA.ChangeNumInteg(iteg);
     switch (temps)
        { case TEMPS_0 : ptret = defA.Position_0(); break;
          case TEMPS_t : ptret = defA.Position_t(); break;
          case TEMPS_tdt : ptret = defA.Position_tdt(); break;
        };
   };

  return ptret; // retour
};

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

// ramene l'element geometrique correspondant au ddl passé en paramètre
ElemGeomC0& ElemMeca::ElementGeometrie(Enum_ddl  ddl) const
   { Enum_ddl enu = PremierDdlFamille(ddl);
     switch (enu)
     { case X1 : return ElementGeometrique(); break;
       case SIG11 : return ElementGeometrique(); break;
       case EPS11 : return ElementGeometrique(); break;
       case DEPS11 : return ElementGeometrique(); break;
       case ERREUR : return ElementGeometrique(); break;
       case TEMP :  return ElementGeometrique(); break;
       case V1 :  return ElementGeometrique(); break;
       case GAMMA1 :  return ElementGeometrique(); break;
       case DELTA_TEMP :  return ElementGeometrique(); break;
       case R_TEMP :  return ElementGeometrique(); break;
       case R_X1 :  return ElementGeometrique(); break;
       case R_V1 :  return ElementGeometrique(); break;
       case R_GAMMA1 :  return ElementGeometrique(); break;
       default :
        {cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl) 
              << "\n ElemMeca::ElementGeometrie(Enum_ddl  ddl) " ;
         Sortie(1);}
      }        
     return ElementGeometrique(); // pour taire le compilo 
    }; 

// ramène le nombre de grandeurs génératrices pour un pt d'integ, correspondant à un type enuméré
// peut-être surchargé pour des éléments particuliers 
int ElemMeca::NbGrandeurGene(Enum_ddl ddl) const
   { Enum_ddl enu = PremierDdlFamille(ddl);
     int nbGG = 0.; // par défaut
     switch (enu)
     { case X1 :  case SIG11 :  case EPS11 :  case DEPS11 : 
         { // cas d'un calcul de mécanique classique
           switch (Type_geom_generique(id_geom))
	        {   case LIGNE : case POINT_G :   nbGG = 1; break; // SIG11
		        case SURFACE : nbGG = 3; break; // SIG11, SIG22, SIG12
		        case VOLUME :  nbGG = 6; break; // les 6 composantes
		        default :
			      cout << "\nErreur : cas non traite, id_geom= :" << Nom_type_geom(Type_geom_generique(id_geom))
			           << "ElemMeca::NbGrandeurGene(.. \n";
			      Sortie(1);
	         };
        	break;	
        	} 
       default :
        {cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl) 
              << "\n ElemMeca::NbGrandeurGene(Enum_ddl  ddl) " ;
         Sortie(1);}
      };        
     return nbGG;  
    }; 
        
// modification de l'orientation de l'élément en fonction de cas_orientation
// =0: inversion simple (sans condition) de l'orientation
// si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration
// 1.  si tous les jacobiens sont négatifs on change d'orientation
// 2.  si tous les jacobiens sont positifs on ne fait rien
// 3.  si certains jacobiens sont positifs et d'autres négatifs message
//     d'erreur et on ne fait rien
// ramène true: s'il y a eu changement effectif, sinon false
bool ElemMeca::Modif_orient_elem(int cas_orientation)
 { // retour: 
   bool retour=false; // par défaut pas de changement
   if (cas_orientation == 0) // cas où on inverse l'orientation sans condition particulière
    { // on change l'orientation de l'élément
      retour = true;
      int nbnoe = tab_noeud.Taille();
      Tableau<Noeud *> tab_inter(tab_noeud); // on crée un tableau intermédiaire
      // on récupère la numérotation locale inverse
      const Tableau<int>  tabi =  ElementGeometrique().InvConnec();
		    // on met à jour le tableau actuel
      for (int n=1;n<=nbnoe;n++)
      	 tab_noeud(n)=tab_inter(tabi(n));
	   }
   else 
	   { // si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration
      // 1.  si tous les jacobiens sont négatifs on change d'orientation
      // 2.  si tous les jacobiens sont positifs on ne fait rien
      // 3.  si certains jacobiens sont positifs et d'autres négatifs message
      //     d'erreur et on ne fait rien
      int cas=1; // a priori tout est ok
      // boucle sur les pt d'integ
      for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg())
        { double jacobien_0 = def->JacobienInitial();
          
          if (jacobien_0 < 0)
          	{if (cas ==1) // on a trouvé un jacobien négatif  
          	  {cas =0; }
          	}// si c'était positif --> négatif
          else // cas positif
           	{if (cas == 0) // on a déjà changé
           	  { cas = 2; break;} // on sort de la boucle
          	};
        };
      // gestion de pb
      if (cas == 2)
      	{ cout << "\n **** Attention **** element nb= "<< this->Num_elt() << " peut-etre trop distordu ?" 
      	       << " pt d'integ meca positif et  negatif !! ";
          if (ParaGlob::NiveauImpression() > 7)
          {// on va essayer de sortir plus d'info
           int pti=1;cout << "\n les jacobiens initiaux \n";
           // on sort les valeurs des jacobiens
           for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg(),pti++)
            { cout << " pti= "<<pti << ", J= "<< def->JacobienInitial();};
           // les coordonnées des noeuds en suivant l'ordre de la table de connection locale
           int nbnoe = tab_noeud.Taille();
           cout << "\n les coordonnees de l'element \n";
           for (int ino=1;ino<=nbnoe;ino++)
            { cout << "\n "<< ino << " " << tab_noeud(ino)->Coord0();
            };
          
          };
      	}
      else if (cas == 0)
      	{ switch (cas_orientation)
           { case 1: // on change l'orientation de l'élément
              {retour = true;
               int nbnoe = tab_noeud.Taille();
               Tableau<Noeud *> tab_inter(tab_noeud); // on crée un tableau intermédiaire
               // on récupère la numérotation locale inverse
               const Tableau<int>  tabi =  ElementGeometrique().InvConnec();
               for (int n=1;n<=nbnoe;n++)
                 tab_noeud(n)=tab_inter(tabi(n));
               break;
              }
           	 case -1: // on sort une information à l'écran
              { cout << "\n element nb= "<< this->Num_elt() << "jacobien negatif " ;
                break;
              }
             default:
              cout << "\n erreur le cas : " << cas_orientation
                   << " n'est pas actuellement pris en compte"
                   << "\n ElemMeca::Modif_orient_elem(...";
              Sortie(1);     
           };
      	};
	 };	
	// retour
	return retour;	
 };

		         
// récupération de  valeurs interpolées pour les grandeur enu ou directement calculées
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
// une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul
Tableau <double> ElemMeca::Valeur_multi_interpoler_ou_calculer
                 (bool absolue, Enum_dure temps,const List_io<Ddl_enum_etendu>& enu
                  ,Deformation & defor
                  ,const Met_abstraite::Impli* ex_impli
                  ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
                  ,const Met_abstraite::Expli* ex_expli
                  )
 {
     // 2) le fait que l'on veut une sortie dans une base ad hoc ou pas
     int dim  = lesPtIntegMecaInterne->DimTens();int dim_sortie_tenseur = dim;
     // dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue
     if (absolue)
       dim_sortie_tenseur = ParaGlob::Dimension();
     // --- pour ne faire qu'un seul test ensuite
     bool prevoir_change_dim_tenseur = false;
     if (absolue)
        // mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire
        // car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue
        // il faut un tenseur intermédiaire  qui a ou non une dimension différente
        prevoir_change_dim_tenseur = true;
     // recup de l'incrément de temps
     double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
     double unSurDeltat=0;
     if (Abs(deltat) >= ConstMath::trespetit)
       {unSurDeltat = 1./deltat;}
     else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
       { // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
         if (unSurDeltat < 0)
          { cout << "\n le pas de temps est négatif !! "; };
         unSurDeltat = ConstMath::tresgrand;
       };

     // -- def des tenseurs locaux
     Coordonnee* Mtdt = NULL; // coordonnées finales éventuelles du point d'intégration considéré
     Coordonnee* Mt=NULL; // coordonnées à t
     Coordonnee* M0 = NULL; // coordonnées initiales éventuelles du point d'intégration considéré
     Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal actuel si c'est adéquate
     Coordonnee* N_surf_t = NULL; // coordonnée d'un vecteur normal à t si c'est adéquate
     Coordonnee* N_surf_t0 = NULL; // coordonnée d'un vecteur normal à t0 si c'est adéquate
     Coordonnee* Vitesse = NULL; // cas des vitesses

     Tableau <double> tab_ret (enu.size());
  
     // éléments de métrique et matrices de passage
     TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH;
     BaseB* giB_0;BaseB* giB_t;
  
     if (ex_impli != NULL)
       {gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt;
        giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt;
        giB_0 = ex_impli->giB_0;giB_t = ex_impli->giB_t;
       }
     else if (ex_expli_tdt != NULL)
       {gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt;
        giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt;
        giB_0 = ex_expli_tdt->giB_0; giB_t = ex_expli_tdt->giB_t;
       }
     else if (ex_expli != NULL)
       {gijHH = ex_expli->gijHH_t;gijBB = ex_expli->gijBB_t;
        giB = giB_t = ex_expli->giB_t; giH_0 = ex_expli->giH_0;giH = ex_expli->giH_t;
        giB_0 = ex_expli->giB_0;
       }
     else
      {
          cout << "\n *** cas non prevu : aucune metrique transmise "
               << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(..." << endl;
          Sortie(1);
      };

     // on définie des indicateurs pour ne pas faire plusieurs fois le même calcul
     List_io<Ddl_enum_etendu>::const_iterator ie,iefin=enu.end();
     bool besoin_coordonnees = false; bool besoin_deplacements = false;
     bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false;
     for (ie=enu.begin(); ie!=iefin;ie++)
       { int posi = (*ie).Position()-NbEnum_ddl();
         switch (posi)
          {case 114: case 115: case 116: // le vecteur normal
              { N_surf = new Coordonnee(ParaGlob::Dimension()); break;}
           case 117: case 118: case 119: // le vecteur normal à t
              { N_surf_t = new Coordonnee(ParaGlob::Dimension()); break;}
           case 120: case 121: case 122: // le vecteur normal à t0
              { N_surf_t0 = new Coordonnee(ParaGlob::Dimension()); break;}

           case 123: case 124: case 125: // la position à t
              { Mt = new Coordonnee(ParaGlob::Dimension());
                besoin_coordonnees_t = true;
                break;
              }
           case 126: case 127: case 128: // la position à t0
              { M0 = new Coordonnee(ParaGlob::Dimension());
                besoin_coordonnees_t0 = true;
                break;
              }
           default:
              break;
          };
       };
  
     // définition des tenseurs si nécessaire
  
  
     // ----- maintenant on calcule les grandeurs nécessaires -----
     if (besoin_coordonnees)
       {Mtdt = new Coordonnee(ParaGlob::Dimension());
        *Mtdt = defor.Position_tdt();
       };
     if (besoin_coordonnees_t )
       {*Mt = defor.Position_tdt();
       };
     if (besoin_deplacements || besoin_coordonnees_t0)
       {if (M0 == NULL)
           M0 = new Coordonnee(ParaGlob::Dimension());
        (*M0) = defor.Position_0();
       };
     if (Vitesse != NULL)
       {Vitesse = new Coordonnee(ParaGlob::Dimension());
        (*Vitesse) = defor.VitesseM_tdt();
       };
  
     // def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D
     if (N_surf != NULL)
      { // on vérifie que la métrique est correcte
        if (giB->NbVecteur() != 2)
         {if (ParaGlob::NiveauImpression() > 0)
            {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
                  << " le vecteur normal ne sera pas correctement calcule";
             if (ParaGlob::NiveauImpression() > 2)
               cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... ";
            };
          cout << endl;
         }
        else // sinon c'est ok
         {// la normale vaut le produit vectoriel des 2 premiers vecteurs
          (*N_surf) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2));
          N_surf->Normer(); // que l'on norme
         };
      };
     // idem à l'instant t
     if (N_surf_t != NULL)
      { // on vérifie que la métrique est correcte
        if (giB_t->NbVecteur() != 2)
         {if (ParaGlob::NiveauImpression() > 0)
            {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
                  << " le vecteur normal ne sera pas correctement calcule";
             if (ParaGlob::NiveauImpression() > 2)
               cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... ";
            };
          cout << endl;
         }
        else // sinon c'est ok
         {// la normale vaut le produit vectoriel des 2 premiers vecteurs
          (*N_surf_t) = Util::ProdVec_coorBN( (*giB_t)(1), (*giB_t)(2));
          N_surf_t->Normer(); // que l'on norme
         };
      };
     // idem à l'instant t0
     if (N_surf_t0 != NULL)
      { // on vérifie que la métrique est correcte
        if (giB_0->NbVecteur() != 2)
         {if (ParaGlob::NiveauImpression() > 0)
            {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
                  << " le vecteur normal ne sera pas correctement calcule";
             if (ParaGlob::NiveauImpression() > 2)
               cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... ";
            };
          cout << endl;
         }
        else // sinon c'est ok
         {// la normale vaut le produit vectoriel des 2 premiers vecteurs
          (*N_surf_t0) = Util::ProdVec_coorBN( (*giB_0)(1), (*giB_0)(2));
          N_surf_t0->Normer(); // que l'on norme
         };
       // on n'arrête pas l'exécution, car on pourrait vouloir sortir les normales pour un ensemble
       // d'éléments contenant des volumes, des surfaces, des lignes: bon... il y aura quand même des
       // pb au niveau des iso par exemple, du au fait que l'on va faire des moyennes sur des éléments
       // de type différents (à moins de grouper par type du coup on n'aura pas le warning
      };
     //----- fin du calcul des grandeurs nécessaires  -----
  
     // on balaie maintenant la liste des grandeurs à sortir
     int it; // it est l'indice dans le tableau de retour
     for (it=1,ie=enu.begin(); ie!=iefin;ie++,it++)
      { if ((Meme_famille((*ie).Enum(),SIG11)) || (Meme_famille((*ie).Enum(),EPS11))
           || (Meme_famille((*ie).Enum(),DEPS11)) || (Meme_famille((*ie).Enum(),X1))
           || (Meme_famille((*ie).Enum(),UX)) )
         { // on recherche en générale une interpolation en fonction des noeuds: il faut donc que la grandeur soit
           // présente aux noeuds !!
          // def du numéro de référence du ddl_enum_etendue
          int posi = (*ie).Position()-NbEnum_ddl();
          // récupération des informations en fonction des différents cas
          // **** 1 >>>>> -- cas des ddl pur, que l'on sort dans le repère global par défaut
          // cas des contraintes
          if ((Meme_famille((*ie).Enum(),SIG11)) && ((*ie).Nom_vide()))
              { tab_ret(it)= defor.DonneeInterpoleeScalaire((*ie).Enum(),temps);
              }
          else if ((Meme_famille((*ie).Enum(),EPS11)) && ((*ie).Nom_vide()))
              { tab_ret(it)= defor.DonneeInterpoleeScalaire((*ie).Enum(),temps);
              }
          else if ((Meme_famille((*ie).Enum(),DEPS11)) && ((*ie).Nom_vide()))
              { tab_ret(it)= defor.DonneeInterpoleeScalaire((*ie).Enum(),temps);
              }
          else if ((Meme_famille((*ie).Enum(),X1)) && ((*ie).Nom_vide()))
              { tab_ret(it)= (*Mtdt)((*ie).Enum() - X1 +1);
              }
          else if ((Meme_famille((*ie).Enum(),UX)) && ((*ie).Nom_vide()))
              { int i_cor = (*ie).Enum() - UX +1; // l'indice de coordonnée
                tab_ret(it)= (*Mtdt)(i_cor) - (*M0)(i_cor);
              }
          else if ((Meme_famille((*ie).Enum(),V1)) && ((*ie).Nom_vide()))
              { int i_cor = (*ie).Enum() - V1 +1; // l'indice de coordonnée
                tab_ret(it)= (*Vitesse)(i_cor);
              }
             // --- a complèter ----
          
          else
             {// **** 2 >>>>> -- cas des grandeurs déduites des ddl pures
              switch (posi)
              { case 114: // le vecteur normal N_surf_1
                 {tab_ret(it)= (*N_surf)(1);break;}
                case 115: // le vecteur normal N_surf_2
                 {tab_ret(it)= (*N_surf)(2);break;}
                case 116: // le vecteur normal N_surf_3
                 {tab_ret(it)= (*N_surf)(3);break;}
                case 117: // le vecteur normal N_surf_1_t
                 {tab_ret(it)= (*N_surf_t)(1);break;}
                case 118: // le vecteur normal N_surf_2_t
                 {tab_ret(it)= (*N_surf_t)(2);break;}
                case 119: // le vecteur normal N_surf_3_t
                 {tab_ret(it)= (*N_surf_t)(3);break;}
                case 120: // le vecteur normal N_surf_1_t0
                 {tab_ret(it)= (*N_surf_t0)(1);break;}
                case 121: // le vecteur normal N_surf_2_t0
                 {tab_ret(it)= (*N_surf_t0)(2);break;}
                case 122: // le vecteur normal N_surf_3_t0
                 {tab_ret(it)= (*N_surf_t0)(3);break;}
                case 123: // la position géométrique Mt
                 {tab_ret(it)= (*Mt)(1);break;}
                case 124: // la position géométrique Mt
                 {tab_ret(it)= (*Mt)(2);break;}
                case 125: // la position géométrique Mt
                 {tab_ret(it)= (*Mt)(3);break;}
                case 126: // la position géométrique M0
                 {tab_ret(it)= (*M0)(1);break;}
                case 127: // la position géométrique M0
                 {tab_ret(it)= (*M0)(2);break;}
                case 128: // la position géométrique M0
                 {tab_ret(it)= (*M0)(3);break;}
                
                default :
                  {cout << "\n cas de ddl actuellement non traite "
                        << "\n pas de ddl " << (*ie).Nom() << "  dans l'element "
                        << "\n ou cas non implante pour l'instant"
                        << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(...";
                   tab_ret(it) = 0.;
                   };
               }  // fin cas **** 2 >>>>>
              }   //   "      "       "
          }  // -- fin
        else if (( (*ie).Enum() == TEMP) && ((*ie).Nom_vide()))
           {// on vérifie que le ddl existe au premier noeud
            if (tab_noeud(1)->Existe_ici(TEMP))
             {tab_ret(it)= def->DonneeInterpoleeScalaire(TEMP,temps);}
            else if (ParaGlob::NiveauImpression()>3)
              {cout << "\n pas de ddl temperature disponible au noeud "
                    << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(...";
              // this->Affiche();
              };
           }
        else
          { tab_ret(it) = 0.;
            cout << "\n cas de ddl actuellement non traite "
                << "\n pas de ddl " << (*ie).Nom() << "  dans l'element "
                << "\n ou cas non implante pour l'instant, on retourne 0"
                << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(...";
           };
      };// -- fin de la boucle sur la liste de Ddl_enum_etendu
  
     // delet e des tenseurs
     if (Mtdt != NULL) delete Mtdt; // coordonnée du point à tdt
     if (Mt != NULL ) delete Mt; // la position à t
     if (M0 != NULL ) delete M0; // coordonnée du point à 0
     if (N_surf != NULL) delete N_surf; // vecteur normal à la surface
     if (N_surf_t != NULL) delete N_surf_t; // vecteur normal à t à la surface
     if (N_surf_t0 != NULL) delete N_surf_t0; // vecteur normal à la surface
     if (Vitesse != NULL) delete Vitesse; // vitesse
     // liberation des tenseurs intermediaires
     LibereTenseur();
     return tab_ret;
    };

// récupération de  valeurs interpolées pour les grandeur enu
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
// une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul
void ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer
                 (bool absolue, Enum_dure temps,List_io<TypeQuelconque>& enu
                  ,Deformation & defor
                  ,const Met_abstraite::Impli* ex_impli
                  ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
                  ,const Met_abstraite::Expli* ex_expli
                 )

{  // ----- def de grandeurs de travail
   // def de la dimension des tenseurs
   // il y a deux pb a gérer: 1) le fait que la dimension absolue peut-être différente de la dimension des tenseurs
   // 2) le fait que l'on veut une sortie dans une base ad hoc ou pas
   int dim  = lesPtIntegMecaInterne->DimTens();int dim_sortie_tenseur = dim;
   // dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue
   int dim_espace = ParaGlob::Dimension();
   if (absolue)
     dim_sortie_tenseur = dim_espace;
   // pour ne faire qu'un seul test ensuite
   bool prevoir_change_dim_tenseur = false;
   // initialement on faisait le test suivant,
   // if  ((absolue) && (dim != dim_sortie_tenseur))
   if (absolue)
      // mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire
      // car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue
      // il faut un tenseur intermédiaire qui a ou non une dimension différente
      prevoir_change_dim_tenseur = true;

     // éléments de métrique et matrices de passage
     TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH;
     BaseB* giB_0;BaseB* giB_t;
 
     if (ex_impli != NULL)
       {gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt;
        giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt;
        giB_t = ex_impli->giB_t; giB_0 = ex_impli->giB_0;
       }
     else if (ex_expli_tdt != NULL)
       {gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt;
        giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt;
        giB_t = ex_expli_tdt->giB_t; giB_0 = ex_expli_tdt->giB_0;
       }
     else if (ex_expli != NULL)
       {gijHH = ex_expli->gijHH_t;gijBB = ex_expli->gijBB_t;
        giB = giB_t = ex_expli->giB_t; giH_0 = ex_expli->giH_0;giH = ex_expli->giH_t;
        giB_0 = ex_expli->giB_0;
       }
     else
      {
          cout << "\n *** cas non prevu : aucune metrique transmise "
               << "\n ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer(..." << endl;
          Sortie(1);
      };

   // def de tenseurs pour la sortie
   // les tenseurs restants en locale
   Coordonnee* Mtdt=NULL; // coordonnées éventuelles du point d'intégration considéré
   Coordonnee* Mt=NULL; // coordonnées à t
   Coordonnee* M0=NULL; // coordonnées à t0
   Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal si c'est adéquate
   Coordonnee* N_surf_t = NULL; // coordonnée d'un vecteur normal à t si c'est adéquate
   Coordonnee* N_surf_t0 = NULL; // coordonnée d'un vecteur normal à t0 si c'est adéquate
   Coordonnee* Vitesse = NULL; // cas des vitesses
   Coordonnee* Acceleration = NULL; // cas des accélérations
   Coordonnee* Deplacement = NULL; // cas du déplacement

   // pour les valeurs propres
   // pour les vecteurs propres
   // pour les bases
   Tableau <Coordonnee *> base_ad_hoc , base_giH , base_giB;
 

   // --- dev d'un ensemble de variable booléenne pour gérer les sorties en une passe -----
   // on se réfère au informations définit dans la méthode: Les_type_evolues_internes()
   bool besoin_coordonnees = false;
   bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false;
   bool besoin_rep_local_ortho=false;
   bool besoin_rep_giH = false;        bool besoin_rep_giB = false;
   // on initialise ces variables booléennes et les conteneurs
   List_io<TypeQuelconque>::iterator ipq,ipqfin=enu.end();
   for (ipq=enu.begin();ipq!=ipqfin;ipq++)
    {switch ((*ipq).EnuTypeQuelconque().EnumTQ())
     { case POSITION_GEOMETRIQUE : {besoin_coordonnees=true;
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         Mtdt = gr.ConteneurCoordonnee();  break;}
       case POSITION_GEOMETRIQUE_t : {besoin_coordonnees_t=true;
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         Mt = gr.ConteneurCoordonnee();  break;}
       case POSITION_GEOMETRIQUE_t0 : {besoin_coordonnees_t0=true;
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         M0 = gr.ConteneurCoordonnee();  break;}
       case REPERE_LOCAL_ORTHO : {besoin_rep_local_ortho=true;
         Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
         base_ad_hoc.Change_taille(dim_espace);
         switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
         switch (dim_espace) {case 3: base_ad_hoc(3) = &(gr(3));
                              case 2: base_ad_hoc(2) = &(gr(2));
                              case 1: base_ad_hoc(1) = &(gr(1));
                              default:break;
                             };
         break;}
       case REPERE_LOCAL_H : {besoin_rep_giH=true;
         Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
         base_giH.Change_taille(dim);
         switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
         switch (dim) {case 3: base_giH(3) = &(gr(3));
                       case 2: base_giH(2) = &(gr(2));
                       case 1: base_giH(1) = &(gr(1));
                       default:break;
                      };
         break;}
       case REPERE_LOCAL_B : {besoin_rep_giB=true;
         Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
         base_giB.Change_taille(dim);
         switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
         switch (dim) {case 3: base_giB(3) = &(gr(3));
                       case 2: base_giB(2) = &(gr(2));
                       case 1: base_giB(1) = &(gr(1));
                       default:break;
                      };
         break;}

       case NN_SURF:{
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         N_surf = gr.ConteneurCoordonnee();  break;}
       case NN_SURF_t:{
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         N_surf_t = gr.ConteneurCoordonnee();  break;}
       case NN_SURF_t0:{
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         N_surf_t0 = gr.ConteneurCoordonnee();  break;}
       case DEPLACEMENT:{besoin_coordonnees=true;
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         Deplacement = gr.ConteneurCoordonnee();  break;}
       case VITESSE:{
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         Vitesse = gr.ConteneurCoordonnee();  break;}
       case ACCELERATION:{
         Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
         Acceleration = gr.ConteneurCoordonnee();
         break;}
       // dans le cas des numéros, traitement direct ici
       case NUM_ELEMENT:
        { *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= num_elt; break;}
       case NUM_MAIL_ELEM:
        { *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= num_maillage; break;}

       default :
         {// on initialise la grandeur pour éviter d'avoir des valeurs aléatoires
          ((*ipq).Grandeur_pointee())->InitParDefaut();
          if (ParaGlob::NiveauImpression() > 0)
           {cout << "\nWarning : attention cas non traite: "
                 <<  (*ipq).EnuTypeQuelconque().NomPlein() << "!\n";
            // on initialise la grandeur pour éviter d'avoir des valeurs aléatoires
            if (ParaGlob::NiveauImpression() > 5)
               cout << "\n ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer(....";
           };
         }
     };
    };

   // récup des bases si besoin
   if (besoin_rep_local_ortho)
     {// on calcule éventuellement la matrice de passage
      Mat_pleine Aa0(dim,dim),Aafin(dim,dim);
      if (   ((!absolue)&&(dim_espace==3)&&(dim==2))
          || ((!absolue)&&(dim_espace>1)&&(dim==1))
         )
       {
         def->BasePassage(absolue,*giB_0,*giB,*giH_0,*giH,Aa0,Aafin);
         //  g^i = Aa^i_{.a} * Ip^a
         // et on a : Ip_a = beta_a^{.j} g_j = [A^j_{.a}]^T g_j
       };
      if ((!absolue)&&(dim_espace==3)&&(dim==2))
        // cas d'éléments 2D, pour lesquels on veut un repère local ad hoc
        // on ramène une base à 3 vecteurs
        { *(base_ad_hoc(1)) = Aafin(1,1) * giB->Coordo(1) + Aafin(2,1) * giB->Coordo(2);
          *(base_ad_hoc(2)) = Aafin(1,2) * giB->Coordo(1) + Aafin(2,2) * giB->Coordo(2);
          *(base_ad_hoc(3)) = Util::ProdVec_coor(*(base_ad_hoc(1)),*(base_ad_hoc(2)));
        }
      else if((!absolue)&&(dim_espace>1)&&(dim==1))
        // cas d'éléments 1D, dans un espace 2D ou 3D
        // on ramène un seul vecteur non nul, les autres ne peuvent être calculé sans info supplémentaire
         { *(base_ad_hoc(1)) = Aafin(1,1) * giB->Coordo(1);
           (base_ad_hoc(2))->Zero(); (base_ad_hoc(3))->Zero(); // init à 0 par défaut
        }
      else // dans tous les autres cas
       { switch (dim_espace) // on se contente de ramener le repère identité
          {case 3: (base_ad_hoc(3))->Zero();(*(base_ad_hoc(3)))(3)=1.;
           case 2: (base_ad_hoc(2))->Zero();(*(base_ad_hoc(2)))(2)=1.;
           case 1: (base_ad_hoc(1))->Zero();(*(base_ad_hoc(1)))(1)=1.;
           default:break;
          };
       };
     };
 
   if (besoin_rep_giH)
     {switch (dim) {case 3: *(base_giH(3)) = giH->Coordo(3);
                    case 2: *(base_giH(2)) = giH->Coordo(2);
                    case 1: *(base_giH(1)) = giH->Coordo(1);
                    default:break;
                   };
     };
   if (besoin_rep_giB)
     {switch (dim) {case 3: *(base_giB(3)) = giB->Coordo(3);
                    case 2: *(base_giB(2)) = giB->Coordo(2);
                    case 1: *(base_giB(1)) = giB->Coordo(1);
                    default:break;
                   };
     };
 
   // ----- calcul des grandeurs à sortir
 
   if (besoin_coordonnees)
        (*Mtdt) = defor.Position_tdt();
   if (besoin_coordonnees_t)
        (*Mt) = defor.Position_t();
   if (besoin_coordonnees_t0)
        (*M0) = defor.Position_0();

   // def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D
   if (N_surf != NULL)
      { // on vérifie que la métrique est correcte
        if (giB->NbVecteur() != 2)
         {if (ParaGlob::NiveauImpression() > 0)
            {cout << "\n *** attention il ne s'agit pas d'un element 2D,"
                  << " le vecteur normal ne sera pas disponible";
             if (ParaGlob::NiveauImpression() > 2)
               cout << "\n ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer(... ";
            };
          cout << endl;
         }
        else // sinon c'est ok
         {// la normale vaut le produit vectoriel des 2 premiers vecteurs
          (*N_surf) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2));
          N_surf->Normer(); // que l'on norme
         };
       // on n'arrête pas l'exécution, car on pourrait vouloir sortir les normales pour un ensemble
       // d'éléments contenant des volumes, des surfaces, des lignes: bon... il y aura quand même des
       // pb au niveau des iso par exemple, du au fait que l'on va faire des moyennes sur des éléments
       // de type différents (à moins de grouper par type du coup on n'aura pas le warning
      };
     // idem à l'instant t
     if (N_surf_t != NULL)
      { // on vérifie que la métrique est correcte
        if (giB_t->NbVecteur() != 2)
         {if (ParaGlob::NiveauImpression() > 0)
            {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
                  << " le vecteur normal ne sera pas correctement calcule";
             if (ParaGlob::NiveauImpression() > 2)
               cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... ";
            };
          cout << endl;
         }
        else // sinon c'est ok
         {// la normale vaut le produit vectoriel des 2 premiers vecteurs
          (*N_surf_t) = Util::ProdVec_coorBN( (*giB_t)(1), (*giB_t)(2));
          N_surf_t->Normer(); // que l'on norme
         };
      };
     // idem à l'instant t0
     if (N_surf_t0 != NULL)
      { // on vérifie que la métrique est correcte
        if (giB_0->NbVecteur() != 2)
         {if (ParaGlob::NiveauImpression() > 0)
            {cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
                  << " le vecteur normal ne sera pas correctement calcule";
             if (ParaGlob::NiveauImpression() > 2)
               cout << "\n ElemMeca::Valeur_multi_interpoler_ou_calculer(... ";
            };
          cout << endl;
         }
        else // sinon c'est ok
         {// la normale vaut le produit vectoriel des 2 premiers vecteurs
          (*N_surf_t0) = Util::ProdVec_coorBN( (*giB_0)(1), (*giB_0)(2));
          N_surf_t0->Normer(); // que l'on norme
         };
      };

   // cas du déplacement et de la vitesse
   if (Deplacement != NULL)
    (*Deplacement) = (*Mtdt) - defor.Position_0();
   if (Vitesse != NULL)
    (*Vitesse) = defor.VitesseM_tdt();
   if (Acceleration != NULL)
//    switch (dim_espace)
//     { case 3: (*Acceleration) = defor.DonneeInterpoleeScalaire(GAMMA3,temps);
//       case 2: (*Acceleration) = defor.DonneeInterpoleeScalaire(GAMMA2,temps);
//       case 1: (*Acceleration) = defor.DonneeInterpoleeScalaire(GAMMA1,temps);
//     };
    (*Acceleration) = defor.AccelerationM_tdt();

   // --- cas des grandeurs de la décomposition polaire

   // delete des tenseurs

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


// mise à jour éventuel de repère d'anisotropie
void ElemMeca::Mise_a_jour_repere_anisotropie
                (BlocGen & bloc,LesFonctions_nD*  lesFonctionsnD)
 { // cas de la définition d'un repère d'anisotropie
   // pour l'instant le repère d'anisotropie est éventuellement
   // utilisable par la loi de comportement
   // 1) tout d'abord on doit définir le repère au niveau de chaque point d'intégration
   // définition d'un repère d'anisotropie à un point d'intégration
   int nb_pti = lesPtIntegMecaInterne->NbPti();
   for (int ni=1;ni<= nb_pti;ni++)
     {BaseH repH = DefRepereAnisotropie(ni,lesFonctionsnD,bloc);
      // on crée les grandeurs de passages pour renseigner les grandeurs particulières
      // de la loi de comportement
      int dim = ParaGlob::Dimension();
      Tableau <Coordonnee> tab_coor(dim);
      for (int i=1;i<=dim;i++)
       tab_coor(i)=repH.Coordo(i);
      // passage des infos: on utilise la méthode Complete_SaveResul
      // qui au premier passage via ElemMeca::Complete_ElemMeca(
      // a créé les conteneurs et au passage courant met à jour seulement
      if (tabSaveDon(ni) != NULL)
         tabSaveDon(ni)->Complete_SaveResul(bloc,tab_coor,loiComp);
     }; 
 };
