// FICHIER : SfeMembT.cc
// CLASSE : SfeMembT


// 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-2021 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 <iostream>
#include <cmath>
using namespace std;  //introduces namespace std
#include <stdlib.h>
#include "Sortie.h"
#include "MathUtil.h"

#include "SfeMembT.h"
#include "TypeConsTens.h"
#include "FrontPointF.h"
#include "Loi_Umat.h"
#include "ExceptionsElemMeca.h"
#include "ExceptionsLoiComp.h"
#include "CharUtil.h"

// =========================== variables communes ===============================
// a priori les pointeurs de fonction pointent sur rien au debut
TenseurHH* SfeMembT::sig_bulk_pourSfe_HH = NULL;  // variable de travail pour le bulk
                    
// =========================== fin variables communes ===========================
// ---- definition du constructeur de la classe conteneur de donnees communes ------------	

SfeMembT::DonnComSfe::DonnComSfe (ElemGeomC0* elCentre,const GeomSeg& seg,const  DdlElement& tab
	                       ,DdlElement&  tabErr,DdlElement&  tab_Err1Sig,DdlElement&  t_ddlXi_C
	                       ,const Met_Sfe1&  met_gene,
	                        Tableau <Vecteur *> & resEr,Mat_pleine& raidEr,
	                        ElemGeomC0* elEr,GeomSeg& segEr,ElemGeomC0* elS,GeomSeg& segmS
	                       ,Vecteur&  residu_int,Mat_pleine&  raideur_int,
	                        Tableau <Vecteur* > & residus_extN,Tableau <Mat_pleine* >&  raideurs_extN,
	                        Tableau <Vecteur* > & residus_extA,Tableau <Mat_pleine* >&  raideurs_extA,
	                        Tableau <Vecteur* >&  residus_extS,Tableau <Mat_pleine* >&  raideurs_extS,
	                        Mat_pleine&  mat_masse,ElemGeomC0* elMas,const GeomSeg& segMas,int nbi,
	                        Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan
	                        ,Tableau <int>& nMetddl,const Met_abstraite&  met_cent ) :
	          eleCentre(elCentre),segment(seg),tab_ddl(tab),tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig)
	          ,tab_ddlXi_C(t_ddlXi_C),met_SfeMembT(met_gene)
	          ,matGeom(tab.NbDdl(),tab.NbDdl()),matInit(tab.NbDdl(),tab.NbDdl())
	          ,d_epsBB(tab.NbDdl()),d_sigHH(tab.NbDdl()),d2_epsBB(nbi)
	          ,resErr(resEr),raidErr(raidEr),eleEr(elEr),segmentEr(segEr),eleS(elS),segS(segmS)
	          ,residu_interne(residu_int),raideur_interne(raideur_int)
	          ,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN)
	          ,residus_externeA(residus_extA),raideurs_externeA(raideurs_extA)
	          ,residus_externeS(residus_extS),raideurs_externeS(raideurs_extS)
	          ,matrice_masse(mat_masse),eleMas(elMas),segmentMas(segMas)
	          ,tabType_rienCL(tabTypeCL),vplan_rien(vplan)
	          ,nMetVTab_ddl(nMetddl),met_surf_cent(met_cent)
           ,sfematD(NULL),sferesD(NULL),noeud_a_prendre_en_compte(NULL)
	{ 
   int nbddl = tab.NbDdl();
   int dim_tens = met_gene.Nbvec_des_bases();
   for (int ni=1;ni<=nbi;ni++)
    {d2_epsBB(ni).Change_taille(nbddl);
     for (int i1=1; i1<= nbddl; i1++)
      for (int i2=1; i2<= nbddl; i2++)
       d2_epsBB(ni)(i1,i2) = NevezTenseurBB (dim_tens);
    };
//      int tailledeps = d_epsBB.Taille();
//	  for (int i=1;i<= tailledeps; i++)
	  for (int i=1;i<= nbddl; i++)
	         d_epsBB(i) = NevezTenseurBB (dim_tens);
//      int tailledsig = d_sigHH.Taille();
//	  for (int j=1;j<= tailledsig; j++)
	  for (int j=1;j<= nbddl; j++)
	         d_sigHH(j) = NevezTenseurHH (dim_tens);
	 };
// Il ne doit exister  qu'un exemplaire de la classe, donc au niveau du constructeur de copie les tableaux
// de pointeurs pointes sur les mêmes entités que la valeur passée en paramètre	              
SfeMembT::DonnComSfe::DonnComSfe(DonnComSfe& a) :
   eleCentre(a.eleCentre),segment(a.segment)
			,tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr),tab_Err1Sig11(a.tab_Err1Sig11)
   ,tab_ddlXi_C(a.tab_ddlXi_C),met_SfeMembT(a.met_SfeMembT),matGeom(a.matGeom)
   ,matInit(a.matInit),d2_epsBB(a.d2_epsBB),resErr(a.resErr),raidErr(a.raidErr)
   ,eleEr(a.eleEr),segmentEr(a.segmentEr),eleS(a.eleS),segS(a.segS)
	  ,d_epsBB(a.d_epsBB),d_sigHH(a.d_sigHH)
	  ,residu_interne(a.residu_interne),raideur_interne(a.raideur_interne)
	  ,residus_externeN(a.residus_externeN),raideurs_externeN(a.raideurs_externeN)
	  ,residus_externeA(a.residus_externeA),raideurs_externeA(a.raideurs_externeA)
	  ,residus_externeS(a.residus_externeS),raideurs_externeS(a.raideurs_externeS)
	  ,matrice_masse(a.matrice_masse),eleMas(a.eleMas),segmentMas(a.segmentMas)
	  ,tabType_rienCL(a.tabType_rienCL),vplan_rien(a.vplan_rien)
	  ,nMetVTab_ddl(a.nMetVTab_ddl),met_surf_cent(a.met_surf_cent)
   ,sfematD(NULL),sferesD(NULL),noeud_a_prendre_en_compte(NULL)

	       { int nbddl = tab_ddl.NbDdl();
	         int nbi=d2_epsBB.Taille();
	         for (int ni=1;ni<=nbi;ni++)
              for (int i1=1; i1<= nbddl; i1++)
               for (int i2=1; i2<= nbddl; i2++)
                 d2_epsBB(ni)(i1,i2) = NevezTenseurBB (*(a.d2_epsBB(ni)(i1,i2)));
//             int tailledeps = d_epsBB.Taille();
//	         for (int i=1;i<= tailledeps; i++)
	         for (int i=1;i<= nbddl; i++)
	            d_epsBB(i) = NevezTenseurBB (*(a.d_epsBB(i)));
//             int tailledsig = d_sigHH.Taille();
//	         for (int j=1;j<= tailledsig; j++)
	         for (int j=1;j<= nbddl; j++)
	            d_sigHH(j) = NevezTenseurHH (*(d_sigHH(j)));
             
          // cas d'une stabilisation éventuelle
          if (a.sfematD != NULL)
            sfematD = new Mat_pleine(*a.sfematD);
          if (a.sferesD != NULL)
            sferesD = new Vecteur (*a.sferesD);
          if (a.noeud_a_prendre_en_compte != NULL)
            noeud_a_prendre_en_compte = new Tableau<int> (*a.noeud_a_prendre_en_compte);
        };
				  
// certaines grandeurs pointées qui n'exitent qu'en un seul exemplaire 
// ne seront supprimées que par l'ordre destruction, 				    
SfeMembT::DonnComSfe::DonnComSfe::~DonnComSfe()
       { int nbddl = tab_ddl.NbDdl();
	 	      int nbi=d2_epsBB.Taille();
	 	      for (int ni=1;ni<=nbi;ni++)
            for (int i1=1; i1<= nbddl; i1++)
               for (int i2=1; i2<= nbddl; i2++)
                  delete d2_epsBB(ni)(i1,i2);
         int tailledeps = d_epsBB.Taille();
	        for (int i=1;i<= tailledeps; i++)
            delete d_epsBB(i);
         int tailledsig = d_sigHH.Taille();
	        for (int j=1;j<= tailledsig; j++)
            delete d_sigHH(j);
         // cas d'une stabilisation éventuelle
         if (sfematD != NULL)
           delete sfematD;
         if (sferesD != NULL)
           delete sferesD;
         if (noeud_a_prendre_en_compte != NULL)
           delete noeud_a_prendre_en_compte;

       };
	             
	            
// ---------- fin definition de la  classe conteneur de donnees communes ------------	
// -+-+ definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
SfeMembT::UneFois::UneFois () : // constructeur par défaut
	doCoMemb(NULL),CalResPrem_t(0),CalResPrem_tdt(0),CalimpPrem(0),dualSortSfe(0)
	,CalSMlin_t(0),CalSMlin_tdt(0),CalSMRlin(0)
	,CalSMsurf_t(0),CalSMsurf_tdt(0),CalSMRsurf(0)
	,CalSMvol_t(0),CalSMvol_tdt(0),CalSMvol(0)
	,CalDynamique(0),CalPt_0_t_tdt(0)
	,nbelem_in_Prog(0)
	        {};
SfeMembT::UneFois::~UneFois ()
	 {	  delete doCoMemb;
	 };   

// -+-+ fin definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
	            
