Herezh_dev/comportement/lois_combinees/LoiContraintesPlanes.cc

// FICHIER : LoiContraintesPlanes.cp
// CLASSE : LoiContraintesPlanes


// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.


//#include "Debug.h"
#include "LesLoisDeComp.h"

# include <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeQuelconqueParticulier.h"
#include "TypeConsTens.h"
#include "NevezTenseurQ.h"
#include "Util.h"
#include "ExceptionsLoiComp.h"
#include "MotCle.h"
#include "MathUtil.h"
#include "CharUtil.h"


#include "LoiContraintesPlanes.h"


//==================== cas de la class de sauvegarde SaveResul ===================

// constructeur  par défaut à ne pas utiliser
LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::SaveResul_LoiContraintesPlanes() :
	le_SaveResul(NULL) ,l_sigoHH(NULL),l_sigoHH_t(NULL),epsInvar(3),epsInvar_t(3),depsInvar(3),depsInvar_t(3)
    ,def_equi(4), def_equi_t(4)
    ,eps_mecaBB(NULL),eps_mecaBB_t(NULL),met_meca()  // -- partie def mécanique éventuelle
    ,eps_P_mecaBB(NULL),eps_P_mecaBB_t(NULL)  // -- suite
    ,l_energ(),l_energ_t(),hsurh0(0.),h_tsurh0(0.),d_hsurh0()
    ,indicateurs_resolution(),indicateurs_resolution_t()
//    ,niveau_sig33(0.),niveau_sig33_t(0.)
  { cout << "\n erreur, le constructeur par defaut ne doit pas etre utilise !"
	        << "\n LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::SaveResul_LoiContraintesPlanes()";
    Sortie(1);
  }; 

// le constructeur courant
LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::SaveResul_LoiContraintesPlanes
                        (SaveResul* l_des_SaveResul):
	le_SaveResul(NULL) ,l_sigoHH(NULL),l_sigoHH_t(NULL),epsInvar(3),epsInvar_t(3),depsInvar(3),depsInvar_t(3)
	,def_equi(4), def_equi_t(4)
 ,eps_mecaBB(NULL),eps_mecaBB_t(NULL),met_meca()  // -- partie def mécanique éventuelle
 ,eps_P_mecaBB(NULL),eps_P_mecaBB_t(NULL)  // -- suite
 ,l_energ(),l_energ_t(),hsurh0(0.),h_tsurh0(0.),d_hsurh0()
 ,indicateurs_resolution(),indicateurs_resolution_t()
// ,niveau_sig33(0.),niveau_sig33_t(0.)
	 { if (l_des_SaveResul != NULL) le_SaveResul = l_des_SaveResul->Nevez_SaveResul();
    // par défaut les contraintes sont en dim 3
    l_sigoHH=NevezTenseurHH(3,0);
    // idem interHH,  par défaut les contraintes sont en dim 3
    l_sigoHH_t =NevezTenseurHH(3,0.);
   
  }; 

// constructeur de copie 
LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::SaveResul_LoiContraintesPlanes
                                         (const LoiContraintesPlanes::SaveResul_LoiContraintesPlanes& sav ):
	le_SaveResul(NULL) ,l_sigoHH(NULL),l_sigoHH_t(NULL),epsInvar(sav.epsInvar),epsInvar_t(sav.epsInvar_t)
    ,depsInvar(sav.depsInvar),depsInvar_t(sav.depsInvar_t),def_equi(sav.def_equi), def_equi_t(sav.def_equi_t)
    ,eps_mecaBB(NULL),eps_mecaBB_t(NULL),met_meca()  // -- partie def mécanique éventuelle
    ,eps_P_mecaBB(NULL),eps_P_mecaBB_t(NULL)  // -- suite
	   ,l_energ(sav.l_energ),l_energ_t(sav.l_energ_t)
    ,hsurh0(sav.hsurh0),h_tsurh0(sav.h_tsurh0),d_hsurh0()
    ,indicateurs_resolution(sav.indicateurs_resolution)
    ,indicateurs_resolution_t(sav.indicateurs_resolution_t)
// ,niveau_sig33(sav.niveau_sig33),niveau_sig33_t(sav.niveau_sig33_t)
	 { if (sav.le_SaveResul != NULL) {le_SaveResul=(sav.le_SaveResul)->Nevez_SaveResul();};
    // automatiquement les tenseur sont non nuls
    l_sigoHH=NevezTenseurHH(*(sav.l_sigoHH));
    // idem interHH, 
    l_sigoHH_t =NevezTenseurHH(*(sav.l_sigoHH_t));
    // -- partie def mécanique éventuelle
    if (sav.eps_mecaBB != NULL)
        {eps_mecaBB = NevezTenseurBB(*sav.eps_mecaBB);
         eps_mecaBB_t = NevezTenseurBB(*sav.eps_mecaBB_t);
         eps_P_mecaBB = NevezTenseurBB(*sav.eps_P_mecaBB);
         eps_P_mecaBB_t = NevezTenseurBB(*sav.eps_P_mecaBB_t);
        };
  }
	        // destructeur
LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::~SaveResul_LoiContraintesPlanes()
  { if (le_SaveResul != NULL) delete le_SaveResul;  
    delete l_sigoHH; 
    delete l_sigoHH_t;
    if (eps_mecaBB != NULL)
      {delete eps_mecaBB;
       delete eps_mecaBB_t;
       delete eps_P_mecaBB;
       delete eps_P_mecaBB_t;
      };
  };
  
// affectation
Loi_comp_abstraite::SaveResul &
     LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::operator = ( const Loi_comp_abstraite::SaveResul & a)
  { LoiContraintesPlanes::SaveResul_LoiContraintesPlanes& sav
                      = *((LoiContraintesPlanes::SaveResul_LoiContraintesPlanes*) &a);
    // on affecte si non dimensionné, sinon on crée à l'identique
    if (sav.le_SaveResul != NULL)
     { if (le_SaveResul == NULL)
         {le_SaveResul = sav.le_SaveResul->Nevez_SaveResul();}
       else
         {(*le_SaveResul) = *(sav.le_SaveResul);};
     };
    // idem pour les tenseurs
    if (sav.l_sigoHH != NULL)
     { if (l_sigoHH == NULL)
         l_sigoHH = NevezTenseurHH(*(sav.l_sigoHH));
       else
         *(l_sigoHH) = *(sav.l_sigoHH);
     };
    if (sav.l_sigoHH_t != NULL)
     { if (l_sigoHH_t == NULL)
         l_sigoHH_t = NevezTenseurHH(*(sav.l_sigoHH_t));
       else
         *(l_sigoHH_t) = *(sav.l_sigoHH_t);
     };
    // puis les grandeurs qui ne posent pas de pb
    epsInvar = sav.epsInvar;
    epsInvar_t = sav.epsInvar_t;
    depsInvar = sav.depsInvar;
    depsInvar_t = sav.depsInvar_t;
    def_equi = sav.def_equi;
    def_equi_t = sav.def_equi_t;
	   l_energ = sav.l_energ;
    l_energ_t = sav.l_energ_t;
    hsurh0 = sav.hsurh0;
    h_tsurh0 = sav.h_tsurh0;
    d_hsurh0 = sav.d_hsurh0;
    indicateurs_resolution = sav.indicateurs_resolution;
    indicateurs_resolution_t = sav.indicateurs_resolution_t;
    // cas conditionnel de eps_mecaBB et eps_mecaBB_t
    if (sav.eps_mecaBB != NULL)
     { if (eps_mecaBB == NULL)
         eps_mecaBB = NevezTenseurBB(*(sav.eps_mecaBB));
       else
         *(eps_mecaBB) = *(sav.eps_mecaBB);
     };
    if (sav.eps_mecaBB_t != NULL)
     { if (eps_mecaBB_t == NULL)
         eps_mecaBB_t = NevezTenseurBB(*(sav.eps_mecaBB_t));
       else
         *(eps_mecaBB_t) = *(sav.eps_mecaBB_t);
     };
    if (sav.eps_P_mecaBB != NULL)
     { if (eps_P_mecaBB == NULL)
         eps_P_mecaBB = NevezTenseurBB(*(sav.eps_P_mecaBB));
       else
         *(eps_P_mecaBB) = *(sav.eps_P_mecaBB);
     };
    if (sav.eps_P_mecaBB_t != NULL)
     { if (eps_P_mecaBB_t == NULL)
         eps_P_mecaBB_t = NevezTenseurBB(*(sav.eps_P_mecaBB_t));
       else
         *(eps_P_mecaBB_t) = *(sav.eps_P_mecaBB_t);
     };

    return *this;
  };

     //============= 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 LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Lecture_base_info 
                                                 (ifstream& ent,const int cas)
  {  // ici toutes les données sont toujours a priori variables
     // ou en tout cas pour les méthodes appelées, elles sont gérées par le paramètre: cas
     string toto; ent >> toto;
     #ifdef MISE_AU_POINT	 	 
      if (toto != "S_C_Plane")
          { cout << "\n erreur en lecture du conteneur pour la loi de contrainte plane generique "
                 << " \n LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Lecture_base_info(..";
            Sortie(1);     
          }
     #endif
     
     if (le_SaveResul != NULL)
	        le_SaveResul->Lecture_base_info(ent,cas);
     // la contrainte sauvegardée est celle stable uniquement
     l_sigoHH_t->Lecture(ent); // lecture du tenseur à t
     (*l_sigoHH) = (*l_sigoHH_t); // init à la même valeur de la contrainte à tdt
     // invariants de déformation
     ent >> toto >> epsInvar_t(1) >> epsInvar_t(2) >> epsInvar_t(3);
     epsInvar = epsInvar_t;
     // invariants de vitesse de déformation
     ent >> toto >> depsInvar_t(1) >> depsInvar_t(2) >> depsInvar_t(3);
     depsInvar = depsInvar_t;
     // les différentes def cumulée
     ent >> toto >> def_equi_t(1) >> def_equi_t(2) >> def_equi_t(3) >> def_equi_t(4);
     def_equi = def_equi_t;
     // lecture de l'énergie à t
     ent >> l_energ_t; 
     l_energ = l_energ_t; // puis init au cas où de la valeur à tdt
     // idem pour l'élongation d'épaisseur à t
     ent >> toto >> h_tsurh0;
     hsurh0 = h_tsurh0; // puis init au cas où de la valeur à tdt
     d_hsurh0.Zero(); // on initialise le vecteur tangent
     // on met à jour éventuellement les indicateurs de résolution
     if (indicateurs_resolution.Taille())
      { indicateurs_resolution_t = indicateurs_resolution;
        // on met les indicateurs à 0
        indicateurs_resolution.Change_taille(2,0.);
      };
//     // le niveau pour sig33
//     ent >> toto >> niveau_sig33_t;
//     niveau_sig33 = niveau_sig33_t;
      // cas de la déformation mécanique éventuelle
      ent >> toto;
      {int dim = ParaGlob::Dimension();
       if (toto == "1")
        {if (eps_mecaBB_t == NULL) eps_mecaBB_t = NevezTenseurBB(dim,0.);
         ent >> toto; // pour passer le mot clé: eps_mecaBB_t:
         #ifdef MISE_AU_POINT
         if ( toto != "eps_mecaBB_t:")
          {  cout << "\n *** erreur en lecture de la deformation mecanique : on attendait "
                  << "le mot cle: eps_mecaBB_t et on a lu  "<< toto
                  << "\n LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Lecture_base_info(...";
             cout << endl;
             Sortie(1);
          };
         #endif
         eps_mecaBB_t->Lecture(ent); // lecture du tenseur à t
         if (eps_mecaBB == NULL) eps_mecaBB = NevezTenseurBB(dim,0.);
         (*eps_mecaBB) = (*eps_mecaBB_t); // init à la même valeur de la contrainte à tdt
         if (eps_P_mecaBB_t == NULL) eps_P_mecaBB_t = NevezTenseurBB(dim,0.);
         ent >> toto; // pour passer le mot clé: eps_P_mecaBB_t:
         #ifdef MISE_AU_POINT
         if ( toto != "eps_P_mecaBB_t:")
          {  cout << "\n *** erreur en lecture de la deformation mecanique : on attendait "
                  << "le mot cle: eps_P_mecaBB_t et on a lu  "<< toto
                  << "\n LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Lecture_base_info(...";
             cout << endl;
             Sortie(1);
          };
         #endif
         eps_P_mecaBB_t->Lecture(ent); // lecture du tenseur à t
         if (eps_P_mecaBB == NULL) eps_P_mecaBB = NevezTenseurBB(dim,0.);
         (*eps_P_mecaBB) = (*eps_P_mecaBB_t); // init à la même valeur de la contrainte à tdt
        }; // sinon on n'a rien à faire a priori
      };
  };

// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Ecriture_base_info
                                                 (ofstream& sort,const int cas)
   {  // ici toutes les données sont toujours a priori variables
      // ou en tout cas pour les méthodes appelées, elles sont gérées par le paramètre: cas
      sort << "\n S_C_Plane  ";
      // données de la loi
      if (le_SaveResul != NULL) 
			       le_SaveResul->Ecriture_base_info(sort,cas);sort << "  ";
      // la contrainte sauvegardée est celle stable uniquement
      l_sigoHH_t->Ecriture(sort); sort << " " ;// écriture du tenseur
      // invariants de déformation
      sort << " invarDef: " << epsInvar_t(1) << " " << epsInvar_t(2) <<" " << epsInvar_t(3) ;
      // invariants de vitesse de déformation
      sort << " invarVitDef: " << depsInvar_t(1) << " " << depsInvar_t(2) <<" " << depsInvar_t(3) ;
      // les différentes def cumulée
      sort << " def_cumulees: " << def_equi_t(1) << " " << def_equi_t(2) <<" " << def_equi_t(3) <<" " << def_equi_t(4);
      sort << " " << l_energ_t; // écriture de l'énergie à t
      sort << " h_tsurh0 " << h_tsurh0; // écriture de l'élongation d'épaisseur à t
//      sort << "\n niveau_sig33 "<<niveau_sig33_t; // le niveau pour sig33
      // cas de la déformation mécanique éventuelle
      if (eps_mecaBB_t != NULL)
         {sort << " 1 eps_mecaBB_t: "; eps_mecaBB_t->Ecriture(sort); sort << " " ;
          sort << "  eps_P_mecaBB_t: "; eps_P_mecaBB_t->Ecriture(sort); sort << " " ;
         }
      else
         {sort << " 0 ";};

	 };
	        
	  // mise à jour des informations transitoires en définitif s'il y a convergence 
	  // par exemple (pour la plasticité par exemple)
void LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::TdtversT() 
  { if (le_SaveResul != NULL)    // mise à jour du conteneur éventuel de la loi interne
      le_SaveResul->TdtversT();
    (*l_sigoHH_t) = (*l_sigoHH);
    epsInvar_t = epsInvar;
    depsInvar_t = depsInvar;
    def_equi_t = def_equi;
    l_energ_t = l_energ;
    h_tsurh0 = hsurh0;
    d_hsurh0.Zero(); // on initialise le vecteur tangent
//    niveau_sig33_t = niveau_sig33;
    // on met à jour éventuellement les indicateurs de résolution
    if (indicateurs_resolution.Taille())
      { indicateurs_resolution_t = indicateurs_resolution;
        // on met les indicateurs à 0
        indicateurs_resolution.Change_taille(2,0.);
      };
    // cas de la déformation mécanique éventuelle
    if (eps_mecaBB_t != NULL)
      {(*eps_mecaBB_t) = (*eps_mecaBB);
       (*eps_P_mecaBB_t) = (*eps_P_mecaBB);
      };
	};

void LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::TversTdt() 
  { if (le_SaveResul != NULL)    // mise à jour du conteneur éventuel de la loi interne
      le_SaveResul->TversTdt();
    (*l_sigoHH) = (*l_sigoHH_t);
    epsInvar = epsInvar_t;
    depsInvar = depsInvar_t;
    def_equi = def_equi_t;
    l_energ = l_energ_t;
    hsurh0 = h_tsurh0;
//    niveau_sig33 = niveau_sig33_t;
    d_hsurh0.Zero(); // on initialise le vecteur tangent
    // cas de la déformation mécanique éventuelle
    if (eps_mecaBB_t != NULL)
      {(*eps_mecaBB) = (*eps_mecaBB_t);
       (*eps_P_mecaBB) = (*eps_P_mecaBB_t);
      };
	};
			  
// affichage à l'écran des infos
void LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Affiche() const 
{ cout << "\n SaveResul_LoiContraintesPlanes: " ; 
  cout << "\n l_sigoHH: ";
  l_sigoHH->Ecriture(cout);
  cout << "\n l_sigoHH_t: ";
  l_sigoHH_t->Ecriture(cout);
  cout << "\n invariants_deformation_t: " << epsInvar_t(1) << " " << epsInvar_t(2) <<" " << epsInvar_t(3) ;
  cout << "\n invariants_vitesse_de_deformation_t: " << depsInvar_t(1) << " " << depsInvar_t(2) <<" " << depsInvar_t(3) ;
  cout << " def_cumulees: " << def_equi_t(1) << " " << def_equi_t(2) <<" " << def_equi_t(3) <<" " << def_equi_t(4);
  cout << "\n l_energ: ";
  cout << l_energ;
  cout << "\n l_energ_t: ";
  cout << l_energ_t;
  cout << "\n h_tsurh0 "<< h_tsurh0;
  // cas de la déformation mécanique éventuelle
  if (eps_mecaBB != NULL)
    {  cout << "\n eps_mecaBB: "; eps_mecaBB->Ecriture(cout);
    };
  // suite
  if (eps_P_mecaBB != NULL)
    {  cout << "\n eps_P_mecaBB: "; eps_P_mecaBB->Ecriture(cout);
    };
  cout << "\n -- partie relative a la loi interne: ";
  le_SaveResul->Affiche();
//  cout << "\n niveau_sig33 "<< niveau_sig33_t;
  cout << "\n .. fin SaveResul_LoiContraintesPlanes:.. \n" ;
};

//changement de base de toutes les grandeurs internes tensorielles stockées
// beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
// gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
// gpH(i) = gamma(i,j) * gH(j)
void LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::ChBase_des_grandeurs(const Mat_pleine& beta,const Mat_pleine& gamma)
{ // on ne s'intéresse qu'aux grandeurs tensorielles
  if (le_SaveResul != NULL) // s'il y a un conteneur affecté
      le_SaveResul->ChBase_des_grandeurs(beta,gamma);
  l_sigoHH->ChBase(gamma);
  l_sigoHH_t->ChBase(gamma);
  // cas de la déformation mécanique éventuelle
  if (eps_mecaBB != NULL)
      {eps_mecaBB->ChBase(beta);
       eps_mecaBB_t->ChBase(beta);
       eps_P_mecaBB->ChBase(beta);
       eps_P_mecaBB_t->ChBase(beta);
      };

  // les vecteurs sont un moyen de stockage mais pas dépendant de la base, donc pas de chgt de coordonnées
};

// procedure permettant de completer éventuellement les données particulières
// de la loi stockées
// au niveau du point d'intégration par exemple: exemple: un repère d'anisotropie
// completer est appelé apres sa creation avec les donnees du bloc transmis
// peut etre appeler plusieurs fois
Loi_comp_abstraite::SaveResul* LoiContraintesPlanes::SaveResul_LoiContraintesPlanes
         ::Complete_SaveResul(const BlocGen & bloc, const Tableau <Coordonnee>& tab_coor
                             ,const Loi_comp_abstraite* loi)
  {// on transmet au conteneur 3D interne
   const LoiContraintesPlanes * loi_CP = (const LoiContraintesPlanes*) loi;
   le_SaveResul->Complete_SaveResul(bloc,tab_coor,loi_CP->lois_interne);
   return this;
  };

// création des conteneurs pour la déformation mécanique
void LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Creation_def_mecanique()
  {// par défaut on crée des tenseurs de la dimension de l'espace de travail
   int dim = ParaGlob::Dimension();
   if (eps_mecaBB == NULL)
     {eps_mecaBB = NevezTenseurBB(dim,0.);
      eps_mecaBB_t = NevezTenseurBB(dim,0.);
      eps_P_mecaBB = NevezTenseurBB(dim,0.);
      eps_P_mecaBB_t = NevezTenseurBB(dim,0.);
     };

  };

	        // ---- récupération d'information: spécifique à certaine classe dérivée
double LoiContraintesPlanes::SaveResul_LoiContraintesPlanes::Deformation_plastique() 
   { cout << "\n pour l'instant cette option n'est pas implante dans le cas d'une loi"
          << "\n contrainte plane quelconque "
          << "\n double Loi_comp_abstraite::SaveResul_LoiContraintesPlanes::Deformation_plastique()";
     Sortie(1);
     return 0.; // pour taire le warning, mais on ne passe jamais là
   };

//==================== fin du cas de la class de sauvegarde SaveResul ============

LoiContraintesPlanes::LoiContraintesPlanes ()  : // Constructeur par defaut
  Loi_comp_abstraite(LOI_CONTRAINTES_PLANES,RIEN_CATEGORIE_LOI_COMP,2)
  ,type_de_contrainte(RIEN_CONTRAINTE_MATHEMATIQUE),fac_penal(30.),prec(0.005)
  ,lois_interne(NULL)
//  ,niveau_mini_sig33(0.)
//  ,niveauF_grandeurGlobale(NULL),niveauF_grandeurConsultable(NULL)
//  ,niveauF_ddlEtendu(NULL),niveauF_temps(NULL)
  ,alg_zero(),maxi_delta_var_eps_sur_iter_pour_Newton(0.1)
  ,sortie_post(0)
  ,fct_tolerance_residu(NULL),fct_tolerance_residu_rel(NULL)

  // ---- réglage par défaut de la résolution avec linéarisation de sig33(hsurh0)=0
  ,val_initiale(1),racine(1),der_at_racine(1,1)
  ,residu(1),derResidu(1,1),mini_hsurh0(0.001),maxi_hsurh0(1000)
  ,giB_normer_3_tdt_3D_sauve(3)
  // -- conteneur des métriques
  ,expli_3D(NULL),impli_3D(NULL),umat_cont_3D(NULL) 
  // -- les variables pointées dans les conteneurs, et leur pointeur associé éventuellement
  ,giB_0_3D(),giH_0_3D(),giB_t_3D(),giH_t_3D(),giB_tdt_3D(),giH_tdt_3D()   
  ,gijBB_0_3D(),gijHH_0_3D(),gijBB_t_3D(),gijHH_t_3D()
  ,gijBB_tdt_3D(),gijHH_tdt_3D()  
  ,gradVmoyBB_t_3D(),gradVmoyBB_tdt_3D(),gradVBB_tdt_3D()
  ,gradVmoyBB_t_3D_P(NULL),gradVmoyBB_tdt_3D_P(NULL),gradVBB_tdt_3D_P(NULL)  
  ,jacobien_tdt_3D(0.),jacobien_0_3D(0.)

       // les vrais tableaux associés aux tableaux de pointeurs
  ,d_giB_tdt_3D(),d_giH_tdt_3D()   
  ,d_gijBB_tdt_3D(),d2_gijBB_tdt_3D(),d_gijHH_tdt_3D()
  ,d_gijBB_tdt_3D_P(NULL),d2_gijBB_tdt_3D_P(NULL),d_gijHH_tdt_3D_P(NULL)  
  ,d_jacobien_tdt_3D(NULL)  
  ,d_gradVmoyBB_t_3D(),d_gradVmoyBB_tdt_3D(),d_gradVBB_t_3D(),d_gradVBB_tdt_3D()
  ,d_gradVmoyBB_t_3D_P(NULL),d_gradVmoyBB_tdt_3D_P(NULL),d_gradVBB_t_3D_P(NULL),d_gradVBB_tdt_3D_P(NULL)
		 
  // puis les grandeurs hors métriques
  ,sig_HH_t_3D(),sig_HH_3D(),Deps_BB_3D(),eps_BB_3D(),delta_eps_BB_3D()
  ,d_eps_BB_3D_P(),d_sig_HH_3D_P()
  ,d_eps_BB_3D(),d_sig_HH_3D(),d_sigma_deps_3D()
  ,d_sigma_deps_2D()

     // un conteneur d'un point d'intégration courant
  ,ptintmeca(3)

    { // conteneurs des métriques
	     // a priori seules les grandeurs principales sont affecté
	     expli_3D = new Met_abstraite::Expli_t_tdt  // constructeur normal
                          (&giB_0_3D,&giH_0_3D,&giB_t_3D,&giH_t_3D,&giB_tdt_3D,&giH_tdt_3D
                          ,&gijBB_0_3D,&gijHH_0_3D,&gijBB_t_3D,&gijHH_t_3D
                          ,&gijBB_tdt_3D,&gijHH_tdt_3D
                          ,gradVmoyBB_t_3D_P,gradVmoyBB_tdt_3D_P,gradVBB_tdt_3D_P // pas affecté par défaut
                          ,&d_gijBB_tdt_3D_P,&jacobien_tdt_3D,&jacobien_t_3D,&jacobien_0_3D);
      impli_3D = new  Met_abstraite::Impli  // constructeur normal 
                       (&giB_0_3D,&giH_0_3D,&giB_t_3D,&giH_t_3D,&giB_tdt_3D,&d_giB_tdt_3D,&giH_tdt_3D,&d_giH_tdt_3D
                        ,&gijBB_0_3D,&gijHH_0_3D,&gijBB_t_3D,&gijHH_t_3D,&gijBB_tdt_3D,&gijHH_tdt_3D
                        ,gradVmoyBB_t_3D_P,gradVmoyBB_tdt_3D_P,gradVBB_tdt_3D_P // pas affecté par défaut
                        ,&d_gijBB_tdt_3D_P
								                ,d2_gijBB_tdt_3D_P // pas affecté par défaut
								                ,&d_gijHH_tdt_3D_P
                        ,&jacobien_tdt_3D,&jacobien_t_3D,&jacobien_0_3D,&d_jacobien_tdt_3D
                        ,d_gradVmoyBB_t_3D_P,d_gradVmoyBB_tdt_3D_P     // pas affecté par défaut
		  	               ,d_gradVBB_t_3D_P,d_gradVBB_tdt_3D_P);	        // pas affecté par défaut
	     umat_cont_3D = new Met_abstraite::Umat_cont  // constructeur normal
                          (&giB_0_3D,&giH_0_3D,&giB_t_3D,&giH_t_3D,&giB_tdt_3D,&giH_tdt_3D
                          ,&gijBB_0_3D,&gijHH_0_3D,&gijBB_t_3D,&gijHH_t_3D
                          ,&gijBB_tdt_3D,&gijHH_tdt_3D
                          ,gradVmoyBB_t_3D_P,gradVmoyBB_tdt_3D_P,gradVBB_tdt_3D_P      // pas affecté par défaut
                          ,&jacobien_tdt_3D,&jacobien_t_3D,&jacobien_0_3D);
		    // on ajoute les invariants au pt integ courant
		    ptintmeca.Change_statut_Invariants_contrainte (true);						  
      // --- initialisation des paramètres de la résolution de newton
      // tolérance absolue et relative sur la résolution de l'équation sig33(hsurh0)=0
      double  tolerance_residu = 5.e-3;
      double  tolerance_residu_rel = 1.e-4;
      int nb_boucle_maxi = 6; // le maximum d'itération permis
      int nb_dichotomie = 4; // le maxi de dichotomie prévu pour l'équation de Newton
      double coef_mini_delta_x = ConstMath::trespetit; // coef pour piloter le mini de delta x dans newton
      alg_zero.Modif_prec_res_abs(tolerance_residu);
      alg_zero.Modif_prec_res_rel(tolerance_residu_rel);
      alg_zero.Modif_iter_max(nb_boucle_maxi);
      alg_zero.Modif_nbMaxiIncre(nb_dichotomie);
      alg_zero.Modif_coef_mini_delta_x(coef_mini_delta_x);
      der_at_racine.Change_Choix_resolution(CRAMER,DIAGONAL);
      derResidu.Change_Choix_resolution(CRAMER,DIAGONAL);
	   };
	
