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

// calcul de la nouvelle épaisseur (sans raideur) avec métrique en explicite
// cas où l'épaisseur est stocké dans l'élément
// mise à jour du volume au pti
void SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(const double& epaisseur0, const ParaAlgoControle &
                                   ,const double & epaisseur_t, PtIntegMecaInterne & ptIntegMeca
                                   ,double & epaisseur_tdt, const Met_abstraite::Expli_t_tdt& ex)
  {  // on commence par calculer la trace du tenseur des contraintes
     double traceSig = (ptIntegMeca.SigHH_const() * (* ex.gijBB_tdt)).Trace();
     double traceSig_t = (ptIntegMeca.SigHH_t_const() * (* ex.gijBB_t)).Trace();
     const double& troisK = 3. * ptIntegMeca.ModuleCompressibilite_const(); // récup de la compressibilité
     // calcul de la nouvelle épaisseur
     
     //--- vérification et gestion des cas particuliers ---
     bool cas_particulier=false;
     {// on suppose un module de compressibilité positif sinon ce n'est pas normal
      if  (troisK < 0.)
       {if (ParaGlob::NiveauImpression()>2)
          cout << "\n erreur*** on a trouve une compressibilite negative = " << troisK/3.
              << "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(...";
        // on ne change pas l'épaisseur
        cas_particulier = true;
       };
     
      // classiquement on suppose que troisK doit-être bien supérieur à traceSig:
      // on considère un facteur 10 arbitraire
      // si ce n'est pas le cas, on risque d'avoir un calcul d'épaisseur erroné
      if (Dabs(traceSig) > 10. * Dabs(troisK/3.))
       {if (ParaGlob::NiveauImpression()>2)
          cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
               << " le coef de compressibilite = "<<troisK/3
               << " est par rapport a delta la trace de sigma = "
               << traceSig << " 10 fois inferieur !! ce n'est pas normal a priori  "
               << " on continue sans changer l'epaisseur  ... "
               << "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(..."
               << endl;
        // on ne change pas l'épaisseur
        cas_particulier = true;
      };
      
      double delta_traceSig = traceSig-traceSig_t;
      // maintenant on regarde si le module peut conduire à une épaisseur négative
      if (delta_traceSig > (troisK/3. - ConstMath::petit))
       //  dans le cas où delta traceSig > troisK, cela va donner une épaisseur négative !!
        {if (ParaGlob::NiveauImpression()>2)
          cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
               << " coef de compressibilite = "<<troisK/3
               << " est inferieur a delta trace de sigma = "
               << delta_traceSig
               << " ce qui va donner une epaisseur finale negative !! "
               << " on continue sans changer l'epaisseur  ... "
               << "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(..."
               << endl;
         // on ne change pas l'épaisseur
         cas_particulier = true;
        };
        
      // maintenant on regarde si le module n'est pas trop petit
      // si oui il y a aussi un pb
      if ((Dabs(troisK) <= ConstMath::unpeupetit)
          && (Dabs(delta_traceSig) > 10. * Dabs(troisK/3.))
         )
         // on ne va pas pouvoir peut-être calculer une nouvelle épaisseur
         {if (ParaGlob::NiveauImpression()>2)
             cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
                  << " coef de compressibilite = "<<troisK/3
                  << " et delta trace de sigma = "
                  << delta_traceSig
                  << " sont de valeurs trop faibles (et un peu etrange) !! "
                  << " on continue sans changer l'epaisseur  ... "
                  << "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(..."
                  << endl;
            // on ne change pas l'épaisseur
            cas_particulier = true;
        };
     };

     if (cas_particulier)
       { epaisseur_tdt = epaisseur_t;
         ptIntegMeca.Volume_pti() *= epaisseur_t;
       }
     else // sinon cas sans pb
      {
//     epaisseur_tdt = (epaisseur0 * (*(ex.jacobien_0))) / (*(ex.jacobien_tdt))
//		              * troisK / (troisK - traceSig);
        epaisseur_tdt = (epaisseur_t * (*(ex.jacobien_t))) / (*(ex.jacobien_tdt))
		                           * troisK / (troisK - traceSig+traceSig_t);
        // calcul du volume au pti
        ptIntegMeca.Volume_pti() *= epaisseur_tdt ;
      };
	};