// CONSTRUCTEURS :
	// Constructeur par defaut
SfeMembT::SfeMembT () :
   ElemMeca(),lesPtMecaInt(),donnee_specif(),areteTypeCL(NULL),vplan(NULL)
	,t_N_centre(),unefois(NULL),nombre(NULL)
      { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca
      };   
// Constructeur fonction  d'un numero
// d'identification , d'identificateur d'interpolation et de geometrie 
SfeMembT::SfeMembT (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,string info) :
     ElemMeca(num_mail,num_id,id_interp_elt,id_geom_elt,info),lesPtMecaInt(),donnee_specif()
     ,areteTypeCL(NULL),vplan(NULL),t_N_centre(),unefois(NULL),nombre(NULL)
      { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca
      };
		
// Constructeur fonction  d'un numero de maillage et d'identification,
// du tableau de connexite des noeuds, d'identificateur d'interpolation et de geometrie
SfeMembT::SfeMembT (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,
                     const Tableau<Noeud *>& tab,string info) :
 ElemMeca(num_mail,num_id,tab,id_interp_elt,id_geom_elt,info),lesPtMecaInt(),donnee_specif()
 ,areteTypeCL(NULL),vplan(NULL),t_N_centre(),unefois(NULL),nombre(NULL)
      { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca
      };
// Constructeur de copie
SfeMembT::SfeMembT (const SfeMembT& SfeM) :
 ElemMeca (SfeM)
 ,lesPtMecaInt(SfeM.lesPtMecaInt)
 ,donnee_specif(SfeM.donnee_specif),unefois(SfeM.unefois),nombre(SfeM.nombre)
 ,areteTypeCL(SfeM.areteTypeCL),vplan(SfeM.vplan),t_N_centre(SfeM.t_N_centre)
   { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca
		   // --- cas des différentes  puissances et ... ---
     SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
     residu = &(CoSfe->residu_interne); // residu local
     raideur = &(CoSfe->raideur_interne); // raideur locale
         // --- cas de la dynamique -----
     mat_masse =  &(CoSfe->matrice_masse);
     	// --- cas des efforts externes concernant noeuds ------
     res_extN = &(CoSfe->residus_externeN); // pour les résidus et second membres
     raid_extN= &(CoSfe->raideurs_externeN);// pour les raideurs
     	// --- cas des efforts externes concernant les aretes ------
     res_extA = &(CoSfe->residus_externeA); // pour les résidus et second membres
     raid_extA= &(CoSfe->raideurs_externeA);// pour les raideurs
     	// --- cas des efforts externes concernant les faces  ------
     res_extS= &(CoSfe->residus_externeS); // pour les résidus et second membres
     raid_extS= &(CoSfe->raideurs_externeS); // pour les raideurs
     // pour les conditions limites on les crée si besoin
     if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL)))
     	{ // si ça ne pointe pas sur la valeur par défaut, alors c'est qu'il y a
     	  // des CL particulières, on les crée donc à l'identique
     	  areteTypeCL = new Tableau <EnuTypeCL> (*(SfeM.areteTypeCL));
     	  vplan = new Tableau <Coordonnee3> (*(SfeM.vplan));
     	};
  };
	
// DESTRUCTEUR :
SfeMembT::~SfeMembT ()
{   /*  // on va libérer les déformations d'arrête et de surface qui sont un peu particulières
  if (defArete.Taille() != 0)
   if (defArete(1) != NULL) 
    { delete defArete(1);
      // on remet à null les pointeurs pour un appel correcte du destructeur d'elemMeca
      for (int ia=1;ia<= 12; ia++)
        defArete(ia) = NULL; 
     }    
  if (defSurf.Taille() != 0)
   if (defSurf(1) != NULL) 
    { delete defSurf(1);
      // on remet à null les pointeurs pour un appel correcte du destructeur d'elemMeca
      for (int ia=1;ia<= 12; ia++)
        defSurf(ia) = NULL; 
     }  */
  // concernant tous les éléments qui sont définit qu'une seule fois, leur destruction est effectuée avec la méthode destruction
  // appelée par les classes descendantes	   
  LibereTenseur();
 };


// Lecture des donnees de la classe sur fichier
void SfeMembT::LectureDonneesParticulieres
    (UtilLecture * entreePrinc,Tableau<Noeud  *> * tabMaillageNoeud)
  { int nb;int nbnte = nombre->nbnte;
    tab_noeud.Change_taille(nbnte);
    for (int i=1; i<= nbnte; i++)
     { *(entreePrinc->entree) >> nb;
       if ((entreePrinc->entree)->rdstate() == 0) 
       // pour mémoire ici on a 
           /*       enum io_state
                    {  badbit   = 1<<0, // -> 1 dans rdstate()
                       eofbit   = 1<<1, // -> 2 
                       failbit  = 1<<2, // -> 4 
                       goodbit  = 0     // -> O 
                     };*/
          tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
     #ifdef ENLINUX
       else  if ((entreePrinc->entree)->fail())
               // on a atteind la fin de la ligne et on appelle un nouvel enregistrement 
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
               *(entreePrinc->entree) >> nb;
              tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
              }
     #else
 /*     #ifdef SYSTEM_MAC_OS_X_unix
       else  if ((entreePrinc->entree)->fail())
               // on a atteind la fin de la ligne et on appelle un nouvel enregistrement 
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
               *(entreePrinc->entree) >> nb;
              tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
              }
      #else */
       else  if ((entreePrinc->entree)->eof())
        // la lecture est bonne mais on a atteind la fin de la ligne
        {  tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture
           // si ce n'est pas la fin de la lecture on appelle un nouvel enregistrement 
           if (i != nbnte) entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
         }
  //    #endif  
     #endif  
       else // cas d'une erreur de lecture
        { cout << "\n erreur de lecture inconnue  ";
          entreePrinc->MessageBuffer("** lecture des données particulières **");
          cout << "SfeMembT::LecDonPart"; 
          Affiche(2);
          Sortie (1);
        }
     };
    // construction du tableau de ddl des noeuds de SfeMembT
    ConstTabDdl(); 
	   // construction du tableau des noeuds du centre
	   t_N_centre.Change_taille(nombre->nbnce);
    for (int i=1;i<= nombre->nbnce;i++)
	    t_N_centre(i)=tab_noeud(i);
  };
  		
// calcul d'un point dans l'élément réel en fonction des coordonnées dans l'élément de référence associé
// temps: indique si l'on veut les coordonnées à t = 0, ou t ou tdt
// 1) cas où l'on utilise la place passée en argument
Coordonnee & SfeMembT::Point_physique(const Coordonnee& ,Coordonnee & co,Enum_dure )
 { cout << "\n methode non encore implantee !!"
        << "\n Coordonnee & SfeMembT::Point_physique(...";
   Sortie(1);         
   // d) retour
   return co;  
  };
  
// 3) cas où l'on veut les coordonnées aux 1, 2 ou trois temps selon la taille du tableau t_co
void SfeMembT::Point_physique(const Coordonnee& ,Tableau <Coordonnee> & )
 { cout << "\n methode non encore implantee !!"
        << "\n void SfeMembT::Point_physique(...";
   Sortie(1);         
  };
 