// Constructeur de copie
LoiContraintesPlanes::LoiContraintesPlanes (const LoiContraintesPlanes& loi) :
   Loi_comp_abstraite(loi)
	 ,type_de_contrainte(loi.type_de_contrainte)
	 ,fac_penal(loi.fac_penal),prec(loi.prec)
  ,lois_interne(NULL)
//  ,niveau_mini_sig33(0.)
//  ,niveauF_grandeurGlobale(NULL),niveauF_grandeurConsultable(NULL)
//  ,niveauF_ddlEtendu(NULL),niveauF_temps(NULL)
	 ,alg_zero(loi.alg_zero),maxi_delta_var_eps_sur_iter_pour_Newton(loi.maxi_delta_var_eps_sur_iter_pour_Newton)
  ,sortie_post(loi.sortie_post)
  ,fct_tolerance_residu(loi.fct_tolerance_residu),fct_tolerance_residu_rel(loi.fct_tolerance_residu_rel)
  // ---- réglage résolution avec linéarisation
  ,val_initiale(1),racine(1),der_at_racine(1,1)
  ,giB_normer_3_tdt_3D_sauve(3)
  ,residu(1),derResidu(1,1),mini_hsurh0(loi.mini_hsurh0),maxi_hsurh0(loi.maxi_hsurh0)
  // -- conteneur des métriques: ce sont des pointeurs, pour l'instant on ne les affecte pas
  ,expli_3D(NULL),impli_3D(NULL),umat_cont_3D(NULL) 
  // -- les variables pointées dans les conteneurs
  ,giB_0_3D(loi.giB_0_3D),giH_0_3D(loi.giH_0_3D),giB_t_3D(loi.giB_t_3D),giH_t_3D(loi.giH_t_3D) 
  ,giB_tdt_3D(loi.giB_tdt_3D),d_giB_tdt_3D(loi.d_giB_tdt_3D),giH_tdt_3D(loi.giH_tdt_3D)
  ,d_giH_tdt_3D(loi.d_giH_tdt_3D),gijBB_0_3D(loi.gijBB_0_3D),gijHH_0_3D(loi.gijHH_0_3D) 
  ,gijBB_t_3D(loi.gijBB_t_3D),gijHH_t_3D(loi.gijHH_t_3D),gijBB_tdt_3D(loi.gijBB_tdt_3D)
  ,gijHH_tdt_3D(loi.gijHH_tdt_3D),gradVmoyBB_t_3D(loi.gradVmoyBB_t_3D)
  ,gradVmoyBB_tdt_3D(loi.gradVmoyBB_tdt_3D),gradVBB_tdt_3D(loi.gradVBB_tdt_3D)
  ,jacobien_0_3D(loi.jacobien_0_3D),d_jacobien_tdt_3D(loi.d_jacobien_tdt_3D)
  ,d_gradVBB_t_3D_P(loi.d_gradVBB_t_3D_P),d_gradVBB_tdt_3D_P(loi.d_gradVBB_tdt_3D_P)
       // avec les vrais tableaux associés aux tableaux de pointeurs: ici uniquement le dimensionnement
  ,d_gijBB_tdt_3D(loi.d_gijBB_tdt_3D),d2_gijBB_tdt_3D(loi.d2_gijBB_tdt_3D)
  ,d_gijHH_tdt_3D(loi.d_gijHH_tdt_3D)
  ,d_gradVmoyBB_t_3D(loi.d_gradVmoyBB_t_3D),d_gradVmoyBB_tdt_3D(loi.d_gradVmoyBB_tdt_3D)
  ,d_gradVBB_t_3D(loi.d_gradVBB_t_3D),d_gradVBB_tdt_3D(loi.d_gradVBB_tdt_3D)
		 
  // puis les grandeurs hors métriques, pour les tableaux de pointeurs, c'est uniquement du dimensionnement
  ,sig_HH_t_3D(loi.sig_HH_t_3D),sig_HH_3D(loi.sig_HH_3D),Deps_BB_3D(loi.Deps_BB_3D)
  ,eps_BB_3D(loi.eps_BB_3D),delta_eps_BB_3D(loi.delta_eps_BB_3D)
  ,d_eps_BB_3D_P(loi.d_eps_BB_3D_P),d_sig_HH_3D_P(loi.d_sig_HH_3D_P)
  ,d_eps_BB_3D(loi.d_eps_BB_3D),d_sig_HH_3D(loi.d_sig_HH_3D),d_sigma_deps_3D(loi.d_sigma_deps_3D)
  ,d_sigma_deps_2D(loi.d_sigma_deps_2D)


     // un conteneur d'un point d'intégration courant
  ,ptintmeca(loi.ptintmeca)
	
     { lois_interne = loi.lois_interne->Nouvelle_loi_identique();
	      // association des pointeurs de grandeurs si nécessaire		 
	      if (loi.gradVmoyBB_t_3D_P != NULL) {gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
	      if (loi.gradVmoyBB_tdt_3D_P != NULL) {gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
	      if (loi.gradVBB_tdt_3D_P != NULL) {gradVBB_tdt_3D_P = &gradVBB_tdt_3D;}; 
		 
	      // def des tableaux de pointeurs pour les conteneurs de métriques
		
       // -- cas des tableaux de pointeurs, ils ont déjà la bonne dimension
	      int ta_d_gijBB_tdt_3D = loi.d_gijBB_tdt_3D.Taille();
	      for (int i=1;i<= ta_d_gijBB_tdt_3D;i++)
	      		d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));

	      if (loi.d2_gijBB_tdt_3D_P!= NULL)
	       { int tai_d2_gijBB_tdt_3D = d2_gijBB_tdt_3D.Taille1();
	         d2_gijBB_tdt_3D_P = new Tableau2<TenseurBB *> (tai_d2_gijBB_tdt_3D);
	         int taj_d2_gijBB_tdt_3D = d2_gijBB_tdt_3D.Taille2();
	         for (int i=1;i<= tai_d2_gijBB_tdt_3D;i++)
	           for (int j=1;j<= taj_d2_gijBB_tdt_3D;j++)
	             (*d2_gijBB_tdt_3D_P)(i,j) = &(d2_gijBB_tdt_3D(i,j));
	       };
			 
	      int ta_d_gijHH_tdt_3D = d_gijHH_tdt_3D.Taille();
	      for (int i=1;i<= ta_d_gijHH_tdt_3D;i++)
	      		d_gijHH_tdt_3D_P(i) = &(d_gijHH_tdt_3D(i));
	     
	      if (loi.d_gradVmoyBB_t_3D_P != NULL)
	       { int ta_d_gradVmoyBB_t_3D = d_gradVmoyBB_t_3D.Taille();
	         d_gradVmoyBB_t_3D_P = new Tableau<TenseurBB *> (ta_d_gradVmoyBB_t_3D);
	      	for (int i=1;i<= ta_d_gradVmoyBB_t_3D;i++)
	      		(*d_gradVmoyBB_t_3D_P)(i) = &(d_gradVmoyBB_t_3D(i));
	       };
			 
	      if (loi.d_gradVmoyBB_tdt_3D_P != NULL)
	       { int ta_d_gradVmoyBB_tdt_3D = d_gradVmoyBB_tdt_3D.Taille();
	         d_gradVmoyBB_tdt_3D_P = new Tableau<TenseurBB *> (ta_d_gradVmoyBB_tdt_3D);
	      	  for (int i=1;i<= ta_d_gradVmoyBB_tdt_3D;i++)
	      		   (*d_gradVmoyBB_tdt_3D_P)(i) = &(d_gradVmoyBB_tdt_3D(i));
	       };		 
		    
	      if (loi.d_gradVBB_t_3D_P != NULL)
	       { int ta_d_gradVBB_t_3D = d_gradVBB_t_3D.Taille();
	         d_gradVBB_t_3D_P = new Tableau<TenseurBB *> (ta_d_gradVBB_t_3D);
	      	  for (int i=1;i<= ta_d_gradVBB_t_3D;i++)
	      		  (*d_gradVBB_t_3D_P)(i) = &(d_gradVBB_t_3D(i));
	       };
			 
	      if (loi.d_gradVBB_tdt_3D_P != NULL)
	       { int ta_d_gradVBB_tdt_3D = d_gradVBB_tdt_3D.Taille();
	         d_gradVBB_tdt_3D_P = new Tableau<TenseurBB *> (ta_d_gradVBB_tdt_3D);
	      	  for (int i=1;i<= ta_d_gradVBB_tdt_3D;i++)
	      		   (*d_gradVBB_tdt_3D_P)(i) = &(d_gradVBB_tdt_3D(i));
	       };		 
		
	      // conteneurs des métriques
	      // a priori seules les grandeurs principales sont affecté
	      expli_3D = new Met_abstraite::Expli_t_tdt  // constructeur normal
                          (&giB_0_3D,&giH_0_3D,&giB_t_3D,&giH_t_3D,&giB_tdt_3D,&giH_tdt_3D
                          ,&gijBB_0_3D,&gijHH_0_3D,&gijBB_t_3D,&gijHH_t_3D
                          ,&gijBB_tdt_3D,&gijHH_tdt_3D
                          ,gradVmoyBB_t_3D_P,gradVmoyBB_tdt_3D_P,gradVBB_tdt_3D_P // pas affecté par défaut
                          ,&d_gijBB_tdt_3D_P,&jacobien_tdt_3D,&jacobien_t_3D,&jacobien_0_3D);
       impli_3D = new  Met_abstraite::Impli  // constructeur normal 
                       (&giB_0_3D,&giH_0_3D,&giB_t_3D,&giH_t_3D,&giB_tdt_3D,&d_giB_tdt_3D,&giH_tdt_3D,&d_giH_tdt_3D
                        ,&gijBB_0_3D,&gijHH_0_3D,&gijBB_t_3D,&gijHH_t_3D,&gijBB_tdt_3D,&gijHH_tdt_3D
                        ,gradVmoyBB_t_3D_P,gradVmoyBB_tdt_3D_P,gradVBB_tdt_3D_P // pas affecté par défaut
                        ,&d_gijBB_tdt_3D_P
								                ,d2_gijBB_tdt_3D_P // pas affecté par défaut
								                ,&d_gijHH_tdt_3D_P
                        ,&jacobien_tdt_3D,&jacobien_t_3D,&jacobien_0_3D,&d_jacobien_tdt_3D
                        ,d_gradVmoyBB_t_3D_P,d_gradVmoyBB_tdt_3D_P     // pas affecté par défaut
		  	               ,d_gradVBB_t_3D_P,d_gradVBB_tdt_3D_P);	        // pas affecté par défaut
	      umat_cont_3D = new Met_abstraite::Umat_cont  // constructeur normal
                          (&giB_0_3D,&giH_0_3D,&giB_t_3D,&giH_t_3D,&giB_tdt_3D,&giH_tdt_3D
                          ,&gijBB_0_3D,&gijHH_0_3D,&gijBB_t_3D,&gijHH_t_3D
                          ,&gijBB_tdt_3D,&gijHH_tdt_3D
                          ,gradVmoyBB_t_3D_P,gradVmoyBB_tdt_3D_P,gradVBB_tdt_3D_P      // pas affecté par défaut
                          ,&jacobien_tdt_3D,&jacobien_t_3D,&jacobien_0_3D);
			 
	      // puis les tableaux de pointeurs de grandeurs hors métriques
	      int ta_d_eps_BB_3D = d_eps_BB_3D.Taille();
	      for (int i=1;i<= ta_d_eps_BB_3D;i++)
	          d_eps_BB_3D_P(i) = &(d_eps_BB_3D(i));
			 
	      int ta_d_sig_HH_3D = d_sig_HH_3D.Taille();
	      for (int i=1;i<= ta_d_sig_HH_3D;i++)
	          d_sig_HH_3D_P(i) = &(d_sig_HH_3D(i));
       // initialisation des paramètres de la résolution de newton
       der_at_racine.Change_Choix_resolution(CRAMER,DIAGONAL);
       derResidu.Change_Choix_resolution(CRAMER,DIAGONAL);
//       // traitement des pondérations
//       if (loi.niveauF_grandeurGlobale != NULL)
//         niveauF_grandeurGlobale = new Ponderation_GGlobal(*(loi.niveauF_grandeurGlobale));
//       if (loi.niveauF_ddlEtendu != NULL)
//         niveauF_ddlEtendu = new Ponderation(*(loi.niveauF_ddlEtendu));
//       if (loi.niveauF_temps != NULL)
//            niveauF_temps = new Ponderation_temps(*(loi.niveauF_temps));;
//       if (loi.niveauF_grandeurConsultable != NULL)
//         niveauF_grandeurConsultable = new Ponderation_Consultable(*(loi.niveauF_grandeurConsultable));
       // on regarde s'il y a un pilotage via des fct nD de la précision
       if (fct_tolerance_residu != NULL)
         { if (fct_tolerance_residu->NomFonction() == "_")
          {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
           string non_fonction("_");
           fct_tolerance_residu = Fonction_nD::New_Fonction_nD(*fct_tolerance_residu);
          };
         };
       if (fct_tolerance_residu_rel != NULL)
         { if (fct_tolerance_residu_rel->NomFonction() == "_")
          {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
           string non_fonction("_");
           fct_tolerance_residu_rel = Fonction_nD::New_Fonction_nD(*fct_tolerance_residu_rel);
          };
         };
   		};

LoiContraintesPlanes::~LoiContraintesPlanes ()
// Destructeur 
     { if (lois_interne != NULL) delete lois_interne;
       // les conteneurs de pointeurs:
       delete expli_3D;delete  impli_3D; delete umat_cont_3D;
		 
       // pour les grandeurs de base, pas de new, donc pas de delete
		 
       // pour les tableaux de pointeurs, idem s'ils ne sont pas addressé par un pointeur
       // sinon il faut détruire le tableau

       if (d2_gijBB_tdt_3D_P!= NULL) delete d2_gijBB_tdt_3D_P;
       if (d_gradVmoyBB_t_3D_P!= NULL) delete d_gradVmoyBB_t_3D_P;
       if (d_gradVmoyBB_tdt_3D_P!= NULL) delete d_gradVmoyBB_tdt_3D_P;
       if (d_gradVBB_t_3D_P!= NULL) delete d_gradVBB_t_3D_P;
       if (d_gradVBB_tdt_3D_P!= NULL) delete d_gradVBB_tdt_3D_P;		 
      
//       // cas des pondérations
//       if (niveauF_grandeurGlobale != NULL)
//          delete niveauF_grandeurGlobale;
//       if (niveauF_ddlEtendu != NULL)
//          delete niveauF_ddlEtendu;
//       if (niveauF_temps != NULL)
//          delete niveauF_temps;
//       if (niveauF_grandeurConsultable != NULL)
//          delete niveauF_grandeurConsultable;
       // cas du pilotage de la précision
       if (fct_tolerance_residu != NULL)
         if (fct_tolerance_residu->NomFonction() == "_") delete fct_tolerance_residu;
       if (fct_tolerance_residu_rel != NULL)
         if (fct_tolerance_residu_rel->NomFonction() == "_") delete fct_tolerance_residu_rel;
	    };

   // def d'une instance de données spécifiques, et initialisation
// valable une fois que les différentes lois internes sont définit
LoiContraintesPlanes::SaveResul * LoiContraintesPlanes::New_et_Initialise()  
   { // on crée éventuellement le conteneur pour la loi
     SaveResul* le_SaveResul = NULL; 
     if (lois_interne != NULL) le_SaveResul = lois_interne->New_et_Initialise(); 
     // on ramène la bonne instance
     LoiContraintesPlanes::SaveResul * retour =  new SaveResul_LoiContraintesPlanes(le_SaveResul);
     // retour
     return retour;	  
	  };