// calcul de la nouvelle épaisseur (sans raideur) avec métrique en implicite
// cas où l'épaisseur est stocké dans l'élément
// mise à jour du volume au pti
void SfeMembT::CalEpaisseurAtdtImp_et_vol_pti(const double& epaisseur0, const ParaAlgoControle &
                                   ,const double & epaisseur_t, PtIntegMecaInterne & ptIntegMeca
                                   ,double & epaisseur_tdt, const Met_abstraite::Impli& ex)
  {  // on commence par calculer la trace du tenseur des contraintes
//     double traceSig = (ptIntegMeca.SigHH_const() * (* ex.gijBB_tdt)).Trace();
//     double traceSig_t = (ptIntegMeca.SigHH_t_const() * (* ex.gijBB_t)).Trace();
     double traceSig = (ptIntegMeca.SigHH_const() && (* ex.gijBB_tdt));
     double traceSig_t = (ptIntegMeca.SigHH_t_const() && (* ex.gijBB_t));
     const double& troisK = 3. * ptIntegMeca.ModuleCompressibilite_const(); // récup de la compressibilité
     // calcul de la nouvelle épaisseur
     
     //--- vérification et gestion des cas particuliers ---
     bool cas_particulier=false;
     {// on suppose un module de compressibilité positif sinon ce n'est pas normal
      if  (troisK < 0.)
       {if (ParaGlob::NiveauImpression()>2)
          cout << "\n erreur*** on a trouve une compressibilite negative = " << troisK
              << "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) ";
        // on ne change pas l'épaisseur
        cas_particulier = true;
       };
     
      // classiquement on suppose que troisK doit-être bien supérieur à traceSig:
      // on considère un facteur 10 arbitraire
      // si ce n'est pas le cas, on risque d'avoir un calcul d'épaisseur erroné
      if (Dabs(traceSig) > 10. * Dabs(troisK/3.))
       {if (ParaGlob::NiveauImpression()>2)
          cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
               << " le coef de compressibilite = "<<troisK/3
               << " est par rapport a delta la trace de sigma = "
               << traceSig << " 10 fois inferieur !! ce n'est pas normal a priori  "
               << " on continue sans changer l'epaisseur  ... "
               << "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) "
               << endl;
        // on ne change pas l'épaisseur
        cas_particulier = true;
      };
      
      double delta_traceSig = traceSig-traceSig_t;
      // maintenant on regarde si le module peut conduire à une épaisseur négative
      if (delta_traceSig > (troisK/3. - ConstMath::petit))
       //  dans le cas où delta traceSig > troisK, cela va donner une épaisseur négative !!
        {if (ParaGlob::NiveauImpression()>2)
          cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
               << " coef de compressibilite = "<<troisK/3
               << " est inferieur a delta trace de sigma = "
               << delta_traceSig
               << " ce qui va donner une epaisseur finale negative !! "
               << " on continue sans changer l'epaisseur  ... "
               << "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) "
               << endl;
         // on ne change pas l'épaisseur
         cas_particulier = true;
        };
        
      // maintenant on regarde si le module n'est pas trop petit
      // si oui il y a aussi un pb
      if ((Dabs(troisK) <= ConstMath::unpeupetit)
          && (Dabs(delta_traceSig) > 10. * Dabs(troisK/3.))
         )
         // on ne va pas pouvoir peut-être calculer une nouvelle épaisseur
         {if (ParaGlob::NiveauImpression()>2)
             cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
                  << " coef de compressibilite = "<<troisK/3
                  << " et delta trace de sigma = "
                  << delta_traceSig
                  << " sont de valeurs trop faibles (et un peu etrange) !! "
                  << " on continue sans changer l'epaisseur  ... "
                  << "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) "
                  << endl;
            // on ne change pas l'épaisseur
            cas_particulier = true;
        };
     };

     if (cas_particulier)
       { epaisseur_tdt = epaisseur_t;
         ptIntegMeca.Volume_pti() *= epaisseur_t;
       }
     else // sinon cas sans pb
      {
        epaisseur_tdt = (epaisseur_t * (*(ex.jacobien_t))) / (*(ex.jacobien_tdt))
                             * troisK / (troisK - traceSig+traceSig_t);
        // calcul du volume au pti
        ptIntegMeca.Volume_pti() *= epaisseur_tdt ;
      };
	};