// vérification que la courbure ne soit pas anormale au temps enu
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1:  si tout est ok,
//        0:  une courbure trop grande a été détecté
// retour = -1 : il y a un pb inconnue dans le calcul
int SfeMembT::Test_courbure_anormale3(Enum_dure enu,double inf_normale)
 { SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
   if (unefois->CalimpPrem == 0)
    { unefois->CalimpPrem = 1;
      Tableau<Enum_variable_metrique> tab(20);
      tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
      tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
      tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
      tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
      tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
      tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
      CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
    };
   // la variable de retour initialisée à 1 par défaut
   int retour = 1;
   // on gère les exceptions éventuelles en mettant le bloc sous surveillance
   try
     { // puis les indices de boucle
       int indice = 1;
       int nbtotsurf = CoSfe->eleCentre->Nbi();
       int nbtotepais =(CoSfe->segment).Nbi();
       for (int nisur =1; nisur<= nbtotsurf; nisur++)  // boucle sur les pt d'integ de surface
        {for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur
          {((DeformationSfe1*) def)->ChangeNumIntegSfe1(nisur,niepais);
           ((DeformationSfe1*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
           retour *= ((DeformationSfe1*) def)->Test_courbure_anormale(enu,inf_normale);
           if (!retour)
             break; // on arrête dès qu'un pb est repéré
          }; // fin boucle sur les points d'intégration
         if (!retour)
          break; // on arrête dès qu'un pb est repéré
        };
     }
   catch (ErrSortieFinale)
        // cas d'une direction voulue vers la sortie
        // on relance l'interuption pour le niveau supérieur
      { ErrSortieFinale toto;
        throw (toto);
      }
   catch ( ... )
      { if (ParaGlob::NiveauImpression() > 3)
          {cout << "\n warning: erreur generee par un element lors du test de verification de la courbure ";
          };
        retour = -1;
      }
   return retour;
 };

// récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu
// ici il s'agit de grandeurs scalaire, 
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
Tableau <double> SfeMembT::Valeur_a_diff_temps(bool absolue,Enum_dure enu_t,const List_io<Ddl_enum_etendu>& enu,int iteg)
    { // appel de la procedure de elem meca
	  int cas;       
      if (!(unefois->dualSortSfe ) && (unefois->CalimpPrem ))
        { cas=1;unefois->dualSortSfe += 1;
         } 
      else if ((unefois->dualSortSfe ) && (unefois->CalimpPrem ))       
         { cas = 11;}
      else if (!(unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))       
        { cas=2;unefois->dualSortSfe += 1;
         }         
      else if ((unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))       
         { cas = 12;}
      // sinon problème
      else
        { cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
               << " pour inialiser les conteneurs des tenseurs  resultats ";
          if (ParaGlob::NiveauImpression() >= 4)     
            cout << "\n cas non prévu, unefois->dualSortSfe= " << unefois->dualSortSfe
               << " unefois->CalimpPrem= " << unefois->CalimpPrem
               << "\n SfeMembT::Valeur_a_diff_temps(Enum_dure enu_t...";
          Sortie(1);
         };       
      return ElemMeca::Valeur_multi(absolue,enu_t,enu,iteg,cas);
    };
    
// récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu
// ici il s'agit de grandeurs tensorielles, le retour s'effectue dans la liste
// de conteneurs quelconque associée
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void SfeMembT::ValTensorielle_a_diff_temps(bool absolue,Enum_dure enu_t,List_io<TypeQuelconque>& enu,int iteg)
    { // appel de la procedure de elem meca
	  int cas;       
      if (!(unefois->dualSortSfe ) && (unefois->CalimpPrem ))
        { cas=1;unefois->dualSortSfe += 1;
         } 
      else if ((unefois->dualSortSfe ) && (unefois->CalimpPrem ))       
         { cas = 11;}
      else if (!(unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))       
        { cas=2;unefois->dualSortSfe += 1;
         }         
      else if ((unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))       
         { cas = 12;}
      // sinon problème
      else
        { cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
               << " pour inialiser les conteneurs des tenseurs  resultats ";
          if (ParaGlob::NiveauImpression() >= 4)     
            cout << "\n cas non prévu, unefois->dualSortSfe= " << unefois->dualSortSfe
               << " unefois->CalimpPrem= " << unefois->CalimpPrem
               << "\n SfeMembT::ValTensorielle_a_diff_temps(Enum_dure enu_t...";
          Sortie(1);
         };       
      ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas);
    };
      
// Calcul du residu local à t ou tdt en fonction du booleen atdt
Vecteur* SfeMembT::CalculResidu (bool atdt,const ParaAlgoControle & pa)
  {   SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
      Tableau <TenseurBB *>& d_epsBB = (CoSfe->d_epsBB);//             "
      energie_totale.Initialisation_differenciee(energie_totale_t); // init de l'énergie totale sur l'élément
      TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
      if (ElemMeca::Bulk_visco_actif())  
          { if (sig_bulk_pourSfe_HH == NULL)  {sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);}
            else if (sig_bulk_pourSfe_HH->Dimension() != sigHH_t_1->Dimension())
               {delete sig_bulk_pourSfe_HH; sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);};           
           };               
      // dimensionnement de la metrique
      if (!atdt)
       {if( !(unefois->CalResPrem_t ))
        { unefois->CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         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 ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
         };}  
      else 
       {if( !(unefois->CalResPrem_tdt ))
        {unefois->CalResPrem_tdt += 1;
         Tableau<Enum_variable_metrique> tab(11);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
         tab(10) = igradVBB_tdt;tab(11) = iVtdt; 
         CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
         };};  
      // dimensionnement du residu
      int nbddl = CoSfe->tab_ddl.NbDdl();
      if ( residu == NULL)
          residu = new Vecteur(nbddl); // cas du premier passage
      else
          for (int i =1; i<= nbddl; i++)  // cas des autres passages
             (*residu)(i) = 0.;             // on initialise a zero
      double jacobien; // def du jacobien
      int indice = 1;
      int nbtotsurf =CoSfe->eleCentre->Nbi();        
      int nbtotepais =(CoSfe->segment).Nbi();        
      // préparation du calcul de l'épaisseur finale (qui dépend de tous les points d'intégration)
      double somme_epaisseurs_finales = 0;        
      // on utilise une épaisseur intermédiaire pour la métrique
      // de manière à utiliser l'épaisseur de l'incrément précédent
      Epai  epais_pour_metrique;

      // on doit boucler à l'extérieur sur les pt de surface et à l'intérieur sur les pt d'épaisseurs
      // pour profiter du fait que une fois un pt calculer, on n'a pas à recalculer la courbure !!!!!
      for (int nisur =1; nisur<= nbtotsurf; nisur++)  // boucle sur les pt d'integ de surface
       for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur


        { if ((donnee_specif.epais != NULL) // on change les épaisseurs si elles sont stockées
              && (!Loi_rien(loiComp->Id_comport()))
             )
            { epais_pour_metrique = *(donnee_specif.epais);
              epais_pour_metrique.epaisseur_tdt = epais_pour_metrique.epaisseur_t;
//              ((DeformationSfe1*) def)->Change_epaisseur(*donnee_specif.epais);  // au niveau de l'élément
              ((DeformationSfe1*) def)->Change_epaisseur(epais_pour_metrique);
            };
          ((DeformationSfe1*) def)->ChangeNumIntegSfe1(nisur,niepais);
          ((DeformationSfe1*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
          
          PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
          TenseurHH & sigHH = *(ptIntegMeca.SigHH());
          TenseurBB & DepsBB_ = *(ptIntegMeca.DepsBB());
          EnergieMeca& energ =  tab_energ(indice);
          EnergieMeca& energ_t =  tab_energ_t(indice);
          double poids = CoSfe->eleCentre->Wi(nisur) * (CoSfe->segment).Wi(niepais) ;
          CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
          if (loiTP != NULL) {sTP = tabSaveTP(indice);};         // au pt d'integ si elles existes 
          if ((loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
              || (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);
          double epaisseur_tdt = 0.; // ne sert que pour les éléments 2D avec épaisseurs stockées à l'élément
          double epaisseur_pris_en_compte = 0.; // celle qui sera réellement prise en compte

          // -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
          try // ce qui permet de donner les numéros d'éléments et de pti
           {if (atdt)
               { const Met_abstraite::Expli_t_tdt& ex = loiComp->Cal_explicit_tdt 
                               (tabSaveDon(indice), *def,CoSfe->tab_ddl,ptIntegMeca,d_epsBB,jacobien
                                 ,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
                 // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
                 if (Bulk_visco_actif()) 
                    ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_ ,(*sig_bulk_pourSfe_HH));
                 // on tient compte de la variation d'épaisseur si l'épaisseur est stockée dans l'élément
                 // et que l'on n'a pas une loi rien
                 if ((donnee_specif.epais != NULL) && (!Loi_rien(loiComp->Id_comport())))
                  { switch (loiComp->Comportement_3D_CP_DP_1D())
                     { case COMP_CONTRAINTES_PLANES :
                        // en contrainte plane on recalcule l'épaisseur
                        { CalEpaisseurAtdtExp_et_vol_pti(donnee_specif.epais->epaisseur0,pa,donnee_specif.epais->epaisseur_t
                                                     ,ptIntegMeca,epaisseur_tdt,ex);
//------ le 7 nov 2016: changement car je ne veux pas utiliser la variation d'épaisseur
// qui est calculée dans les lois. C'est sans doute possible, bien qu'avec les lois additives et mélange
// ce n'est pas sûr que cela veuille dire quelque chose. Donc on essaie de sans passer
// De plus c'est cohérent avec les éléments membranes ...

// ---- remodif le 1 avril 2017 car la méthode précédente fait que l'élément converge très mal !!
// je ne sais pas trop pourquoi (avec une loi de déformation plane, cela converge immédiatement !)
//                   // en fait, on récupère la variation d'épaisseur calculée via la condition de contrainte plane
//                   // on utilise uniquement la variation d'épaisseur calculée à la surface médiane ??? pour l'instant
//                   epaisseur_tdt = loiComp->HsurH0(tabSaveDon(nisur)) * donnee_specif.epais->epaisseur0;

// il faudrait peut-être utiliser finalement la variation d'épaisseur calculée par la loi, comme pour les éléments
// membranes qui ont été modifiées pour tenir compte de la différence de def entre géométrique et la def méca
// Par exemple lors de critère de rupture, on aura le même pb qu'avec les plis
//*** à faire
                          epaisseur_pris_en_compte = donnee_specif.epais->epaisseur_tdt;
                          break;
                        }
                      case COMP_DEFORMATIONS_PLANES :
                        // en deformation plane, l'épaisseur ne change pas
                        {  epaisseur_tdt = donnee_specif.epais->epaisseur0;
                           ptIntegMeca.Volume_pti() *= donnee_specif.epais->epaisseur0;
                           epaisseur_pris_en_compte = donnee_specif.epais->epaisseur0;
                           break;
                        }
                      default :
                      cout << "\nErreur : cas de loi de comportement : " << Nom_comp((loiComp->Id_comport())) 
                         << " non prevu pour le calcul le l'epaisseur avec l'element SFE ";
                        cout << "\n SfeMembT::Calcul_implicit(.... \n";
                        Sortie(1);
                     };
                  };
                 somme_epaisseurs_finales += epaisseur_tdt; //pour le calcul final
                }               
            else   
               { const Met_abstraite::Expli& ex = loiComp->Cal_explicit_t
                              (tabSaveDon(indice), *def,CoSfe->tab_ddl,ptIntegMeca,d_epsBB,jacobien
                               ,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
                 // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
                 if (Bulk_visco_actif())
                   ModifContrainteBulk(TEMPS_t,*ex.gijHH_t,sigHH,DepsBB_ ,(*sig_bulk_pourSfe_HH));
                 // pour l'instant on expédie un message, car il y a des chances que cette partie disparaisse !!
                 cout << "\n attention *** calcul non operationnelle en explicite t, on ne tient pas compte de la"
                     << " variation d'epaisseur !! (a modifier ) "; 
               };
           }
          catch (ErrNonConvergence_loiDeComportement excep)
           // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
           // d'une loi de comportement incrémentale
            { cout << "\n erreur de loi de comportement element= " << this->Num_elt()
                   << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw excep;
            }
          catch (ErrSortieFinale)
               // cas d'une direction voulue vers la sortie
               // on relance l'interuption pour le niveau supérieur
             { ErrSortieFinale toto;
               throw (toto);
             }
          catch ( ... )
            { cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
                   << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw (Err_inconnue_ElemMeca());
            };

          if ((jacobien <= 0.)|| (std::isinf(jacobien))||( std::isnan(jacobien))) // vérif du jacobien
            { if ((std::isinf(jacobien))||( std::isnan(jacobien)))
                // on met un message quelque soit le niveau d'impression
                { cout << "\n **********  attention on a trouve un jacobien infini ou nan !! = ("<<jacobien<<")********* " << endl; };
              if (ParaGlob::NiveauImpression() >= 1)
               { cout << "\n **********  attention on a trouve un jacobien  negatif!! ("<<jacobien<<")********* " << endl; };
            };
          if ((jacobien <= 0.) && (pa.JabobienNegatif() != 0)) // vérif du jacobien  et traitement adéquate si besoin
              { switch (pa.JabobienNegatif())
                { case 1:
                  { cout << "\n on annulle sa contribution. Element nb: " <<  Num_elt() << " ni_surf: " << nisur
                         << " ni_ep: " << niepais;
                    break;
                   }
                // dans les autres cas on ne fait rien   
                };       
              }
          else
           { // on prend en compte toutes les  variations
            // dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
            // épaisseur_tdt ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
            // pour éviter des tests et une recopie du code on laisse epaisseur_tdt, avec une valeur 1
            // ce qui fait qu'elle n'intervient pas 
            if (donnee_specif.epais == NULL) epaisseur_pris_en_compte = 1.;
            // il ne faut pas interpoler l'énergie à t, car elle était associée à un volume qui a changé
            // donc il ne faut considérer que l'accroissement d'énergie
            energie_totale.Ajout_differenciee(energ,energ_t,(poids * jacobien * epaisseur_pris_en_compte));
//            energie_totale += (energ-energ_t) * (poids * jacobien * epaisseur_tdt);
//            energie_totale += energ * (poids * jacobien * epaisseur_tdt);

            //test si le jacobien due aux gi finaux est très différent du jacobien de la facette
            //test si le jacobien due aux gi finaux est très différent du jacobien de la facette
            Delta_Jacobien_anormal(tabSaveDefDon(indice),nisur,niepais);
            
            for (int j =1; j<= nbddl; j++)
				          {int jloc =  CoSfe->nMetVTab_ddl(j); // adressage indirecte
                 (*residu)(jloc) -=  (sigHH && (*(d_epsBB(j)))) * jacobien * poids * epaisseur_pris_en_compte;
                 volume += jacobien * poids  * epaisseur_tdt; // là on prend la vrai épaisseur
				          };
           };
          if (ElemMeca::Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
           { sigHH -= (*sig_bulk_pourSfe_HH); }  // ceci pour la sauvegarde  ou autre utilisation       
 
        }; // fin boucle sur les points d'intégration
      // on calcul éventuellement l'épaisseur moyenne finale si c'est un élément 2D
      if ((donnee_specif.epais != NULL)&&(!Loi_rien(loiComp->Id_comport()))) // cas 2D
         donnee_specif.epais->epaisseur_tdt = somme_epaisseurs_finales / (indice - 1);  
      // liberation des tenseurs intermediaires
      LibereTenseur();   
      return residu;  
     };

// Calcul du residu local et de la raideur locale,
//  pour le schema implicite
Element::ResRaid  SfeMembT::Calcul_implicit (const ParaAlgoControle & pa)
 { SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
   Tableau <TenseurBB *>& d_epsBB = (CoSfe->d_epsBB);//             "
   Tableau <TenseurHH *>& d_sigHH = (CoSfe->d_sigHH);//             "
   energie_totale.Initialisation_differenciee(energie_totale_t); // init de l'énergie totale sur l'élément
   E_elem_bulk_tdt = E_elem_bulk_t; P_elem_bulk = 0.; // init pour l'énergie et la puissance associées au bulk
   // init éventuel de la contrainte de bulk viscosity
   TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
   if (ElemMeca::Bulk_visco_actif())  
    { if (sig_bulk_pourSfe_HH == NULL)  {sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);}
      else if (sig_bulk_pourSfe_HH->Dimension() != sigHH_t_1->Dimension())
               {delete sig_bulk_pourSfe_HH; sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);};           
    };
    
   // init éventuel des intégrales de volume et temps
   Init_Integ_vol_et_temps();
   // --- init pour l'init de la stabilisation de membrane-biel éventuelle
   const Met_abstraite::Impli* ex_final=NULL; // contiendra après les boucles la dernière métrique
   Mise_a_jour_A_calculer_force_stab(); // permettra ensuite de savoir où le calcul doit être fait

   bool cald_Dvirtuelle = false;
   if (unefois->CalimpPrem == 0)
    { unefois->CalimpPrem = 1;
      Tableau<Enum_variable_metrique> tab(20);
      tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
      tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
      tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
      tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
      tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
      tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
      CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
      // on ne calcul la dérivée de la déformation virtuelle qu'une fois
      // car elle est constante dans le temps et indépendante des coordonnées
      cald_Dvirtuelle=true; 
    };
   // dimensionnement du residu
   int nbddl = CoSfe->tab_ddl.NbDdl();
   if ( residu == NULL)
       residu = new Vecteur(nbddl); // cas du premier passage
   else
       for (int i =1; i<= nbddl; i++)  // cas des autres passages
          (*residu)(i) = 0.;             // on initialise a zero
   // dimensionnement de la raideur
   if ( raideur == NULL)
       raideur = new Mat_pleine(nbddl,nbddl); // cas du premier passage
   else
       for (int i =1; i<= nbddl; i++)  // cas des autres passages
         for (int j=1; j<= nbddl; j++) //
           (*raideur)(i,j) = 0.;       // on initialise a zero
   double jacobien; // def du jacobien
   Vecteur d_jacobien_tdt; 
   int indice = 1;
   int nbtotsurf = CoSfe->eleCentre->Nbi();        
   int nbtotepais =(CoSfe->segment).Nbi(); 
   // préparation du calcul de l'épaisseur finale (qui dépend de tous les points d'intégration)
   double somme_epaisseurs_finales = 0;
   // on utilise une épaisseur intermédiaire pour la métrique
   // de manière à utiliser l'épaisseur de l'incrément précédent
   Epai  epais_pour_metrique;
//--debug
//if (num_elt ==  89)
// { cout << "\n SfeMembT::Calcul_implicit debug ";
//   cout << endl;
// };

   // on doit boucler à l'extérieur sur les pt de surface et à l'intérieur sur les pt d'épaisseurs
   // pour profiter du fait que une fois un pt calculer, on n'a pas à recalculer la courbure !!!!!
// ne sert plus**   int nbint = nbtotsurf*nbtotepais; // nb total de pti
   for (int nisur =1; nisur<= nbtotsurf; nisur++)  // boucle sur les pt d'integ de surface
    for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur
     {if ((donnee_specif.epais != NULL) // on change les épaisseurs si elles sont stockées
          && (!Loi_rien(loiComp->Id_comport()))
          )
            { epais_pour_metrique = *(donnee_specif.epais);
              epais_pour_metrique.epaisseur_tdt = epais_pour_metrique.epaisseur_t;
//              ((DeformationSfe1*) def)->Change_epaisseur(*donnee_specif.epais);  // au niveau de l'élément
              ((DeformationSfe1*) def)->Change_epaisseur(epais_pour_metrique);
            };
      ((DeformationSfe1*) def)->ChangeNumIntegSfe1(nisur,niepais);
      ((DeformationSfe1*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
      
      
      PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
      TenseurHH & sigHH = *(ptIntegMeca.SigHH());
      TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
      EnergieMeca& energ =  tab_energ(indice);
      EnergieMeca& energ_t =  tab_energ_t(indice);
      double poids = (CoSfe->eleCentre->Wi(nisur) * (CoSfe->segment).Wi(niepais)) * 0.5 ;
      CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
      if (loiTP != NULL) {sTP = tabSaveTP(indice);};         // au pt d'integ si elles existes 
      if ((loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
         || (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
         ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);

      // -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
      const Met_abstraite::Impli* ex_pointe = NULL; // c'est pour construire ex
      try // ce qui permet de donner les numéros d'éléments et de pti
       { ex_pointe = &loiComp->Cal_implicit(tabSaveDon(indice), *def,CoSfe->tab_ddl
                                               ,ptIntegMeca,d_epsBB,jacobien,d_jacobien_tdt
                                               ,d_sigHH,pa,sTP,loiTP,dilatation,energ,energ_t
                                               ,premier_calcul_meca_impli_expli);
       }
      catch (ErrNonConvergence_loiDeComportement excep)
       // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
       // d'une loi de comportement incrémentale
        { cout << "\n erreur de loi de comportement element= " << this->Num_elt()
                 << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
          throw excep;
        }
      catch (ErrSortieFinale)
           // cas d'une direction voulue vers la sortie
           // on relance l'interuption pour le niveau supérieur
         { ErrSortieFinale toto;
           throw (toto);
         }
      catch ( ... )
        { cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
                 << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
          throw (Err_inconnue_ElemMeca());
        };
      const Met_abstraite::Impli& ex=(*ex_pointe);
      ex_final = ex_pointe; // stockage pour la stabilisation éventuelle de membraneBiel

      // on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
      Tableau2 <TenseurBB *>& d2_epsBB =  (CoSfe->d2_epsBB)(indice); // pour simplifier
      
      if (cald_Dvirtuelle)  
           def->Cal_var_def_virtuelle(false,d2_epsBB);
      // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
      if (Bulk_visco_actif()) 
           ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_tdt,(*sig_bulk_pourSfe_HH) );
      double rap=1.; // pour le calcul éventuel du rapport de jacobien actuel
      if (pa.MaxiVarJacobien() > 1.) // cas de la demande du calcul du rapport de jacobien 
         {if (Dabs(jacobien) > (*(ex.jacobien_0))) {rap = Dabs(jacobien) / (*(ex.jacobien_0));}
          else                                     {rap = (*(ex.jacobien_0)) / Dabs(jacobien) ;}
          };  
      // on tient compte de la variation d'épaisseur pour ce point d'integ
      double epaisseur_tdt = 0.; // ne sert que pour les éléments 2D avec épaisseurs stockées à l'élément
      double epaisseur_pris_en_compte = 0.; // celle qui sera réellement prise en compte
      if ((donnee_specif.epais != NULL) && (!Loi_rien(loiComp->Id_comport())) )
 		     { switch (loiComp->Comportement_3D_CP_DP_1D()) 
            { case COMP_CONTRAINTES_PLANES :
					            {
                   // en contrainte plane on recalcule l'épaisseur
                   CalEpaisseurAtdtImp_et_vol_pti(donnee_specif.epais->epaisseur0
                                     ,pa,donnee_specif.epais->epaisseur_t
                                     ,ptIntegMeca,epaisseur_tdt,ex);
//------ le 7 nov 2016: changement car je ne veux pas utiliser la variation d'épaisseur
// qui est calculée dans les lois. C'est sans doute possible, bien qu'avec les lois additives et mélange
// ce n'est pas sûr que cela veuille dire quelque chose. Donc on essaie de sans passer
// De plus c'est cohérent avec les éléments membranes ...

// ---- remodif le 1 avril 2017 car la méthode précédente fait que l'élément converge très mal !!
// je ne sais pas trop pourquoi (avec une loi de déformation plane, cela converge immédiatement !)
//                   // en fait, on récupère la variation d'épaisseur calculée via la condition de contrainte plane
//                   // on utilise uniquement la variation d'épaisseur calculée à la surface médiane ??? pour l'instant
//                   epaisseur_tdt = loiComp->HsurH0(tabSaveDon(nisur)) * donnee_specif.epais->epaisseur0;
// ************$la ligne qui suit est uniquement pour faire passer qq calcul ---> mais c'est faux !!!
// ************$la ligne qui suit est uniquement pour faire passer qq calcul ---> mais c'est faux !!!
// ************$la ligne qui suit est uniquement pour faire passer qq calcul ---> mais c'est faux !!!
                   
//     epaisseur_tdt = donnee_specif.epais->epaisseur0;
 // on utilise l'épaisseur de l'incrément précédent
                   epaisseur_pris_en_compte = donnee_specif.epais->epaisseur_tdt;
					              break;
                 }
              case COMP_DEFORMATIONS_PLANES :
                 // en deformation plane, l'épaisseur ne change pas
                 {  epaisseur_tdt = donnee_specif.epais->epaisseur0;
                    ptIntegMeca.Volume_pti() *= donnee_specif.epais->epaisseur0;
                    epaisseur_pris_en_compte = epaisseur_tdt;
                    break;
                 }
              default :
                 cout << "\nErreur : cas de loi de comportement : " << Nom_comp((loiComp->Id_comport())) 
                    << " non prevu pour le calcul le l'epaisseur avec l'element SFE ";
                 cout << "\n SfeMembT::Calcul_implicit(.... \n";
                 Sortie(1);
            };
        };
		    somme_epaisseurs_finales += epaisseur_tdt; //pour le calcul final

//--------- pour le débug --------
//      Deformation & defS = *defSurf(1); // pour simplifier l'écriture
//      int nbddls = ddlS.NbDdl();


//------- pour le débug fin -------


      if ((jacobien <= 0.)|| (std::isinf(jacobien))||( std::isnan(jacobien))) // vérif du jacobien
          { if ((std::isinf(jacobien))||( std::isnan(jacobien)))
                 // on met un message quelque soit le niveau d'impression
                 { cout << "\n **********  attention on a trouve un jacobien infini ou nan !! = ("<<jacobien<<")********* " << endl; };
               if (ParaGlob::NiveauImpression() >= 1)
                 { cout << "\n **********  attention on a trouve un jacobien  negatif!! ("<<jacobien<<")********* " << endl; };
          };
      if ((jacobien <= 0.) && (pa.JabobienNegatif()!= 0)) // vérif du jacobien  et traitement adéquate si besoin
          { cout << "\n ParaGlob::NiveauImpression() = "<<ParaGlob::NiveauImpression() << " jacobien= " << jacobien << endl;
               switch (pa.JabobienNegatif())
               { case 1:   
                   { cout << "\n on annulle sa contribution. Element nb: " <<  Num_elt()  << " ni_surf: " << nisur
                          << " ni_ep: " << niepais;
                     break;
                    }
                 case 2:
                   { // on génère une exception
                     throw ErrJacobienNegatif_ElemMeca();break;
                    }
                 // dans les autres cas on ne fait rien   
               };       
          }
      else if ((pa.MaxiVarJacobien() > 1.) && ( rap > pa.MaxiVarJacobien()))
          { if (ParaGlob::NiveauImpression() >= 3)
                { cout << "\n ***  attention la variation maximale du jacobien est atteinte  !! *** " 
                       << "\n ele= "<<  Num_elt() << " ni_surf: " << nisur
                       << " ni_ep: " << niepais << " jacobien= " << jacobien << " jacobien_0= " 
                       << (*(ex.jacobien_0)) << endl;
                };
               // on génère une exception
               throw ErrVarJacobienMini_ElemMeca();break;
          }
      else 
       { // on prend en compte toutes les  variations
       
         // dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
         // épaisseur_tdt ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
         // pour éviter des tests et une recopie du code on laisse epaisseur_tdt, avec une valeur 1
         // ce qui fait qu'elle n'intervient pas
         if (donnee_specif.epais == NULL) epaisseur_pris_en_compte = 1.;
         
         // on continue que s'il s'agit d'une loi différente de rien
         if (!Loi_rien(loiComp->Id_comport()))
           {// il ne faut pas interpoler l'énergie à t, car elle était associée à un volume qui a changé
            // donc il ne faut considérer que l'accroissement d'énergie
            energie_totale.Ajout_differenciee(energ,energ_t,(poids * jacobien * epaisseur_pris_en_compte));
//         energie_totale += (energ-energ_t) * (poids * jacobien * epaisseur_tdt);
//         energie_totale += energ * (poids * jacobien * epaisseur_tdt);
         // l'ordre des ddl vient de la métrique: tous le xi, puis les épaisseurs, par contre dans le stockage
         // local on doit avoir noeud1: ddlxi puis epais, noeud2: idem etc pour les noeuds centraux, puis ddl xi
         // pour les noeuds externes. Si pas de ddl epais, tous les noeuds sont traités de la même manière
         // on utilise un adressage indirecte
////--- debug
//if (Dabs(epaisseur_pris_en_compte - donnee_specif.epais->epaisseur0) > 0.01)
// {cout << "\n debug SfeMembT::Calcul_implicit: epai= "
//       << epaisseur_pris_en_compte << " epai_0= "<< donnee_specif.epais->epaisseur0 << endl;
// };
//// --- fin debug

            for (int j =1; j<= nbddl; j++)  // 1ere boucle sur les ddl
              { int jloc =  CoSfe->nMetVTab_ddl(j); // adressage indirecte           
               (*residu)(jloc) -= ((*(d_epsBB(j))) && sigHH ) * poids * jacobien* epaisseur_pris_en_compte;
/*//--debug
if ((*residu)(jloc) > 10.)
 { cout << "\n SfeMembT::Calcul_implicit debug ";
   cout << "\n (*ex.giB_0)(1): " ;(*ex.giB_0)(1).Affiche();
   cout << " (*ex.giB_0)(2): " ;(*ex.giB_0)(2).Affiche();
   cout << " (*ex.giH_0)(1): " ;(*ex.giH_0)(1).Affiche();
   cout << " (*ex.giH_0)(2): " ;(*ex.giH_0)(2).Affiche();
   cout << "\n (*ex.giBB_0): " ;(*ex.gijBB_0).Ecriture(cout);
   cout << "\n (*ex.giHH_0): " ;(*ex.gijHH_0).Ecriture(cout);
  
   cout << "\n (*(d_epsBB(j))): " ;(*(d_epsBB(j))).Ecriture(cout);
   cout << "\n sigHH: " ;sigHH.Ecriture(cout);
   cout << endl;

   ex_pointe = &loiComp->Cal_implicit(tabSaveDon(indice), *def,CoSfe->tab_ddl
   ,ptIntegMeca,d_epsBB,jacobien,d_jacobien_tdt
   ,d_sigHH,pa,sTP,loiTP,dilatation,energ,energ_t
   ,premier_calcul_meca_impli_expli);
  
   Sortie(1);
 };
//--fin debug */
               //test si le jacobien due aux gi finaux est très différent du jacobien de la facette
               //test si le jacobien due aux gi finaux est très différent du jacobien de la facette
               Delta_Jacobien_anormal(tabSaveDefDon(indice),nisur,niepais);

               volume += jacobien * poids * epaisseur_tdt; // là on prend la vraie épaisseur
               for (int i =1; i<= nbddl; i++)  // 2ere boucle sur les ddl
                  { int iloc =  CoSfe->nMetVTab_ddl(i); // adressage indirecte
                    (*raideur)(jloc,iloc) +=  ( ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien 
                                       +  (sigHH && (*((CoSfe->d2_epsBB)(indice)(j,i)))) * jacobien
                                       +  (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i)
                                       ) * poids * epaisseur_pris_en_compte;

// on utilise toujour la variation de D, donc c'est maintenant directement intégré dans l'expression précédente, et la suite est commentée
//                   // on intègre la variation seconde si c'est demandé
//                   if (pa.Var_D())
//                      (*raideur)(jloc,iloc) +=  (sigHH && (*((CoSfe->d2_epsBB)(indice)(j,i)))) 
//                                            * jacobien * poids * epaisseur_tdt;                    
// ---- pour test débug
//                    if (i < 10) 
//                       (*raideur)(j,i) +=  (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i) * poids * donnee_specif.epaisseur_tdt;                 
                  };
              };
            };
           // --- cas éventuel d'une stabilisation membrane-biel ---
           // ici il s'agit de la contribution  précise à chaque pti
           if (pt_StabMembBiel != NULL)
             if (!(pt_StabMembBiel->Aa_calculer()))
               { double poids_vol = poids * epaisseur_pris_en_compte;
                 Cal_implicit_StabMembBiel
                    (nisur,*ex_final,nbtotsurf,poids_vol,unefois->doCoMemb->noeud_a_prendre_en_compte);
               };
          };// fin du if sur la valeur non négative du jacobien
////--debug
//for (int i = 1; i<= 18; i++)
// { if ((*residu)(i) != 0. )
//     { cout << "\n SfeMembT::Calcul_implicit debug ";
//       cout << endl;
//     };
////if (num_elt ==  176)
//// { cout << "\n SfeMembT::Calcul_implicit debug ";
////   cout << endl;
// };
////--fin debug
      if (ElemMeca::Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
        { sigHH -= (*sig_bulk_pourSfe_HH); }  // ceci pour la sauvegarde  ou autre utilisation
 
     }; // fin boucle sur les points d'intégration
  #ifdef MISE_AU_POINT
  if (ParaGlob::NiveauImpression() > 9)
   { cout << "\n Raideur et second membre locaux: ? (o/n) "; string rep(" ");
     // procédure de lecture avec prise en charge d'un retour chariot
     rep = lect_return_defaut(false,"n");

     if ((rep == "0")||(rep == "o"))
       {cout << "\n Raideur et second membre locaux: elem= " << Num_elt_const()
             << ", maillage= " << Num_maillage();
        cout << "\n raideur: ";
        raideur->Affiche();
        cout << "\n second membre: " << (*residu);
       };
   };
  #endif
  // --- intervention dans le cas d'une stabilisation d'hourglass
  // stabilisation pour un calcul implicit
  if (type_stabHourglass)
      Cal_implicit_hourglass();
  // --- cas éventuel d'une stabilisation membrane-biel ---
  // ici il s'agit soit du calcul approché d'initialisation, soit de la fin du calcul après la boucle
  // modif: 10 janvier 2019 non c'est le calcul correct une fois la raideur calculée
  if (pt_StabMembBiel != NULL)
   {if (pt_StabMembBiel->Aa_calculer())
          {Cal_implicit_StabMembBiel
             (0,*ex_final, nbtotsurf,volume,unefois->doCoMemb->noeud_a_prendre_en_compte);}
    else  {Cal_implicit_StabMembBiel
             (-1,*ex_final, nbtotsurf,volume,unefois->doCoMemb->noeud_a_prendre_en_compte);} ;
   };
  #ifdef MISE_AU_POINT
  if (ParaGlob::NiveauImpression() > 9)
   { if ((type_stabHourglass)||(pt_StabMembBiel != NULL))
      {cout << "\n apres stabilisation: Raideur et second membre locaux: ? (o/n) "; string rep(" ");
       // procédure de lecture avec prise en charge d'un retour chariot
       rep = lect_return_defaut(false,"n");

       if ((rep == "0")||(rep == "o"))
         {cout << "\n Raideur et second membre locaux: elem= " << Num_elt_const()
               << ", maillage= " << Num_maillage();
          cout << "\n raideur: ";
          raideur->Affiche();
          cout << "\n second membre: " << (*residu);
         };
     };
   };
  #endif

////--debug
//for (int i = 1; i<= 18; i++)
// { if ((*residu)(i) != 0. )
//     { cout << "\n SfeMembT::Calcul_implicit debug ";
//       cout << endl;
//     };
////if (num_elt ==  176)
//// { cout << "\n SfeMembT::Calcul_implicit debug ";
////   cout << endl;
// };
////--fin debug


   // on calcul éventuellement l'épaisseur moyenne finale si c'est un élément 2D
   if ((donnee_specif.epais != NULL) && (!Loi_rien(loiComp->Id_comport())) )// cas 2D
      donnee_specif.epais->epaisseur_tdt = somme_epaisseurs_finales / (indice - 1);  
   // liberation des tenseurs intermediaires
   LibereTenseur();   
   Element::ResRaid el;
   el.res = residu; 
   el.raid = raideur;       
   return el;
  };
		
// Calcul de la matrice masse pour l'élément
// celle-ci est calculée uniquement pour l'élément centrale et en considérant
// que la masse volumique est constante dans l'épaisseur
Mat_pleine * SfeMembT::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse) 
 {    SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
      // dimensionement de la métrique si nécessaire
      if (!(unefois->CalDynamique ))
       { unefois->CalDynamique += 1;
         Tableau<Enum_variable_metrique> tab(5);
         tab(1) = igiB_0;  tab(2) = igiB_tdt; tab(3) = igijBB_0;
         tab(4) = igijBB_tdt; tab(5) = igradVmoyBB_t; 
         CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
         // on vérifie le bon dimensionnement de la matrice
         if (type_calcul_masse == MASSE_CONSISTANTE)
           // dans le cas où la masse est consistante il faut la redimensionner
          { int nbddl = CoSfe->tab_ddl.NbDdl();
            (CoSfe->matrice_masse).Initialise (nbddl,nbddl,0.);
           } 
        };
      if ((donnee_specif.epais != NULL) // on change les épaisseurs si elles sont stockées
          && (!Loi_rien(loiComp->Id_comport()))
          )
              ((DeformationSfe1*) def)->Change_epaisseur(*donnee_specif.epais);  // au niveau de l'élément
       // appel de la routine générale 
////************ non, ce n'est pas correcte pour les éléments 3D, il faut revoir le calcul de la masse
//// en utilisant d(OM)/d_ddl au lieu de phi^r       
       ElemMeca::Cal_Mat_masse (CoSfe->tab_ddl,type_calcul_masse,
                            nombre->nbisMas,CoSfe->eleMas->TaPhi(),nombre->nbnce
                            ,CoSfe->eleMas->TaWi());
       if ((donnee_specif.epais != NULL) // prise en compte de l'épaisseur finale, si elle est stockée
           && (!Loi_rien(loiComp->Id_comport()))
          )
           (*mat_masse) *= donnee_specif.epais->epaisseur_tdt; // au niveau de l'élément, sinon on est en 3D, donc c'est
                                                        // intégré dans le jacobien                    
       return mat_masse;                     
  };
		
//============= lecture écriture dans base info ==========
	
// 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 SfeMembT::Lecture_base_info
    (ifstream& ent,const Tableau<Noeud  *> * tabMaillageNoeud,const int cas) 
{// tout d'abord appel de la lecture de la classe elem_meca
 ElemMeca::Lecture_bas_inf(ent,tabMaillageNoeud,cas);
 // traitement du cas particulier de l'élément centrale
 string nom;
 switch (cas)
  { case 1 : // ------- on récupère tout -------------------------
     { // construction du tableau de ddl des noeuds de l'élément centrale
       ConstTabDdl(); 
       // récup contraintes et déformation
       lesPtMecaInt.Lecture_base_info(ent,cas);
       
       // les données spécifiques
       ent >> nom;
       if (nom == "epaisStockeDansElement")
         { // cas où les épaisseurs sont stockées dans l'élément
           if (donnee_specif.epais == NULL) donnee_specif.epais = new Epai;
           ent >>  nom >> donnee_specif.epais->epaisseur0
               >>  nom >> donnee_specif.epais->epaisseur_t
               >>  nom >> donnee_specif.epais->epaisseur_tdt;
         }
       else 
       	 // cas d'un stockage au niveau des noeuds
       	 { if (donnee_specif.epais != NULL) delete donnee_specif.epais;};
       // CL éventuelles : 
       ent >> nom ;
       bool existe_CL; ent >> existe_CL;
       if (existe_CL)
	        {if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL)))
           { // si les tableaux n'existent pas on les crée
             if (areteTypeCL == (&(unefois->doCoMemb->tabType_rienCL)))
                 areteTypeCL = new Tableau <EnuTypeCL> (3);
             if (vplan == (&(unefois->doCoMemb->vplan_rien)))
                 vplan = new Tableau <Coordonnee3> (3);
             // et on récupère les infos
             ent >> nom; (*areteTypeCL)(1) = Id_nomTypeCL(nom); ent >> (*vplan)(1);
             ent >> nom; (*areteTypeCL)(2) = Id_nomTypeCL(nom); ent >> (*vplan)(2);
             ent >> nom; (*areteTypeCL)(3) = Id_nomTypeCL(nom); ent >> (*vplan)(3);
           }
         }
       else
        {// on supprime les tableaux locaux éventuels et on pointe vers les valeurs par défaut
         if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL))) delete areteTypeCL;
         areteTypeCL = (&(unefois->doCoMemb->tabType_rienCL));
         if (vplan != (&(unefois->doCoMemb->vplan_rien))) delete vplan;
         vplan = (&(unefois->doCoMemb->vplan_rien));
        };	
       break;
	    }
    case  2 : // ----------- lecture uniquement de se qui varie --------------------
     { // récup contraintes et déformation
       lesPtMecaInt.Lecture_base_info(ent,cas);
       // les données spécifiques
       bool avecEpaisDansElement = true;
       ent >> avecEpaisDansElement;
       if (avecEpaisDansElement)
       	{ if (donnee_specif.epais == NULL) donnee_specif.epais = new Epai;
       	  ent >> nom >> donnee_specif.epais->epaisseur_tdt;	
       	  donnee_specif.epais->epaisseur_t = donnee_specif.epais->epaisseur_tdt;
       	};

       // CL éventuelles : on ne les lit pas car elles sont généré éventuellement 
       // à chaque fois que l'on met les CL, sinon elles ne changent pas 
	      break;
     }
    default :
     { cout << "\nErreur : valeur incorrecte du type de lecture !\n";
       cout << "SfeMembT::::Lecture_base_info(ifstream& ent,const "
            << "Tableau<Noeud  *> * tabMaillageNoeud,const int cas)"
            << " cas= " << cas << endl;
       Sortie(1);
     }
   }
  };
  
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void SfeMembT::Ecriture_base_info(ofstream& sort,const int cas) 
{// tout d'abord appel de l'écriture de la classe elem_meca
 ElemMeca::Ecriture_bas_inf(sort,cas);
 // traitement du cas particulier de l'élément central
 switch (cas)
  { case 1 : // ------- on sauvegarde tout -------------------------
     { 
       // des tenseurs déformation et contrainte,
       lesPtMecaInt.Ecriture_base_info(sort,cas);
       // les données spécifiques
       if (donnee_specif.epais != NULL)
         { sort << "\n epaisStockeDansElement " << " epaisseur0= " << donnee_specif.epais->epaisseur0
                << " epaisseur_t= " << donnee_specif.epais->epaisseur_t
                << " epaisseur_tdt= " << donnee_specif.epais->epaisseur_tdt;
         }
       else { sort << "\n epaisStockeAuxNoeuds ";	};
       // CL éventuelles : 
       sort << " CL= ";
	      if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL)))
        { // si ça ne pointe pas sur la valeur par défaut, alors c'est qu'il y a
          // des CL particulières,
          sort << " 1 ";
          sort << NomTypeCL((*areteTypeCL)(1)) << " " << (*vplan)(1)
               << NomTypeCL((*areteTypeCL)(2)) << " " << (*vplan)(2)
               << NomTypeCL((*areteTypeCL)(3)) << " " << (*vplan)(3);
        }
	      else
   	   	{sort << " 0 ";
	   	   };
            
	      break;
	    }
    case  2 : // ----------- sauvegarde uniquement de se qui varie --------------------
     { // des tenseurs déformation et contrainte, 
       lesPtMecaInt.Ecriture_base_info(sort,cas);
       // les données spécifiques
       if (donnee_specif.epais != NULL)
        { sort << "\n 1 epaisseur_tdt= " << donnee_specif.epais->epaisseur_tdt;}
       else
        { sort << "\n 0 ";}; 
       // CL éventuelles : on ne les sort pas car elles sont généré éventuellement 
       // à chaque fois que l'on met les CL, sinon elles ne changent pas 
	      break;
	    }
    default :
	    { cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
	      cout << "SfeMembT::::Ecriture_base_info(ofstream& sort,const int cas)"
			         << " cas= " << cas << endl;
	      Sortie(1);
	    }
  };
};