// Lecture des donnees de la classe sur fichier
void LoiContraintesPlanes::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
  { // vérification que la dimension de l'espace est 3D
    if (ParaGlob::Dimension() != 3)
	     { cout << "\n *** erreur 1 : la dimension de l'espace doit etre 3, pour pouvoir utiliser une loi de contrainte plane ";
        throw (UtilLecture::ErrNouvelleDonnee(-1));
	       if (ParaGlob::NiveauImpression() > 5)
		        cout << "\n LoiContraintesPlanes::LectureDonneesParticulieres (... ";
	       Sortie(1);
		    };
   
    // on regarde s'il y a une indication de niveau mini pour la contrainte sig33
    string nom_class_methode("LoiContraintesPlanes::LectureDonneesParticulieres");
{
//    // ---  niveau_mini_sig33
////    bool passer_une_ligne = false;
//    if(strstr(entreePrinc->tablcar,"plus_gestion_valeur_sig33_")!=0)
//     { entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
//       string mot_cle="valeur_mini_pour_sig33_=";string st1,st2;
//       bool retour = entreePrinc->Lecture_un_parametre_double
//          (0.,nom_class_methode,-ConstMath::grand,ConstMath::grand,mot_cle,niveau_mini_sig33);
//       entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
//       // on regarde maintenant s'il y a un contrôle du niveau du mini via des pondérations
//       if(strstr(entreePrinc->tablcar,"les_grandeurs_de_controle_=")!=0)
//        { Ponderation ponder; // un élément courant
//          // on initialise
//          string mot_cle, st1,st2;
//          int compteur = 1; // pour éviter une boucle infinie
//          List_io < Ddl_enum_etendu > listddlenum;
//          List_io <bool> listbool;
//         
//          // on lit les grandeurs pour l'instant en string
//          List_io <string> list_id_grand; //liste de travail de toutes les grandeurs
//          List_io <string> list_id_grand_glob; //liste des grandeurs globales
//          List_io <string> list_id_ddl_etendu; //liste des ddl étendus
//         
//          int nb_grandeurs_globales = 0; // init
//          while (st1 != "fin_grandeurs_")
//           {*(entreePrinc->entree)  >> st1 ;
//            if (st1 == "fin_grandeurs_") break; // pas terrible mais c'est ce qui me vient à l'idée !
//            // on vérifie que le mot clé lu est exploitable
//            string minus = Minuscules(st1); // la forme minuscules
//            if ((st1=="TEMPS") || (Ddl_enum_etendu::VerifExistence(st1))
//                || EstUneGrandeurGlobale(minus)
//               )
//              {list_id_grand.push_back(st1);
//               if (EstUneGrandeurGlobale(minus)) // si c'est une grandeur globale on cumule le compteur
//                  nb_grandeurs_globales++;
//              }
//            else
//              { cout << "\n erreur en lecture, le type de grandeur lu" << st1
//                     << " n'est pas acceptable  "
//                     << "\n LoiContraintesPlanes::LectureDonneesParticulieres (...";
//                entreePrinc->MessageBuffer("**erreur05**");
//                throw (UtilLecture::ErrNouvelleDonnee(-1));
//                Sortie(1);     
//              };
//           // on rempli les listes spécifiques
//           if (Ddl_enum_etendu::VerifExistence(st1))
//             list_id_ddl_etendu.push_back(st1);
//           if (EstUneGrandeurGlobale(minus))
//             list_id_grand_glob.push_back(st1);
//           };
//         
//          // --- maintenant on va lire les fonctions
//          entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
//          mot_cle="deb_fonct_="; // le début des fonctions
//          entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
//          List_io <string>::iterator ili,ilifin = list_id_grand.end();
//          List_io <Courbe1D*> li_courbe_globale; // inter
//          List_io <bool> val_aux_noeuds; // inter
//          int nb_GGlob=0;int nb_enum_etendu = 0; // compteur pour le dimensionnement
//          // premier passage on crée les fonctions
//          for (ili=list_id_grand.begin();ili != ilifin;ili++)
//            {string minus = Minuscules(*ili); // la forme minuscules
//             if ((*ili)=="TEMPS") // cas de la fonction du temps
//               { // lecture de la loi d'évolution
//                 val_aux_noeuds.push_back(false); // sert à rien, mais permet
//                      // de garder le même ordre de rangement que list_id_grand
//                 if (niveauF_temps == NULL)
//                    niveauF_temps = new Ponderation_temps;
//                 // lecture de la pondération
//                 niveauF_temps->LectureDonneesPonderation(entreePrinc, lesCourbes1D);
//               }
//              else if (EstUneGrandeurGlobale(minus))
//               { val_aux_noeuds.push_back(false); // sert à rien, mais permet
//                 nb_GGlob++;             // de garder le même ordre de rangement que list_id_grand
//                 // pour le premier élément de la liste on lit la fonction
//                 // pour les autres, on ne fait rien (la fct est déjà lue)
//                 if (nb_GGlob == 1)
//                  {// création de la pondération si nécessaire
//                   if (niveauF_grandeurGlobale == NULL)
//                      niveauF_grandeurGlobale = new Ponderation_GGlobal;
//                   // on lit la fonction
//                   niveauF_grandeurGlobale->LecturePonderation
//                             (list_id_grand_glob,entreePrinc,lesFonctionsnD);
//                  };
//               }
//              else if (Ddl_enum_etendu::VerifExistence(*ili))
//               { nb_enum_etendu++;
//                 // pour le premier élément de la liste on lit la pondération
//                 // pour les autres, on ne fait rien (la pondération est déjà lue)
//                 if (nb_enum_etendu == 1)
//                  { if (niveauF_ddlEtendu==NULL) // on le crée
//                       niveauF_ddlEtendu = new Ponderation;
//                    // lecture
//                    niveauF_ddlEtendu->LectureDonneesPonderation
//                                   (list_id_ddl_etendu,entreePrinc,lesCourbes1D);
//                  };
//               };
//            };
//          // on passe le mot clé de fin
//          mot_cle="fin_fonct_"; // la fin des fonctions
//          entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
//             
//             
//   // *** à changer en fonction de l'utilisation
//              // deuxième passage on met à jour les conteneurs
//   //           niveauF_grandeurGlobale->type_grandeur_GGlob.Change_taille(nb_GGlob);
//   //           niveauF_grandeurGlobale->c_proport.Change_taille(nb_GGlob);
//             
//          entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
////          passer_une_ligne = false;
//        };
////       if (passer_une_ligne)
////         entreePrinc->NouvelleDonnee();
//     };
}
  
    // on lit maintenant la méthode pour prendre en compte la contrainte
	   *(entreePrinc->entree)  >> type_de_contrainte;
   
	   // on met tout d'abord les valeurs par défaut
	   prec = 0.005; fac_penal = 30.;
    string nom;
    // on traite en fonction du type de contrainte
    switch (type_de_contrainte)
	     {
		      case PERTURBATION :
         { if (strstr(entreePrinc->tablcar,"deformation_epaisseur")!=NULL)
            { *(entreePrinc->entree)  >> nom;
              if (nom != "deformation_epaisseur")
                { cout << "\n erreur en lecture de la loi Contraintes planes, on attendait le mot cle deformation_epaisseur"
                        << " alors que l'on a lue " << nom;
                  entreePrinc->MessageBuffer("**erreur 2 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
                  throw (UtilLecture::ErrNouvelleDonnee(-1));
                  Sortie(1);
                };
            }
           else // sinon les autres cas de perturbation ne sont actuellement pas pris en compte 
            { cout << "\n erreur en lecture du type de perturbation : pour l'instant seule "
                << " le type: deformation_epaisseur , est pris en compte " ;
              entreePrinc->MessageBuffer("**erreur 3 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);
            };
           break;
         }
		      case PENALISATION : case MULTIPLICATEUR_DE_LAGRANGE :
         { // ---- on regarde s'il faut lire une précision
           if (strstr(entreePrinc->tablcar,"prec=")!=NULL)
            { *(entreePrinc->entree)  >> nom;
              if (nom != "prec=")
                { cout << "\n erreur en lecture de la loi Contraintes planes, on attendait le mot cle prec="
                       << " alors que l'on a lue " << nom;
                  entreePrinc->MessageBuffer("**erreur 4 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
                  throw (UtilLecture::ErrNouvelleDonnee(-1));
                  Sortie(1);
                };
              *(entreePrinc->entree)  >> prec;
            };

           // --- on regarde s'il faut lire la pénalisation
           if (strstr(entreePrinc->tablcar,"fac=")!=NULL)
            { *(entreePrinc->entree)  >> nom;
              if (nom != "fac=")
                { cout << "\n erreur en lecture de la loi Contraintes planes, on attendait le mot cle fac="
                       << " alors que l'on a lue " << nom;
                  entreePrinc->MessageBuffer("**erreur 5 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
                  throw (UtilLecture::ErrNouvelleDonnee(-1));
                  Sortie(1);
                };
               *(entreePrinc->entree)  >> fac_penal;
            };
           // --- on regarde s'il faut lire le niveau local d'affichage
           if (strstr(entreePrinc->tablcar,"permet_affichage_")!=NULL)
            { *(entreePrinc->entree)  >> nom;
              if (nom != "permet_affichage_")
                { cout << "\n erreur en lecture de la loi Contraintes planes, on attendait le mot cle permet_affichage_ "
                       << " alors que l'on a lue " << nom;
                  entreePrinc->MessageBuffer("**erreur 51 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
                  throw (UtilLecture::ErrNouvelleDonnee(-1));
                  Sortie(1);
                };
               Lecture_permet_affichage(entreePrinc,lesFonctionsnD);
               alg_zero.Modif_affichage(Permet_affichage()); // on met à jour l'algo de newton
            };
           break;
         }
        
		      case NEWTON_LOCAL :
         { // ---- on regarde s'il faut lire des paramètres de réglage
        
            // --- lecture éventuelle des paramètres de réglage ----
            //     de l'algo de résolution de l'équation d'avancement temporel
            if(strstr(entreePrinc->tablcar,"avec_parametres_de_reglage_")!=0)
             {entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
              // on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_reglage_Algo_Newton_"
              // ou un nouveau mot clé global auquel cas il y a pb !!
              MotCle motCle; // ref aux mots cle
              while (strstr(entreePrinc->tablcar,"fin_parametres_reglage_Algo_Newton_")==0)
               {
                // si on a  un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur
                if ( motCle.SimotCle(entreePrinc->tablcar))
                 { cout << "\n erreur de lecture des parametre de reglage de l'algorithme de Newton: on n'a pas trouve le mot cle "
                        << " fin_parametres_reglage_Algo_Newton_ et par contre la ligne courante contient un mot cle global  ";
                   entreePrinc->MessageBuffer("** erreur5 des parametres de reglage de la loi de comportement de contraintes planes**");
                   throw (UtilLecture::ErrNouvelleDonnee(-1));
                   Sortie(1);     
                 };
                 
                // lecture d'un mot clé
                *(entreePrinc->entree) >> nom;
                
                if ((entreePrinc->entree)->rdstate() == 0) 
                  {} // 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) >>nom;
                   }
                #else
                else  if ((entreePrinc->entree)->eof())
                  // la lecture est bonne mais on a atteind la fin de la ligne
                  { if(nom != "fin_parametres_reglage_Algo_Newton_")
                      {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom;};
                  } 
                #endif  
                else // cas d'une erreur de lecture
                 { cout << "\n erreur de lecture inconnue  ";
                   entreePrinc->MessageBuffer("** erreur4 des parametres de reglage de l'algoritheme de Newton de la loi de comportement de contraintes planes**");
                   throw (UtilLecture::ErrNouvelleDonnee(-1));
                   Sortie(1);     
                 };
                
                // nombre d'itération maxi
                if (nom == "nb_iteration_maxi_")
                 {int nb_boucle_maxi; // le maximum d'itération permis
                  *(entreePrinc->entree) >> nb_boucle_maxi;
                  alg_zero.Modif_iter_max(nb_boucle_maxi);
                 }
                // nombre de dichotomie maxi 
                else if (nom == "nb_dichotomie_maxi_")
                 {int nb_dichotomie; // le maxi de dichotomie prévu pour l'équation de Newton
                  *(entreePrinc->entree) >> nb_dichotomie;
                  alg_zero.Modif_nbMaxiIncre(nb_dichotomie);
                 }
                // tolérance absolue sur le résidu
                else if (nom == "tolerance_residu_")
                 { string mot_cle("=fonction_nD:");
                   if(strstr(entreePrinc->tablcar,mot_cle.c_str())==0)
                     {// lecture du paramètre
                      double  tolerance_residu; // tolérance absolue sur la résolution de l'équation
                      *(entreePrinc->entree) >> tolerance_residu;
                      alg_zero.Modif_prec_res_abs(tolerance_residu);
                     }
                   else // on lit une fonction
                     {// on lit le nom de la fonction
                      string nom_fonct;
                      bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_fonct);
                      if (!lec )
                         { entreePrinc->MessageBuffer("**erreur en lecture**  tolerance_residu_ via une fonction nD");
                           throw (UtilLecture::ErrNouvelleDonnee(-1));
                           Sortie(1);
                         };
                      // maintenant on définit la fonction
                      if (lesFonctionsnD.Existe(nom_fonct))
                       {fct_tolerance_residu = lesFonctionsnD.Trouve(nom_fonct);
                       }
                      else
                       {// sinon il faut la lire maintenant
                        string non("_");
                        fct_tolerance_residu = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
                        // lecture de la courbe
                        fct_tolerance_residu->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
                        // maintenant on vérifie que la fonction est utilisable
                        if (fct_tolerance_residu->NbComposante() != 1 )
                         { cout << "\n erreur en lecture, la fonction " << nom_fonct
                                << " est une fonction vectorielle a   " << fct_tolerance_residu->NbComposante()
                                << " composantes alors qu'elle devrait etre scalaire ! "
                                << " elle n'est donc pas utilisable !! ";
                           string message("\n**erreur** \n"+nom_class_methode+"(...");
                           entreePrinc->MessageBuffer(message);
                           throw (UtilLecture::ErrNouvelleDonnee(-1));
                           Sortie(1);
                         };
                       };
                      // on mettra à jour pendant le calcul, la valeur de la précision
                     };
                 }
                // tolérance relative sur le résidu
                else if (nom == "tolerance_residu_rel_")
                 { string mot_cle("=fonction_nD:");
                   if(strstr(entreePrinc->tablcar,mot_cle.c_str())==0)
                     {// lecture du paramètre
                      double  tolerance_residu_rel; // tolérance absolue sur la résolution de l'équation
                      *(entreePrinc->entree) >> tolerance_residu_rel;
                      alg_zero.Modif_prec_res_rel(tolerance_residu_rel);
                     }
                   else // on lit une fonction
                     {// on lit le nom de la fonction
                      string nom_fonct;
                      bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_fonct);
                      if (!lec )
                         { entreePrinc->MessageBuffer("**erreur en lecture**  tolerance_residu_ via une fonction nD");
                           throw (UtilLecture::ErrNouvelleDonnee(-1));
                           Sortie(1);
                         };
                      // maintenant on définit la fonction
                      if (lesFonctionsnD.Existe(nom_fonct))
                       {fct_tolerance_residu_rel = lesFonctionsnD.Trouve(nom_fonct);
                       }
                      else
                       {// sinon il faut la lire maintenant
                        string non("_");
                        fct_tolerance_residu_rel = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
                        // lecture de la courbe
                        fct_tolerance_residu_rel->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
                        // maintenant on vérifie que la fonction est utilisable
                        if (fct_tolerance_residu_rel->NbComposante() != 1 )
                         { cout << "\n erreur en lecture, la fonction " << nom_fonct
                                << " est une fonction vectorielle a   " << fct_tolerance_residu_rel->NbComposante()
                                << " composantes alors qu'elle devrait etre scalaire ! "
                                << " elle n'est donc pas utilisable !! ";
                           string message("\n**erreur** \n"+nom_class_methode+"(...");
                           entreePrinc->MessageBuffer(message);
                           throw (UtilLecture::ErrNouvelleDonnee(-1));
                           Sortie(1);
                         };
                       };
                      // on mettra à jour pendant le calcul, la valeur de la précision relative
                     };
                 }
                // maxi_delta_var_eps_sur_iter_pour_Newton
                else if (nom == "maxi_delta_var_eps_sur_iter_pour_Newton_")
                 { *(entreePrinc->entree) >> maxi_delta_var_eps_sur_iter_pour_Newton;
                 }
                // le mini de variation de x
                else if (nom == "mini_delta_x_")
                 {double  mini_delta_x;
                  *(entreePrinc->entree) >> mini_delta_x;
                  alg_zero.Modif_mini_delta_x(mini_delta_x);
                 }
                // le coef_mini_delta_x
                else if (nom == "coef_mini_delta_x_")
                 {double  coef_mini_delta_x;
                  *(entreePrinc->entree) >> coef_mini_delta_x;
                  alg_zero.Modif_coef_mini_delta_x(coef_mini_delta_x);
                 }
                // le minimum de hsurh0
                else if (nom == "mini_hsurh0_")
                 {*(entreePrinc->entree) >> mini_hsurh0;
                 }
                // le maximum de hsurh0
                else if (nom == "maxi_hsurh0_")
                 {*(entreePrinc->entree) >> maxi_hsurh0;
                 }
                // forcer un stockage des indicateurs de la résolution
                else if (nom == "sortie_post_")
                 {*(entreePrinc->entree) >> sortie_post;
                 }
                // forcer un affichage particulier pour les méthodes
                else if (nom == "permet_affichage_")
                 {Lecture_permet_affichage(entreePrinc,lesFonctionsnD);
                  alg_zero.Modif_affichage(Permet_affichage()); // on met à jour l'algo de newton
                 }
                // forcer un affichage particulier pour newton
                else if (nom == "permet_affichage_specifique_newton_")
                 {double permet_affichage_specifique_newton=0;
                  *(entreePrinc->entree) >> permet_affichage_specifique_newton;
                  alg_zero.Modif_affichage(permet_affichage_specifique_newton); // on met à jour l'algo de newton
                 }
 
                // sinon ce n'est pas un mot clé connu, on le signale
                else if (nom != "fin_parametres_reglage_Algo_Newton_")
                 { cout << "\n erreur en lecture d'un parametre, le mot cle est inconnu "
                        << " on a lu : " << nom << endl;
                   if (ParaGlob::NiveauImpression()>3)
                        cout << "\n LoiContraintesPlanes::LectureDonneesParticulieres(UtilLecture * entreePrinc) " << endl ;
                   throw (UtilLecture::ErrNouvelleDonnee(-1));
                   Sortie(1);          
                 }
               }; //-- fin du while
             }; //-- fin de la lecture des paramètres de réglage

           break;
         }
        
        
		      default :
			           cout << "\nErreur : valeur incorrecte du type de contrainte lue !: " << type_de_contrainte << " \n";
              entreePrinc->MessageBuffer("**erreur 6 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);
	     };
	// -- on lit maintenant au choix soit rien, soit un facteur de pénalisation et/ou une précision
											 
	 // maintenant lecture de la loi  
    entreePrinc->NouvelleDonnee();
    // lecture du nom de la loi
    string st2,nom3;
    *(entreePrinc->entree)  >> st2;
    // definition de la loi
    LoiAbstraiteGeneral * pt = LesLoisDeComp::Def_loi(st2);
    lois_interne = (Loi_comp_abstraite*) LesLoisDeComp::Def_loi(st2);
    pt = lois_interne;

    // --- lecture des informations particulières propres à la  loi   
    entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
    pt->LectureDonneesParticulieres (entreePrinc,lesCourbes1D,lesFonctionsnD);
    // on s'occupe de la catégorie après la lecture des informations particulières (variable def dans LoiAbstraiteGeneral)
	   LoiAbstraiteGeneral::categorie_loi_comp = pt->Id_categorie();
    if (!GroupeMecanique(categorie_loi_comp))
      { cout << "\n erreur1 en lecture des lois constitutives elementaire d'une loi LoiContraintesPlanes"
             << "\n la loi lue: " << pt->Nom_comport() << " n'est pas une loi mecanique, elle fait partie "
             << " de la categorie: "<< Nom_categorie_loi_comp(categorie_loi_comp);
        entreePrinc->MessageBuffer("**erreur 7 LoiContraintesPlanes::LectureDonneesParticulieres (...**");
        throw (UtilLecture::ErrNouvelleDonnee(-1));
        Sortie(1);     
       };
    
    // si la loi est thermo dépendante on indique que la loi de contrainte plane l'est aussi
    if (((Loi_comp_abstraite*)pt)->ThermoDependante()) this->thermo_dependant = true;
       
//    entreePrinc->NouvelleDonnee(); // prepa du flot de lecture pour d'autre loi éventuellement 
    // appel au niveau de la classe mère
    Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
                 (*entreePrinc,lesFonctionsnD,true);
  };

// affichage de la loi
void LoiContraintesPlanes::Affiche() const 
  {  cout << "\n .......  loi de comportement LoiContraintesPlanes ........";
     cout << "\n type_de_contrainte: " << Nom_contrainte_mathematique(type_de_contrainte)
	         << " prec= " << prec << " ";
	    if (type_de_contrainte == PENALISATION) 
	      { cout << " fac_penal= " << fac_penal << " ";}
     else if (type_de_contrainte == NEWTON_LOCAL)
       { // --- paramètre de réglage
         cout << "\n reglage_algo_newton_equadiff: " ;
         alg_zero.Affiche();
         cout << " maxi_delta_var_eps_sur_iter_pour_Newton_ "<<maxi_delta_var_eps_sur_iter_pour_Newton << " ";
         if (fct_tolerance_residu != NULL)
           {cout << "\n pilotage tol_residu: ";
            if (fct_tolerance_residu->NomFonction() != "_")
              cout << fct_tolerance_residu->NomFonction();
            else
              fct_tolerance_residu->Affiche();
            cout << "\n";
           };
         if (fct_tolerance_residu_rel != NULL)
           {cout << "\n pilotage tol_residu_rel: ";
            if (fct_tolerance_residu_rel->NomFonction() != "_")
              cout << fct_tolerance_residu_rel->NomFonction();
            else
              fct_tolerance_residu_rel->Affiche();
            cout << "\n";
           };
         // niveau d'affichage
         Affiche_niveau_affichage();
         cout << " sortie_post "<< sortie_post
              << " ";
        };
//     // pondération éventuelle
//     if (niveauF_grandeurGlobale != NULL)
//       niveauF_grandeurGlobale->Affiche();
//     if (niveauF_ddlEtendu!= NULL)
//       niveauF_ddlEtendu->Affiche();
//     if (niveauF_temps!= NULL)
//       niveauF_temps->Affiche();

     lois_interne->Affiche();
     cout << "\n .......  fin de la loi de comportement LoiContraintesPlanes ........";     
  };
            
// affichage et definition interactive des commandes particulières à chaques lois
void LoiContraintesPlanes::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
 {ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
	 cout << "\n definition standart pour CP (rep o) ou exemples exhaustifs (rep n'importe quoi) ? ";
	 string rep = "_";     // procédure de lecture avec prise en charge d'un retour chariot
  rep = lect_return_defaut(true,"o");

	 if ((rep == "o") || (rep != "O" ))
    sort  << "\n#   --- exemple de declaration pour une loi elastique --- "
          << "\n#   "
          << "\n   toto LOI_CONTRAINTES_PLANES    "
          << "\n          NEWTON_LOCAL 	avec_parametres_de_reglage_   "
          << "\n                    nb_iteration_maxi_  20   "
          << "\n                    nb_dichotomie_maxi_  20   "
          << "\n                    tolerance_residu_ 1.e-3   "
          << "\n                    tolerance_residu_rel_ 1.e-4   "
          << "\n                 fin_parametres_reglage_Algo_Newton_   "
          << "\n   "
          << "\n     ISOELAS   "
          << "\n       200000  0.3   "
          << "\n   fin_loi_contrainte_plane  # --- fin de la loi de contrainte plane  "
          << "\n\n \n ";
  
	 if ((rep != "o") && (rep != "O" ) && (rep != "0") )
   {sort << "\n# .......  loi de comportement LoiContraintesPlanes ........"
         << "\n#  a)sur la premiere ligne on indique  : "
         << "\n#     Obligatoirement : la methode utilisee pour imposer les contraintes plane: par defaut par multiplicateur de Lagrange "
         << "\n#     les differents choix sont:  "
         << "\n#										   "
         << "\n#	PERTURBATION deformation_epaisseur :  a chaque iteration (en implicite) ou increment (en explicite) la deformation d'epaisseur	   "
         << "\n#                                          est mise a jour, c'est a dire la deformation eps33  "
         << "\n#                                          pour cela on se sert du module de compressibilite et de la condition de "
         << "\n#                                          de contrainte plane "
         << "\n#										   "
         << "\n#										   "
/*
         << "\n#     MULTIPLICATEUR_LAGRANGE : utilisation d'un multiplicateur de Lagrange, l'equilibre n'est cependant pas exact "
         << "\n#										  il depend de la valeur des ddl calcules globalement. La precision de la condition "
         << "\n#										  sig33=0 est consultee lors de la resolution globale, on peut ainsi indiquer une precision  "
         << "\n#										  en dessous de laquelle on considerera que la condition est ok. Par defaut cette precision "
         << "\n#										  est relative et vaut 0.005 * Max|sig_ij|. On peut indiquer a la suite du mot cle prec= valeur"
         << "\n#										  ou prec= est un mot cle facultatif, valeur indique la precision que l'on desire "
         << "\n#										   "
         << "\n#     PENALISATION            : utilisation d'un facteur de penalisation dont la valeur est par defaut "
         << "\n#                               30 fois le maximum de |d_sig_ij/d_ddl| , a la suite du  "
         << "\n#                               mot cle on peut indiquer fac= facteur , ou fac= est un mot cle et facteur "
         << "\n#                               est un nombre qui multiplie par le maxi de |d_sig_ij/d_ddl| sera la penalisation  "
         << "\n#										  comme pour la methode avec multiplicateur de Lagrange, la condition est approchee aussi  "
         << "\n#										  la precision de la contrainte sig33=0 est consultee lors de la resolution globale. On "
         << "\n#										  peut donc aussi indiquer une precision differente de celle par defaut de la meme maniere "
         << "\n#										  que pour le multiplicateu de Lagrange "
         << "\n#										   "
*/
         << "\n#     NEWTON_LOCAL     : utilisation d'une methode interne de Newton pour imposer precisemment la condition "
         << "\n#										  La methode itere localement en chaque point d'integration pour imposer la condition "
         << "\n#										  Par defaut, on considere que la convergence est ok lorsque la condition est satisfaite  "
         << "\n#										  pour une precision relative de 0.005 * Max|sig_ij|. "
         << "\n#			---	exemple de declaration:					---	   "
         << "\n#				NEWTON_LOCAL 				"
         << "\n#										   "
         << "\n# ** il est egalement possible (mais pas obligatoire) de definir des parametres de reglage "
         << "\n# de la resolution. Dans ce cas, a la suite du mot cle NEWTON_LOCAL "
         << "\n# on indique le mot cle: avec_parametres_de_reglage_ "
         << "\n# ensuite on defini (dans un ordre quelconque) les parametres que l'on souhaites "
         << "\n# chaque parametre est precede d'un mot cle, on a donc une succession de mot cle suivi d'une grandeur "
         << "\n# on peut avoir un ou plusieur couple parametre-grandeur sur chaque ligne "
         << "\n# par contre la derniere ligne doit comporter uniquement le mot cle: "
         << "\n#  fin_parametres_reglage_Algo_Newton_ "
         << "\n# les differents parametres sont: "
         << "\n# le nombre d'iteration  ex: nb_iteration_maxi_  20 "
         << "\n# le nombre de dichotomie ex: nb_dichotomie_maxi_  20 "
         << "\n# la tolerance absolue sur le residu ex: tolerance_residu_ 5.e-3 "
         << "\n# la tolerance relative sur le residu ex: tolerance_residu_rel_ 1.e-4 "
         << "\n# le minimum de variation de h sur h0 (par defaut 0.001) ex: mini_hsurh0_ 1.e-4 "
         << "\n# le maximum de variation de h sur h0 (par defaut 1000) ex: maxi_hsurh0_ 4 "
         << "\n# la valeur absolue maximale du delta deformation qui est permise a chaque iteration de Newton  "
         << "\n#   par defaut = 0.1 , si on veut une autre valeur: exe: "
         << "\n# maxi_delta_var_eps_sur_iter_pour_Newton_ 0.2           "
         << "\n#   si on donne une valeur negative, il n'y a plus de limite           "
         << "\n#      "
         << "\n#  le mot cle  sortie_post_ , par defaut il vaut 0, dans ce cas aucun indicateur n'est stoke"
         << "\n#  s'il est different de 0, on peut acceder aux indicateurs en post-traitement (nombre d'iteration, dichotomie ... "
         << "\n#  seules les indicateurs en cours sont disponibles, il n'y a pas de stockage sur plusieurs increment "
         << "\n# "
         << "\n# ex: sortie_post_ 1 "
         << "\n# "
         << "\n# -------------- affichage des erreurs et des warning ---------- "
         << "\n# - l'affichage normale est fonction du parametre global d'affichage gerer par le niveau d'affichage"
         << "\n# cependant pour des raisons par exemple de mise au point, il est possible de permettre l'affichage "
         << "\n# a un niveau particulier (mot cle : permet_affichage_ suivi d'un nombre entier) en plus de l'affichage normal. "
         << "\n# l'affichage s'effectuera donc en fonction de l'affichage normale et de l'affichage particulier."
         << "\n# Le fonctionnement de l'affichage particulier suit les mêmes règles que l'affichage globale"
         << "\n# soit permet_affichage_ est nulle (cas par defaut), dans ce cas l'affichage est fonction du niveau global"
         << "\n# soit permet_affichage_ vaut n par exemple, dans ce cas l'affichage est fonction uniquement de n "
         << "\n# "
         << "\n# ex: permet_affichage_ 5  "
         << "\n#										   "
         << "\n#				-- exemple de declaration:				--		   "
         << "\n#				NEWTON_LOCAL 				avec_parametres_de_reglage_"
         << "\n#                  nb_iteration_maxi_  20 "
         << "\n#                  nb_dichotomie_maxi_  20 "
         << "\n#                  tolerance_residu_ 5.e-3 "
         << "\n#                  tolerance_residu_rel_ 1.e-4 "
         << "\n#               fin_parametres_reglage_Algo_Newton_ "
         << "\n#										   "
         << "\n#             "
         << "\n#  NB: il est possible de piloter les tolerances de l'algo de Newton avec une fonction nD  "
         << "\n#     ce qui permet par exemple d'avoir des tolerances qui varient en fct de la precision globlae "
         << "\n#             "
         << "\n#    exemples:  avec fc1 et fc2, 2 fct nD  "
         << "\n#      tolerance_residu_  =fonction_nD: fc1 "
         << "\n#      tolerance_residu_rel_  =fonction_nD: fc2 "
         << "\n#             "
         << "\n#             "
         << "\n#-----------------------------------"
         << "\n#	b)puis sur la ligne suivante: 									   "
         << "\n#     Obligatoirement:  le nom de la loi de comportement 3D sur laquelle on "
         << "\n#     veut imposer une condition de contrainte plane "
         << "\n#								La suite des informations est relative a la loi 3D, voir donc la syntaxe associee"
         << "\n#										   "
         << "\n#										   "
         << "\n#	exemple 0: 		   "
         << "\n		PERTURBATION deformation_epaisseur  "
         << "\n		  HYSTERESIS_3D	"
         << "\n#		.... # partie relative aux parametres specifiques a la loi d'Hysteresis								   "
         << "\n#										   "
/*
         << "\n#	exemple 1: 		   "
         << "\n#		MULTIPLICATEUR_LAGRANGE prec= 1.e-3 "
         << "\n#		  HYSTERESIS_3D	"
         << "\n#		.... # partie relative aux parametres specifiques a la loi d'Hysteresis								   "
         << "\n#										   "
         << "\n#	exemple 2: 	prec et fac sont facultatifs, mais prec doit etre avant fac	   "
         << "\n#		PENALISATION prec= 1.e-3 fac= 40. 	"
         << "\n#		  HYSTERESIS_3D	"
         << "\n#		.... # partie relative aux parametres specifiques a la loi d'Hysteresis								   "
         << "\n#										   "
         << "\n#	exemple 3: 		   "
         << "\n#		NEWTON_LOCAL_LOCAL prec= 1.e-3  	"
         << "\n#		  HYSTERESIS_3D	"
         << "\n#		.... # partie relative aux parametres specifiques a la loi d'Hysteresis								   "
*/
         << "\n#										   "
         << "\n#	a la fin de fin de la loi, on indique un mot cle de fin		";
    sort << "\n          fin_loi_contrainte_plane   # ----- fin de Loi de contrainte plane " << endl;
   };
  };  		  	  

// test si la loi est complete
int LoiContraintesPlanes::TestComplet()
    { int ret = LoiAbstraiteGeneral::TestComplet();
//      // cas éventuelle des pondérations
//      if (niveauF_grandeurGlobale != NULL)
//        niveauF_grandeurGlobale->Verif_complet();
//      if (niveauF_ddlEtendu!= NULL)
//        niveauF_ddlEtendu->Verif_complet();
//      if (niveauF_temps!= NULL)
//        niveauF_temps->Verif_complet();

      ret *=lois_interne->TestComplet();
      return ret; 
    }; 
		          
// calcul d'un module d'young équivalent à la loi, ceci pour un
// chargement nul
double LoiContraintesPlanes::Module_young_equivalent(Enum_dure temps,const Deformation & def ,SaveResul * saveResul_ex)
  {
    SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul_ex);
    double E =lois_interne->Module_young_equivalent(temps,def,save_resul.le_SaveResul);
    return E;
  };
		          
// récupération d'un module de compressibilité équivalent à la loi, ceci pour un chargement nul
// il s'agit ici de la relation -pression = sigma_trace/3. = module de compressibilité * I_eps
double LoiContraintesPlanes::Module_compressibilite_equivalent(Enum_dure temps,const Deformation & def ,SaveResul * saveResul_ex)
  {
    SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul_ex);
    double module_compressibilite =lois_interne->Module_compressibilite_equivalent(temps,def,save_resul.le_SaveResul);
    return module_compressibilite;
  };
		  	
  
 
// récupération de la dernière déformation d'épaisseur calculée: cette déformaion n'est utile que pour des lois en contraintes planes
// - pour les lois 3D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// - pour les lois 2D def planes: retour de 0
// les infos nécessaires à la récupération de la def, sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double LoiContraintesPlanes::Eps33BH(SaveResul * saveResul) const
  { // retour de la def
    return eps_BB_3D(3,3);
  };
   
// récupération de la variation relative d'épaisseur calculée: h/h0
//  cette variation n'est utile que pour des lois en contraintes planes
// - pour les lois 3D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// - pour les lois 2D def planes: retour de 0
// les infos nécessaires à la récupération , sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double LoiContraintesPlanes::HsurH0(SaveResul * saveResul) const
 { // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
   // retour de la variation relative
   return save_resul.hsurh0;
 };


// activation des données des noeuds et/ou elements nécessaires au fonctionnement de la loi
// exemple: mise en service des ddl de température aux noeuds
// ici la grandeur qui sert de proportion entre la première loi et la seconde
void LoiContraintesPlanes::Activation_donnees(Tableau<Noeud *>& tabnoeud,bool dilatation,LesPtIntegMecaInterne& lesPtMecaInt)
  { // pondération éventuelle
//    if (niveauF_ddlEtendu != NULL)
//        niveauF_ddlEtendu->Activation_donnees(tabnoeud, dilatation, lesPtMecaInt);
    // appel relatif à la lois associée
    lois_interne->Activation_donnees(tabnoeud,dilatation,lesPtMecaInt);    
    // appel de la méthode de la classe mère
    Loi_comp_abstraite::Activ_donnees(tabnoeud,dilatation,lesPtMecaInt);     
  };
   
// récupération des grandeurs particulière (hors ddl )
// correspondant à liTQ
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void LoiContraintesPlanes::Grandeur_particuliere
              (bool absolue,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list<int>& decal) const
 { // tout d'abord on récupère le conteneur
   SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveDon);
   int dim = ParaGlob::Dimension();
   // maintenant on s'occupe des grandeurs de la loi  elle même,
   List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end();
   list<int>::iterator idecal=decal.begin();
   for (itq=liTQ.begin();itq!=itqfin;itq++,idecal++)
     { 
      SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveDon);
      TypeQuelconque& tipParticu = (*itq); // pour simplifier
      if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
       switch (tipParticu.EnuTypeQuelconque().EnumTQ())
       {  case CONTRAINTE_INDIVIDUELLE_A_CHAQUE_LOI_A_T:
           // 2) -----cas des contraintes individuelles à chaque loi à t uniquement
            { Tab_Grandeur_TenseurHH& tyTQ= *((Tab_Grandeur_TenseurHH*) (*itq).Grandeur_pointee()); // pour simplifier
              // en fait on utilise systématiquement un tenseur d'ordre le + élevé, car c'est le conteneur le plus générique
              // et Tab_Grandeur_TenseurHH ne supporte que des tenseurs du même ordre donc s'il y a un tenseur élevé
              // interne il faut que tous les tenseurs soient du même ordre
              TenseurHH* sigHH = (save_resul.l_sigoHH_t); // pour simplifier
              if (Dabs(sigHH->Dimension()) != dim)
                {tyTQ(1+(*idecal)).Affectation_trans_dimension(*sigHH,true);
                }
              else // cas même dimension
                {tyTQ(1+(*idecal)) = (*(save_resul.l_sigoHH_t));
                };
              (*idecal)++;
              break;
            }
         case DEF_MECANIQUE:
           // 2) -----cas des déformations dites mécaniques, éventuelle
            { Tab_Grandeur_TenseurBB& tyTQ= *((Tab_Grandeur_TenseurBB*) (*itq).Grandeur_pointee()); // pour simplifier
              // en fait on utilise systématiquement un tenseur d'ordre le + élevé, car c'est le conteneur le plus générique
              // et Tab_Grandeur_TenseurBB ne supporte que des tenseurs du même ordre donc s'il y a un tenseur élevé
              // interne il faut que tous les tenseurs soient du même ordre
              TenseurBB* eps_P_mecaBB_t = (save_resul.eps_P_mecaBB_t); // pour simplifier
              if (eps_P_mecaBB_t != NULL)
               {if (Dabs(eps_P_mecaBB_t->Dimension()) != dim)
                  {tyTQ(1+(*idecal)).Affectation_trans_dimension(*eps_P_mecaBB_t,true);
                  }
                else // cas même dimension
                  {tyTQ(1+(*idecal)) = (*(save_resul.eps_P_mecaBB_t));
                  };
               }
              else // sinon on met à 0
                {tyTQ(1+(*idecal)).Inita(0.);};
              (*idecal)++;
              break;
            }
          case ENERGIE_ELASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T:
           // 3) -----cas de l'énergie élastique individuelles à chaque loi à t
            { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
              tyTQ(1+(*idecal)) = save_resul.l_energ_t.EnergieElastique();  (*idecal)++;
              break;
            }
//          case SIG_EPAISSEUR:
//           // 3) -----cas de la contrainte d'épaisseur  à la loi à t: en fait il s'agit de sig33
//            { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
//              tyTQ(1+(*idecal)) = (*(save_resul.l_sigoHH))(3,3);  (*idecal)++;
//              break;
//            }
          case ENERGIE_PLASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T:
          // 4) -----cas de l'énergie plastique individuelles à chaque loi à t
          if ((*itq).EnuTypeQuelconque() == ENERGIE_PLASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T)
           { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
             tyTQ(1+(*idecal)) = save_resul.l_energ_t.DissipationPlastique();  (*idecal)++;
              break;
           }
          case ENERGIE_VISQUEUSE_INDIVIDUELLE_A_CHAQUE_LOI_A_T:
          // 5) -----cas de l'énergie visqueuse individuelles à chaque loi à t
          if ((*itq).EnuTypeQuelconque() == ENERGIE_VISQUEUSE_INDIVIDUELLE_A_CHAQUE_LOI_A_T)
           { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
             tyTQ(1+(*idecal)) = save_resul.l_energ_t.DissipationVisqueuse();  (*idecal)++;
              break;
           }
          case DEF_EPAISSEUR:
           // 6) -----cas de la déformation d'épaisseur à t
           { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
             // -- on calcul en fonction de l'élongation d'épaisseur
             double hsurh0=1. ; // init par défaut: on utilise une variable inter, car s'il n'y a pas eu
                 // de calcul préalable,save_resul.hsurh0 == 0 ce qui conduit à une def infinie
             if (save_resul.hsurh0 != 0.) // cas où il y a  eu un calcul
                hsurh0=save_resul.hsurh0;
             
             switch (type_de_deformation)
                {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
                      // cas d'une déformation d'Almansi
                      {  // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1.
                         // donc dans le repère global : epsBB33 = 1/2 * (1. - 1./(h/h0)^2)
                         tyTQ(1+(*idecal)) = 0.5 * (1. - 1./(hsurh0 * hsurh0));(*idecal)++;
                      };
                      break;
                 case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
                     // cas d'une def logarithmique ou une approximation
                      {  tyTQ(1+(*idecal)) = log(hsurh0);(*idecal)++;
                      };
                      break;
                 default :
                   cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
                        << Nom_type_deformation(type_de_deformation);
                   cout << "\n LoiContraintesPlanes::Grandeur_particuliere \n";
                   Sortie(1);
                };
              break;
           
           }
          case NB_INCRE_TOTAL_RESIDU:
           { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
             if ((save_resul.indicateurs_resolution_t.Taille()))
                  tyTQ(1+(*idecal)) =  save_resul.indicateurs_resolution_t(1);
             else tyTQ(1+(*idecal)) = 0.;
             (*idecal)++; break;
           }
          case NB_ITER_TOTAL_RESIDU:
           { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
             if ((save_resul.indicateurs_resolution_t.Taille()))
                  tyTQ(1+(*idecal)) =  save_resul.indicateurs_resolution_t(2);
             else tyTQ(1+(*idecal)) = 0.;
             (*idecal)++; break;
           }
         
          default: ;// on ne fait rien
       };
     };
   // puis appel pour la lois associée
   lois_interne->Grandeur_particuliere(absolue,liTQ,save_resul.le_SaveResul,decal);
 };