// calcul de la nouvelle épaisseur moyenne finale (sans raideur) 
// mise à jour des volumes aux pti
// ramène l'épaisseur moyenne calculée à atdt
const double& SfeMembT::CalEpaisseurMoyenne_et_vol_pti(bool atdt)
  {  //     ..  en bouclant sur les pt d'integ enregistré ..
	  // -- on récupère et on calcule les jacobiens moyens à 0 et final	  
	  double jacobien_moy_fin = 0.; // init
   double jacobien_moy_ini = 0.;
	  // -- de même on récupère et on calcul la trace moyenne de la contrainte
	  double traceSig_moy = 0.;double traceSig_moy_ini = 0.;
	  // -- de même on récupère et on calcul le module de compressibilité moyen
	  double troisK_moy = 0.;
   if (donnee_specif.epais != NULL)
	   {Epai&  epais = *(donnee_specif.epais); // pour simplifier
     int indice = 1;
     SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb;         // pour simplifier l'écriture
     int nbtotsurf =CoSfe->eleCentre->Nbi();        
     int nbtotepais =(CoSfe->segment).Nbi();
     int nbi = nbtotsurf * nbtotepais;
     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
//	    for (int i=1;i<= nombre->nbi;i++)
        { // cas de la compressibilité
          const double& troisK = 3. * (*lesPtIntegMecaInterne)(indice).ModuleCompressibilite_const();
          troisK_moy += troisK;
          // cas des jacobiens
          const double& jacobien_0 = *(tabSaveDefDon(indice)->Meti_00().jacobien_);
          if (atdt)
            { const double& jacobien_ini = *(tabSaveDefDon(indice)->Meti_t().jacobien_);
              jacobien_moy_ini += jacobien_ini;
              const double& jacobien_fin = *(tabSaveDefDon(indice)->Meti_tdt().jacobien_);
              jacobien_moy_fin += jacobien_fin;
              // cas de la trace de sigma
              const double traceSig = (*lesPtIntegMecaInterne)(indice).SigHH_const() && (*(tabSaveDefDon(indice)->Meti_tdt().gijBB_));
              traceSig_moy += traceSig;
              const double traceSig_ini = (*lesPtIntegMecaInterne)(indice).SigHH_t_const() && (*(tabSaveDefDon(indice)->Meti_t().gijBB_));
              traceSig_moy_ini += traceSig_ini;
             
              (*lesPtIntegMecaInterne)(indice).Volume_pti() *= (epais.epaisseur_t * jacobien_ini) / jacobien_fin
                                * troisK / (troisK - traceSig+traceSig_ini); // la pression = - traceSig/3 !!!
            }
          else
            { const double& jacobien_ini = *(tabSaveDefDon(indice)->Meti_00().jacobien_);
              jacobien_moy_ini += jacobien_ini;
              const double& jacobien_fin = *(tabSaveDefDon(indice)->Meti_t().jacobien_);
              jacobien_moy_fin += jacobien_fin;
              // cas de la trace de sigma
              double traceSig = (*lesPtIntegMecaInterne)(indice).SigHH_const() && (*(tabSaveDefDon(indice)->Meti_t().gijBB_));
              traceSig_moy += traceSig;
              (*lesPtIntegMecaInterne)(indice).Volume_pti() *= (epais.epaisseur0 * jacobien_ini) / jacobien_fin
                                  * troisK / (troisK - traceSig);
            };
        };
     jacobien_moy_ini /= nbi;
     jacobien_moy_fin /= nbi;
     traceSig_moy_ini /= nbi;
     traceSig_moy /= nbi;
     troisK_moy /= nbi;
     // d'où le calcul de la nouvelle épaisseur en utilisant la relation:
	    // (V-V_0)/V = trace(sigma)/3 /K_moy
   
     //--- vérification et gestion des cas particuliers ---
     bool cas_particulier=false;
     {// on suppose un module de compressibilité positif sinon ce n'est pas normal
      if  (troisK_moy < 0.)
       {if (ParaGlob::NiveauImpression()>2)
          cout << "\n erreur*** on a trouve une compressibilite moyenne negative = "
               << troisK_moy/3.
               << "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(...";
        // on ne change pas l'épaisseur
        cas_particulier = true;
       };
     
      // classiquement on suppose que troisK_moy doit-être bien supérieur à traceSig_moy:
      // on considère un facteur 10 arbitraire
      // si ce n'est pas le cas, on risque d'avoir un calcul d'épaisseur erroné
      double delta_traceSig = traceSig_moy-traceSig_moy_ini;
      if (Dabs(delta_traceSig) > 10. * Dabs(troisK_moy))
       {if (ParaGlob::NiveauImpression()>2)
          cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
               << " le coef de compressibilite moyen = "<<troisK_moy/3.
               << " est par rapport a delta la trace de sigma = "
               << delta_traceSig << " 10 fois inferieur !! ce n'est pas normal a priori  "
               << " on continue sans changer l'epaisseur  ... "
               << "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(... "
               << endl;
        // on ne change pas l'épaisseur
        cas_particulier = true;
      };
      
      // maintenant on regarde si le module peut conduire à une épaisseur négative
      if (delta_traceSig > (troisK_moy/3. - ConstMath::petit))
       //  dans le cas où delta_traceSig > troisK_moy, cela va donner une épaisseur négative !!
        {if (ParaGlob::NiveauImpression()>2)
          cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
               << " coef de compressibilite moyen = "<<troisK_moy/3.
               << " est inferieur a delta la trace de sigma = "
               << delta_traceSig
               << " ce qui va donner une epaisseur finale negative !! "
               << " on continue sans changer l'epaisseur  ... "
               << "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(... "
               << endl;
         // on ne change pas l'épaisseur
         cas_particulier = true;
        };
        
      // maintenant on regarde si le module n'est pas trop petit
      // si oui il y a aussi un pb
      if ((Dabs(troisK_moy) <= ConstMath::unpeupetit)
          && (Dabs(traceSig_moy) > 10. * Dabs(troisK_moy/3.))
         )
         // on ne va pas pouvoir peut-être calculer une nouvelle épaisseur
         {if (ParaGlob::NiveauImpression()>2)
             cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
                  << " coef de compressibilite moyen = "<<troisK_moy/3
                  << " et delta trace de sigma = "
                  << delta_traceSig
                  << " sont de valeurs trop faibles (et un peu etrange) !! "
                  << " on continue sans changer l'epaisseur  ... "
                  << "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(... "
                  << endl;
            // on ne change pas l'épaisseur
            cas_particulier = true;
        };
     };
     
     // si cas particulier: on ne peut pas calculer, on reste sur les valeurs précédentes
     if (cas_particulier)
      {if (atdt)
         { // on met à jour l'épaisseur
           epais.epaisseur_tdt = epais.epaisseur_t;
           // puis retour
           return epais.epaisseur_tdt;
         }
       else
         {// retour
          return epais.epaisseur_t;
         };
      }
     else // sinon c'est ok
       {if (atdt)
        {epais.epaisseur_tdt = (epais.epaisseur_t * jacobien_moy_ini) / jacobien_moy_fin
                               * troisK_moy / (troisK_moy - traceSig_moy+traceSig_moy_ini);
      
      //epais.epaisseur_tdt = (epais.epaisseur0 * jacobien_moy_0) / jacobien_moy_fin
		    //                       * troisK_moy / (troisK_moy + traceSig_moy);
//--debug
//if (num_elt==1)
// { Noeud* noe = tab_noeud(1);
//   double R_0 = noe->Coord0().Norme(); double R = noe->Coord2().Norme();
//cout << "\n e0= " << epais.epaisseur_tdt<< " troisK_moy=" << troisK_moy << " traceSig_moy=" << traceSig_moy
//     << " J0= " << jacobien_moy_0 << " J= " << jacobien_moy_fin << " R_0 " << R_0 << " R= " << R;
//	};
//-- fin debug									
		       return epais.epaisseur_tdt;
		      }
	       else
        {epais.epaisseur_t = (epais.epaisseur0 * jacobien_moy_ini) / jacobien_moy_fin
                             * troisK_moy / (troisK_moy - traceSig_moy);
         return epais.epaisseur_t;
        };
       };
   };
   // sinon c'est une erreur
   cout << "\n *** erreur dans SfeMembT::CalEpaisseurMoyenne_et_vol_pti( ";
   Sortie(1);
	};