//------- calcul d'erreur, remontée des contraintes -------------------

// calcul du résidu et de la matrice de raideur pour le calcul d'erreur
// en fait on ne calculera l'erreur que sur la contrainte sur toute la surface
// mais sur le premier pt d'integ d'épaisseur
Element::Er_ResRaid    SfeMembT::ContrainteAuNoeud_ResRaid()
    {  SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
       // dimensionnement de la metrique
       if(!( unefois->CalResPrem_t ))
          {unefois->CalResPrem_t += 1;
           Tableau<Enum_variable_metrique> tab(10);
           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 ;
           tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
           CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
          };  
       // appel du programme général
       int tabn_taille = tab_noeud.Taille();
       ElemMeca:: SigmaAuNoeud_ResRaid(tabn_taille,
                                        CoSfe->eleCentre->TaPhi(),
                                        CoSfe->eleCentre->TaWi(),
                                        CoSfe-> resErr,CoSfe->raidErr,
                                        CoSfe->eleEr->TaPhi(),
                                        CoSfe->eleEr->TaWi());
     
       // --- on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
       // dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
       // épaisseur_t ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
       // pour éviter des tests et une recopie du code on laisse epaisseur_t, avec une valeur 1
       // ce qui fait qu'elle n'intervient pas
       double epaisseur = 0.; // init
       if (donnee_specif.epais == NULL) {epaisseur = 1.;}
       else {epaisseur = donnee_specif.epais->epaisseur_t;};
       // on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
       int taille_resErr=CoSfe->resErr.Taille();
       for (int i=1;i<= taille_resErr;i++)
          (*CoSfe-> resErr(i)) *= epaisseur;
       CoSfe->raidErr *= epaisseur;

       return    (Element::Er_ResRaid( &(CoSfe-> resErr),&(CoSfe->raidErr)));
    } ;