// récupération de la liste de tous les grandeurs particulières
// ces grandeurs sont ajoutées à la liste passées en paramètres
void LoiContraintesPlanes::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
 { // tout d'abord on passe en revue les grandeurs des lois associées
   // ** au niveau de l'exécution ce sera l'inverse de cette ordre: on s'occupera d'abord de this puis les lois internes
	  // ** mais a priori cela n'a pas d'importance 
   // appel de la loi 3D
	  int nb_loi = 1; // pour être identique pour la loi additive !!
   lois_interne->ListeGrandeurs_particulieres(absolue,liTQ);
	
   // ... maintenant on s'occupe des grandeurs de la loi  elle même,   
   // 1) -----cas des contraintes individuelles à la loi à t uniquement
   // ici il s'agit du tenseur des contraintes du pas précédent ou du pas actuel
   //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
   {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
    for (itq=liTQ.begin();itq!=itqfin;itq++)
      if ((*itq).EnuTypeQuelconque() == CONTRAINTE_INDIVIDUELLE_A_CHAQUE_LOI_A_T)
        { Tab_Grandeur_TenseurHH& tyTQ= *((Tab_Grandeur_TenseurHH*) (*itq).Grandeur_pointee()); // pour simplifier
		        int taille = tyTQ.Taille()+nb_loi;
			       tyTQ.Change_taille(taille); nexistePas = false;
		      };
    if (nexistePas)
	      { TenseurHH* tens = NevezTenseurHH(ParaGlob::Dimension()); // un tenseur typique
                                               // qui fonctionnera en absolue ou non
         Tab_Grandeur_TenseurHH gtHH(*tens,nb_loi);
         // def d'un type quelconque représentatif 
         TypeQuelconque typQ(CONTRAINTE_INDIVIDUELLE_A_CHAQUE_LOI_A_T,SIG11,gtHH);
         liTQ.push_back(typQ);
         delete tens; // car on n'en a plus besoin
		     };
   };
   // cas de la déformation mécanique éventuelle: comme on ne peut pas savoir si elle existe
   // sans avoir accès au save result, on définit un conteneur possible
   {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
    for (itq=liTQ.begin();itq!=itqfin;itq++)
      if ((*itq).EnuTypeQuelconque() == DEF_MECANIQUE)
        { Tab_Grandeur_TenseurBB& tyTQ= *((Tab_Grandeur_TenseurBB*) (*itq).Grandeur_pointee()); // pour simplifier
          int taille = tyTQ.Taille()+nb_loi;
          tyTQ.Change_taille(taille); nexistePas = false;
        };
    if (nexistePas)
       { TenseurBB* tens = NevezTenseurBB(ParaGlob::Dimension()); // un tenseur typique
                                               // qui fonctionnera en absolue ou non
         Tab_Grandeur_TenseurBB gtBB(*tens,nb_loi);
         // def d'un type quelconque représentatif
         TypeQuelconque typQ(DEF_MECANIQUE,SIG11,gtBB);
         liTQ.push_back(typQ);
         delete tens; // car on n'en a plus besoin
       };
   };
   // 2) -----cas de la contrainte sig33 spécifiquement
   Tableau <double> tab_double(nb_loi);   
   Tab_Grandeur_scalaire_double grand_courant(tab_double); 
//   //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
//   {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
//    for (itq=liTQ.begin();itq!=itqfin;itq++)
//      if ((*itq).EnuTypeQuelconque() == SIG_EPAISSEUR)
//        { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
//		        int taille = tyTQ.Taille()+nb_loi;
//			       tyTQ.Change_taille(taille); nexistePas = false;
//		      };
//    if (nexistePas)
//	     {TypeQuelconque typQ1(SIG_EPAISSEUR,SIG11,grand_courant);
//       liTQ.push_back(typQ1);
//		    };
//   };
   // pour toutes les énergies
   // 3) -----cas de l'énergie élastique individuelles à chaque loi à t uniquement
   //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
	  {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
	   for (itq=liTQ.begin();itq!=itqfin;itq++)
      if ((*itq).EnuTypeQuelconque() == ENERGIE_ELASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T)
        { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
		        int taille = tyTQ.Taille()+nb_loi;
			       tyTQ.Change_taille(taille); nexistePas = false;
        };
	   if (nexistePas)
	     {TypeQuelconque typQ1(ENERGIE_ELASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T,SIG11,grand_courant);
       liTQ.push_back(typQ1);
		    };
   };
   // 4) -----cas de l'énergie plastique individuelles à chaque loi à t uniquement
   //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
	  {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
	   for (itq=liTQ.begin();itq!=itqfin;itq++)
      if ((*itq).EnuTypeQuelconque() == ENERGIE_PLASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T)
        { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
		        int taille = tyTQ.Taille()+nb_loi;
			       tyTQ.Change_taille(taille); nexistePas = false;
		      };
	   if (nexistePas)
      {TypeQuelconque typQ1(ENERGIE_PLASTIQUE_INDIVIDUELLE_A_CHAQUE_LOI_A_T,SIG11,grand_courant);
       liTQ.push_back(typQ1);
		    };
   };
   // 5) -----cas de l'énergie visqueuse individuelles à chaque loi à t uniquement
   //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
	  {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
	   for (itq=liTQ.begin();itq!=itqfin;itq++)
      if ((*itq).EnuTypeQuelconque() == ENERGIE_VISQUEUSE_INDIVIDUELLE_A_CHAQUE_LOI_A_T)
        { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
		        int taille = tyTQ.Taille()+nb_loi;
			       tyTQ.Change_taille(taille); nexistePas = false;
		      };
	   if (nexistePas)
      {TypeQuelconque typQ1(ENERGIE_VISQUEUSE_INDIVIDUELLE_A_CHAQUE_LOI_A_T,SIG11,grand_courant);
       liTQ.push_back(typQ1);
		    };
   };
   // 6) -----cas de la déformation d'épaisseur à t uniquement
   //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
	  {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
	   for (itq=liTQ.begin();itq!=itqfin;itq++)
      if ((*itq).EnuTypeQuelconque() == DEF_EPAISSEUR)
        { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
		        int taille = tyTQ.Taille()+nb_loi;
			       tyTQ.Change_taille(taille); nexistePas = false;
		      };
	   if (nexistePas)
	     {TypeQuelconque typQ1(DEF_EPAISSEUR,EPS11,grand_courant);
       liTQ.push_back(typQ1);
		    };
   };


   // ---- la suite dépend de l'indicateur : sortie_post
   if (sortie_post)
    {// j) ----- NB_INCRE_TOTAL_RESIDU
     {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
      for (itq=liTQ.begin();itq!=itqfin;itq++)
        if ((*itq).EnuTypeQuelconque() == NB_INCRE_TOTAL_RESIDU)
         {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
          int taille = tyTQ.Taille()+1;
          tyTQ.Change_taille(taille); nexistePas = false;
         };
      if (nexistePas)
        {TypeQuelconque typQ1(NB_INCRE_TOTAL_RESIDU,SIG11,grand_courant);
         liTQ.push_back(typQ1);
        };
     };
     // k) ----- NB_ITER_TOTAL_RESIDU
     {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
      for (itq=liTQ.begin();itq!=itqfin;itq++)
        if ((*itq).EnuTypeQuelconque() == NB_ITER_TOTAL_RESIDU)
         {Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
          int taille = tyTQ.Taille()+1;
          tyTQ.Change_taille(taille); nexistePas = false;
         };
      if (nexistePas)
        {TypeQuelconque typQ1(NB_ITER_TOTAL_RESIDU,SIG11,grand_courant);
         liTQ.push_back(typQ1);
        };
     };
    }; // fin du cas ou sortie_post est actif, c-a-d que l'on veut des infos sur les indicateurs  de résolution




 };
  
	   //----- lecture écriture de restart -----
	   // 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 LoiContraintesPlanes::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
  { if (cas == 1) 
     { string st1; string nom;
       ent >> st1 >> nom >> type_de_contrainte >> nom >> fac_penal >> nom >> prec;
//       // niveau mini sig33
//       ent >> st1  >> niveau_mini_sig33;
//       // pondération éventuelle du niveau mini de sig33
//       ent >> nom;
//       if (nom == "avec_Ponderation_GGlobal")
//        {if (niveauF_grandeurGlobale==NULL)
//           niveauF_grandeurGlobale = new Ponderation_GGlobal;
//         niveauF_grandeurGlobale->Lecture_base_info(ent,cas, lesFonctionsnD);
//        }
//       else
//        { if (niveauF_grandeurGlobale!=NULL)
//            delete niveauF_grandeurGlobale;
//        };
//       // les pondérations avec les ddl
//       ent >> nom;
//       if (nom == "avec_Ponderation_ddl")
//        {if (niveauF_ddlEtendu==NULL)
//           niveauF_ddlEtendu = new Ponderation;
//         niveauF_ddlEtendu->Lecture_base_info(ent,cas, lesCourbes1D);
//        }
//       else
//        { if (niveauF_ddlEtendu!=NULL)
//            delete niveauF_ddlEtendu;
//        };
//       // la pondération temps
//       ent >> nom;
//       if (nom == "avec_Ponderation_temps")
//        {if (niveauF_temps==NULL)
//           niveauF_temps = new Ponderation_temps;
//         niveauF_temps->Lecture_base_info(ent,cas, lesCourbes1D);
//        }
//       else
//        { if (niveauF_temps!=NULL)
//            delete niveauF_temps;
//        };
      
       // --- paramètre de réglage
       ent >> nom ; // lecture de "parametre_algo_newton_equadiff:"
       alg_zero.Lecture_base_info(ent,cas);
       // le paramètre maxi_delta_var_eps_sur_iter_pour_Newton
       ent >> nom >> maxi_delta_var_eps_sur_iter_pour_Newton ;
       // pilotage éventuel de la précision
       ent >> nom;
       if (nom != "non_pilot_tol_residu")
        { ent >> nom; // lecture fct locale ou globale
          if (nom == "fct_globale")
            {ent >> nom; // lecture du nom de la fonction
             if (fct_tolerance_residu != NULL)
               delete fct_tolerance_residu;
             // maintenant on définit la fonction
             if (lesFonctionsnD.Existe(nom))
                  {fct_tolerance_residu = lesFonctionsnD.Trouve(nom);}
             else { cout << "\n *** erreur en lecture de la fonction nD " << nom
                         << " on ne la trouve pas !! ";
                    Sortie(1);
                  };
            }
          else // sinon c'est une fonction locale
           { if (nom != "fct_locale")
               { cout << "\n *** erreur en lecture de la fonction nD "
                      << " on attendait le mot cle fct_locale !! ";
                 Sortie(1);
               }
             else
              {// on lit le type de fonction
               ent >> nom; EnumFonction_nD enu = Id_Nom_Fonction_nD(nom);
               nom = "_";
               fct_tolerance_residu = Fonction_nD::New_Fonction_nD(nom, enu);
               fct_tolerance_residu->Lecture_base_info(ent,cas);
              };
           };
        };
      
       if (nom != "non_pilot_tol_residu")
        { ent >> nom; // lecture fct locale ou globale
          if (nom == "fct_globale")
            {ent >> nom; // lecture du nom de la fonction
             if (fct_tolerance_residu_rel != NULL)
               delete fct_tolerance_residu_rel;
             // maintenant on définit la fonction
             if (lesFonctionsnD.Existe(nom))
                  {fct_tolerance_residu_rel = lesFonctionsnD.Trouve(nom);}
             else { cout << "\n *** erreur en lecture de la fonction nD " << nom
                         << " on ne la trouve pas !! ";
                    Sortie(1);
                  };
            }
          else // sinon c'est une fonction locale
           { if (nom != "fct_locale")
               { cout << "\n *** erreur en lecture de la fonction nD "
                      << " on attendait le mot cle fct_locale !! ";
                 Sortie(1);
               }
             else
              {// on lit le type de fonction
               ent >> nom; EnumFonction_nD enu = Id_Nom_Fonction_nD(nom);
               nom = "_";
               fct_tolerance_residu_rel = Fonction_nD::New_Fonction_nD(nom, enu);
               fct_tolerance_residu_rel->Lecture_base_info(ent,cas);
              };
           };
        };
      
       // les autres paramètres
       ent >>  nom; // entête
       ent >>  nom >>  mini_hsurh0;
       ent >>  nom >>  maxi_hsurh0;
       // le niveau d'affichage
       Lecture_permet_affichage(ent,cas,lesFonctionsnD);
       ent >>  nom >> sortie_post;
       alg_zero.Modif_affichage(Permet_affichage()); // on met à jour l'algo de newton
       // --- la loi associée
       ent >> nom ;
 	     if (st1 != "LOI_CONTRAINTES_PLANES")
	       { cout << "\n erreur en lecture de la loi : LOI_CONTRAINTES_PLANES, on attendait le mot cle : LOI_CONTRAINTES_PLANES "
	              << " et on a lue: " << st1
					          << "\n LoiContraintesPlanes::Lecture_base_info_loi(...";
	         Sortie(1);
        };

       // on lit la loi 3D
       ent >> st1; lois_interne = (Loi_comp_abstraite *) LesLoisDeComp::Def_loi(st1);
       lois_interne->Lecture_base_info_loi(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
     }
    else
     { // on utilise directement la  loi déjà défini
       lois_interne->Lecture_base_info_loi(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
     };
  }; 

       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void LoiContraintesPlanes::Ecriture_base_info_loi(ofstream& sort,const int cas)
  { if (cas == 1) 
     { sort << "\n  LOI_CONTRAINTES_PLANES " << " type_de_contrainte= " << type_de_contrainte 
            << " fac_penal= " << fac_penal << " prec= " << prec ;
      
//       // niveau mini sig33
//       sort << "\n niveau_mini_sig33= "<< niveau_mini_sig33;
//       // pondération éventuelle du niveau mini de sig33
//       bool sans_courbe = true;
//       if (niveauF_grandeurGlobale != NULL)
//         {sort << " avec_Ponderation_GGlobal ";
//          niveauF_grandeurGlobale->Ecriture_base_info(sort,cas,sans_courbe);
//         }
//       else {sort << " sans_Ponderation_GGlobal ";};
//       if (niveauF_ddlEtendu!= NULL)
//         {sort << " avec_Ponderation_ddl ";
//          niveauF_ddlEtendu->Ecriture_base_info(sort,cas);
//         }
//       else {sort << " sans_Ponderation_ddl ";};
//       if (niveauF_temps!= NULL)
//         {sort << " avec_Ponderation_temps ";
//          niveauF_temps->Ecriture_base_info(sort,cas);
//         }
//       else {sort << " sans_Ponderation_temps ";};
//      
      
       // --- paramètre de réglage
       sort << "\n parametre_algo_newton_equadiff: ";
       alg_zero.Ecriture_base_info(sort,cas);
       // le paramètre maxi_delta_var_eps_sur_iter_pour_Newton
       sort << " maxi_delta_var_eps_sur_iter_pour_Newton_ "<<maxi_delta_var_eps_sur_iter_pour_Newton << " ";
       // pilotage éventuel de la précision
       if (fct_tolerance_residu == NULL)
          sort << " non_pilot_tol_residu ";
       else
          {sort << " pilot_tol_residu ";
           if (fct_tolerance_residu->NomFonction() != "_")
             sort << " fct_globale " << fct_tolerance_residu->NomFonction();
           else
             {sort << " fct_locale " << Nom_Fonction_nD(fct_tolerance_residu->Type_Fonction());
              fct_tolerance_residu->Ecriture_base_info(sort,cas);
             }
           sort << "\n";
          };
       if (fct_tolerance_residu_rel == NULL)
          sort << " non_pilot_tol_resi_rel ";
       else
          {sort << " pilot_tol_resi_rel ";
           if (fct_tolerance_residu_rel->NomFonction() != "_")
             sort << " fct_globale " << fct_tolerance_residu_rel->NomFonction();
           else
             {sort << " fct_locale "<< Nom_Fonction_nD(fct_tolerance_residu->Type_Fonction());
              fct_tolerance_residu_rel->Ecriture_base_info(sort,cas);
             }
           sort << "\n";
          };
      
       // les autres paramètres
       sort << "\n autres-parametres_loiCP:--> ";
       sort << "\n mini_hsurh0: "<<mini_hsurh0;
       sort << "\n maxi_hsurh0: "<<maxi_hsurh0;
       // niveau d'affichage
       Affiche_niveau_affichage(sort,cas);
       sort << " sortie_post " << sortie_post
            << " ";
       // --- la loi associée
       sort <<"\n -->loi_associee_CP: ";

       sort << lois_interne->Nom_comport() << "  "; lois_interne->Ecriture_base_info_loi(sort,cas);
     }
    else
     { lois_interne->Ecriture_base_info_loi(sort,cas);
     };
  }; 
            

 // ========== codage des METHODES VIRTUELLES  protegees:================
        // calcul des contraintes a t+dt
void LoiContraintesPlanes::Calcul_SigmaHH (TenseurHH& sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl,
             TenseurBB & gijBB_t,TenseurHH & gijHH_t,BaseB& giB,BaseH& gi_H, TenseurBB & epsBB_,
             TenseurBB & delta_epsBB,
             TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,
             double& jacobien_0,double& jacobien,TenseurHH & sigHH
		  	,EnergieMeca & energ,const EnergieMeca & 
		  	,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Expli_t_tdt& ex)
{ // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
  module_compressibilite=module_cisaillement=0.; // init
  energ.Inita(0.); // initialisation des énergies mises en jeux
  
  bool affichage = (Permet_affichage() > 3);
  if (affichage)
   {cout << "\n  --- LoiContraintesPlanes::Calcul_SigmaHH --- ";
    Signature_pti_encours(cout);
   };
  if (Permet_affichage() > 3)
   {cout << "\n ===  donnees d'entree ";
    cout << "\n epsBB_tdt(local)= "; epsBB_.Ecriture(cout);
    cout << "\n DepsBB(local)= "; DepsBB.Ecriture(cout);
    cout << "\n delta_epsBB(local)= "; delta_epsBB.Ecriture(cout);
    cout << "\n sigHH_t(local)= "; sigHH_t.Ecriture(cout);
    // en absolue
    Tenseur3BB tiutiu;
    epsBB_.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n eps en absolu :";tiutiu.Ecriture(cout);
    DepsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n Deps en absolu :";tiutiu.Ecriture(cout);
    delta_epsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n delta_eps en absolu :";tiutiu.Ecriture(cout);
    Tenseur3HH titi;
    sigHH_t.BaseAbsolue(titi,*(ex.giB_tdt));
    cout << "\n sig_t en absolu :";titi.Ecriture(cout);
   };

  // initialisation du tenseurs contrainte 
//  sig_HH_3D.Inita(0.);
  // pour les contraintes
  bool plusZero = true;
  sig_HH_t_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_t),plusZero);
//  sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH),plusZero);
  // passage des informations spécifique à la loi le_SaveResul
  lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
  lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
  lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
  // on sauvegarde la nouvelle épaisseur, car elle va être modifiée et si on doit
  // repartir au départ, ce ne sera plus possible sans cette sauvegarde
  double sauve_hsurh0=save_resul.hsurh0;

  // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
  Vecteur* d_jacobien_tdt = NULL;      // ne sert pas ici
  Tableau <TenseurBB *>* d_epsBB=NULL; //   "    "
//  Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_,d_epsBB,delta_epsBB
//                                               ,jacobien_0,jacobien,d_jacobien_tdt);
	
  // choix entre calcul avec une boucle locale de Newton on non. Dans le premier cas il nous faut la version d_sig/d_eps
  // au lieu de d_sig/d_ddl
  bool mauvaise_convergence = false; // init
 
//cout << "\n boucle 1000 pour debug LoiContraintesPlanes::Calcul_SigmaHH "<<endl;
//for (int jj=1;jj<=1000;jj++)
 
  switch (type_de_contrainte)
   {case NEWTON_LOCAL: // cas de la version d_sig/d_eps
     {
       // initialisation du tenseurs contrainte
       sig_HH_3D.Inita(0.);
       sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH),plusZero);
       // passage des métriques de l'ordre 2 vers 3 
       Passage_metrique_ordre2_vers_3(ex);
       // prise en compte de h sur métriques
       Prise_en_compte_h_sur_metrique();
       // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
       Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_,d_epsBB,delta_epsBB
                                                    ,jacobien_0,jacobien,d_jacobien_tdt);
     
       // on appel la procédure de résolution de sig33(hsurh0)=0
       if ((sortie_post)&&(save_resul.indicateurs_resolution.Taille()!= 2)) // dimensionnement éventuelle de la sortie d'indicateurs
          save_resul.indicateurs_resolution.Change_taille(2);
       // ----- pour ce faire on appelle une méthode de recherche de zero
       val_initiale(1)=save_resul.h_tsurh0; // on démarre la recherche à la valeur à t
       // on impose que les grandeurs soient dans les limites admises
       if (Limitation_h(val_initiale(1)))
         {if (Permet_affichage() > 3)
            cout << "\n loi CP: (newton local) val initiales:  " << flush;
         };
       racine.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
       der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
       try // on met le bloc sous surveillance
        {
         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs, pas utilisées ici
         
         // dans le cas où on utilise une précision qui est pilotée
         {// opération de transmission de la métrique: encapsulé ici
          const Met_abstraite::Impli* ex_impli = NULL;
          const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &ex;
          const Met_abstraite::Umat_cont* ex_expli = NULL;

          if (fct_tolerance_residu != NULL)
           {// ici on utilise les variables connues aux pti, ou calculées à partir de
            // on commence par récupérer les conteneurs des grandeurs à fournir
            List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu->Li_enu_etendu_scalaire();
            List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu->Li_equi_Quel_evolue();
            bool absolue = true; // on se place systématiquement en absolu
            // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
            // pour les grandeurs strictement scalaire
             Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
                         (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
                                           );
            // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
            // pour les Coordonnees et Tenseur
            Valeurs_Tensorielles_interpoler_ou_calculer
                          (absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
            // calcul de la valeur et retour dans tab_ret
            Tableau <double> & tab_val = fct_tolerance_residu->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
            #ifdef MISE_AU_POINT
            if (tab_val.Taille() != 1)
               { cout << "\n loi CP: Erreur : la fonction nD de pilotage de la precision "
                      << " doit calculer un scalaire or le tableau de retour est de taille "
                      << tab_val.Taille() << " ce n'est pas normal !\n";
                 cout << " LoiContraintesPlanesDouble::Calcul_SigmaHH\n";
                 Sortie(1);
               };
            #endif
            // on récupère le premier élément du tableau uniquement
            double tol = tab_val(1);
              #ifdef MISE_AU_POINT
              if (Permet_affichage() > 5)
                 cout << "\n Newton_prec: tol_= "<< tol;
              #endif
            alg_zero.Modif_prec_res_abs(tol);
          };

          if (fct_tolerance_residu_rel != NULL)
           {// ici on utilise les variables connues aux pti, ou calculées à partir de
            // on commence par récupérer les conteneurs des grandeurs à fournir
            List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu_rel->Li_enu_etendu_scalaire();
            List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu_rel->Li_equi_Quel_evolue();
            bool absolue = true; // on se place systématiquement en absolu
            // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
            // pour les grandeurs strictement scalaire
             Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
                         (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
                                           );
            // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
            // pour les Coordonnees et Tenseur
            Valeurs_Tensorielles_interpoler_ou_calculer
                          (absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
            // calcul de la valeur et retour dans tab_ret
            Tableau <double> & tab_val = fct_tolerance_residu_rel->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
            #ifdef MISE_AU_POINT
            if (tab_val.Taille() != 1)
               { cout << "\n loi CP: Erreur : la fonction nD de pilotage de la precision relative "
                      << " doit calculer un scalaire or le tableau de retour est de taille "
                      << tab_val.Taille() << " ce n'est pas normal !\n";
                 cout << " LoiContraintesPlanesDouble::Calcul_SigmaHH\n";
                 Sortie(1);
               };
            #endif
            // on récupère le premier élément du tableau uniquement
            double tol_rel = tab_val(1);
            #ifdef MISE_AU_POINT
            if (Permet_affichage() > 5)
               cout << ", tol_relative= "<< tol_rel;
            #endif
            alg_zero.Modif_prec_res_rel(tol_rel);
          };
         };
         
         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
         bool conver=alg_zero.Newton_raphson
                  (*this,&LoiContraintesPlanes::Residu_constitutif,&LoiContraintesPlanes::Mat_tangente_constitutif
                   ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                   ,maxi_delta_var_eps_sur_iter_pour_Newton);
         if(sortie_post) // sauvegarde éventuelle des indicateurs
              {save_resul.indicateurs_resolution(1)+=nb_incr_total;
               save_resul.indicateurs_resolution(2)+=nb_iter_total;
              };
         // on vérifie qu'il n'y a pas de pb de convergence
         double absracinemax=racine.Max_val_abs();
         if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
           { if (Permet_affichage() > 0)
               { cout << "\n loi CP: non convergence sur l'algo de la resolution de sig33(hsurh0)=0 "
                      << " h/h_0(t+dt)= " << racine(1) << " a t on avait: " << save_resul.h_tsurh0
                      << "\n nb_incr_total=" << nb_incr_total
                      << " nb_iter_total=" << nb_iter_total;
               };
             // on garde en mémoire
             mauvaise_convergence=true;
           };
         // on vérifie que les nouvelles dimensions transversales ne sont pas négatives
         if (racine(1) < 0.)
             { if (Permet_affichage() > 0)
                 { cout << "\n loi CP: **** erreur dans la resolution : on obtient une "
                        << " epaisseur  negative "
                        << ", h/h_0(t+dt)= " << racine(1) << " a t on avait: " << save_resul.h_tsurh0;
                 };
               // on garde en mémoire
               mauvaise_convergence=true;
             };
        }
       catch (ErrNonConvergence_Newton erreur)
         {
          if (Permet_affichage() > 0)
               { cout << "\n loi CP: erreur de non convergence avec l'algo de la resolution de sig33(hsurh0)=0 "
                      << " on obtient une valeur infinie  ou NaN ";
               };
          // on garde en mémoire
          mauvaise_convergence=true;
         }
       catch (ErrSortieFinale)
            // cas d'une direction voulue vers la sortie
            // on relance l'interuption pour le niveau supérieur
          { ErrSortieFinale toto;
            throw (toto);
          }
       catch ( ... )
        { // dans le cas d'une erreur inconnue, on génère également une exception
          if (Permet_affichage() > 0)
               { cout << "\n loi CP: erreur non identifiee sur l'algo de la resolution de sig33(hsurh0)=0 ";
//// ----- debug
//cout << "\n debug :  LoiContraintesPlanes::Calcul_SigmaHH ( ";
//       racine.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
//       der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
//         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs, pas utilisées ici
//         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
//         bool conver=alg_zero.Newton_raphson
//                  (*this,&LoiContraintesPlanes::Residu_constitutif,&LoiContraintesPlanes::Mat_tangente_constitutif
//                   ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total);
//// ---- fin debug
               };
          // on garde en mémoire
          mauvaise_convergence=true;
        };
       if (!mauvaise_convergence)
         { // arrivée ici, cela veut dire que tout à bien fonctionné
           save_resul.hsurh0 = racine(1); // récup de la solution
           // on met à jour les modules
           module_compressibilite= module_compressibilite_3D;
           module_cisaillement= module_cisaillement_3D;
           #ifdef MISE_AU_POINT
            if  (Permet_affichage() > 7)
            { cout << "\n -- loi CP: bonne convergence de Newton ";
              TenseurHH* ptHH = NevezTenseurHH(sig_HH_3D);
              sig_HH_3D.BaseAbsolue(*ptHH,giB_tdt_3D);
              cout << "\n sigma apres CP: ";
              ptHH->Ecriture(cout);
              delete ptHH;
            };
           #endif

           // passage des tenseurs résultats à l'ordre 2
           ((Tenseur2HH*) &sigHH)->Affectation_3D_a_2D(sig_HH_3D);
           break;
         };
       // sinon on ne fait pas de break, donc on continue avec la méthode de perturbation
     };

    default: // cas de la version d_sig/d_ddl
     {
      // initialisation du tenseurs contrainte
      sig_HH_3D.Inita(0.);
      sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH),plusZero);
      // récup de l'épaisseur de départ
      // on impose que les grandeurs sauvegardées soient dans les limites admises
      if (Limitation_h(sauve_hsurh0))
       {if (Permet_affichage() > 3)
          cout << "loi CP: calcul direct sans Newton " << flush;
       };
      save_resul.hsurh0 = sauve_hsurh0;

      // passage des métriques de l'ordre 2 vers 3 
      Passage_metrique_ordre2_vers_3(ex);
      // prise en compte de h  sur métriques
      Prise_en_compte_h_sur_metrique();
      // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
      Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_,d_epsBB,delta_epsBB
                                                   ,jacobien_0,jacobien,d_jacobien_tdt);
     
      lois_interne->Calcul_SigmaHH(sig_HH_t_3D,Deps_BB_3D,tab_ddl,gijBB_t_3D,gijHH_t_3D,giB_tdt_3D,giH_tdt_3D
                 ,eps_BB_3D,delta_eps_BB_3D
                 ,gijBB_tdt_3D,gijHH_tdt_3D,d_gijBB_tdt_3D_P,jacobien_0_3D,jacobien_tdt_3D,sig_HH_3D
                 ,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite,module_cisaillement,*expli_3D
                 );

      // passage des tenseurs résultats à l'ordre 2
      ((Tenseur2HH*) &sigHH)->Affectation_3D_a_2D(sig_HH_3D);
      if (Permet_affichage() > 0)
         { cout << "\n loi CP: valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0 << flush;
         };
     }
    break;
   };
 
 // sig_HH_3D.Ecriture(cout);
			       
	(*save_resul.l_sigoHH) = sig_HH_3D; // sauvegarde en locale
	// ---calcul des invariants de déformation et de vitesse de déformation
 Calcul_invariants_et_def_cumul();
  
 // dans le cas d'une contrainte par perturbation