//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
bool SfeMembT::Delta_Jacobien_anormal(const Deformation::SaveDefResul* don
                                      ,int nisur,int niepais)
  {// on récupère le jacobien de la surface médiane sans courbure
   const DeformationSfe1::SaveDefResulSfe1* donsfe = ((DeformationSfe1::SaveDefResulSfe1*) don);
   #ifdef MISE_AU_POINT
   if (donsfe->Meti_a_tdt().jacobien_ == NULL)
    {cout << "\n*** erreur, le jacobien de la facette n'est pas attribue ! "
          << " pti de surface "<<nisur
          << "\n element SFE: num= " << Num_elt_const()
          << ", num maillage= " << Num_maillage()
          << "\n SfeMembT::Delta_Jacobien_anormal(..."<<flush;
     Sortie(1);
    };
   if (donsfe->Meti_tdt().giB_ == NULL)
    {cout << "\n*** erreur, la bas giB_tdt n'est pas attribue ! "
          << " pti de surface "<<nisur << " et au pti d'epaisseur " << niepais
          << "\n element SFE: num= " << Num_elt_const()
          << ", num maillage= " << Num_maillage()
          << "\n SfeMembT::Delta_Jacobien_anormal(..."<<flush;
     Sortie(1);
    };
   #endif

   double jaco_init = *(donsfe->Meti_a_tdt().jacobien_);
   // on récupère la base au pti courant
   const BaseB * giB_ = donsfe->Meti_tdt().giB_;
   
   // on calcule le produit mixte pour les gi finaux
   double jacobien_final = Util::ProduitMixte(*giB_);
   double ratio = jacobien_final/(jaco_init+ConstMath::pasmalpetit);
   double limite = ParaGlob::param->ParaAlgoControleActifs().Ratio_maxi_jacoMembrane_jacoPti();
   if ((ratio > limite) || (ratio < 1./(limite+ConstMath::pasmalpetit)))
    cout << "\n *** attention *** le rapport entre le jacobien initial ("<< jaco_init
         << ") et le jacobien ("<< jacobien_final
         << ") calcule au pti de surface "<<nisur << " et au pti d'epaisseur " << niepais
         << " est en dehors de l'intervalle ["<<1./(limite+ConstMath::pasmalpetit) <<","<<limite<<"] "
         << "( peut-etre que les courbure-epaisseur sont trop importantes localement ? )"
         << "\n element SFE: num= " << Num_elt_const()
         << ", num maillage= " << Num_maillage()
         << flush;
  };