// calcul du résidu  pour le calcul d'erreur ramenée aux noeuds
Element::Er_ResRaid   SfeMembT::ErreurAuNoeud_ResRaid()
    {  SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
       // dimensionnement de la metrique
       if(!( unefois->CalResPrem_t ))
       { unefois->CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         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 ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
        };  
       // appel du programme général
       int tabn_taille = tab_noeud.Taille();
       ElemMeca::Cal_ErrAuxNoeuds(tabn_taille,
                                   CoSfe->eleCentre->TaPhi(),CoSfe->eleCentre->TaWi()
                                   ,CoSfe-> resErr );
     
       // --- on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
       // dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
       // épaisseur_t ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
       // pour éviter des tests et une recopie du code on laisse epaisseur_t, avec une valeur 1
       // ce qui fait qu'elle n'intervient pas
       double epaisseur = 0.; // init
       if (donnee_specif.epais == NULL) {epaisseur = 1.;}
       else {epaisseur = donnee_specif.epais->epaisseur_t;};
       int taille_resErr=CoSfe->resErr.Taille();
       for (int i=1;i<= taille_resErr;i++)
          (*CoSfe-> resErr(i)) *= epaisseur;
       CoSfe->raidErr *= epaisseur;

       return    (Element::Er_ResRaid( &(CoSfe-> resErr),&(CoSfe->raidErr)));
     } ;
     