//  if ((type_de_contrainte == PERTURBATION )||mauvaise_convergence) // non, la def d'épaisseur ne devrait pas changer !
// dans le cas d'une mauvaise convergence comme c'est traité plus haut, donc cela ne sert à rien de la recalculer
  if (type_de_contrainte == PERTURBATION )
    // calcul de la déformation d'épaisseur correspondant à la condition de contraintes planes
    {Tenseur3BH sigBH = gijBB_tdt_3D * sig_HH_3D;
     Calcul_eps33_parVarVolume(jacobien_0,module_compressibilite,jacobien,sigBH);
     if  ((Permet_affichage() > 0) && mauvaise_convergence)
         { cout << "\n loi CP: valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0 << flush;
         };
    };
 if (Permet_affichage() > 3)
  {cout << "\n -- loi CP: resultats ";
   cout << "\n sigHH_tdt= "; sigHH.Ecriture(cout);
   Tenseur3HH titi;
   sigHH.BaseAbsolue(titi,*(ex.giB_tdt));
   cout << "\n sig_tdt en absolu :";titi.Ecriture(cout);
   cout << "\n energ= " << energ;
   cout << "\n module_compressibilite= " << module_compressibilite
        << ", module_cisaillement= " << module_cisaillement;
  };

	energ = save_resul.l_energ; // récup des énergies
	// passage des tenseurs résultats à l'ordre 2
	((Tenseur2HH*) &sigHH)->Affectation_3D_a_2D(sig_HH_3D);
 //	sigHH_tdt.Ecriture(cout);
	LibereTenseur();
 LibereTenseurQ();
 if (Permet_affichage() > 5)
    cout << "\n fin: LoiContraintesPlanes::Calcul_SigmaHH(... " << flush;

 };
        
        // calcul des contraintes a t+dt et de ses variations  
void LoiContraintesPlanes::Calcul_DsigmaHH_tdt (TenseurHH& sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl
            ,BaseB& giB_t,TenseurBB & gijBB_t,TenseurHH & gijHH_t
            ,BaseB& giB_tdt,Tableau <BaseB> & d_giB_tdt,BaseH& giH_tdt,Tableau <BaseH> & d_giH_tdt
            ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB
            ,TenseurBB & delta_epsBB,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt
		          ,Tableau <TenseurBB *>& d_gijBB_tdt
		          ,Tableau <TenseurHH *>& d_gijHH_tdt,double& jacobien_0,double& jacobien
		          ,Vecteur& d_jacobien_tdt,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
		          ,EnergieMeca & energ,const EnergieMeca &
		          ,double& module_compressibilite,double&  module_cisaillement
		          ,const Met_abstraite::Impli& ex)
{ // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
  module_compressibilite=module_cisaillement=0.; // init
  energ.Inita(0.); // initialisation des énergies mises en jeux
  
  bool affichage = (Permet_affichage() > 3);
  if (affichage)
   {cout << "\n  --- LoiContraintesPlanes::Calcul_DsigmaHH_tdt --- ";
    Signature_pti_encours(cout);
   };
  if (Permet_affichage() > 3)
   {cout << "\n ===  donnees d'entree ";
    cout << "\n epsBB_tdt(local)= "; epsBB_tdt.Ecriture(cout);
    cout << "\n DepsBB(local)= "; DepsBB.Ecriture(cout);
    cout << "\n delta_epsBB(local)= "; delta_epsBB.Ecriture(cout);
    cout << "\n sigHH_t(local)= "; sigHH_t.Ecriture(cout);
    // en absolue
    Tenseur3BB tiutiu;
    epsBB_tdt.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n eps en absolu :";tiutiu.Ecriture(cout);
    DepsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n Deps en absolu :";tiutiu.Ecriture(cout);
    delta_epsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n delta_eps en absolu :";tiutiu.Ecriture(cout);
    Tenseur3HH titi;
    sigHH_t.BaseAbsolue(titi,*(ex.giB_tdt));
    cout << "\n sig_t en absolu :";titi.Ecriture(cout);
   };

  // on vérifie que le tableau de travail intermédiaire statique est correctement
  // dimensionné sinon on le modifie   
  int taille = d_sigHH.Taille();

  // on affecte et/ou on redimensionne éventuellement les tableaux contraintes-déformation fonction du nombre de ddl 
  // pour le passage 3D: on considère que tous les tableaux doivent avoir la même dimension: la même que dans le cas 2D
  int ta_d_sig_HH_3D = d_sig_HH_3D.Taille();
  if (ta_d_sig_HH_3D != taille)
    { // cela veut dire que tous les tableaux sont mal dimensionnés
      ta_d_sig_HH_3D = d_sigHH.Taille(); // nouvelle taille
      d_sig_HH_3D.Change_taille(ta_d_sig_HH_3D); d_sig_HH_3D_P.Change_taille(ta_d_sig_HH_3D);
      for (int i=1;i<=ta_d_sig_HH_3D;i++) {d_sig_HH_3D_P(i) = &(d_sig_HH_3D(i));}
    };
//  for (int i=1;i<=ta_d_sig_HH_3D;i++) {d_sig_HH_3D_P(i)->Inita(0.);};
        
  // initialisation du tenseurs contrainte 
//  sig_HH_3D.Inita(0.);
  // --- pour les contraintes passage en 3D
  bool plusZero = true;
  sig_HH_t_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_t),plusZero);
//  sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_tdt),plusZero);
  // --- pour la cinématique
  // on sauvegarde la nouvelle épaisseur, car elle va être modifiée et si on doit
  // repartir au départ, ce ne sera plus possible sans cette sauvegarde
  double sauve_hsurh0=save_resul.hsurh0;
  // passage des informations spécifique à la loi le_SaveResul
  lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
  lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
  lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours

 
  // choix entre calcul avec une boucle locale de Newton on non. Dans le premier cas il nous faut la version d_sig/d_eps
  // au lieu de d_sig/d_ddl
  bool mauvaise_convergence = false; // init
  switch (type_de_contrainte)
   {case NEWTON_LOCAL: // cas de la version d_sig/d_eps
     {
       // initialisation du tenseurs contrainte
       sig_HH_3D.Inita(0.);
       sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_tdt),plusZero);
       // passage des métriques de l'ordre 2 vers 3
       Passage_metrique_ordre2_vers_3(ex);
       // prise en compte de h et variation sur métriques
       Prise_en_compte_h_et_variation_sur_metrique(ex);
       // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
       Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_tdt,&d_epsBB,delta_epsBB
                                             ,jacobien_0,jacobien,&d_jacobien_tdt);
       for (int i=1;i<=ta_d_sig_HH_3D;i++) {d_sig_HH_3D_P(i)->Inita(0.);};
       // on appel la procédure de résolution de sig33(hsurh0)=0
       if ((sortie_post)&&(save_resul.indicateurs_resolution.Taille()!= 2)) // dimensionnement éventuelle de la sortie d'indicateurs
          save_resul.indicateurs_resolution.Change_taille(2);
       // ----- pour ce faire on appelle une méthode de recherche de zero
       val_initiale(1)=save_resul.h_tsurh0; // on démarre la recherche à la valeur à t
       // on impose que les grandeurs soient dans les limites admises
       if (Limitation_h(val_initiale(1)))
         {if (Permet_affichage() > 3)
            cout << "\n  loi CP: (newton local) val initiales: " << flush;
         };
       racine.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
       der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
       try // on met le bloc sous surveillance
        {
         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs, pas utilisées ici

         // dans le cas où on utilise une précision qui est pilotée
         { // opération de transmission de la métrique: encapsulé ici
           const Met_abstraite::Impli* ex_impli = &ex;
           const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
           const Met_abstraite::Umat_cont* ex_expli = NULL;

           if (fct_tolerance_residu != NULL)
             {// ici on utilise les variables connues aux pti, ou calculées à partir de
              // on commence par récupérer les conteneurs des grandeurs à fournir
              List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu->Li_enu_etendu_scalaire();
              List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu->Li_equi_Quel_evolue();
              bool absolue = true; // on se place systématiquement en absolu
              // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
              // pour les grandeurs strictement scalaire
               Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
                           (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
                                             );
              // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
              // pour les Coordonnees et Tenseur
              Valeurs_Tensorielles_interpoler_ou_calculer
                            (absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
              // calcul de la valeur et retour dans tab_ret
              Tableau <double> & tab_val = fct_tolerance_residu->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
              #ifdef MISE_AU_POINT
              if (tab_val.Taille() != 1)
                 { cout << "\n loi CP: Erreur : la fonction nD de pilotage de la precision "
                        << " doit calculer un scalaire or le tableau de retour est de taille "
                        << tab_val.Taille() << " ce n'est pas normal !\n";
                   cout << " LoiContraintesPlanesDouble::Calcul_DsigmaHH_tdt\n";
                   Sortie(1);
                 };
              #endif
              // on récupère le premier élément du tableau uniquement
              double tol = tab_val(1);
              #ifdef MISE_AU_POINT
              if (Permet_affichage() > 5)
                 cout << "\n loi CP: Newton_prec: tol_= "<< tol;
              #endif
              alg_zero.Modif_prec_res_abs(tol);
            };

           if (fct_tolerance_residu_rel != NULL)
             {// ici on utilise les variables connues aux pti, ou calculées à partir de
              // on commence par récupérer les conteneurs des grandeurs à fournir
              List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu_rel->Li_enu_etendu_scalaire();
              List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu_rel->Li_equi_Quel_evolue();
              bool absolue = true; // on se place systématiquement en absolu
              // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
              // pour les grandeurs strictement scalaire
               Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
                           (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
                                             );
              // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
              // pour les Coordonnees et Tenseur
              Valeurs_Tensorielles_interpoler_ou_calculer
                            (absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
              // calcul de la valeur et retour dans tab_ret
              Tableau <double> & tab_val = fct_tolerance_residu_rel->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
              #ifdef MISE_AU_POINT
              if (tab_val.Taille() != 1)
                 { cout << "\n loi CP: Erreur : la fonction nD de pilotage de la precision relative "
                        << " doit calculer un scalaire or le tableau de retour est de taille "
                        << tab_val.Taille() << " ce n'est pas normal !\n";
                   cout << " LoiContraintesPlanesDouble::Calcul_DsigmaHH_tdt\n";
                   Sortie(1);
                 };
              #endif
              // on récupère le premier élément du tableau uniquement
              double tol_rel = tab_val(1);
              #ifdef MISE_AU_POINT
              if (Permet_affichage() > 5)
                 cout << ", tol_relative= "<< tol_rel;
              #endif
              alg_zero.Modif_prec_res_rel(tol_rel);
            };
         };
       
         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
         bool conver=alg_zero.Newton_raphson
                  (*this,&LoiContraintesPlanes::Residu_constitutif,&LoiContraintesPlanes::Mat_tangente_constitutif
                   ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                   ,maxi_delta_var_eps_sur_iter_pour_Newton);
         if(sortie_post) // sauvegarde éventuelle des indicateurs
              {save_resul.indicateurs_resolution(1)+=nb_incr_total;
               save_resul.indicateurs_resolution(2)+=nb_iter_total;
              };
         // on vérifie qu'il n'y a pas de pb de convergence
         double absracinemax=racine.Max_val_abs();
         if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
           { if (Permet_affichage() > 0)
               { cout << "\n loi CP: non convergence sur l'algo de la resolution de sig33(hsurh0)=0 "
                      << " h/h_0(t+dt)= " << racine(1) << " a t on avait: " << save_resul.h_tsurh0
                      << "\n nb_incr_total=" << nb_incr_total
                      << " nb_iter_total=" << nb_iter_total;
               };
             // on garde en mémoire
             mauvaise_convergence=true;
           };
         // on vérifie que les nouvelles dimensions transversales ne sont pas négatives
         if (racine(1) < 0.)
             { if (Permet_affichage() > 0)
                 { cout << "\n loi CP: **** erreur dans la resolution : on obtient une "
                        << " epaisseur  negative "
                        << ", h/h_0(t+dt)= " << racine(1) << " a t on avait: " << save_resul.h_tsurh0;
                 };
               // on garde en mémoire
               mauvaise_convergence=true;
             };
        }
       catch (ErrNonConvergence_Newton erreur)
         {
          if (Permet_affichage() > 0)
               { cout << "\n loi CP: erreur de non convergence avec l'algo de la resolution de sig33(hsurh0)=0 "
                      << " on obtient une valeur infinie  ou NaN ";
               };
          // on garde en mémoire
          mauvaise_convergence=true;
         }
       catch (ErrSortieFinale)
            // cas d'une direction voulue vers la sortie
            // on relance l'interuption pour le niveau supérieur
          { ErrSortieFinale toto;
            throw (toto);
          }
       catch ( ... )
        { // dans le cas d'une erreur inconnue, on génère également une exception
          if (Permet_affichage() > 0)
               { cout << "\n loi CP: erreur non identifiee sur l'algo de la resolution de sig33(hsurh0)=0 ";
//// ----- debug
//cout << "\n debug :  LoiContraintesPlanes::Calcul_DSigmaHH ( ";
//       racine.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
//       der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
//         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs, pas utilisées ici
//         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
//         bool conver=alg_zero.Newton_raphson
//                  (*this,&LoiContraintesPlanes::Residu_constitutif,&LoiContraintesPlanes::Mat_tangente_constitutif
//                   ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total);
//// ---- fin debug
               };
          // on garde en mémoire
          mauvaise_convergence=true;
        };
       if (!mauvaise_convergence)
         { // arrivée ici, cela veut dire que tout à bien fonctionné
           save_resul.hsurh0 = racine(1); // récup de la solution
           // on met à jour les modules
           module_compressibilite= module_compressibilite_3D;
           module_cisaillement= module_cisaillement_3D;
          
           // calcul de la variation de la déformation d'épaisseur en fonction des déformations planes
           Tenseur2HH d_EPS33_EPS_2D_BB; // init à 0
           double inv_dsig33_eps33 = 0.;
           if (Dabs(d_sigma_deps_3D(3,3,3,3)) > ConstMath::trespetit)
            {inv_dsig33_eps33 = -1 / d_sigma_deps_3D(3,3,3,3);
                //- 1./ MaX(ConstMath::pasmalpetit,d_sigma_deps_3D(3,3,3,3));
             d_EPS33_EPS_2D_BB.Coor(1,1) = inv_dsig33_eps33 * d_sigma_deps_3D(3,3,1,1);
             d_EPS33_EPS_2D_BB.Coor(1,2) = inv_dsig33_eps33 * d_sigma_deps_3D(3,3,1,2);
             d_EPS33_EPS_2D_BB.Coor(2,2) = inv_dsig33_eps33 * d_sigma_deps_3D(3,3,2,2);
            };
////----- debug
//{cout << "\n debug :  LoiContraintesPlanes::Calcul_DSigmaHH ( ";
// cout << "\n inv_dsig33_eps33= "<<inv_dsig33_eps33
//      << " d_sigma_deps_3D(3,3,3,3)= " << d_sigma_deps_3D(3,3,3,3);
// cout << "\n d_sigma_deps_3D(3,3,1,1)= "<< d_sigma_deps_3D(3,3,1,1)
//      << " d_sigma_deps_3D(3,3,1,2)= "<< d_sigma_deps_3D(3,3,1,2)
//      << " d_sigma_deps_3D(3,3,2,2)= "<< d_sigma_deps_3D(3,3,2,2);
//
//}
////--- fin debug
          
           // récup de la variation de d sig^{alpha,beta}/ d eps_33
           Tenseur2HH d_sigAlphaBeta3D_deps33_HH(d_sigma_deps_3D(1,1,3,3),d_sigma_deps_3D(2,2,3,3),d_sigma_deps_3D(1,2,3,3));
           // construction de la variation de d sig^{alpha,beta} / d eps_{gamma eta) du uniquement à la variation de eps_33
          
           // passage des tenseurs résultats à l'ordre 2
           ((Tenseur2HH*) &sigHH_tdt)->Affectation_3D_a_2D(sig_HH_3D);
           // on calcul la matrice tangente en contrainte plane
           // on a : d_sigHH_tdt/d_ddl = d_sigHH_tdt/d_eps * d_eps/d_ddl
           Tenseur3HH dsigHH_inter; // une variable de travail
          
           if (Permet_affichage() > 7)
              { cout << "\n ContraintesPlanes: dans le repere locale d_sigma_deps_3D=  \n";
                int e=1;
                for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
                   { cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<d_sigma_deps_3D(i,j,k,l) << " ; ";
                     if (e>6) {cout << "\n"; e=1;}
                   };
                Tenseur3HHHH  inter_HHHH;
                d_sigma_deps_3D.ChangeBase(inter_HHHH,giB_tdt_3D);
                cout << "\n dans le repere orthonormee d_sigma_deps_3D=  \n";
                e=1;
                for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
                   { cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
                     if (e>6) {cout << "\n"; e=1;}
                   };
              };
            #ifdef MISE_AU_POINT
            if (Permet_affichage() > 7)
              {  cout << "\n d_EPS33_EPS_2D_BB";
                 d_EPS33_EPS_2D_BB.Ecriture(cout);
              };
              
            #endif

           for (int i = 1; i<= taille; i++)
            { // on fait uniquement une égalité d'adresse et de ne pas utiliser
              // le constructeur d'ou la profusion d'* et de ()
              Tenseur2HH & dsigHH = *((Tenseur2HH*) (d_sigHH(i)));  // passage en dim 2
              const Tenseur2BB & depsBB = *((Tenseur2BB *) (d_epsBB(i))); //    "
              dsigHH_inter = d_sigma_deps_3D && d_eps_BB_3D(i);
              dsigHH.Affectation_3D_a_2D(dsigHH_inter);
              // on s'occupe maintenant de la variation de d sig^{alpha,beta} / d eps_{gamma eta)
              // du uniquement à la variation de eps_33
              dsigHH += d_sigAlphaBeta3D_deps33_HH * (d_EPS33_EPS_2D_BB && depsBB);
            };
           #ifdef MISE_AU_POINT
           if (Permet_affichage() > 4)
              { for (int i = 1; i<= taille; i++)
                 {  cout << "\n et au final: d_sigHH("<<i<<") ";
                    d_sigHH(i)->Ecriture(cout);
                 };
              };
           #endif
           break;
         };
       // sinon on ne fait pas de break, donc on continue avec la méthode de perturbation
     };

    default: // cas de la version d_sig/d_ddl, on encore du cas où la méthode de Newton n'a pas convergé
     {
       // initialisation du tenseurs contrainte
       sig_HH_3D.Inita(0.);
       sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_tdt),plusZero);
       // récup de l'épaisseur de départ
       // on impose que les grandeurs sauvegardées soient dans les limites admises
       if (Limitation_h(sauve_hsurh0))
        {if (Permet_affichage() > 3)
           cout << "calcul direct sans Newton " << flush;
        };
       save_resul.hsurh0 = sauve_hsurh0;
       // passage des métriques de l'ordre 2 vers 3
       Passage_metrique_ordre2_vers_3(ex);
       // prise en compte de h et variation sur métriques
       Prise_en_compte_h_et_variation_sur_metrique(ex);
       // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
       Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_tdt,&d_epsBB,delta_epsBB
                                             ,jacobien_0,jacobien,&d_jacobien_tdt);
       for (int i=1;i<=ta_d_sig_HH_3D;i++) {d_sig_HH_3D_P(i)->Inita(0.);};
     
       lois_interne->Calcul_DsigmaHH_tdt(sig_HH_t_3D,Deps_BB_3D,tab_ddl,giB_t_3D,gijBB_t_3D,gijHH_t_3D
                                      ,giB_tdt_3D,d_giB_tdt_3D,giH_tdt_3D,d_giH_tdt_3D
                                      ,eps_BB_3D,d_eps_BB_3D_P,delta_eps_BB_3D,gijBB_tdt_3D,gijHH_tdt_3D,d_gijBB_tdt_3D_P
                                      ,d_gijHH_tdt_3D_P,jacobien_0_3D,jacobien_tdt_3D,d_jacobien_tdt_3D
                                      ,sig_HH_3D,d_sig_HH_3D_P
                                      ,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite,module_cisaillement,*impli_3D
                                      );

      // passage des tenseurs résultats à l'ordre 2
      ((Tenseur2HH*) &sigHH_tdt)->Affectation_3D_a_2D(sig_HH_3D);
     
      if (mauvaise_convergence)
       {if (Permet_affichage() > 0)
         { cout << "\n loi CP: valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0;
         };
        // on annulle les dérivées des épaisseurs
        save_resul.d_hsurh0.Zero();
       };

      // récup de l'opérateur tangent
      for (int k=1;k<=taille;k++)
            ((Tenseur2HH*) d_sigHH(k))->Affectation_3D_a_2D(*d_sig_HH_3D_P(k));
     }
    break;
   };
  // informations éventuelles
  #ifdef MISE_AU_POINT
  if  (Permet_affichage() > 7)
    { TenseurHH* ptHH = NevezTenseurHH(sig_HH_3D);
      sig_HH_3D.BaseAbsolue(*ptHH,giB_tdt_3D);
      cout << "\n sigma apres CP: ";
      ptHH->Ecriture(cout);
      delete ptHH;
    };
  #endif
  // sig_HH_3D.Ecriture(cout);
           
  (*save_resul.l_sigoHH) = sig_HH_3D; // sauvegarde en locale
  // ---calcul des invariants de déformation et de vitesse de déformation
  Calcul_invariants_et_def_cumul();
  
  energ = save_resul.l_energ; // récup des énergies
  
  // --- on remet à jour éventuellement l'épaisseur
  // dans le cas d'une contrainte par perturbation