// retourne la liste abond้e de tous les  donn้es particuli่res interne actuellement utilis้s
// par l'้l้ment (actif ou non), sont exclu de cette liste les donn้es particuli่res des noeuds
// reli้s เ l'้l้ment
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
List_io <TypeQuelconque> SfeMembT::Les_types_particuliers_internes(bool absolue) const
  { // on commence par récupérer la liste général provenant d'ElemMeca
    List_io <TypeQuelconque> ret = ElemMeca::Les_types_particuliers_internes(absolue);
    // ensuite on va ajouter les données particulières aux sfe
    Grandeur_scalaire_double grand_courant; // def d'une grandeur courante 
    if (donnee_specif.epais != NULL)
      { //  $$$ cas de l'้paisseur initiale 
        TypeQuelconque typQ1(EPAISSEUR_MOY_INITIALE,SIG11,grand_courant);
        ret.push_back(typQ1);
        //  $$$ cas de l'้paisseur finale 
        TypeQuelconque typQ2(EPAISSEUR_MOY_FINALE,SIG11,grand_courant);
        ret.push_back(typQ2);
      };
    //  $$$ cas des efforts normaux classiques, linéique
    {TypeQuelconque typQ1(NN_11,SIG11,grand_courant); ret.push_back(typQ1);};
    {TypeQuelconque typQ1(NN_22,SIG11,grand_courant); ret.push_back(typQ1);};
    {TypeQuelconque typQ1(NN_12,SIG11,grand_courant); ret.push_back(typQ1);};
    //  $$$ cas des moments linéiques
    {TypeQuelconque typQ1(MM_11,SIG11,grand_courant); ret.push_back(typQ1);};
    {TypeQuelconque typQ1(MM_22,SIG11,grand_courant); ret.push_back(typQ1);};
    {TypeQuelconque typQ1(MM_12,SIG11,grand_courant); ret.push_back(typQ1);};
    
    //  $$$ cas du tenseur de courbure
    // on ne l'exprime jamais (pour l'instant) dans le repère absolu)
    // mais on sort toutes les infos: repère local et courbures principales
    Tenseur2BB tens; // un tenseur typique
    Grandeur_TenseurBB grand2(tens); // def d'une grandeur courante TenseurBB
    TypeQuelconque typQ3(TENSEUR_COURBURE,SIG11,grand2);
    ret.push_back(typQ3);
    //  $$$ cas du vecteur contenant les courbures principales
    Tableau <double> t_d(2); // un tableau qui sert d'indicateur de type
    Tab_Grandeur_scalaire_double grand3(t_d); // def d'une grandeur courante tab de scalaire
    TypeQuelconque typQ4(COURBURES_PRINCIPALES,SIG11,grand3);
    ret.push_back(typQ4);
    //  $$$ cas des directions principales de courbure (dans le rep่re local
    Coordonnee3 v_c_princ; // un vecteur qui sert d'indicateur de type 
    Tab_Grandeur_Coordonnee grand4(v_c_princ,2); // def d'une grandeur courante coordonn้e
    TypeQuelconque typQ5(DIRECTIONS_PRINC_COURBURE,SIG11,grand4);
    ret.push_back(typQ5);
    //  $$$ cas du repère local orthonormé d'expression
    Coordonnee3 v_rep;
    Tab_Grandeur_Coordonnee grand5(v_rep,3); // def d'une grandeur courante coordonn้e
    TypeQuelconque typQ6(REPERE_LOCAL_ORTHO,SIG11,grand5);
    ret.push_back(typQ6);
    
    // liberation des tenseurs intermediaires
    LibereTenseur(); 
    return ret;
    }; 
        