// actualisation des ddl et des grandeurs actives de t+dt vers t      
void SfeMembT::TdtversT()
 { lesPtMecaInt.TdtversT(); // contrainte
   if (donnee_specif.epais != NULL)
      donnee_specif.epais->epaisseur_t = donnee_specif.epais->epaisseur_tdt;
   int taille = lesPtMecaInt.NbPti();
   for (int ni=1;ni<= taille; ni++)
     { if (tabSaveDon(ni) != NULL)  tabSaveDon(ni)->TdtversT();
       if (tabSaveTP(ni) != NULL)  tabSaveTP(ni)->TdtversT();
       if (tabSaveDefDon(ni) != NULL)  tabSaveDefDon(ni)->TdtversT();
       }
   ElemMeca::TdtversT_(); // appel de la procédure mère
  };
// actualisation des ddl et des grandeurs actives de t vers tdt      
void SfeMembT::TversTdt()
 { lesPtMecaInt.TversTdt(); // contrainte
   if (donnee_specif.epais != NULL)
      donnee_specif.epais->epaisseur_tdt = donnee_specif.epais->epaisseur_t;
   int taille = lesPtMecaInt.NbPti();
   for (int ni=1;ni<= taille; ni++)
    { if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->TversTdt();
      if (tabSaveTP(ni) != NULL) tabSaveTP(ni)->TversTdt();
      if (tabSaveDefDon(ni) != NULL) tabSaveDefDon(ni)->TversTdt();
      }
   ElemMeca::TversTdt_(); // appel de la procédure mère
  };