//  if ((type_de_contrainte == PERTURBATION )||mauvaise_convergence)
  if (type_de_contrainte == PERTURBATION )
     // calcul de la déformation d'épaisseur correspondant à la condition de contraintes planes
     {Tenseur2BH sigBH = gijBB_tdt * sigHH_tdt;
      Calcul_d_eps33_parVarVolume(jacobien_0,module_compressibilite,jacobien,sigHH_tdt,d_jacobien_tdt
                     ,d_sigHH,d_gijBB_tdt,gijBB_tdt);
      if  ((Permet_affichage() > 0) && mauvaise_convergence)
          { cout << "\n loi CP: valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0;
          };
     };
  if (Permet_affichage() > 3)
   {cout << "\n --  resultats ";
    cout << "\n sigHH_tdt= "; sigHH_tdt.Ecriture(cout);
    Tenseur3HH titi;
    sigHH_tdt.BaseAbsolue(titi,*(ex.giB_tdt));
    cout << "\n sig_tdt en absolu :";titi.Ecriture(cout);
    cout << "\n energ= " << energ;
    cout << "\n module_compressibilite= " << module_compressibilite
         << ", module_cisaillement= " << module_cisaillement;
   };

  LibereTenseur();
  LibereTenseurQ();
  if (Permet_affichage() > 5)
    cout << "\n fin: LoiContraintesPlanes::Calcul_DsigmaHH_tdt(... " << flush;
};
		  	
        // calcul des contraintes et ses variations  par rapport aux déformations a t+dt
        // en_base_orthonormee:
        //    true : le tenseur de contrainte en entrée est  en orthonormee
        //                  le tenseur de déformation et son incrémentsont également en orthonormees
        //    si = false: les bases transmises sont utilisées
        // ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void LoiContraintesPlanes::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB& DepsBB
            ,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
            ,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps
            ,EnergieMeca & energ,const EnergieMeca & 
            ,double& module_compressibilite,double&  module_cisaillement
            ,const Met_abstraite::Umat_cont& ex)
{ // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
  module_compressibilite=module_cisaillement=0.; // init
  energ.Inita(0.); // initialisation des énergies mises en jeux
  
  bool affichage = (Permet_affichage() > 3);
  if (affichage)
   {cout << "\n  --- LoiContraintesPlanes::Calcul_dsigma_deps --- ";
    Signature_pti_encours(cout);
   };

  if (Permet_affichage() > 3)
   {cout << "\n ===  donnees d'entree ";
    cout << "\n epsBB_tdt(local)= "; epsBB_tdt.Ecriture(cout);
    cout << "\n DepsBB(local)= "; DepsBB.Ecriture(cout);
    cout << "\n delta_epsBB(local)= "; delta_epsBB.Ecriture(cout);
    cout << "\n sigHH_t(local)= "; sigHH_t.Ecriture(cout);
    // en absolue
    Tenseur3BB tiutiu;
    epsBB_tdt.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n eps en absolu :";tiutiu.Ecriture(cout);
    DepsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n Deps en absolu :";tiutiu.Ecriture(cout);
    delta_epsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
    cout << "\n delta_eps en absolu :";tiutiu.Ecriture(cout);
    Tenseur3HH titi;
    sigHH_t.BaseAbsolue(titi,*(ex.giB_tdt));
    cout << "\n sig_t en absolu :";titi.Ecriture(cout);
   };

  // initialisation du tenseurs contrainte 
//  sig_HH_3D.Inita(0.);
  // --- pour les contraintes passage en 3D
  bool plusZero = true;
  sig_HH_t_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_t),plusZero);
//  sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_tdt),plusZero);
  // --- pour la cinématique
  // on sauvegarde la nouvelle épaisseur, car elle va être modifiée et si on doit
  // repartir au départ, ce ne sera plus possible sans cette sauvegarde
  double sauve_hsurh0=save_resul.hsurh0;
  // passage des informations spécifique à la loi le_SaveResul
  lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
  lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
  lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours

  // choix entre calcul avec une boucle locale de Newton on non. Dans le premier cas il nous faut la version d_sig/d_eps
  // au lieu de d_sig/d_ddl
  bool mauvaise_convergence = false; // init
  switch (type_de_contrainte)
   {case NEWTON_LOCAL: // cas de la version d_sig/d_eps
     {
       // initialisation du tenseurs contrainte
       sig_HH_3D.Inita(0.);
       sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_tdt),plusZero);
       // passage des métriques de l'ordre 2 vers 3
       Passage_metrique_ordre2_vers_3(ex);
       // prise en compte de h sur métriques
       Prise_en_compte_h_sur_metrique();
       // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
       Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_tdt,delta_epsBB,jacobien_0,jacobien);
       // on appel la procédure de résolution de sig33(hsurh0)=0
       if ((sortie_post)&&(save_resul.indicateurs_resolution.Taille()!= 2)) // dimensionnement éventuelle de la sortie d'indicateurs
          save_resul.indicateurs_resolution.Change_taille(2);
       // ----- pour ce faire on appelle une méthode de recherche de zero
       val_initiale(1)=save_resul.h_tsurh0; // on démarre la recherche à la valeur à t
       // on impose que les grandeurs soient dans les limites admises
       if (Limitation_h(val_initiale(1)))
         {if (Permet_affichage() > 3)
            cout << "\n  (newton local) val initiales: " << flush;
         };
       racine.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
       der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
       try // on met le bloc sous surveillance
        {
         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs, pas utilisées ici

         // dans le cas où on utilise une précision qui est pilotée
         { // opération de transmission de la métrique: encapsulé ici
           const Met_abstraite::Impli* ex_impli = NULL;
           const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
           const Met_abstraite::Umat_cont* ex_expli = &ex;

           if (fct_tolerance_residu != NULL)
             {// ici on utilise les variables connues aux pti, ou calculées à partir de
              // on commence par récupérer les conteneurs des grandeurs à fournir
              List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu->Li_enu_etendu_scalaire();
              List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu->Li_equi_Quel_evolue();
              bool absolue = true; // on se place systématiquement en absolu
              // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
              // pour les grandeurs strictement scalaire
               Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
                           (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
                                             );
              // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
              // pour les Coordonnees et Tenseur
              Valeurs_Tensorielles_interpoler_ou_calculer
                            (absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
              // calcul de la valeur et retour dans tab_ret
              Tableau <double> & tab_val = fct_tolerance_residu->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
              #ifdef MISE_AU_POINT
              if (tab_val.Taille() != 1)
                 { cout << "\n loi CP: Erreur : la fonction nD de pilotage de la precision "
                        << " doit calculer un scalaire or le tableau de retour est de taille "
                        << tab_val.Taille() << " ce n'est pas normal !\n";
                   cout << " LoiContraintesPlanesDouble::Calcul_dsigma_deps\n";
                   Sortie(1);
                 };
              #endif
              // on récupère le premier élément du tableau uniquement
              double tol = tab_val(1);
              #ifdef MISE_AU_POINT
              if (Permet_affichage() > 5)
                 cout << "\n loi CP: Newton_prec: tol_= "<< tol;
              #endif
              alg_zero.Modif_prec_res_abs(tol);
            };

           if (fct_tolerance_residu_rel != NULL)
             {// ici on utilise les variables connues aux pti, ou calculées à partir de
              // on commence par récupérer les conteneurs des grandeurs à fournir
              List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu_rel->Li_enu_etendu_scalaire();
              List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu_rel->Li_equi_Quel_evolue();
              bool absolue = true; // on se place systématiquement en absolu
              // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
              // pour les grandeurs strictement scalaire
               Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
                           (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
                                             );
              // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
              // pour les Coordonnees et Tenseur
              Valeurs_Tensorielles_interpoler_ou_calculer
                            (absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
              // calcul de la valeur et retour dans tab_ret
              Tableau <double> & tab_val = fct_tolerance_residu_rel->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
              #ifdef MISE_AU_POINT
              if (tab_val.Taille() != 1)
                 { cout << "\n loi CP: Erreur : la fonction nD de pilotage de la precision relative "
                        << " doit calculer un scalaire or le tableau de retour est de taille "
                        << tab_val.Taille() << " ce n'est pas normal !\n";
                   cout << " LoiContraintesPlanesDouble::Calcul_dsigma_deps\n";
                   Sortie(1);
                 };
              #endif
              // on récupère le premier élément du tableau uniquement
              double tol_rel = tab_val(1);
              #ifdef MISE_AU_POINT
              if (Permet_affichage() > 5)
                 cout << ", tol_relative= "<< tol_rel;
              #endif
              alg_zero.Modif_prec_res_rel(tol_rel);
            };
         };
       
         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
         bool conver=alg_zero.Newton_raphson
                  (*this,&LoiContraintesPlanes::Residu_constitutif,&LoiContraintesPlanes::Mat_tangente_constitutif
                   ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                   ,maxi_delta_var_eps_sur_iter_pour_Newton);
         if(sortie_post) // sauvegarde éventuelle des indicateurs
              {save_resul.indicateurs_resolution(1)+=nb_incr_total;
               save_resul.indicateurs_resolution(2)+=nb_iter_total;
              };
         // on vérifie qu'il n'y a pas de pb de convergence
         double absracinemax=racine.Max_val_abs();
         if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
           { if (Permet_affichage() > 0)
               { cout << "\n loi CP: non convergence sur l'algo de la resolution de sig33(hsurh0)=0 "
                      << " h/h_0(t+dt)= " << racine(1) << " a t on avait: " << save_resul.h_tsurh0
                      << "\n nb_incr_total=" << nb_incr_total
                      << " nb_iter_total=" << nb_iter_total;
               };
             // on garde en mémoire
             mauvaise_convergence=true;
           };
         // on vérifie que les nouvelles dimensions transversales ne sont pas négatives
         if (racine(1) < 0.)
             { if (Permet_affichage() > 0)
                 { cout << "\n loi CP: **** erreur dans la resolution : on obtient une "
                        << " epaisseur  negative "
                        << ", h/h_0(t+dt)= " << racine(1) << " a t on avait: " << save_resul.h_tsurh0;
                 };
               // on garde en mémoire
               mauvaise_convergence=true;
             };
        }
       catch (ErrNonConvergence_Newton erreur)
         {
          if (Permet_affichage() > 0)
               { cout << "\n loi CP: erreur de non convergence avec l'algo de la resolution de sig33(hsurh0)=0 "
                      << " on obtient une valeur infinie  ou NaN ";
               };
          // on garde en mémoire
          mauvaise_convergence=true;
         }
       catch (ErrSortieFinale)
            // cas d'une direction voulue vers la sortie
            // on relance l'interuption pour le niveau supérieur
          { ErrSortieFinale toto;
            throw (toto);
          }
       catch ( ... )
        { // dans le cas d'une erreur inconnue, on génère également une exception
          if (Permet_affichage() > 0)
               { cout << "\n loi CP: erreur non identifiee sur l'algo de la resolution de sig33(hsurh0)=0 ";
//// ----- debug
//cout << "\n debug :  LoiContraintesPlanes::Calcul_DSigma_deps ( ";
//       racine.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
//       der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
//         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs, pas utilisées ici
//         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
//         bool conver=alg_zero.Newton_raphson
//                  (*this,&LoiContraintesPlanes::Residu_constitutif,&LoiContraintesPlanes::Mat_tangente_constitutif
//                   ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total);
//// ---- fin debug
               };
          // on garde en mémoire
          mauvaise_convergence=true;
        };
       if (!mauvaise_convergence)
         { // arrivée ici, cela veut dire que tout à bien fonctionné
           save_resul.hsurh0 = racine(1); // récup de la solution
           // on met à jour les modules
           module_compressibilite= module_compressibilite_3D;
           module_cisaillement= module_cisaillement_3D;
          
           // calcul de la variation de la déformation d'épaisseur en fonction des déformations planes
           Tenseur2HH d_EPS33_EPS_2D_BB; // init
           double inv_dsig33_eps33 = - 1./ MaX(ConstMath::pasmalpetit,d_sigma_deps_3D(3,3,3,3));
           d_EPS33_EPS_2D_BB.Coor(1,1) = inv_dsig33_eps33 * d_sigma_deps_3D(3,3,1,1);
           d_EPS33_EPS_2D_BB.Coor(1,2) = inv_dsig33_eps33 * d_sigma_deps_3D(3,3,1,2);
           d_EPS33_EPS_2D_BB.Coor(2,2) = inv_dsig33_eps33 * d_sigma_deps_3D(3,3,2,2);
          
           // récup de la variation de d sig^{alpha,beta}/ d eps_33
           Tenseur2HH d_sigAlphaBeta3D_deps33_HH(d_sigma_deps_3D(1,1,3,3),d_sigma_deps_3D(2,2,3,3),d_sigma_deps_3D(1,2,3,3));
           // passage des tenseurs résultats à l'ordre 2
           ((Tenseur2HH*) &sigHH_tdt)->Affectation_3D_a_2D(sig_HH_3D);
            #ifdef MISE_AU_POINT
            if  (Permet_affichage() > 7)
             { cout << "\n --- bonne convergence de Newton ";
               TenseurHH* ptHH = NevezTenseurHH(sig_HH_3D);
               sig_HH_3D.BaseAbsolue(*ptHH,giB_tdt_3D);
               cout << "\n sigma apres CP: ";
               ptHH->Ecriture(cout);
               delete ptHH;
             };
            if (Permet_affichage() > 8)
              { cout << "\n dans le repere locale d_sigma_deps_3D=  \n";
                int e=1;
                for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
                   { cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<d_sigma_deps_3D(i,j,k,l) << " ; ";
                     if (e>6) {cout << "\n"; e=1;}
                   };
                Tenseur3HHHH  inter_HHHH;
                d_sigma_deps_3D.ChangeBase(inter_HHHH,giB_tdt_3D);
                cout << "\n dans le repere orthonormee d_sigma_deps_3D=  \n";
                e=1;
                for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
                   { cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
                     if (e>6) {cout << "\n"; e=1;}
                   };
              };
           #endif
           // on calcul la matrice tangente en contrainte plane
           // on commence par transférer la partie 3D en la restriction 2D
           d_sigma_deps_2D.TransfertDunTenseurGeneral(d_sigma_deps_3D);
           // pour on ajoute la variation due à eps33
           d_sigma_deps_2D += Tenseur2HHHH::Prod_tensoriel(d_sigAlphaBeta3D_deps33_HH,d_EPS33_EPS_2D_BB);
           // maintenant on renseigne le tenseur de sortie
           bool pluszero = true;
           d_sigma_deps.Affectation_trans_dimension(d_sigma_deps_2D, pluszero);
           break;
         };
       // sinon on ne fait pas de break, donc on continue avec la méthode de perturbation
     };

    default: // cas de la version d_sig/d_ddl, on encore du cas où la méthode de Newton n'a pas convergé
     {bool en_base_orthonormee = false; // ici les tenseurs ne sont pas forcément en orthonormee
      // initialisation du tenseurs contrainte
      sig_HH_3D.Inita(0.);
      sig_HH_3D.Affectation_2D_a_3D(*((Tenseur2HH*) &sigHH_tdt),plusZero);
      // récup de l'épaisseur de départ
      // on impose que les grandeurs sauvegardées soient dans les limites admises
      if (Limitation_h(sauve_hsurh0))
       {if (Permet_affichage() > 3)
          cout << "calcul direct sans Newton " << flush;
       };
      save_resul.hsurh0 = sauve_hsurh0;
      // passage des métriques de l'ordre 2 vers 3
      Passage_metrique_ordre2_vers_3(ex);
      // prise en compte de h sur métriques
      Prise_en_compte_h_sur_metrique();
      // passage des informations liées à la déformation  de 2 vers 3, et variation de volume éventuelle
      Passage_deformation_volume_ordre2_vers_3(DepsBB,epsBB_tdt,delta_epsBB,jacobien_0,jacobien);
      lois_interne->Calcul_dsigma_deps(en_base_orthonormee,sig_HH_t_3D,Deps_BB_3D
                                      ,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
                                      ,sig_HH_3D,d_sigma_deps_3D
                                      ,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite,module_cisaillement
                                      ,*umat_cont_3D);

      // passage des tenseurs résultats à l'ordre 2
      ((Tenseur2HH*) &sigHH_tdt)->Affectation_3D_a_2D(sig_HH_3D);
     
      if (mauvaise_convergence)
       {if (Permet_affichage() > 0)
         { cout << "\n loi CP: valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0;
         };
        // on annulle les dérivées des épaisseurs
        save_resul.d_hsurh0.Zero();
       };

      // maintenant on renseigne le tenseur de sortie
      bool pluszero = true;
      d_sigma_deps.Affectation_trans_dimension(d_sigma_deps_3D, pluszero);
  //debug
  //cout << "\n debug LoiContraintesPlanes::Calcul_DsigmaHH_tdt sig33_3D= "<<sig_HH_3D(3,3)<<endl;*/
  // fin debug
     }
    break;
   };
 
  // sig_HH_3D.Ecriture(cout);
           
  (*save_resul.l_sigoHH) = sig_HH_3D; // sauvegarde en locale
  // ---calcul des invariants de déformation et de vitesse de déformation
  Calcul_invariants_et_def_cumul();
  
  energ = save_resul.l_energ; // récup des énergies
  
  // --- on remet à jour éventuellement l'épaisseur
  // dans le cas d'une contrainte par perturbation