// récupération de grandeurs particulières au numéro d'ordre = iteg
// celles-ci peuvent être quelconques
// en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void SfeMembT::Grandeur_particuliere (bool absolue,List_io<TypeQuelconque>& liTQ,int iteg)
   { // *** pour l'instant on ne sort rien de particulier dans le repère absolue !! cf. les grandeurs
     // on balaie la liste transmise pour les grandeurs propres
     List_io<TypeQuelconque>::iterator il,ilfin = liTQ.end();
     // on fait un premier passage pour voir s'il y a une demande d'infos relatives à la courbure
     bool besoin_tenseur_courbure=false;TenseurBB * curbBB_tdt = NULL;
     bool besoin_courbures_principales=false;Coordonnee2 val_princ;
     bool besoin_dir_princ_courb=false; Tableau <Coordonnee3 > t_vec_princ;
     bool besoin_rep_local_ortho=false; Tableau <Coordonnee3 > t_ortho_local;
     for (il=liTQ.begin();il!=ilfin;il++)
      {TypeQuelconque& tipParticu = (*il); // pour simplifier
       if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
        switch (tipParticu.EnuTypeQuelconque().EnumTQ())
         { case EPAISSEUR_MOY_INITIALE: (*il).Inactive(); break; // on inactive la grandeur quelconque
           case EPAISSEUR_MOY_FINALE: (*il).Inactive(); break; // on inactive la grandeur quelconque
           // --- cas des éléments liés  à la courbure on inactive et on crée des booléen adoc
           case TENSEUR_COURBURE:  
             { (*il).Inactive();besoin_tenseur_courbure=true;break;}
           case COURBURES_PRINCIPALES:   
             { (*il).Inactive();besoin_courbures_principales=true;break;}
           case DIRECTIONS_PRINC_COURBURE:   
             { (*il).Inactive();besoin_dir_princ_courb=true;break;}
           case REPERE_LOCAL_ORTHO:   
             { (*il).Inactive();besoin_rep_local_ortho=true;break;}
           default: break; // on ne fait rien
         };
      };
     // dans le cas où on demande des infos relatives à la courbure on les récupères    
     if ((besoin_tenseur_courbure)||(besoin_courbures_principales)||(besoin_dir_princ_courb)
         ||(besoin_rep_local_ortho))
      {// ----- def de grandeurs de travail 
       // def de la dimension des tenseurs  
       int dim = lesPtIntegMecaInterne->DimTens();
       PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(iteg);
       // on recupère le tableau pour la lecture des coordonnées des tenseurs
       int nbcompo = ParaGlob::NbCompTens();
       const Tableau2 <int> & OrdreCont = OrdreContrainteR(nbcompo);
       // définition des grandeurs qui sont indépendante de la boucle sur les ddl_enum_etendue
       def->ChangeNumInteg(iteg);  // on change le num้ro de point d'int้gration courant
       // def du cas de sortie
       int cas=0;       
       if ((unefois->dualSortSfe == 0) && (unefois->CalimpPrem != 0))
         { cas=1;unefois->dualSortSfe = 1;
          } 
       else if ((unefois->dualSortSfe == 1) && (unefois->CalimpPrem != 0))       
          { cas = 11;}
       else if ((unefois->dualSortSfe == 0) && (unefois->CalResPrem_tdt != 0))       
         { cas=2;unefois->dualSortSfe = 1;
          }         
       else if ((unefois->dualSortSfe == 1) && (unefois->CalResPrem_tdt != 0))       
          { 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::Grandeur_particuliere(List_io<TypeQuelconque>&...";
           Sortie(1);
          };
       // choix d'init en fonction de la valeur de "cas"          
       if (((cas == 1) || (cas == 2)))
         { // cas d'une premiere initialisation
           Tableau<Enum_variable_metrique> tab(5);
           tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt;
           tab(4) = igiH_0; tab(5) = igijHH_0;
           met->PlusInitVariables(tab) ;
          };
       //él้ments de m้trique et matrices de passage   
       TenseurHH* gijHH;TenseurBB* gijBB;
       Mat_pleine jB0(dim,dim),jBfin(dim,dim);
       if ((cas==1) || (cas==11))
                // calcul d'une  ortho interessance de visualisation des tenseurs
                // cas de tenseur 3D -> Ia, cas 1D on prend un vecteur norme collineaire
                // avec g1, dans le cas 2D
                // la nouvelle base jB est calculee dans def par projection de "Ia" sur Galpha
                // le resultat est une matrice de passage utilisable pour le changement de base
                // jB0 a t=0, B pour les tenseurs BB, jH0 idem pour les tenseurs HH
                // resultat a t+dt
            {const Met_abstraite::InfoImp& ex = def->RemontImp(absolue,jB0,jBfin);
             gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt;
             }
       else if ((cas==2) || (cas==12)) 
                // resultat a t
            {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(absolue,jB0,jBfin);
             gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt;
             };
       // maintenant on récupère le tenseur de courbure
       const Met_Sfe1::Courbure_t_tdt& curb = ((DeformationSfe1*)def)->RecupCourbure0_t_tdt();
       // passage en absolu si nécessaire
       if (besoin_tenseur_courbure)
        {Tenseur_ns2BB& cur_ns_bBB_tdt = *((Tenseur_ns2BB*) curb.curbBB_tdt) ;
         // on le symétrise
         curbBB_tdt = NevezTenseurBB(dim);
         (*curbBB_tdt).Coor(1,1) = cur_ns_bBB_tdt(1,1);(*curbBB_tdt).Coor(2,2) = cur_ns_bBB_tdt(2,2);
         (*curbBB_tdt).Coor(1,2) = 0.5*(cur_ns_bBB_tdt(1,2) + cur_ns_bBB_tdt(2,1));
         curbBB_tdt->ChBase(jBfin); // le tenseur de courbure dans le rep่re local
        };
       // --- bases ortho locales (du plan de la facette ---
       if (besoin_dir_princ_courb || besoin_rep_local_ortho) 
       	{ t_ortho_local.Change_taille(3);
          t_ortho_local(1) = jBfin(1,1) * curb.aiB_tdt->Coordo(1) + jBfin(2,1) * curb.aiB_tdt->Coordo(2);
          t_ortho_local(2) = jBfin(1,2) * curb.aiB_tdt->Coordo(1) + jBfin(2,2) * curb.aiB_tdt->Coordo(2);
          t_ortho_local(3) = Util::ProdVec_coor(t_ortho_local(1),t_ortho_local(2)); 
       	};
       // --- les directions principales et/ou les valeurs principales
       int ret=0;  
       if (besoin_dir_princ_courb)
       	{ Tenseur_ns2BB& cur_ns_bBB_tdt = *((Tenseur_ns2BB*) curb.curbBB_tdt) ;
          // on le sym้trise
          Tenseur2BB curb_symBB_tdt;
          curb_symBB_tdt.Coor(1,1) = cur_ns_bBB_tdt(1,1);curb_symBB_tdt.Coor(2,2) = cur_ns_bBB_tdt(2,2);
          curb_symBB_tdt.Coor(1,2) = 0.5*(cur_ns_bBB_tdt(1,2) + cur_ns_bBB_tdt(2,1));
          // on cr้e le BH correspondant
       	  Tenseur2BH curbBH = curb_symBB_tdt  * (*gijHH);
          Mat_pleine mat(2,2);
          val_princ = curbBH.ValPropre(ret,mat);
//          // calcul des vecteurs principaux dans la base locale
//          Coordonnee3 V_1=mat(1,1) * curb.aiB_tdt->Coordo(1) + mat(2,1) * curb.aiB_tdt->Coordo(2);
//          Coordonnee3 V_2=mat(1,2) * curb.aiB_tdt->Coordo(1) + mat(2,2) * curb.aiB_tdt->Coordo(2);
//          t_vec_princ.Change_taille(2);
//          t_vec_princ(1)(1)= V_1 * t_ortho_local(1); t_vec_princ(1)(2)= V_1 * t_ortho_local(2);
//          t_vec_princ(2)(1)= V_2 * t_ortho_local(1); t_vec_princ(2)(2)= V_2 * t_ortho_local(2);
          // non finalement on les exprimes dans la base globale sinon on ne peut pas visualiser sans manip
          t_vec_princ.Change_taille(2);
          // s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
          if (ret != -1)
           {t_vec_princ(1)=mat(1,1) * curb.aiH_tdt->Coordo(1) + mat(2,1) * curb.aiH_tdt->Coordo(2);
            t_vec_princ(2)=mat(1,2) * curb.aiB_tdt->Coordo(1) + mat(2,2) * curb.aiB_tdt->Coordo(2);
            t_vec_princ(1).Normer();t_vec_princ(2).Normer();
           };
       	}
       else if (besoin_courbures_principales)
        // si on arrive ici c'est que l'on n'a pas besoin des directions principales mais
        // seulement des valeurs propres de courbure	
       	{ Tenseur_ns2BB& cur_ns_bBB_tdt = *((Tenseur_ns2BB*) curb.curbBB_tdt) ;
          // on le sym้trise
          Tenseur2BB curb_symBB_tdt;
          curb_symBB_tdt.Coor(1,1) = cur_ns_bBB_tdt(1,1);curb_symBB_tdt.Coor(2,2) = cur_ns_bBB_tdt(2,2);
          curb_symBB_tdt.Coor(1,2) = 0.5*(cur_ns_bBB_tdt(1,2) + cur_ns_bBB_tdt(2,1));
          // on cr้e le BH correspondant
       	  Tenseur2BH curbBH = curb_symBB_tdt  * (*gijHH);
          val_princ = curbBH.ValPropre(ret);
       	};
       if (ret == -1)
       	{ cout << "\n pb dans le calcul de valeurs ou vecteurs propres du tenseur de courbure !";
          if (ParaGlob::NiveauImpression() > 6 )
            cout << "\n SfeMembT::Grandeur_particuliere (...";
       	};
      };
     
     // ------ maintenant on fait la sortie réelle  -----
     // on commence par appeler la fonction de la classe mère
     //  il n'y aura pas de calcul des grandeurs inactivées
     // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
     ElemMeca::Grandeur_particuliere (absolue,liTQ,iteg);
     
     // puis les grandeurs sp้cifiques 
     for (il=liTQ.begin();il!=ilfin;il++)
      {TypeQuelconque& tipParticu = (*il); // pour simplifier
       if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
        switch (tipParticu.EnuTypeQuelconque().EnumTQ())
         { 
        // 1) -----cas de l'้paisseur moyenne initiale, ici elle ne d้pend pas du point d'int้gration
           case EPAISSEUR_MOY_INITIALE: 
            {if (donnee_specif.epais != NULL)
                 { *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=donnee_specif.epais->epaisseur0;
                   (*il).Active(); 
                 }; break;
            }
        // 2) -----cas de l'้paisseur moyenne finale, ici elle ne d้pend pas du point d'int้gration
           case EPAISSEUR_MOY_FINALE: (*il).Inactive(); break; // on inactive la grandeur quelconque
            { if (donnee_specif.epais != NULL)
                 { *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=donnee_specif.epais->epaisseur_tdt;
                   (*il).Active();
                 }; break;
             }
        // 3) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas du tenseur de courbure ฐฐฐฐฐฐฐฐฐฐ
           case TENSEUR_COURBURE:  
            { ((Grandeur_TenseurBB*) tipParticu.Grandeur_pointee())->RefConteneurTenseur() = *curbBB_tdt;
              (*il).Active();break;
             }
        // 4) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas du vecteur contenant les courbures principales ฐฐฐฐฐฐฐฐฐฐ
           case COURBURES_PRINCIPALES:   
            { (*((Tab_Grandeur_scalaire_double*) tipParticu.Grandeur_pointee()))(1) = val_princ(1);
              (*((Tab_Grandeur_scalaire_double*) tipParticu.Grandeur_pointee()))(2) = val_princ(2);
              (*il).Active();break;
             }
        // 5) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas des directions principales de courbure (dans le rep่re local ฐฐฐฐฐฐฐฐฐฐ
           case DIRECTIONS_PRINC_COURBURE:   
            { (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(1) = t_vec_princ(1);
              (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(2) = t_vec_princ(2);
              (*il).Active();break;
             }
        // 6) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas du rep่re local orthonorm้ d'expression ฐฐฐฐฐฐฐฐฐฐ
           case REPERE_LOCAL_ORTHO:   
            { (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(1) = t_ortho_local(1);
              (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(2) = t_ortho_local(2);
              (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(3) = t_ortho_local(3);
              (*il).Active();break;
             };
           default: break; // on ne fait rien
         };
       };
      // liberation des tenseurs intermediaires
      LibereTenseur(); 
      if (curbBB_tdt != NULL) delete curbBB_tdt;
    };                             	

 
// calcul éventuel de la normale à un noeud
// ce calcul existe pour les éléments 2D, 1D axi, et aussi pour les éléments 1D
// qui possède un repère d'orientation
// en retour coor = la normale si coor.Dimension() est = à la dimension de l'espace
// si le calcul n'existe pas -->  coor.Dimension() = 0
// ramène un entier :
//    == 1 : calcul normal
//    == 0 : problème de calcul -> coor.Dimension() = 0
//    == 2 : indique que le calcul n'est pas licite pour le noeud passé en paramètre
//           c'est le cas par exemple des noeuds exterieurs pour les éléments SFE
//           mais il n'y a pas d'erreur, c'est seulement que l'élément n'est pas ad hoc pour
//           calculer la normale à ce noeud là
// temps: indique  à quel moment on veut le calcul
// pour des éléments particulier (ex: SFE) la méthode est surchargée
int SfeMembT::CalculNormale_noeud(Enum_dure temps,const Noeud& noe,Coordonnee& coor)
  {  // === par rapport à la méthode d'ElemMeca, ici on ne s'occupe des noeuds de l'élément central
     int retour = 1; // on part avec l'idée que tout va être ok
     Enum_type_geom enutygeom =  Type_geom_generique(this->Id_geometrie());
//     if ((enutygeom == LIGNE) || (enutygeom == SURFACE))
     if (enutygeom != SURFACE) // dans une première étape on ne s'occupe que des surfaces
        {coor.Libere(); // pas de calcul possible
         retour = 0;
        }
     else // sinon le calcul est possible
       { // on commence par repérer le noeud dans la numérotation locale
         int nuoe=0;
//         int borne_nb_noeud=nombre->nbnce+1; // ici on ne prend en compte que les noeuds centraux
         int borne_nb_noeud=tab_noeud.Taille()+1;
         for (int i=1;i< borne_nb_noeud;i++)
           {Noeud& noeloc = *tab_noeud(i);
            if (   (noe.Num_noeud() == noeloc.Num_noeud())
                && (noe.Num_Mail() == noeloc.Num_Mail())
               )
              {nuoe = i; break;
              };
            };
          // on ne continue que si on a trouvé le noeud
          if (nuoe != 0)
           { // il y a deux cas: soit le noeud fait partie des noeuds centraux ou non
             if (nuoe <= nombre->nbnce) // cas où le noeud fait partie du centre
               {ElemGeomC0& elemgeom = ElementGeometrique(); // récup de la géométrie
                // récup des coordonnées locales du noeuds
                const Coordonnee& theta_noeud = elemgeom.PtelemRef()(nuoe);
                // récup des phi et dphi au noeud
                const Vecteur & phi = elemgeom.Phi(theta_noeud);
                const Mat_pleine& dphi = elemgeom.Dphi(theta_noeud);
                switch (temps)
                 {case TEMPS_0 :
                    {const BaseB& baseB = met->BaseNat_0(tab_noeud,dphi,phi);
                     coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
                     coor.Normer();
                     break;
                    }
                  case TEMPS_t :
                   {const BaseB& baseB = met->BaseNat_t(tab_noeud,dphi,phi);
                    coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
                    coor.Normer();
                    break;
                   }
                  case TEMPS_tdt :
                   {const BaseB& baseB = met->BaseNat_tdt(tab_noeud,dphi,phi);
                    coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
                    coor.Normer();
                    break;
                   }
                  default :
                   cout << "\nErreur : valeur incorrecte du  temps demande = "
                        << Nom_dure(temps) << " !\n";
                   cout << "\n SfeMembT::CalculNormale_noeud(Enum_dure temps... \n";
                   retour = 0;
                   Sortie(1);
                 };
               }
             else // sinon le noeud fait partie des noeuds externes
                  // on considère que la normale n'est pas à calculer
               {coor.Libere(); // pas de calcul possible
                retour = 2;
               };
           }
          else
           {cout << "\n *** erreur le noeud demande num= "<<noe.Num_noeud()
                 << " du maillage "<< noe.Num_Mail()
                 << " ne fait pas parti de l'element "
                 << num_elt << " du maillage "<< noe.Num_Mail()
                 << " on ne peut pas calculer la normale au noeud !!"
                 << "\n SfeMembT::CalculNormale_noeud(...";
            retour = 0;
            Sortie(1);
           }
       };
     // retour
     return retour;

  };