// calcul de l'erreur sur l'élément.
// en fait on ne calculera l'erreur que sur la contrainte sur toute la surface
// mais sur le premier pt d'integ d'épaisseur
void SfeMembT::ErreurElement(int type,double& errElemRelative
                            ,double& numerateur, double& denominateur) 
    {  SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
       // dimensionnement de la metrique
       if(!( unefois->CalResPrem_t ))
       { unefois->CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         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 ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
        };  
       // appel du programme général
        ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur,
                              tab_noeud.Taille(),CoSfe->eleCentre->TaPhi(),
                              CoSfe->eleCentre->TaWi(),
                              CoSfe->eleEr->TaPhi(),CoSfe->eleEr->TaWi());
     } ;

// mise à jour de la boite d'encombrement de l'élément, suivant les axes I_a globales
// en retour coordonnées du point mini dans retour.Premier() et du point maxi dans .Second()
// la méthode est différente de la méthode générale car il faut prendre en compte l'épaisseur de l'élément
const DeuxCoordonnees& SfeMembT::Boite_encombre_element(Enum_dure temps)
  { int tab_taille= tab_noeud.Taille();
    // on commence par calculer la boite d'encombrement pour l'élément médian
    Element::Boite_encombre_element( temps);
    // ensuite on augment sytématiquement dans toutes directions d'une valeur h/2 majorée
    // h étant l'épaisseur maxi constatée
    double epais_maxi=0.; // init
    if (donnee_specif.epais != NULL)
      { // cas où les épaisseurs sont stockées sur l'élément
        epais_maxi = donnee_specif.epais->epaisseur_tdt;
       }
    else
      { // cas où les épaisseurs sont stockées aux noeuds
        for (int ine=1;ine<=nombre->nbnce;ine++) // balayage des  noeuds
           epais_maxi = MaX(epais_maxi,tab_noeud(ine)->Valeur_tdt(EPAIS)); 
      };
    
    double hSur2maj = epais_maxi * 0.5
                         * (ParaGlob::param->ParaAlgoControleActifs().Extra_boite_prelocalisation()-1.);
    // ajout d'un extra dans toutes les directions
    hSur2maj += ParaGlob::param->ParaAlgoControleActifs().Ajout_extra_boite_prelocalisation();
    // mise à jour
    boite_encombre.Premier().Ajout_meme_valeur(-hSur2maj);	// le min
	   boite_encombre.Second().Ajout_meme_valeur(hSur2maj);    // le max
    // retour
    return boite_encombre;
   };