//  if ((type_de_contrainte == PERTURBATION )||mauvaise_convergence)
  if (type_de_contrainte == PERTURBATION )
     // calcul de la déformation d'épaisseur correspondant à la condition de contraintes planes
     {Tenseur2BH sigBH = *(ex.gijBB_tdt) * sigHH_tdt;
      Calcul_eps33_parVarVolume(jacobien_0,module_compressibilite,jacobien,sigBH);
      if  ((Permet_affichage() > 0) && mauvaise_convergence)
          { cout << "\n loi CP: valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0;
          };
     };

  if (Permet_affichage() > 3)
   {cout << "\n --  resultats ";
    cout << "\n sigHH_tdt= "; sigHH_tdt.Ecriture(cout);
    Tenseur3HH titi;
    sigHH_tdt.BaseAbsolue(titi,*(ex.giB_tdt));
    cout << "\n sig_tdt en absolu :";titi.Ecriture(cout);
    if (Permet_affichage() > 4)
      cout << "\n d_sigma_deps= "; d_sigma_deps.Ecriture(cout);
    cout << "\n energ= " << energ;
    cout << "\n module_compressibilite= " << module_compressibilite
         << ", module_cisaillement= " << module_cisaillement;
   };

 // fin debug
  LibereTenseur();
  LibereTenseurQ();
 if (Permet_affichage() > 5)
   cout << "\n fin: LoiContraintesPlanes::Calcul_dsigma_deps(... " << flush;
 };


   
       // fonction interne utilisée par les classes dérivées de Loi_comp_abstraite
       // pour répercuter les modifications de la température
       // ici utiliser pour modifier la température des lois élémentaires
       // l'Enum_dure: indique quel est la température courante : 0 t ou tdt
 void LoiContraintesPlanes::RepercuteChangeTemperature(Enum_dure temps)
  { lois_interne->temperature_0 = this->temperature_0;
    lois_interne->temperature_t = this->temperature_t;
    lois_interne->temperature_tdt = this->temperature_tdt;
    lois_interne->dilatation=dilatation;

    // on répercute également les déformations thermiques, qui ne sont utilisées
    // telles quelles que pour certaines lois: ex: loi hyper-élastique
    if (dilatation)
     {//  a- dimensionnement des tenseurs intermédiaires 
      int dim_tens = epsBB_therm->Dimension();
          // -- cas de la déformation
      if (lois_interne->epsBB_therm == NULL) { lois_interne->epsBB_therm = NevezTenseurBB(dim_tens);}
      else if (lois_interne->epsBB_therm->Dimension() != dim_tens)
         { delete lois_interne->epsBB_therm;lois_interne->epsBB_therm = NevezTenseurBB(dim_tens);};
          // -- cas de la vitesse de déformation
      if (lois_interne->DepsBB_therm == NULL) { lois_interne->DepsBB_therm = NevezTenseurBB(dim_tens);}
      else if (lois_interne->DepsBB_therm->Dimension() != dim_tens)
       { delete lois_interne->DepsBB_therm;lois_interne->DepsBB_totale = NevezTenseurBB(dim_tens);};
      // b- affectation des tenseurs
      (*lois_interne->epsBB_therm)=(*epsBB_therm);
      (*lois_interne->DepsBB_therm)=(*DepsBB_therm);
     };
    // puis loi interne
    lois_interne->RepercuteChangeTemperature(temps);
    switch (temps)
      { case TEMPS_0: 
         {lois_interne->temperature = &lois_interne->temperature_0;
          break;
          }
        case TEMPS_t: 
         {lois_interne->temperature = &lois_interne->temperature_t;
          break;
          }
        case TEMPS_tdt: 
         {lois_interne->temperature = &lois_interne->temperature_tdt;
          break;
          }
        default:
             { cout << "\n erreur, cas de temps non prevu !! "
                    << "\n LoiContraintesPlanes::RepercuteChangeTemperature(...";
               Sortie(1);      
             };
      }; 

  };
	
void LoiContraintesPlanes::CalculGrandeurTravail
               (const PtIntegMecaInterne& ptintmec
               ,const Deformation & def,Enum_dure temps,const ThermoDonnee& dTP
               ,const Met_abstraite::Impli* ex_impli
               ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
               ,const Met_abstraite::Umat_cont* ex_umat
               ,const List_io<Ddl_etendu>* exclure_dd_etend
               ,const List_io<const TypeQuelconque *>* exclure_Q
               )
{ // récup des infos spécifiques  
  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
  // par défaut on active les invariants
 
  // cas du niveau mini de sig33
  bool existence_ponderation = false; // init par défaut
   {
//     double proportion = 1.; // initialisation
//     if (niveauF_grandeurGlobale != NULL)
//       { // il faut que l'on récupère les grandeurs globales
//         List_io <Enum_GrandeurGlobale>::iterator il,ilfin
//                = niveauF_grandeurGlobale->Type_grandeur_GGlob().end();
//         // on parcours la liste
//         int i = 1; // indice du tableau d'argument
//         for (il = niveauF_grandeurGlobale->Type_grandeur_GGlob().begin()
//              ; il != ilfin; il++,i++)
//          { // on récupère le pointeur correspondant à la grandeur
//            const void* pointe =  (ParaGlob::param->GrandeurGlobal(*il));
//            #ifdef MISE_AU_POINT
//            if (pointe == NULL)
//              { cout << "\n *** pb dans la loi critere !! "
//                     << " la variable globale  "<< Nom_GrandeurGlobale(*il)
//                     << ", n'est pas disponible, on ne peut pas continuer  "
//                     << "\n LoiContraintesPlanes::CalculGrandeurTravail(..."<<endl;
//                Sortie(1);
//              };
//            #endif
//            TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
//            switch(gr_quelc->Grandeur_pointee()->Type_structure_grandeurAssocie())
//             { case  TYPE_SIMPLE:
//                 { switch(gr_quelc->Grandeur_pointee()->Type_enumGrandeurParticuliere())
//                    {case PARTICULIER_SCALAIRE_ENTIER:
//                       {Grandeur_scalaire_entier& gr
//                             = *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
//                        niveauF_grandeurGlobale->Tab_argument()(i) = *(gr.ConteneurEntier());
//                        break;
//                       }
//                     case PARTICULIER_SCALAIRE_DOUBLE:
//                       {Grandeur_scalaire_double& gr
//                             = *((Grandeur_scalaire_double*) gr_quelc->Grandeur_pointee()); // pour simplifier
//                        niveauF_grandeurGlobale->Tab_argument()(i) = *(gr.ConteneurDouble());
//                        break;
//                       }
//                     default:
//                       { cout << "\n *** pb dans la loi des melanges !! "
//                              << " la variable globale  "<< Nom_GrandeurGlobale(*il)
//                              << ", n'est pas prise en compte actuellement , on ne peut pas continuer  "
//                              << "\n LoiContraintesPlanes::CalculGrandeurTravail(..."<<endl;
//                         Sortie(1);
//                       };
//                    }
//                 
//                  break;
//                 }
//               default:
//                 { cout << "\n *** pb dans la loi des melanges !! "
//                        << " la variable globale  "<< Nom_GrandeurGlobale(*il)
//                        << ", n'est pas prise en compte actuellement , on ne peut pas continuer  "
//                        << "\n LoiContraintesPlanes::CalculGrandeurTravail(..."<<endl;
//                   Sortie(1);
//                 };
//             };
//          };
//         // maintenant on appelle la fonction
////// ------ debug
//// cout << "\n debug LoiContraintesPlanes::CalculGrandeurTravail( ";
//// cout << "niveauF_grandeurGlobale->Tab_argument()= "<< niveauF_grandeurGlobale->Tab_argument()
////      << endl;
////// ------ fin debug
////          
//          
//         double fonc = niveauF_grandeurGlobale->C_proport()
//                       ->Valeur(niveauF_grandeurGlobale->Tab_argument())(1);
//         proportion *= fonc; // mise à jour de la proportion
//         existence_ponderation = true;
//       };
//     if (niveauF_temps != NULL) // cas d'une dépendance au temps
//        {proportion *= niveauF_temps->CalculPonder();
//         existence_ponderation = true;
//        };
//     // avec une dépendance via éventuellement des ddl étendu
//     if (niveauF_ddlEtendu != NULL)
//        {proportion *= niveauF_ddlEtendu->CalculPonderMultiplicatif(ptintmeca, def, temps, dTP);
//         existence_ponderation= true;
//        }
//     if (existence_ponderation) // si une des pondérations est actives
//          {save_resul.niveau_sig33 = proportion*niveau_mini_sig33;}
//     else {save_resul.niveau_sig33 = niveau_mini_sig33;};

   };
 
 
  
  // ptintmec concerne des tenseurs  d'ordres 2 donc il faut transformer en 3D 
  bool plusZero=true; // on ajoute des 0 pour les grandeurs manquantes dans une première étape
  // on commence par transférer le tout
  ptintmeca.Affectation_2D_a_3D(ptintmec,plusZero);
  // ensuite on traite les points particuliers
  // pour tout ce qui est contrainte: l'ajout de 0 est ok
 
  // pour les déformations on met à jour en fonction des grandeurs sauvegardées
  double eps33=0.; double eps33_t=0.; //init par défaut
 
  // -- maintenant on tient compte de l'élongation d'épaisseur
  switch (type_de_deformation)
	  {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
         // cas d'une déformation d'Almansi
         {  // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
            eps33 = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
            eps33_t = 0.5 * (save_resul.h_tsurh0 * save_resul.h_tsurh0 - 1.);
         };
         break;
		  case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
        // cas d'une def logarithmique ou une approximation
         {  eps33 = log(save_resul.hsurh0);
            eps33_t = log(save_resul.h_tsurh0);
         };
         break;
		  default :
			   cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
           << Nom_type_deformation(type_de_deformation);
			   cout << "\n LoiContraintesPlanes::CalculGrandeurTravail \n";
			   Sortie(1);
	  };
 
  (*ptintmeca.EpsBB()).Coor(3,3) = eps33;
  (*ptintmeca.DeltaEpsBB()).Coor(3,3) = eps33 - eps33_t;
  
  // recup de l'incrément de temps
  double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
  double unSurDeltat=0;    
  if (Abs(deltat) >= ConstMath::trespetit)
   {unSurDeltat = 1./deltat;}
  else
  // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
   { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
     if (unSurDeltat < 0)
      { cout << "\n le pas de temps est négatif !! "; };
     unSurDeltat = ConstMath::tresgrand;
   };
  (*ptintmeca.DepsBB()).Coor(3,3) = (eps33 - eps33_t) * unSurDeltat;
  // on met à jour les invariants 3D, on ne peut pas les recalculer, donc on récupère ce qui a été sauvegardé
  if (ptintmeca.Statut_Invariants_deformation()) // le conteneur des invariants a pu être effacé par Affectation_2D_a_3D
     ptintmeca.EpsInvar() =  save_resul.epsInvar;
  if (ptintmeca.Statut_Invariants_vitesseDeformation()) // le conteneur des invariants a pu être effacé par Affectation_2D_a_3D
  ptintmeca.DepsInvar() =  save_resul.depsInvar;
  // idem au niveau des déformations cumulées
  ptintmeca.Deformation_equi() = save_resul.def_equi;
  ptintmeca.Deformation_equi_t() = save_resul.def_equi_t;

  // passage des informations spécifique à la loi le_SaveResul
  lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
  lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
  lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours

  lois_interne->CalculGrandeurTravail(ptintmeca,def,temps,dTP
                             ,ex_impli,ex_expli_tdt,ex_umat,exclure_dd_etend,exclure_Q); // passage  à la loi 3D
 }; 

// passage des grandeurs métriques de l'ordre 2 à 3
// ramène un conteneur dont les éléments sont gérés par la loi CP, et ne peuvent être modifié que par elle
const Met_abstraite::Impli*  LoiContraintesPlanes::Passage_metrique_ordre2_vers_3(const Met_abstraite::Impli& ex)
{// on s'occupe du redimensionnement éventuel
	// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
	// la partie dépendant des vitesses: entre accolades pour pouvoir fermer 
	{if (ex.gradVmoyBB_t != NULL) {impli_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;}; 
	 if (ex.gradVmoyBB_tdt != NULL) {impli_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;}; 
	 if (ex.gradVBB_tdt != NULL) {impli_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
	
	 // on s'occupe tout d'abord des tableaux de vitesses qui peuvent ne pas exister
	 if (ex.d_gradVmoyBB_t != NULL) // cas où il existe
	  { int tail=ex.d_gradVmoyBB_t->Taille();
	    if (d_gradVmoyBB_t_3D_P == NULL) {d_gradVmoyBB_t_3D_P = new Tableau<TenseurBB *> (tail);
	                                      impli_3D->d_gradVmoyBB_t = d_gradVmoyBB_t_3D_P;}
	    else                             {d_gradVmoyBB_t_3D_P->Change_taille(tail);};
		   d_gradVmoyBB_t_3D.Change_taille(tail);
		   for (int i=1;i<= tail;i++) 
			   (*d_gradVmoyBB_t_3D_P)(i) = &(d_gradVmoyBB_t_3D(i)); 
	  };		 
	 if (ex.d_gradVmoyBB_tdt != NULL) // cas où il existe
	  { int tail=ex.d_gradVmoyBB_tdt->Taille();
	    if (d_gradVmoyBB_tdt_3D_P == NULL) {d_gradVmoyBB_tdt_3D_P = new Tableau<TenseurBB *> (tail);
	                                        impli_3D->d_gradVmoyBB_tdt = d_gradVmoyBB_tdt_3D_P;}
	    else									  {d_gradVmoyBB_tdt_3D_P->Change_taille(tail);};
		   d_gradVmoyBB_tdt_3D.Change_taille(tail);
		   for (int i=1;i<= tail;i++) 
			    (*d_gradVmoyBB_tdt_3D_P)(i) = &(d_gradVmoyBB_tdt_3D(i)); 
	  };		 
	 if (ex.d_gradVBB_t != NULL) // cas où il existe
	  { int tail=ex.d_gradVBB_t->Taille();
	    if (d_gradVBB_t_3D_P == NULL) {d_gradVBB_t_3D_P = new Tableau<TenseurBB *> (tail);
	                                   impli_3D->d_gradVBB_t = d_gradVBB_t_3D_P;}
	    else									  {d_gradVBB_t_3D_P->Change_taille(tail);};
		   d_gradVBB_t_3D.Change_taille(tail);
		   for (int i=1;i<= tail;i++) 
			   (*d_gradVBB_t_3D_P)(i) = &(d_gradVBB_t_3D(i)); 
	  };		 
	 if (ex.d_gradVBB_tdt != NULL) // cas où il existe
	  { int tail=ex.d_gradVBB_tdt->Taille();
	    if (d_gradVBB_tdt_3D_P == NULL) {d_gradVBB_tdt_3D_P = new Tableau<TenseurBB *> (tail);
	                                     impli_3D->d_gradVBB_tdt = d_gradVBB_tdt_3D_P;}
	    else									  {d_gradVBB_tdt_3D_P->Change_taille(tail);};
		   d_gradVBB_tdt_3D.Change_taille(tail);
		   for (int i=1;i<= tail;i++) 
			    (*d_gradVBB_tdt_3D_P)(i) = &(d_gradVBB_tdt_3D(i)); 
   };		 

	}; // fin de la partie dédiée à la vitesse
	// maintenant on s'occupe uniquement du redimensionnement des tableaux restants
	   // -- on s'occupe des tableaux nécessaire à la métrique
	int ta_ex_d_giB_tdt = ex.d_giB_tdt->Taille(); // la dimension
 if (d_giB_tdt_3D.Taille() != ta_ex_d_giB_tdt)
	   { 	d_giB_tdt_3D.Change_taille(ta_ex_d_giB_tdt);};
	int ta_ex_d_giH_tdt = ex.d_giH_tdt->Taille(); // la dimension
 if (d_giH_tdt_3D.Taille() != ta_ex_d_giH_tdt)
	   { 	d_giH_tdt_3D.Change_taille(ta_ex_d_giH_tdt);};
	int ta_ex_d_gijBB_tdt = ex.d_gijBB_tdt->Taille(); // la dimension
 if (d_gijBB_tdt_3D.Taille() != ta_ex_d_gijBB_tdt)
	   { d_gijBB_tdt_3D.Change_taille(ta_ex_d_gijBB_tdt); d_gijBB_tdt_3D_P.Change_taille(ta_ex_d_gijBB_tdt);
	     for (int i=1;i<=ta_ex_d_gijBB_tdt;i++) {d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));};
		  };
	int ta_ex_d_gijHH_tdt = ex.d_gijHH_tdt->Taille();
 if (d_gijHH_tdt_3D.Taille() != ta_ex_d_gijHH_tdt)
	   { d_gijHH_tdt_3D.Change_taille(ta_ex_d_gijHH_tdt); d_gijHH_tdt_3D_P.Change_taille(ta_ex_d_gijHH_tdt);
	     for (int i=1;i<=ta_ex_d_gijHH_tdt;i++) {d_gijHH_tdt_3D_P(i) = &(d_gijHH_tdt_3D(i));};
		  };
	if(	d_jacobien_tdt_3D.Taille() != ex.d_jacobien_tdt->Taille())
		{d_jacobien_tdt_3D.Change_taille(ex.d_jacobien_tdt->Taille());}		
																				
		       // le tableau à double dimension 				 
	if (ex.d2_gijBB_tdt != NULL) 	    
		 { int taille1=ex.d2_gijBB_tdt->Taille1(); int taille2=ex.d2_gijBB_tdt->Taille2();
			  if (d2_gijBB_tdt_3D_P == NULL) 
			   { d2_gijBB_tdt_3D_P = impli_3D->d2_gijBB_tdt =  new Tableau2<TenseurBB *> (taille1,taille2);
				    d2_gijBB_tdt_3D.Change_taille(taille1,taille2);
				    for (int i=1;i<=taille1;i++) for (int j=1;j<= taille2;j++)
					     { (*d2_gijBB_tdt_3D_P)(i,j) = &d2_gijBB_tdt_3D(i,j);};
			   }
			  else if ((taille1 != d2_gijBB_tdt_3D.Taille1()) || (taille2 != d2_gijBB_tdt_3D.Taille2()))
			   { d2_gijBB_tdt_3D_P->Change_taille(taille1,taille2); 
				    d2_gijBB_tdt_3D.Change_taille(taille1,taille2);
				    for (int i=1;i<=taille1;i++) for (int j=1;j<= taille2;j++)
					     { (*d2_gijBB_tdt_3D_P)(i,j) = &d2_gijBB_tdt_3D(i,j);};
			   };
		 };
 // on commence par recopier les grandeurs de l'ordre 2 à 3
 bool plusZero = true; // on complète avec des 0 dans un premier temps
 int type_recopie=0; //   = 0 -> on transfert les grandeurs à 0, t et tdt
	impli_3D->Passage_de_Ordre2_vers_Ordre3(ex,plusZero,type_recopie);
		 	
	// maintenant on s'occupe de mettre à jour les grandeurs manquantes
	// - les bases naturelles: le vecteur normal est normé et est identique pour les bases naturelles et duales
	giB_0_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_0_3D(1),giB_0_3D(2))).Normer();
	giH_0_3D.CoordoH(3) = giB_0_3D(3).Bas_haut();
	giB_t_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_t_3D(1),giB_t_3D(2))).Normer();
	giH_t_3D.CoordoH(3) = giB_t_3D(3).Bas_haut();
	   // cas particulier du vecteur tdt
	giB_tdt_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_tdt_3D(1),giB_tdt_3D(2))); // calcul du vecteur normal non normé
	double norme_N_tdt = giB_tdt_3D(3).Norme(); // calcul de la norme qui nous servira pour les variations
 giB_tdt_3D.CoordoB(3) /= norme_N_tdt;
 giB_normer_3_tdt_3D_sauve = giB_tdt_3D(3); // on sauvegarde la valeur du vecteur giB(3) normée
	giH_tdt_3D.CoordoH(3) = giB_tdt_3D(3).Bas_haut();
   //  - les tenseurs métriques: au début 1 pour la direction 3
	gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)
        =gijHH_t_3D.Coor(3,3)=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
	//  - les variations du vecteur normal	
	CoordonneeB d_NN; // un vecteur de travail  
	#ifdef MISE_AU_POINT
	if (ex.d_giB_tdt->Taille() != ex.d_giH_tdt->Taille())
	   { cout << "\n **** sans doute une erreur: (ex.d_giB_tdt->Taille() "<< ex.d_giB_tdt->Taille()
		         << " != ex.d_giH_tdt->Taille()) " << ex.d_giH_tdt->Taille()
				       << "\n LoiDeformationsPlanes::Passage_metrique_ordre2_vers_3(const Met_abstraite::Impli& ex)";
		    Sortie(1);
	   };
 #endif
	for (int iddl=1;iddl<= ta_ex_d_giB_tdt;iddl++)
      { // tout d'abord on calcul la variation de la normale non normée
        d_NN =   (Util::ProdVec_coorB(d_giB_tdt_3D(iddl)(1),giB_tdt_3D(2)))
               + (Util::ProdVec_coorB(giB_tdt_3D(1),d_giB_tdt_3D(iddl)(2)));
        // maintenant on calcul la variation du vecteur normal unitaire     
        d_giB_tdt_3D(iddl).CoordoB(3) =  Util::VarUnVect_coorB(giB_tdt_3D(3),d_NN,norme_N_tdt);
      };
 // retour de la nouvelle métrique
 return impli_3D;

};

// passage des grandeurs métriques de l'ordre 2 à 3: cas explicite
// ramène un conteneur dont les éléments sont gérés par la loi CP, et ne peuvent être modifié que par elle
const Met_abstraite::Expli_t_tdt* LoiContraintesPlanes::Passage_metrique_ordre2_vers_3(const Met_abstraite::Expli_t_tdt& ex)
{ // on s'occupe du redimensionnement éventuel
	 // on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
	 // la partie dépendant des vitesses: entre accolades pour pouvoir fermer
	 {if (ex.gradVmoyBB_t != NULL) {expli_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;}; 
	  if (ex.gradVmoyBB_tdt != NULL) {expli_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;}; 
	  if (ex.gradVBB_tdt != NULL) {expli_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
	 }; // fin de la partie dédiée à la vitesse
	 // la dimension des tableaux : on considère que tous les tableaux doivent avoir la même dimension
	 int ta_ex_d_gijBB_tdt = ex.d_gijBB_tdt->Taille(); // la dimension
	 // maintenant on s'occupe uniquement du redimensionnement des tableaux restants
	 int ta_d_gijBB_tdt_3D = d_gijBB_tdt_3D.Taille();
	 if (ta_d_gijBB_tdt_3D != ta_ex_d_gijBB_tdt)
	  {// cela veut dire que tous les tableaux sont mal dimensionnés
	   // -- on s'occupe des tableaux nécessaire à la métrique
	   d_gijBB_tdt_3D.Change_taille(ta_ex_d_gijBB_tdt); d_gijBB_tdt_3D_P.Change_taille(ta_ex_d_gijBB_tdt);
	   for (int i=1;i<=ta_ex_d_gijBB_tdt;i++) {d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));};
	  };	 
  // on commence par recopier les grandeurs de l'ordre 2 à 3
  bool plusZero = true; // on complète avec des 0 dans un premier temps
  int type_recopie=0; //   = 0 -> on transfert les grandeurs à 0, t et tdt
  expli_3D->Passage_de_Ordre2_vers_Ordre3(ex,plusZero,type_recopie);
	
	 // maintenant on s'occupe de mettre à jour les grandeurs manquantes
	 // - les bases naturelles: le vecteur normal est normé et est identique pour les bases naturelles et duales
	 giB_0_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_0_3D(1),giB_0_3D(2))).Normer();
	 giH_0_3D.CoordoH(3) = giB_0_3D(3).Bas_haut();
	 giB_t_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_t_3D(1),giB_t_3D(2))).Normer();
	 giH_t_3D.CoordoH(3) = giB_t_3D(3).Bas_haut();
	   // cas particulier du vecteur tdt
	 giB_tdt_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_tdt_3D(1),giB_tdt_3D(2))); // calcul du vecteur normal non normé
	 double norme_N_tdt = giB_tdt_3D(3).Norme(); // calcul de la norme qui nous servira pour les variations
  giB_tdt_3D.CoordoB(3) /= norme_N_tdt;
	 giH_tdt_3D.CoordoH(3) = giB_tdt_3D(3).Bas_haut();
  giB_normer_3_tdt_3D_sauve = giB_tdt_3D(3); // on sauvegarde la valeur du vecteur giB(3) normée
    //  - les tenseurs métriques: a priori 1 pour la direction 3, initialement 
//	 gijBB_0_3D(3,3)=gijHH_0_3D(3,3)=1.;
	 gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
                      =gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
  // retour de la nouvelle métrique
  return expli_3D;

};