// ramène le nombre de points d'intégration de surface correspondant à un type énuméré
// ramène 0 si l'élément n'est pas une plaque ou coque
int SfeMembT::NbPtIntegSurface(Enum_ddl ddl ) const 
 { switch (ddl)
 	{ case ERREUR : return nombre->nbisEr; break;
 	  default : return nombre->nbis; break;
 	};
 };
// ramène le nombre de points d'intégration en épaisseur correspondant à un type énuméré
// ramène 0 si l'élément n'est pas une plaque ou coque
int SfeMembT::NbPtIntegEpaiss(Enum_ddl ddl) const 
 { switch (ddl)
 	{ case ERREUR : return nombre->nbieEr; break;
 	  default : return nombre->nbie; break;
 	};
 };
// ramene l'element geometrique de surface correspondant au ddl passé en paramètre
// ou null si ce n'est pas définie, dans ce cas si l'élément géométrique de surface est 2D
// cela signifie qu'il faut se référer à l'élément générique: ElementGeometrie(...
ElemGeomC0* SfeMembT::ElementGeometrieSurface(Enum_ddl ddl ) const 
 { switch (ddl)
 	{ case ERREUR : return (unefois->doCoMemb->eleEr); break;
 	  default : return (unefois->doCoMemb->eleCentre); break;
 	};
 };
// ramene l'element geometrique d'épaisseur correspondant au ddl passé en paramètre
// ou null si ce n'est pas définie, dans ce cas si l'élément géométrique de surface est 2D
// cela signifie que tout est constant (identique) dans l'épaisseur
ElemGeomC0* SfeMembT::ElementGeometrieEpaiss(Enum_ddl ddl ) const 
 { switch (ddl)
 	{ case ERREUR : return &(unefois->doCoMemb->segmentEr); break;
 	  default : return &(unefois->doCoMemb->segment); break;
 	};
 };

// retourne un numero d'ordre d'un point le plus près ou est exprimé la grandeur enum
// par exemple un point d'intégration, mais n'est utilisable qu'avec des méthodes particulières
// par exemple CoordPtInteg, ou Valeur_a_diff_temps
// car le numéro d'ordre peut-être différent du numéro d'intégration au sens classique 
// temps: dit si c'est à 0 ou t ou tdt
int SfeMembT::PointLePlusPres(Enum_dure temps,Enum_ddl enu, const Coordonnee& M) 
  { int iret;
    // récupération des éléments géométriques correspondant à Enu
	ElemGeomC0* elSur = this->ElementGeometrieSurface(enu);
	ElemGeomC0* elEpa = this->ElementGeometrieEpaiss(enu);	
	// récupération des tableaux des fonctions d'interpolations
	const Tableau <Vecteur> & tabphiS =   elSur->TaPhi();
	const Tableau < Mat_pleine > & tabDphiS =   elSur->TaDphi();
	const Tableau <Vecteur> & tabphiH =   elEpa->TaPhi();
	
	// on boucle sur les pt d'integ pour avoir le point le plus près
	int nbptSur = tabphiS.Taille();
	int nbptEpa = tabphiH.Taille();
	
	Coordonnee P; iret=1; double dist= ConstMath::tresgrand;
	
	int ipt=1;
	for (int iptH = 1;iptH <= nbptEpa;iptH++)
	 for (int iptS = 1; iptS <= nbptSur;iptS++,ipt++)
	 { switch (temps)
	   { case TEMPS_0 : P = ((Met_Sfe1*) met)->PointSfe1M_0
	              (tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
	     case TEMPS_t : P = ((Met_Sfe1*) met)->PointSfe1M_t
	              (tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
	     case TEMPS_tdt : P = ((Met_Sfe1*) met)->PointSfe1M_tdt
	              (tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
	    }
	   double di=(M-P).Norme();
	   if (di < dist) { dist = di; iret = ipt;}; 
     }
    return iret; // retour
   }; 
		         
// 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
Coordonnee SfeMembT::CoordPtInteg(Enum_dure temps,Enum_ddl enu,int iteg,bool& erreur)
  { Coordonnee ptret;  // le retour
    // récupération des éléments géométriques correspondant à Enu
	ElemGeomC0* elSur = this->ElementGeometrieSurface(enu);
	ElemGeomC0* elEpa = this->ElementGeometrieEpaiss(enu);	
	// récupération des tableaux des fonctions d'interpolations
	const Tableau <Vecteur> & tabphiS =   elSur->TaPhi();
	const Tableau < Mat_pleine > & tabDphiS =   elSur->TaDphi();
	const Tableau <Vecteur> & tabphiH =   elEpa->TaPhi();
	// def du numéro d'integ en surface et épaisseur
	int nbptSur = tabphiS.Taille();
	int nbptEpa = tabphiH.Taille();
	int iptS = iteg / nbptEpa +1; // division entière
	int iptH = iteg - (iptS-1) * nbptEpa ; 
    // calcul du point
	if ((iteg < 0) || (iteg > (nbptSur*nbptEpa)))
	 { erreur = false;}
	else
	 { switch (temps)
	   { case TEMPS_0 : ptret = ((Met_Sfe1*) met)->PointSfe1M_0
	              (tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
	     case TEMPS_t : ptret = ((Met_Sfe1*) met)->PointSfe1M_t
	              (tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
	     case TEMPS_tdt : ptret = ((Met_Sfe1*) met)->PointSfe1M_tdt
	              (tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
	    }
     }
    return ptret; // retour
   };