// prise en compte de h  sur métriques
void LoiContraintesPlanes::Prise_en_compte_h_sur_metrique()
 {
  //  -- prise en compte de la variation de la déformation d'épaisseur : on considère une déformation de type logarithmique
  // eps33 = log (h/h0) -> h = h0 * exp(eps33) , pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = 1
  // donc à t et tdt il faut multiplier par exp(eps33) à t et tdt
     // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
  // lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
  // ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
  // ensuite elles évolueront
  bool premier_passage = false;
  if (save_resul.h_tsurh0 == 0.)
    { save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.; premier_passage = true;};
  if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
    { save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise

  // on calcul le vecteur 3 naturel et dual
  giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
  giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
  //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
  //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
  double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
  gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;

  double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
  gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;

 };
// prise en compte de h et variation  sur métriques
void LoiContraintesPlanes::Prise_en_compte_h_et_variation_sur_metrique(const Met_abstraite::Impli& ex)
{ //  -- prise en compte de la variation de la déformation d'épaisseur : on considère une déformation de type logarithmique
  // eps33 = log (h/h0) -> h = h0 * exp(eps33) , pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = 1
  // donc à t et tdt il faut multiplier par exp(eps33) à t et tdt
     // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
  // lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
  // ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
  // ensuite elles évolueront
  int ta_ex_d_giB_tdt = ex.d_giB_tdt->Taille(); // la dimension
  if (save_resul.d_hsurh0.Taille() != ta_ex_d_giB_tdt)
     save_resul.d_hsurh0.Change_taille(ta_ex_d_giB_tdt); // donc à 0 la première fois
  bool premier_passage = false;
  if (save_resul.h_tsurh0 == 0.)
    { save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.; premier_passage = true;
    };
  if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
    { save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
  // on calcul le vecteur 3 naturel et dual
  giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
  giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
  //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
  //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
  double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
  gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;

  double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
  gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
 
  // maintenant on s'occupe des variations
  if (!premier_passage)
   {Vecteur& d_hsurh0 = save_resul.d_hsurh0;
    for (int iddl=1;iddl<= ta_ex_d_giB_tdt;iddl++)
       { double d_hsurh0_iddl = d_hsurh0(iddl);
         // gib(3) = h/h0  * N -> d_gib(3) = h/h0  * N = d_(h/h0)  * N + h/h0 * d_N
         // actuellement, c'est la partie d_N qui est stockée dans d_gib(3), on rajoute donc la partie restante
         d_giB_tdt_3D(iddl).CoordoB(3) *= save_resul.hsurh0;
         d_giB_tdt_3D(iddl).CoordoB(3) += d_hsurh0_iddl * giB_normer_3_tdt_3D_sauve; // base naturelle
         // gih(3) = 1/giB(3) -> d_gih(3)(3) = - d_giB(3)(3) / (giB(3)*giB(3))
         d_giH_tdt_3D(iddl).CoordoH(3) = (d_giB_tdt_3D(iddl)(3) / (-hsurh0_carre)).Bas_haut() ; // base duale
         // gijBB(3)= giB(3) * giB(3) -> d_gijBB(3)= 2. * d_giB(3) * giB(3)
         d_gijBB_tdt_3D(iddl).Coor(3,3) = 2. * (d_giB_tdt_3D(iddl)(3).ScalBB(giB_tdt_3D(3)));
         // gijHH(3)= giH(3) * giH(3) -> d_gijHH(3)= 2. * d_giH(3) * giH(3)
         d_gijHH_tdt_3D (iddl).Coor(3,3) = 2. * (d_giH_tdt_3D(iddl)(3).ScalHH(giH_tdt_3D(3)));
       };
   };

};

// passage des grandeurs métriques de l'ordre 2 à 3
// ramène un conteneur dont les éléments sont gérés par la loi CP, et ne peuvent être modifié que par elle
const Met_abstraite::Umat_cont* LoiContraintesPlanes::Passage_metrique_ordre2_vers_3(const Met_abstraite::Umat_cont& ex)
{// on s'occupe du redimensionnement éventuel
	// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
	{if (ex.gradVmoyBB_t != NULL) {umat_cont_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
	 if (ex.gradVmoyBB_tdt != NULL) {umat_cont_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
	 if (ex.gradVBB_tdt != NULL) {umat_cont_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
	}; // fin de la partie dédiée à la vitesse
 // on commence par recopier les grandeurs de l'ordre 2 à 3
 bool plusZero = true; // on complète avec des 0 dans un premier temps
 int type_recopie=0; //   = 0 -> on transfert les grandeurs à 0, t et tdt
	umat_cont_3D->Passage_de_Ordre2_vers_Ordre3(ex,plusZero,type_recopie);
		 	
	// maintenant on s'occupe de mettre à jour les grandeurs manquantes
	// - les bases naturelles: le vecteur normal est normé et est identique pour les bases naturelles et duales
	giB_0_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_0_3D(1),giB_0_3D(2))).Normer();
	giH_0_3D.CoordoH(3) = giB_0_3D(3).Bas_haut();
	giB_t_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_t_3D(1),giB_t_3D(2))).Normer();
	giH_t_3D.CoordoH(3) = giB_t_3D(3).Bas_haut();
	   // cas particulier du vecteur tdt
	giB_tdt_3D.CoordoB(3) = (Util::ProdVec_coorB(giB_tdt_3D(1),giB_tdt_3D(2))); // calcul du vecteur normal non normé
	double norme_N_tdt = giB_tdt_3D(3).Norme(); // calcul de la norme qui nous servira pour les variations
 giB_tdt_3D.CoordoB(3) /= norme_N_tdt;
 giB_normer_3_tdt_3D_sauve = giB_tdt_3D(3); // on sauvegarde la valeur du vecteur giB(3) normée
	giH_tdt_3D.CoordoH(3) = giB_tdt_3D(3).Bas_haut();
   //  - les tenseurs métriques: au début 1 pour la direction 3
	gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
                     =gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
 
 return umat_cont_3D;

// //   -- pour les tenseurs métriques, à ce niveau, la valeur de la norme du vecteur g3 est fixe donc pas de variation
// // ensuite quand on intègrera la variation d'épaisseur, on introduira également les dérivées
//
//
//  //  -- prise en compte de la variation de la déformation d'épaisseur : on considère une déformation de type logarithmique
//  // eps33 = log (h/h0) -> h = h0 * exp(eps33) , pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = 1
//  // donc à t et tdt il faut multiplier par exp(eps33) à t et tdt
//     // récup du conteneur spécifique
//  SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
//  // lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
//  // ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
//  // ensuite elles évolueront
//  bool premier_passage = false;
//  if (save_resul.h_tsurh0 == 0.)
//    { save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.; premier_passage = true;
//    };
//  if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
//    { save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
//  // on calcul le vecteur 3 naturel et dual
//  giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
//  giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
//  //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
//  //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
//  double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
//  gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
//
//  double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
//  gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
 
};

  // passage des informations liées à la déformation  de 2 vers 3, et variation de volume

void LoiContraintesPlanes::Passage_deformation_volume_ordre2_vers_3
								(const TenseurBB& DepsBB,const TenseurBB & epsBB_tdt,const Tableau <TenseurBB *>* d_epsBB
								,const TenseurBB & delta_epsBB,const double& jacobien_0,const double& jacobien
								,const Vecteur* d_jacobien_tdt)

 { // au début on complète avec des 0, puis dans un second temps on tient compte de la déformation d'épaisseur eps33
   bool plusZero = true;
   Deps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &DepsBB),plusZero);
   eps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &epsBB_tdt),plusZero);
   delta_eps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &delta_epsBB),plusZero);
   // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
   if(d_epsBB != NULL)
    {  int taille = d_epsBB->Taille();
       // redimensionnement éventuel
       int ta_d_eps_BB_3D = d_eps_BB_3D.Taille();
       if (ta_d_eps_BB_3D != taille)
        {  // cela veut dire que tous les tableaux sont mal dimensionnés
           d_eps_BB_3D.Change_taille(taille); d_eps_BB_3D_P.Change_taille(taille);
           for (int i=1;i<=taille;i++) {d_eps_BB_3D_P(i) = &(d_eps_BB_3D(i));}			
        };
       // passage 2D 3D
       for (int i=1;i<=taille;i++)
        { d_eps_BB_3D(i).Affectation_2D_a_3D(*((Tenseur2BB*) (*d_epsBB)(i)),plusZero);}
       if (save_resul.d_hsurh0.Taille() != taille)
           save_resul.d_hsurh0.Change_taille(taille); // donc à 0 la première fois

    };
   // -- calcul du jacobien 3D 
   // première phase on se contente de recopier, ce qui suppose une épaisseur = 1.
   jacobien_0_3D = jacobien_0;
   jacobien_tdt_3D = sauve_jacobien2D_tdt = jacobien;
   // ici, l'épaisseur initiale h0 est toujours 1., donc h/h0 = également h
   jacobien_tdt_3D *= save_resul.hsurh0;
   // -- maintenant on tient compte de la déformation d'épaisseur
   switch (type_de_deformation)
    {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
          // cas d'une déformation d'Almansi
          {  // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
             // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
             eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
             delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
           
             // dans le cas où on calcule les variations, on intègre également les variations de ep33
             if (d_epsBB != NULL)
              { Vecteur& d_hsurh0 = save_resul.d_hsurh0;
                int taille = d_epsBB->Taille();
                for (int i=1;i<=taille;i++)
                  {d_eps_BB_3D(i).Coor(3,3)= save_resul.hsurh0 * d_hsurh0(i);}
              };
          };
          break;
     case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
         // cas d'une def logarithmique ou une approximation
          {  eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
             delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
           
             // dans le cas où on calcule les variations, on intègre également les variations de ep33
             if (d_epsBB != NULL)
              { Vecteur& d_hsurh0 = save_resul.d_hsurh0;
                int taille = d_epsBB->Taille();
                // d_eps_BB_3D(3,3) = d_hsurh0 * exp(save_resul.hsurh0)
                for (int i=1;i<=taille;i++)
                  {d_eps_BB_3D(i).Coor(3,3)=  d_hsurh0(i) * eps_BB_3D(3,3);}
              };
          };
          break;
     default :
       cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
            << Nom_type_deformation(type_de_deformation);
       cout << "\n LoiContraintesPlanes::Passage_deformation_volume_ordre2_vers_3 \n";
       Sortie(1);
    };
  
   if(d_jacobien_tdt != NULL)
     { Vecteur& d_hsurh0 = save_resul.d_hsurh0;
       // jacobien_tdt_3D = hsurh0 * jacobien_2D
       // --> d_jacobien_tdt_3D = d_hsurh0 * jacobien_2D + hsurh0 * d_jacobien_2D
       d_jacobien_tdt_3D = d_hsurh0 * jacobien + save_resul.hsurh0 * (*d_jacobien_tdt);
     };
  
 };

   // idem mais sans variation: passage des informations liées à la déformation  de 2 vers 3
void LoiContraintesPlanes::Passage_deformation_volume_ordre2_vers_3
                       (const TenseurBB& DepsBB,const TenseurBB & epsBB_tdt
                        ,const TenseurBB & delta_epsBB,const double& jacobien_0,const double& jacobien)
 { // au début on complète avec des 0, puis dans un second temps on tient compte de la déformation d'épaisseur eps33
   bool plusZero = true;
   Deps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &DepsBB),plusZero);
   eps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &epsBB_tdt),plusZero);
   delta_eps_BB_3D.Affectation_2D_a_3D(*((Tenseur2BB*) &delta_epsBB),plusZero);
   // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
   // -- calcul du jacobien 3D
   // première phase on se contente de recopier, ce qui suppose une épaisseur = 1.
   jacobien_0_3D = jacobien_0;
   jacobien_tdt_3D = sauve_jacobien2D_tdt = jacobien;
   // ici, l'épaisseur initiale h0 est toujours 1., donc h/h0 = également h
   jacobien_tdt_3D *= save_resul.hsurh0;
   // -- maintenant on tient compte de la déformation d'épaisseur
   switch (type_de_deformation)
    {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
          // cas d'une déformation d'Almansi
          {  // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
             // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
             eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
             delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
          };
          break;
     case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
         // cas d'une def logarithmique ou une approximation
          {  eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
             delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
          };
          break;
     default :
       cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
            << Nom_type_deformation(type_de_deformation);
       cout << "\n LoiContraintesPlanes::Passage_deformation_volume_ordre2_vers_3 \n";
       Sortie(1);
    };
  
 };



// mise à jour des informations liées à la déformation  de 2 vers 3
void LoiContraintesPlanes::Mise_a_jour_deformations_et_Jacobien_en_3D
								()

 { // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
   // -- calcul du jacobien 3D
   // première phase on se contente de recopier, ce qui suppose une épaisseur = 1.
   jacobien_tdt_3D = sauve_jacobien2D_tdt * save_resul.hsurh0;
   // -- maintenant on tient compte de la déformation d'épaisseur
   switch (type_de_deformation)
    {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
          // cas d'une déformation d'Almansi
          {  // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
             // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
             eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
             delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
          };
          break;
     case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
         // cas d'une def logarithmique ou une approximation
          {  eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
             delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
          };
          break;
     default :
       cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
            << Nom_type_deformation(type_de_deformation);
       cout << "\n LoiContraintesPlanes::Passage_deformation_volume_ordre2_vers_3 \n";
       Sortie(1);
    };
   // -- cas de la vitesse de déformation
   // recup de l'incrément de temps
   double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
   double unSurDeltat=0;    
   if (Abs(deltat) >= ConstMath::trespetit)
    {unSurDeltat = 1./deltat;}
   else
   // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
    { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
      if (unSurDeltat < 0)
       { cout << "\n le pas de temps est négatif !! "; };
      unSurDeltat = ConstMath::tresgrand;
    };
   Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3) * unSurDeltat;
  
  
 };

   
// calcul de la déformation d'épaisseur
void LoiContraintesPlanes::Calcul_eps33_parVarVolume(double& jaco_2D_0,const double& module_compressibilite,double& jaco_2D_tdt
                                        ,TenseurBH& sigBH)
 { // récup des infos spécifiques
   SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
   // on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
   // avec trace(sigma) en 2D = idem en 3D car on est en contraintes planes
   // jacobien 2D = la surface donc le volume = jacobien 2D * h -> (jaco2D_tdt * h-jaco2D_0*h_0) = delta_V
   // en appelant B= trace(sigma)/3 /K_moy on a:
   // h/h_0 = jaco2D_0 / jaco2D_tdt * 1./(1.-B) d'où la déformation logarithmique
   if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
     { //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
       double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
       double exp_log_VsurV0 = exp(log_VsurV0);
       double hSurh_0 = (jaco_2D_0 / jaco_2D_tdt) * exp_log_VsurV0 ;

       // on test les mini-maxi
       if (Limitation_h(hSurh_0))
         {if (Permet_affichage() > 3)
           {cout << "\n LoiContraintesPlanes::Calcul_eps_trans_parVarVolume(... " << flush;
           };
         };
       save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
     }
   else
     { save_resul.hsurh0 = 1.;};
 };

// calcul de la déformation d'épaisseur et de sa variation 
void LoiContraintesPlanes::Calcul_d_eps33_parVarVolume(double& jaco_2D_0,const double& module_compressibilite,double& jaco_2D_tdt
                    ,TenseurHH& sigHH_,Vecteur& d_jaco_2D
                    ,Tableau <TenseurHH *>& d_sigHH,Tableau <TenseurBB *>& d_gijBB_tdt,TenseurBB & gijBB_)
  { // récup des infos spécifiques
    SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
    // on utilise la relation: log(V/V0) = trace(sigma)/3 /K_moy
    // avec trace(sigma) en 2D = idem en 3D car on est en contraintes planes
    // jacobien 2D = la surface, donc le volume: V = jacobien 2D * h -> (jaco2D_tdt * h)/(jaco2D_0*h_0) = V/V0
    // en appelant log_VsurV0= trace(sigma)/3 /K_moy on a:
    // h/h_0 = jaco2D_0 / jaco2D_tdt * exp(log_VsurV0)
    if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
      { const Tenseur2BB & gijBB = *((Tenseur2BB*) &gijBB_); //   pour simplifier
        const Tenseur2HH & sigHH = *((Tenseur2HH*) &sigHH_); //   pour simplifier
        Tenseur2BH sigBH = gijBB * sigHH;
        double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
        double exp_log_VsurV0 = exp(log_VsurV0);
        double hSurh_0 = (jaco_2D_0 / jaco_2D_tdt) * exp_log_VsurV0 ;

        // on test les mini-maxi
        if (Limitation_h(hSurh_0))
         {if (Permet_affichage() > 3)
           {cout << "\n LoiContraintesPlanes::Calcul_d_eps33_parVarVolume(... " << flush;
           };
         };
        save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
       
        // maintenant le cas de la variation de h/h0
        // on tient compte de la variation relative à la dépendance du jacobien et à la trace de sigma
        // a) vérification et modif éventuelle de la taille
        Vecteur& d_hsurh0 = save_resul.d_hsurh0;
        // b) calcul de la variation
        // d h/h0 = - djaco_2D * (jaco_2D_0 / ((jaco_2D_tdt)^2)) * exp(log_VsurV0)
        //          + d_(sigBH.Trace()) * (jaco_2D_0 / jaco_2D_tdt) * exp(log_VsurV0)  / (3. * module_compressibilite )
        int taille = d_hsurh0.Taille();
        double fact1 = - jaco_2D_0  / (jaco_2D_tdt*jaco_2D_tdt) * exp_log_VsurV0;
        double fact2 = (jaco_2D_0 / jaco_2D_tdt) * exp_log_VsurV0 / (3. * module_compressibilite );
       
        for (int i=1;i<=taille;i++)
         {
          Tenseur2HH & dsigHH = *((Tenseur2HH*) (d_sigHH(i)));  // passage en dim 2
          const Tenseur2BB & dgijBB = *((Tenseur2BB*)(d_gijBB_tdt(i))) ; // pour simplifier l'ecriture
          Tenseur2BH d_sigBH = gijBB * dsigHH + dgijBB * sigHH;
          d_hsurh0(i) = fact1 * d_jaco_2D(i) + fact2 * (d_sigBH.Trace());
         };
       }
    else
      { save_resul.hsurh0 = 1.;save_resul.d_hsurh0.Zero();};
  };

// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D
void LoiContraintesPlanes::Calcul_invariants_et_def_cumul()
  { // récup des infos spécifiques
    SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
    // recup de l'incrément de temps
    double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
    double unSurDeltat=0;    
    if (Abs(deltat) >= ConstMath::trespetit)
        {unSurDeltat = 1./deltat;}
    else
        // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
        { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
            if (unSurDeltat < 0)
                { cout << "\n le pas de temps est négatif !! "; };
            unSurDeltat = ConstMath::tresgrand;
        }; 

    // ---calcul des invariants de déformation et de vitesse de déformation
    Tenseur3BH espBH =  eps_BB_3D * gijHH_tdt_3D;
    Tenseur3BH delta_espBH =  delta_eps_BB_3D * gijHH_tdt_3D;
    save_resul.epsInvar(1) =  espBH.Trace();
    save_resul.epsInvar(2) =  espBH.II();
    save_resul.epsInvar(3) =  espBH.Det();
    save_resul.depsInvar(1) =  unSurDeltat * delta_espBH.Trace();
    save_resul.depsInvar(2) =  unSurDeltat * delta_espBH.II();
    save_resul.depsInvar(3) =  unSurDeltat * delta_espBH.Det();
    
    // cas des grandeurs cumulées (cf. la classe Deformation)
    
    // tableau relatif aux différentes grandeurs de type def scalaires équivalentes
    // def_equi(1) = deformation cumulée = somme des sqrt(2./3. * (delta_eps_barre_BH && delta_eps_barre_BH)) ;
    // def_equi(2) = deformation duale de la contrainte de mises = sqrt(2./3. * (eps_barre_BH && eps_barre_BH)) ;
    // def_equi(3) = niveau maxi atteind par def_equi(2) 
    // def_equi(4) = delta def cumulée = sqrt(2./3. * (delta_eps_barre_BH && delta_eps_barre_BH));

    double delta_eps_equi = sqrt(2./3. * ( (delta_espBH && delta_espBH) - Sqr(delta_espBH.Trace()) /3. ));
    save_resul.def_equi(1) = save_resul.def_equi_t(1) + delta_eps_equi ;
    save_resul.def_equi(2) = sqrt(2./3. * (espBH && espBH)) ;
    if (save_resul.def_equi(2) > save_resul.def_equi_t(3))
       save_resul.def_equi(3) = save_resul.def_equi(2);
    save_resul.def_equi(4) = delta_eps_equi;
    
  };

// calcul de la fonction résidu de la résolution de l'équation constitutive: sig33(h/h0) = 0
// h/h0 est l'inconnue du problème et est l'élément d'entrée: x(1)
// l'argument test ramène
//     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
//       fatal, qui invalide le calcul du résidu
Vecteur& LoiContraintesPlanes::Residu_constitutif (const double & alpha,const Vecteur & x, int& test)
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
    test = 1; // init par défaut
   
    // %%% on met à jour les conteneurs locaux en fonction des paramètres d'entrée %%%
    //  -- prise en compte de la variation de la déformation d'épaisseur
    // pour nous ici on travaille avec giB_03D qui est normé et = 1 c-a-d h0 = 1
    // donc à tdt il faut multiplier par h/h0 à tdt
    double hsurh0 = x(1);
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    // vérif des grandeurs
   
    if (Limitation_h(hsurh0))
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanes::Residu_constitutif(... " << flush;
       erreur_sortie = true;
       test = -1; // l'erreur n'est pas fatale mais elle signifie que ce n'est pas très normal !!
      };

    if (erreur_sortie && (Permet_affichage() > 4))
     {cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
      cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
      cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
     };
   
    save_resul.hsurh0 = hsurh0; // sauvegarde de l'élongation d'épaisseur
    giB_tdt_3D.CoordoB(3) = giB_normer_3_tdt_3D_sauve * hsurh0;
    giH_tdt_3D.CoordoH(3) = (giB_normer_3_tdt_3D_sauve / hsurh0).Bas_haut();
    //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
    //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
    double hsurh0_2 = hsurh0 * hsurh0;
    gijBB_tdt_3D.Coor(3,3) = hsurh0_2; gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
    // mise à jour des informations liées à la déformation  de 2 vers 3
    Mise_a_jour_deformations_et_Jacobien_en_3D();
   
    // calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
    // correspondant aux cas 3D
    Calcul_invariants_et_def_cumul();
   
    // initialisation du comportement tangent / au def
    d_sigma_deps_3D.Inita(0.);
    // appel du calcul de sig et dsig
    bool en_base_orthonormee = false; // ici les tenseurs ne sont pas forcément en orthonormee
    lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
                  ,sig_HH_3D,d_sigma_deps_3D
                  ,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D,module_cisaillement_3D,*umat_cont_3D);
    // retour des grandeurs voulues
    residu(1) = sig_HH_3D(3,3); // - save_resul.niveau_sig33;
//    double sig33 = sig_HH_3D(3,3); // pour simplifier
//    double val_orthonormee_sig33 = Dabs(sig33);
//    // soit on est inférieur à la limite et dans ce cas c'est ok -> on met à 0 le résidu
//    if (val_orthonormee_sig33 > save_resul.niveau_sig33)
//      {residu(1) = sig_HH_3D(3,3) - Signe(sig33) * save_resul.niveau_sig33;}
//    else
//      {residu(1) = 0.;}
   
    return residu;

    };

// calcul de la matrice tangente de la résolution de l'équation constitutive: sig33(h/h0) = 0
// h/h0 est l'inconnue du problème et est l'élément d'entrée: x(1)
// l'argument test ramène
//     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
//       fatal, qui invalide le calcul du résidu et de la dérivée

Mat_abstraite& LoiContraintesPlanes::Mat_tangente_constitutif
                              (const double & alphap,const Vecteur & x, Vecteur& residu, int& test)
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanes & save_resul = *((SaveResul_LoiContraintesPlanes*) saveResul);
    test = 1; // init par défaut

    // %%% on met à jour les conteneurs locaux en fonction des paramètres d'entrée %%%
    //  -- prise en compte de la variation de la déformation d'épaisseur
    // pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = 1
    // donc à tdt il faut multiplier par h/h0 à tdt
    double hsurh0 = x(1);
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    // vérif des grandeurs
    if (Limitation_h(hsurh0))
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanes::Mat_tangente_constitutif(... " << flush;
       erreur_sortie = true;
       test = -1; // l'erreur n'est pas fatale mais elle signifie que ce n'est pas très normal !!
      };

    if (erreur_sortie && (Permet_affichage() > 4))
     {cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
      cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
      cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
     };
    save_resul.hsurh0 = hsurh0; // sauvegarde de l'élongation d'épaisseur
    giB_tdt_3D.CoordoB(3) = giB_normer_3_tdt_3D_sauve * hsurh0;
    giH_tdt_3D.CoordoH(3) = (giB_normer_3_tdt_3D_sauve / hsurh0).Bas_haut();
    //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
    //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
    double hsurh0_2 = hsurh0 * hsurh0;
    gijBB_tdt_3D.Coor(3,3) = hsurh0_2;
    gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
    // mise à jour des informations liées à la déformation  de 2 vers 3
    Mise_a_jour_deformations_et_Jacobien_en_3D();
   
   
    // calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
    // correspondant aux cas 3D
    Calcul_invariants_et_def_cumul();
   
    // initialisation du comportement tangent / au def
    d_sigma_deps_3D.Inita(0.);
    // appel du calcul de sig et dsig
    bool en_base_orthonormee = false; // ici les tenseurs ne sont pas a priori en orthonormee
    lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
                  ,sig_HH_3D,d_sigma_deps_3D
                  ,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D,module_cisaillement_3D,*umat_cont_3D);

    // retour des grandeurs voulues

    residu(1) = sig_HH_3D(3,3);// - save_resul.niveau_sig33;
//    double sig33 = sig_HH_3D(3,3); // pour simplifier
//    double val_orthonormee_sig33 = Dabs(sig33);
//    // soit on est inférieur à la limite et dans ce cas c'est ok -> on met à 0 le résidu
//    if (val_orthonormee_sig33 > save_resul.niveau_sig33)
//      {residu(1) = sig_HH_3D(3,3) - Signe(sig33) * save_resul.niveau_sig33;}
//    else
//      {residu(1) = 0.;}

    // la matrice tangente dépend du type de mesure de déformation, avec la relation
    // d_sig33(h/h0) / d_h/h0 = sig33(h/h0)/d_eps33 * d_eps33/d_h/h0
    double d_eps33_sur_d_hSurh0 = 0.; // init
   
    // -- choix en fonction du type de déformation
    switch (type_de_deformation)
     {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
           // cas d'une déformation d'Almansi
           {  // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
              // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
              d_eps33_sur_d_hSurh0 = hsurh0;
           };
           break;
      case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
          // cas d'une def logarithmique ou une approximation: eps_BB_3D(3,3) = log(hsurh0);
           {  d_eps33_sur_d_hSurh0 = log(hsurh0);
           };
           break;
      default :
        cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
             << Nom_type_deformation(type_de_deformation);
        cout << "\n LoiContraintesPlanes::Mat_tangente_constitutif \n";
        Sortie(1);
     };
   
    // d'où la dérivée du résidu
    derResidu(1,1) = d_sigma_deps_3D(3,3,3,3) * d_eps33_sur_d_hSurh0;
    return derResidu;
   };

// limitation des variations d'épaisseurs
// ramène true s'il y a eu une modif
bool LoiContraintesPlanes::Limitation_h(double& hsurh0)
 {  bool retour = false;
    // il ne faut pas que la nouvelle épaisseur soit inf au mini
    if (hsurh0 < mini_hsurh0)
      {if (Permet_affichage() > 3)
         {cout << "\n *** attention hsurh0 " << hsurh0 << ", est inferieur a la limite fixee "<<mini_hsurh0
               << " on le ramene a la limite ";
         };
       hsurh0 = mini_hsurh0;retour = true;
      };
    // on regarde également pour les maxi
    if (hsurh0 > maxi_hsurh0)
     {if (Permet_affichage() > 3)
        {cout << "\n *** attention hsurh0 " << hsurh0 << ", est superieur a la limite fixee "<<maxi_hsurh0
              << " on le ramene a la limite ";
        };
      hsurh0 = maxi_hsurh0;retour = true;
     };
    // retour
    return retour;
 };