// FICHIER : LoiContraintesPlanesDouble.cp
// CLASSE : LoiContraintesPlanesDouble


// 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 "MathUtil2.h"

#include "LoiContraintesPlanesDouble.h"

            

// calcul des contraintes et ses variations  par rapport aux déformations a t+dt
// ceci dans le cas où l'axe de traction simple est quelconque
// ici on suppose que:
//   - tous les calculs sont en dim 3
//   - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
//
// ViB et ViH : correspond à la base de traction: les vecteurs sont identiques
//              par contre leurs composantes sont différentes car exprimés dans g^i et g_i
//   telle que \vec ViB(1) = \vec ViH(1) = l'axe de traction
//              La base Vi est supposée orthonormée
// en_base_orthonormee:  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 !!** non pas tout le temps a revoir
void LoiContraintesPlanesDouble::Calcul_dsigma_deps
     (bool en_base_orthonormee,const BaseB * ViB,const BaseH * ViH
      ,const TenseurHH & sigHH_t,const TenseurBB& DepsBB
      ,const TenseurBB & epsBB_tdt,const 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_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   module_compressibilite=module_cisaillement=0.; // init
   energ.Inita(0.); // initialisation des énergies mises en jeux
   //récup de eps_mecaBB_t
   // 1) dans le repère gi
   Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
   Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
   // 2) dans le dernier repère V_a: a utiliser en sortie post, pas pour le calcul !!
   Tenseur3BB& eps_P_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_P_mecaBB_t));
   Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_resul.eps_P_mecaBB));

   // on récupère la dimension de l'espace de travail qui doit-être de dimension 3
   // ainsi que le nombre de vecteur: pour l'instant 3
   #ifdef MISE_AU_POINT
   {if (Permet_affichage() > 5)
       cout << "\n  LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
    if  ((ViB->Dimension() != 3) || (ViH->Dimension() !=  3))
        { cout << "\n***Erreur : la dimension de la base ViB = "<< ViB->Dimension()
               << " ou la dimension de la base ViH = " << ViH->Dimension()
               << " est different de 3, on ne peut pas continuer ! "
               << "\n  LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
          Sortie(1);
        };
    if  ((ViB->NbVecteur() != 3) || (ViH->NbVecteur() !=  3))
        { cout << "\n***Erreur : le nombre de vecteur de la base ViB = "
               << ViB->NbVecteur()
               << " ou le nombre de vecteur de la base ViH = " << ViH->NbVecteur()
               << " est different de 3,  on ne peut pas continuer ! "
               << "\n  LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
          Sortie(1);
        };
    // idem pour l'Umat: là on ne test que les  vecteurs à 0, pour alléger le test
    if  (ex.giB_0->Dimension() != 3)
        { cout << "\n***Erreur : la dimension de la base ex.giB_0 = "
               << ex.giB_0->Dimension()
               << " du conteneur Umat, est different de 3, on ne peut pas continuer ! "
               << "\n  LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
          Sortie(1);
        };
    if  (ex.giB_0->NbVecteur() != 3)
        { cout << "\n***Erreur : le nombre de vecteur de la base ex.giB_0 = "
               << ex.giB_0->NbVecteur()
               << " du conteneur Umat est different de 3,  on ne peut pas continuer ! "
               << "\n  LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
          Sortie(1);
        };
   }
   #endif
  
   // les coordonnées des vecteurs \vec V_e sont exprimées
   // dans la base de travail g_i finale on a donc
   //    Vi_B(i) = {\beta'}_i^{.j} \hat{\vec g}_j = betaP(i,j) * giB_tdt(j)
   //   i indice de ligne, j indice de colonne
   //   \vec V_i = \vec V^i : les vecteurs sont égaux mais pas les composantes !

   // 1 -- passage de la base \vec V_e à la base \hat{\vec g}_i
   //                      et \vec V^e à la base \hat{\vec g}^i
   // calcul des grandeurs gamma et beta:
   // \hat{\vec g}^i  =   \gamma^i_{~.e}~\vec V^e
   // \hat{\vec g}_i  =   \beta_i^{~.e}~\vec V_e
  
  
   // 2 -- passage de la base \hat{\vec g}_i à la base \vec V_e
   //                      et \hat{\vec g}^i à la base \vec V^e
   //
   //   Vi_B(i) = {\beta'}_i^{.j} \hat{\vec g}_j
   //   Vi_H(i) = {\gamma'}^i_{.j} \hat{\vec g}^j
  
   // gamma3D(i,e)    =   \gamma^i_{~.e} = \hat{\vec g}^i . \vec V_e
   // ce qui donne que [gamma3D(i,e)] = [{\beta'}_i^{.j}]^T
   for (int i=1;i<4;i++) // calcul de la matrice beta'
    for (int e=1;e<4;e++) // au lieu de passer par un produit scalaire on récupère
      {betaP_3D(i,e) = (*ViH)(i)(e); // directement les coordonnées
       gammaP_3D(i,e) = (*ViB)(i)(e); // "   "
      };

   // \beta_i^{~.e} = \hat{\vec g}_i . \vec V^e
   // c'est aussi l'inverse de la transposé de gamma3D
   // on choisit de calculer avec la formule
   gamma3D = betaP_3D.Transpose(); // ok
   beta3D = gammaP_3D.Transpose(); // betaP_3D.Inverse();
 //  gammaP_3D = beta3D.Transpose();
 
 
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 6)
    {cout << "\n Vi_B(1)= "<<(*ViB)(1) << " Vi_B(2)= "<<(*ViB)(2) << " Vi_B(3)= "<<(*ViB)(3);
     cout << "\n Vi_H(1)= "<<(*ViH)(1) << " Vi_H(2)= "<<(*ViH)(2) << " Vi_H(3)= "<<(*ViH)(3);
     cout << "\n Vi_B(i) = {beta'}_i^{.j} hat{vec g}_j --> beta'= ";betaP_3D.Affiche();
     cout << "\n Vi_H(i) = {gamma'}^i_{.j} hat{vec g}^j--> gamma'= ";gammaP_3D.Affiche();
     cout << "\n hat{vec g}^i  =   gamma^i_{~.e}~vec V^e --> gamma= ";gamma3D.Affiche();
     cout << "\n hat{vec g}_i  =   beta_i^{~.e}~vec V_e--> beta= ";beta3D.Affiche();
    };
   #endif


   // on alimente les pointeurs de tenseurs d'entrée de la cinématique
   epsBB_tdt_cin = ((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
   DepsBB_cin = ((Tenseur3BB*) &DepsBB); // passage en dim 3
   delta_epsBB_cin = ((Tenseur3BB*) &delta_epsBB); // passage en dim 3
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 6)
     {cout << "\n epsBB_tdt_cin "<< (*epsBB_tdt_cin) << "\n DepsBB_cin "<< (*DepsBB_cin)
           << "\n delta_epsBB_cin "<< (*delta_epsBB_cin) << flush;
      // en base orthonormee
      Tenseur3BB epsBB_tdt_cin_inter;
      epsBB_tdt_cin->BaseAbsolue(epsBB_tdt_cin_inter,(* ex.giH_tdt));
      cout << "\n epsBB_tdt_cin= en absolue 3D : "<<epsBB_tdt_cin_inter;
     };
   #endif

   // --- initialisation des déformations mécaniques dans le repère V_a
   eps_P_BB = eps_mecaBB_t; // init: là on est dans le repère g^i
   // on passe ensuite dans le repère V_a qui est donc le nouveau repère, pour ce passage
   eps_P_BB.ChBase(betaP_3D);
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 6)
     {cout << "\n eps_mecaBB_t= "; eps_mecaBB_t.Ecriture(cout);
      cout << "\n eps_P_BB apres change base donc en V_a "<< eps_P_BB;
     };
   #endif

   // on calcul la déformation cinématique dans la direction V_1
   // \epsilon'_{11}(t+\Delta t) = \epsilon_{kl}(t+\Delta t) ~\gamma^k_{~.1} ~\gamma^l_{~.1}
   // idem pour l'accroissement de déformation et la vitesse de déformation
   eps_P_BB.Coor(1,1) = 0.; Deps_P_BB.Coor(1,1)  = 0.; delta_eps_P_BB.Coor(1,1) = 0.;
   for (int k=1;k<4;k++) // calcul de la matrice gamma
    for (int l=1;l<4;l++)
      {eps_P_BB.Coor(1,1) += (*epsBB_tdt_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
       Deps_P_BB.Coor(1,1) += (*DepsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
       delta_eps_P_BB.Coor(1,1) += (*delta_epsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
      };
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 6)
     {cout << "\n eps_P_BB avec la cinematique "<< eps_P_BB;
      cout << "\n Deps_P_BB avec la cinematique "<< Deps_P_BB;
      cout << "\n delta_eps_P_BB avec la cinematique "<< delta_eps_P_BB << flush;
     };
   #endif
   // on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
   // sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
   sig_HH_t_3D = sigHH_t;
   // concernant la métrique relative au repère de travail, on utilise celle passée en paramètre
   (*umat_cont_3D)= ex;
   // concernant le jacobien, le jacobien du conteneur est supposé être celui 1D
   sauve_jacobien1D_tdt = *(ex.jacobien_tdt);
  
   // --- pour les contraintes passage en 3D
//   bool plusZero = true;
   // --- pour la cinématique
   // 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 coursVhVg
   // la base de traction
   ViB_3D = ViB; ViH_3D = ViH;
   // calcul d'élements intermédiaires
   V1V1_BB = Tenseur3BB::Prod_tensoriel((*ViB_3D)(1),(*ViB_3D)(1));
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 7)
     {  cout << "\n V1V1_BB: "<< V1V1_BB;
     }
   #endif
   // pour (m,n) = (2,2); (3,3); (1,2),  et pour (g,h) = (2,2); (3,3); et (1,2)
   // on va utiliser des matrices intermédiaires tels que:
   // in = 1 pour (2,2) , 2 pour (3,3) et 3 pour (1,2)

   for (int gh=1;gh<4;gh++)
       {//if (gh <3)
         {VhVg_BB(gh) = Tenseur3BB::Prod_tensoriel((*ViB_3D)(jnd(gh)),(*ViB_3D)(ind(gh)));}
//        else // donne des tenseurs non symétriques pour (1,2) -> on symétrise
//         {VhVg_BB(gh) = 0.5* Tenseur3BB::Prod_tensoriel((*ViB_3D)(ind(gh)),(*ViB_3D)(jnd(gh)))
//                       + 0.5* Tenseur3BB::Prod_tensoriel((*ViB_3D)(jnd(gh)),(*ViB_3D)(ind(gh)));
//         }
        #ifdef MISE_AU_POINT
        if (Permet_affichage() > 7)
          {  cout << "\n VhVg_BB("<<gh<<")= "<< VhVg_BB(gh);
          }
        #endif
       };
  
   // on sauvegarde les dimensions transverse, car elles vont être modifiées et si on doit
   // repartir au départ, ce ne sera plus possible sans cette sauvegarde
   double sauve_hsurh0=save_resul.hsurh0;
   double sauve_bsurb0=save_resul.bsurb0;
   // on initialise les jacobiens à 0 et t
   jacobien_0_3D = jacobien_0; //  à t=0 -> save_resul.hsurh0 * save_resul.bsurb0 == 1;
   jacobien_t_3D = (*ex.jacobien_t) * save_resul.h_tsurh0 * save_resul.b_tsurb0;
  
   // 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
   bool newtonlocal = false;
   switch (type_de_contrainte)
   {case NEWTON_LOCAL: // cas de la version d_sig/d_eps
     { newtonlocal=true; // pour la gestion du cas par défaut
       // initialisation du tenseur final de contrainte
       sig_HH_3D = sigHH_tdt; // c-a-d de la solution: peut-être pas nécessaire ***
       // passage des informations liées à la déformation  de 1 vers 3, et variation de volume
       Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
     
       // on appel la procédure de résolution de sig33(..)=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
       // on initialise les inconnues à 0. On pourrait partir des def_P sauvegardés
       // mais si on change d'orientation dans la pratique la convergence est moins bonne ??
       val_initiale_3D.Zero(); // = save_resul.def_P_t;
//       val_initiale_3D = save_resul.def_P_t;
       // on impose que les grandeurs soient dans les limites admises
       double& hsurh0 = save_resul.h_tsurh0;double& bsurb0 = save_resul.b_tsurb0; // init
          // la méthode va mettre à jour les variations bsurb0 et hsurh0 et le jacobien
       if (Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D))
         {if (Permet_affichage() > 3)
            cout << "\n pb limitation val initiales: LoiContraintesPlanesDouble::Calcul_dsigma_deps(... "
                 << flush;
         };

//// ---- debut
//         // on vérifie que les anciennes transversales ne sont pas négatives
//         if ((val_initiale(1) < 0.) || (val_initiale(2) < 0.) )
//             { cout << "\n debug LoiContraintesPlanesDouble::Calcul_DsigmaHH_tdt :  "
//                    << " initialement "
//                    << "\n save_resul.b_tsurb0 = " << save_resul.b_tsurb0
//                    << ", save_resul.h_tsurh0 " << save_resul.h_tsurh0 << endl;
//               Sortie(1);
//             };
////------ fin debug

       residu_3D.Zero();            // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
       derResidu_3D.Initialise(0.); // init de la matrice dérivée à 0.
       try // on met le bloc sous surveillance
        {
         //  résolution de l'équation constitutive d'avancement discrétisée en euler implicite
         int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs

         // 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 << "\nErreur : 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 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 << "\nErreur : 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);
            };
         };
       
         bool conver=alg_zero.Newton_raphson
                             (*this,&LoiContraintesPlanesDouble::Residu_constitutif_3D
                             ,&LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D
                             ,val_initiale_3D,racine_3D,derResidu_3D,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_3D.Max_val_abs();
         if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
             { if (Permet_affichage() > 0)
                 { cout << "\n**>> non convergence sur l'algo de la resolution de (sigma.V_2=0, sigma.V_3=0),\n";
                   if (ptintmeca_en_cours != NULL) ptintmeca_en_cours->Signature();
                   cout  << ", nb_incr_total=" << nb_incr_total
                         << ", nb_iter_total=" << nb_iter_total<<", val obtenues:"
                         << "\n eps_P(2,2)= " << racine_3D(1) << ", "
                         << " eps_P(3,3)= " << racine_3D(2) << ", "<< " eps_P(1,2)= " << racine_3D(3)
                         << "\n precedentes: "<<eps_mecaBB_t(2,2)<<", "<<eps_mecaBB_t(3,3)<<", "<<eps_mecaBB_t(1,2);
                   cout  << "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
                         << "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
                 };
               // on garde en mémoire
               mauvaise_convergence=true;
             };
        }
       catch (ErrNonConvergence_Newton erreur)
           {
            if (Permet_affichage() > 0)
                 { cout << "\n erreur de non convergence avec l'algo de la resolution de"
                        << " (sig'22, sig'33, sig'12)=(0,0,0) "
                        << " on obtient une valeur infinie  ou NaN "
                        << "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
                        << "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
                 };
            // 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 erreur non identifiee sur l'algo de la resolution de "
                        << " (sig'22, sig'33, sig'12)=(0,0,0) "
                        << "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
                        << "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
                 };
            // on garde en mémoire
            mauvaise_convergence=true;
          };
       if (!mauvaise_convergence)
          { // arrivée ici, cela veut dire que tout à bien fonctionné
            save_resul.def_P = racine_3D; // récup de la solution dans le repère V_a
            // on doit sauvegarder les def méca dans le repère g^i, ce sont directement
            // les dernières composantes utilisées pour l'appel de la loi 3D
            eps_mecaBB = eps_BB_3D; // dans g^i
            eps_P_mecaBB = eps_P_BB; // dans V_a pour le post traitement
            // on met à jour les modules
            module_compressibilite= module_compressibilite_3D;
            module_cisaillement= module_cisaillement_3D;
           
            // retour du tenseurs résultat
            sigHH_tdt = sig_HH_3D;

            // calcul de l'opérateur tangent final: ici en 3D
            // d sigma^{ij} / d epsilon_{kl}
            Calcul2_dsigma_deps(d_sigma_deps);
           
            #ifdef MISE_AU_POINT
            if (Permet_affichage() > 6)
             { cout << "\nLoiContraintesPlanesDouble::Calcul_dsigma_deps: bonne convergence de Newton ";
               TenseurHH* ptHH = NevezTenseurHH(sig_HH_3D);
               sig_HH_3D.BaseAbsolue(*ptHH,giB_tdt_3D);
               cout << "\n sigma apres CP2: ";
               ptHH->Ecriture(cout);
               delete ptHH;
             };
            if (Permet_affichage() > 5)
              {  cout << "\n def_P(2,2) " << save_resul.def_P(1)
                      << " def_P(3,3) " << save_resul.def_P(2)
                      << " def_P(1,2) " << save_resul.def_P(3)
                      << "\n en g^i: eps_mecaBB.Coor(2,2) " << eps_mecaBB.Coor(2,2)
                      << " eps_mecaBB.Coor(3,3) " << eps_mecaBB.Coor(3,3)
                      << " eps_mecaBB.Coor(1,2) " << eps_mecaBB.Coor(1,2)
                      << "\n en V_a: eps_P_mecaBB.Coor(2,2) " << eps_P_mecaBB.Coor(2,2)
                      << " eps_P_mecaBB.Coor(3,3) " << eps_P_mecaBB.Coor(3,3)
                      << " eps_P_mecaBB.Coor(1,2) " << eps_P_mecaBB.Coor(1,2)
                      << "\n module_compressibilite= " << module_compressibilite
                      << " module_cisaillement= " << module_cisaillement;
                 // on affiche les épaisseur et largeur
                 double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
                 cout << "\n hsurh0= "<< hsurh0 << " bsurb0= " << bsurb0 << flush;
              };
            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
            break;
           };
         // sinon on ne fait pas de break, donc on continue avec la méthode de perturbation
     };


    default: // cas de la version avec perturbation, on encore du cas où la méthode de Newton n'a pas convergé
       {
        // initialisation du tenseurs contrainte
        sig_HH_3D.Inita(0.);
        // dans le cas ou on est passé par la méthode de newton et que cela n'a pas
        // convergé, on doit réinitialiser les variables de départ
        if ((newtonlocal) && (Permet_affichage() > 3))
           cout << "\n calcul direct sans Newton ";
        save_resul.hsurh0 = sauve_hsurh0;
        save_resul.bsurb0 = sauve_bsurb0;
        
        // appel du calcul en perturbation
        Cal_avec_perturbation2();
        // on met à jour les modules
        module_compressibilite= module_compressibilite_3D;
        module_cisaillement= module_cisaillement_3D;

        if (mauvaise_convergence)
         {if (Permet_affichage() > 0)
           { cout << ": valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0
                  << " et de  b/b_0(t+dt)= " << save_resul.bsurb0;
           };
          // on annulle les dérivées des épaisseurs
          save_resul.d_hsurh0.Zero();save_resul.d_bsurb0.Zero();
         };

     }
    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
    { //Tenseur1BH sigBH = gijBB_tdt * sigHH_tdt;
      Tenseur3BH sigBH = gijBB_tdt_3D * sig_HH_3D;
      Tenseur3BH sigBH_t = gijBB_t_3D * sig_HH_t_3D;
      Calcul_eps_trans_parVarVolume2(sigBH_t,jacobien_0,module_compressibilite,jacobien,sigBH,*(ex.jacobien_t));
      if  ((Permet_affichage() > 0) && mauvaise_convergence)
          { cout << "\n valeur final de  h/h_0(t+dt)= " << save_resul.hsurh0
                 << " et de  b/b_0(t+dt)= " << save_resul.bsurb0;
          };
    };

   LibereTenseur();
   LibereTenseurQ();
};

   
       // 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 LoiContraintesPlanesDouble::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 en 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 LoiContraintesPlanesDouble::RepercuteChangeTemperature(...";
               Sortie(1);      
             };
      }; 

  };
	
void LoiContraintesPlanesDouble::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_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) 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;
   };

  // par défaut on active les invariants et on intègre les infos def  en  3D
  if (!calcul_en_3D_via_direction_quelconque)
   { //dans le cas classique, on dispose de la variation d'épaisseur et de largeur pour déterminer les def
    
     // ptintmec concerne des tenseurs  d'ordres 1 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_1D_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
     double eps22=0.; double eps22_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.);
               eps22 = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
               eps22_t = 0.5 * (save_resul.b_tsurb0 * save_resul.b_tsurb0 - 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);
               eps22 = log(save_resul.bsurb0);
               eps22_t = log(save_resul.b_tsurb0);
            };
            break;
       default :
         cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
              << Nom_type_deformation(type_de_deformation);
         cout << "\n LoiContraintesPlanesDouble::CalculGrandeurTravail \n";
         Sortie(1);
      };
    
     (*ptintmeca.EpsBB()).Coor(3,3) = eps33;
     (*ptintmeca.DeltaEpsBB()).Coor(3,3) = eps33 - eps33_t;
     (*ptintmeca.EpsBB()).Coor(2,2) = eps22;
     (*ptintmeca.DeltaEpsBB()).Coor(2,2) = eps22 - eps22_t;
    
     (*ptintmeca.DepsBB()).Coor(3,3) = (eps33 - eps33_t) * unSurDeltat;
     (*ptintmeca.DepsBB()).Coor(2,2) = (eps22 - eps22_t) * unSurDeltat;
    }
  else
    // dans le cas où on travaille dans le repère gi en 3D on dispose directement des infos de déformation
    {//récup de eps_mecaBB
     Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
     Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
     // on sauve ce qu'il y a dans la direction 11 pour éviter les erreurs d'arrondis
     double eps11= (*ptintmeca.EpsBB())(1,1);
     double Deps11= (*ptintmeca.DepsBB())(1,1);
     double delta_eps11= (*ptintmeca.DeltaEpsBB())(1,1);
     // les manips entre tenseurs
     (*ptintmeca.EpsBB()) = eps_mecaBB;
     (*ptintmeca.DeltaEpsBB()) = eps_mecaBB - eps_mecaBB_t;
     (*ptintmeca.DepsBB()) = (*ptintmeca.DeltaEpsBB()) * unSurDeltat;
     // récup des valeurs sauvegardées (normalement ne devrait pas changer...)
     (*ptintmeca.EpsBB()).Coor(1,1) = eps11;
     (*ptintmeca.DeltaEpsBB()).Coor(1,1) = delta_eps11;
     (*ptintmeca.DepsBB()).Coor(1,1) = Deps11;
    };

  // on met à jour les invariants 3D,  on récupère ce qui a été sauvegardé
  if (ptintmeca.Statut_Invariants_deformation()) // le conteneur des invariants a pu être effacé par Affectation_1D_a_3D
     ptintmeca.EpsInvar() =  save_resul.epsInvar;
  if (ptintmeca.Statut_Invariants_vitesseDeformation()) // le conteneur des invariants a pu être effacé par Affectation_1D_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 1 à 3, cas implicite
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_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);};
			   };
		 };
 // récup du conteneur spécifique
 SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);

 // si les grandeurs à 0 ont déjà été calculées, on les sauvegarde d'abord car elles vont
 // être écrasées 
 // on commence par recopier les grandeurs de l'ordre 1 à 3
 bool plusZero = true; // on complète avec des 0 dans un premier temps
 int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
 if ( save_resul.meti_ref_00.giB_ != NULL)
    type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
	impli_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
 #ifdef MISE_AU_POINT
 if (Permet_affichage() > 8)
  { cout << "\n LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3("
         << "\n 1) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(2)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
  };
 #endif

	// maintenant on s'occupe de mettre à jour les grandeurs manquantes
	// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
 //   et sont identiques pour les bases naturelles et duales
 // Dans le cas d'un comportement isotrope, on choisit un repère par défaut
 
 // si les repères associés ne sont pas définis on les définit
 if ( save_resul.meti_ref_00.giB_ == NULL)
  { Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
    save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
    // on s'occupe de la création du premier repère
    CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
    giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
    giB0_1.Normer(); // on le norme
    // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
    CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
    CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
    // deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
    if (ex.giB_0->NbVecteur() > 1)
      { // a priori le second vecteur de la base est un vecteur d'un plan de ref
        giB0_3  = (Util::ProdVec_coorB(giB_0_3D(1), (*(ex.giB_0))(2))).Normer();
        giB0_2 =  Util::ProdVec_coorB(giB0_3,giB0_1);
      }
    else
      { // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
        MathUtil2::Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
      };
////--- debug
//     cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
//     cout << "\n giB0_1: "; giB0_1.Affiche();
//     cout << "\n giB0_2: "; giB0_2.Affiche();
//     cout << "\n giB0_3: "; giB0_3.Affiche();
//// --- fin debug
    giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
    giH_0_3D.CoordoH(2) = giB_0_3D(2).Bas_haut();
    giH_0_3D.CoordoH(3) = giB_0_3D(3).Bas_haut();
  };
////--- debug
//     cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
//     cout << "\n giH_0(2): "; giH_0_3D(2).Affiche();
//     cout << "\n giH_0(3): "; giH_0_3D(3).Affiche();
//// --- fin debug
 // on met à jour les deux autres directions
 // *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
 // *** , mais il faudra surement améliorer
 CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
 giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
 giBt_1.Normer(); // on le norme
 // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
 CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
 CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
 // deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
 if (ex.giB_t->NbVecteur() > 1)
   { // a priori le second vecteur de la base est un vecteur d'un plan de ref
     giBt_3  = (Util::ProdVec_coorB(giB_t_3D(1), (*(ex.giB_t))(2))).Normer();
     giBt_2 =  Util::ProdVec_coorB(giBt_3,giBt_1);
   }
 else
   { // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
     MathUtil2::Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
   };
 giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
 giH_t_3D.CoordoH(2) = giB_t_3D(2).Bas_haut();
 giH_t_3D.CoordoH(3) = giB_t_3D(3).Bas_haut();
 #ifdef MISE_AU_POINT
 if (Permet_affichage() > 8)
  { cout << "\n 2) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(2)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
  };
 #endif

 
 CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
 giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
 giBtdt_1.Normer(); // on le norme
 // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
 CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
 CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
 // deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
 if (ex.giB_tdt->NbVecteur() > 1)
   { // a priori le second vecteur de la base est un vecteur d'un plan de ref
     giBtdt_3  = (Util::ProdVec_coorB(giB_tdt_3D(1), (*(ex.giB_tdt))(2))).Normer();
     giBtdt_2 =  Util::ProdVec_coorB(giBtdt_3,giBtdt_1);
   }
 else
   { // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
     MathUtil2::Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
   };
 giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
 giH_tdt_3D.CoordoH(2) = giB_tdt_3D(2).Bas_haut();
 giH_tdt_3D.CoordoH(3) = giB_tdt_3D(3).Bas_haut();

 #ifdef MISE_AU_POINT
 if (Permet_affichage() > 8)
  { cout << "\n 3) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(2)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
  };
 #endif
 // ici on calcul les vecteurs
 // >>>> et la première partie de leurs variations, fonction du premier vecteur uniquement
 
 //  - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
	gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
                     =gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=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.;

 // lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0 et idem pour l'épaisseur
 // ne sont pas affectés car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
 // ensuite elles évolueront,  idem pour la largeur
 bool premier_passage = false;
 if (save_resul.h_tsurh0 == 0.)
   { save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
     save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 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
  if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
    { save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
 // on calcul les vecteur 2 et 3 naturels et duals
 giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
 giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;

 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;

 #ifdef MISE_AU_POINT
 if (Permet_affichage() > 8)
  { cout << "\n 4) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(2)
         << "\n    base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
  };
 #endif
 //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
 //    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;
   
 double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
	gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;

 double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
 gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;

 // -- cas des variations

 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
      save_resul.d_bsurb0.Change_taille(ta_ex_d_giB_tdt);
     };
	
 //  - les variations du vecteur normal
	#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_ordre1_vers_3(const Met_abstraite::Impli& ex)";
		  Sortie(1);
	   };
 #endif

 
  // maintenant on s'occupe de la seconde partie des variations
  // >>>> fonction du changement d'épaisseur et de largeur

  if (!premier_passage)
   {Vecteur& d_hsurh0 = save_resul.d_hsurh0;
    Vecteur& d_bsurb0 = save_resul.d_bsurb0;
    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) = d (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)));
        
         double d_bsurb0_iddl = d_bsurb0(iddl);
         // gib(2) = b/b0  * N -> d_gib(2) = d (b/b0  * N) = d_(b/b0)  * N + b/b0 * 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(2) *= save_resul.bsurb0;
         d_giB_tdt_3D(iddl).CoordoB(2) += d_bsurb0_iddl * giB_normer_3_tdt_3D_sauve; // base naturelle
         // gih(2) = 1/giB(2) -> d_gih(2)(2) = - d_giB(3)(3) / (giB(3)*giB(3))
         d_giH_tdt_3D(iddl).CoordoH(2) = (d_giB_tdt_3D(iddl)(2) / (-bsurb0_carre)).Bas_haut() ; // base duale
         // gijBB(2)= giB(2) * giB(2) -> d_gijBB(2)= 2. * d_giB(2) * giB(2)
         d_gijBB_tdt_3D(iddl).Coor(2,2) = 2. * (d_giB_tdt_3D(iddl)(2).ScalBB(giB_tdt_3D(2)));
         // gijHH(2)= giH(2) * giH(2) -> d_gijHH(2)= 2. * d_giH(2) * giH(2)
         d_gijHH_tdt_3D (iddl).Coor(2,2) = 2. * (d_giH_tdt_3D(iddl)(2).ScalHH(giH_tdt_3D(2)));
        
       };
   };

};

// passage des grandeurs métriques de l'ordre 1 à 3: cas explicite
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_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));};
	  };	 
 
  // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
  
  // on commence par recopier les grandeurs de l'ordre 1 à 3
  bool plusZero = true; // on complète avec des 0 dans un premier temps
  int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
  if ( save_resul.meti_ref_00.giB_ != NULL)
    type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
  expli_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
	
	 // maintenant on s'occupe de mettre à jour les grandeurs manquantes
	 // - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
  //   et sont identiques pour les bases naturelles et duales
  // Dans le cas d'un comportement isotrope, on choisit un repère par défaut

  // si les repères associés ne sont définis on les définit
  if ( save_resul.meti_ref_00.giB_ == NULL)
   { Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
     save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
     // on s'occupe de la création du premier repère
     CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
     giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
     giB0_1.Normer(); // on le norme
     // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
     CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
     CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
//     Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
     giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
   };
  // on met à jour les deux autres directions
  // *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
  // *** , mais il faudra surement améliorer
  CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
  giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
  giBt_1.Normer(); // on le norme
  // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
  CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
  CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
//  Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
  giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;

  
  CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
  giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
  giBtdt_1.Normer(); // on le norme
  // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
  CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
  CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
//  Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
  giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
 
  //  - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
	 gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
                =gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=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.;
 
  // 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,  idem pour la largeur
  bool premier_passage = false;
  if (save_resul.h_tsurh0 == 0.)
    { save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
      save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 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
  if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
    { save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
  // on calcul les vecteur 2 et 3 naturels et duals
  giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
  giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;

  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 et largeurs
  //    sont normées initialement, seule l'intensité
  //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
  double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
	 gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;

  double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
  gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;

  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 grandeurs métriques de l'ordre 1 à 3: cas Umat
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Umat_cont& 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) {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
  // récup du conteneur spécifique
  SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);

  // on commence par recopier les grandeurs de l'ordre 1 à 3
  bool plusZero = true; // on complète avec des 0 dans un premier temps
  int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
  if ( save_resul.meti_ref_00.giB_ != NULL)
    type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
  umat_cont_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
	
	 // maintenant on s'occupe de mettre à jour les grandeurs manquantes
	 // - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
  //   et sont identiques pour les bases naturelles et duales
  // Dans le cas d'un comportement isotrope, on choisit un repère par défaut
 
  // si les repères associés ne sont définis on les définit
  if ( save_resul.meti_ref_00.giB_ == NULL)
   { Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
     save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
     // on s'occupe de la création du premier repère
     CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
     giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
     giB0_1.Normer(); // on le norme
     // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
     CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
     CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
     MathUtil2::Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
     // comme la base est orthonormée, les giH sont identiques aux giB
     (*(save_resul.meti_ref_00.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_00.giB_)));
     // on affecte les vecteurs qui seront transmis via l'umat
     giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
     BaseH& a0H = (*(save_resul.meti_ref_00.giH_)); // pour simplifier
     giH_0_3D.CoordoH(1)=a0H(1);giH_0_3D.CoordoH(2)=a0H(2);giH_0_3D.CoordoH(3)=a0H(3);
   };
////--- debug
//     cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
//     cout << "\n giH_0(2): "; giH_0_3D(2).Affiche();
//     cout << "\n giH_0(3): "; giH_0_3D(3).Affiche();
//// --- fin debug
  // on met à jour les deux autres directions
  // *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
  // *** , mais il faudra surement améliorer
  CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
  giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
  giBt_1.Normer(); // on le norme
  // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
  CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
  CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
  MathUtil2::Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
  // comme la base est orthonormée, les giH sont identiques aux giB
  // ensuite plus bas, les normes seront changées pour tenir compte des variations d'épaisseur et de largeur
  (*(save_resul.meti_ref_t.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_t.giB_)));
  // on affecte les vecteurs qui seront transmis via l'umat
  giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
  BaseH& a_t_H = (*(save_resul.meti_ref_t.giH_)); // pour simplifier
  giH_t_3D.CoordoH(1)=a_t_H(1);giH_t_3D.CoordoH(2)=a_t_H(2);giH_t_3D.CoordoH(3)=a_t_H(3);

  //--- cas à tdt
  CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
  giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
  giBtdt_1.Normer(); // on le norme
  // on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
  CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
  CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
  MathUtil2::Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
  giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
  // comme la base est orthonormée ici, les giH sont identiques aux giB
  // ensuite plus bas, les normes seront changées pour tenir compte des variations d'épaisseur et de largeur
  (*(save_resul.meti_ref_tdt.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_tdt.giB_)));
  // on affecte les vecteurs qui seront transmis via l'umat
  giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
  BaseH& a_tdt_H = (*(save_resul.meti_ref_tdt.giH_)); // pour simplifier
  giH_tdt_3D.CoordoH(1)=a_tdt_H(1);giH_tdt_3D.CoordoH(2)=a_tdt_H(2);giH_tdt_3D.CoordoH(3)=a_tdt_H(3);

  //  - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
	 gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
                      =gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=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.;
 
  // 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,  idem pour la largeur
  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 peut avoir la même chose en b, indépendamment, lorsque l'on passe de contrainte plane
  // à doublement plane
  if (save_resul.b_tsurb0 == 0.)
    { save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
      premier_passage = true;
    };
  if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
    { save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
  // on calcul les vecteur 2 et 3 naturels et duals
  giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
  giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;

  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 et largeurs
  //    sont normées initialement, seule l'intensité
  //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs
  double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
	 gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;

  double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
  gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;

  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 1 vers 3, et variation de volume
//  utilisée pour les méthodes 1 et 2
// entre autre: calcul du jacobien 3D en fonction de l'élongation 1D et des variations h/h0 et b/b0

void LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3
//								(int cas, TenseurBB& DepsBB,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>* d_epsBB
        (TenseurBB& DepsBB,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>* d_epsBB
								,TenseurBB & delta_epsBB,const double& jacobien_0,double& jacobien
								,Vecteur* d_jacobien_tdt,const double& jacobien_t)

 { // 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_trans_dimension(*((Tenseur1BB*) &DepsBB),plusZero);
   eps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &epsBB_tdt),plusZero);
   delta_eps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &delta_epsBB),plusZero);
   // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) 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 1D 3D
       for (int i=1;i<=taille;i++)
        { d_eps_BB_3D(i).Affectation_trans_dimension(*((Tenseur1BB*) (*d_epsBB)(i)),plusZero);}
       if (save_resul.d_hsurh0.Taille() != taille)
           save_resul.d_hsurh0.Change_taille(taille); // donc à 0 la première fois
       if (save_resul.d_bsurb0.Taille() != taille)
           save_resul.d_bsurb0.Change_taille(taille); // donc à 0 la première fois

    };
   // -- calcul du jacobien 3D 
  
/* la partie cas = 2: ne fonctionne pas bien le choix est fait
d'utiliser les variations d'épaisseur et de largeur
   // dans le cas d'un calcul dans une direction quelconque
   // on utilise la déformation mécanique sauvegardée pour calculer les jacobiens
   if (cas == 2)
    { //récup de eps_mecaBB_t
      Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
      // --- on s'occupe des jacobiens
      // à t le jacobien est calculée à partir de la déformation supposée d'Almansi
      // ici il s'agit du jacobien relatif au comportement mécanique
      // |J| = sqrt(det(2 * eps_ij_t + g_ij(0)))
      // calcul de la matrice correspondante à (2 * eps_ij_t + g_ij(0))
      gij_meca_BB_t = (2. * eps_mecaBB_t + (*(umat_cont_3D->gijBB_0)));
      gij_meca_BB_t.Matrice_composante(mat_inter);
      // calcul du jacobien
      jacobien_t_3D = sqrt(Dabs(mat_inter.Determinant()));
     
      // on fait de même pour le jacobien initial ce qui permet d'être cohérent
      // avec la métrique initiale et la déformation au sens d'Almansi
      (umat_cont_3D->gijBB_0)->Matrice_composante(mat_inter);
      jacobien_0_3D = sqrt(Dabs(mat_inter.Determinant()));
      // le jacobien à tdt sera calculé pendant la résolution
    }
   else if (cas == 1)
*/
    {// première phase on se contente de recopier, ce qui suppose une épaisseur = 1.
     jacobien_0_3D = jacobien_0; // -> partie 1D
     jacobien_tdt_3D = sauve_jacobien1D_tdt = jacobien;
     // ici, l'épaisseur initiale h0 est toujours 1., donc h/h0 = également h et idem pour b/b0
     jacobien_tdt_3D *= save_resul.hsurh0 * save_resul.bsurb0;
     // même travail pour t
     jacobien_t_3D = jacobien_t;
     jacobien_t_3D *= save_resul.h_tsurh0 * save_resul.b_tsurb0;
    };
  
 
   // 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;
    };

   // cas des def, delta def et D
   if (!calcul_en_3D_via_direction_quelconque)
    { //dans le cas classique, on dispose de la variation d'épaisseur et de largeur pour déterminer les def
      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);
                // idem pour la largeur
                eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
                delta_eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - save_resul.b_tsurb0 * save_resul.b_tsurb0);
                Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3)* unSurDeltat;
                Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2)* unSurDeltat;


                // dans le cas où on calcule les variations, on intègre également les variations de eps33 et eps22
                if (d_epsBB != NULL)
                 { Vecteur& d_hsurh0 = save_resul.d_hsurh0;
                   Vecteur& d_bsurb0 = save_resul.d_bsurb0;
                   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);
                      d_eps_BB_3D(i).Coor(2,2)= save_resul.bsurb0 * d_bsurb0(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);
                // idem pour la largeur
                eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0);
                delta_eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0) - log(save_resul.b_tsurb0);
                Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3)* unSurDeltat;
                Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2)* unSurDeltat;

                // dans le cas où on calcule les variations, on intègre également les variations de eps33 et eps22
                if (d_epsBB != NULL)
                 { Vecteur& d_hsurh0 = save_resul.d_hsurh0;
                   Vecteur& d_bsurb0 = save_resul.d_bsurb0;
                   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);
                      d_eps_BB_3D(i).Coor(2,2)=  d_bsurb0(i) * eps_BB_3D(2,2);
                     }
                 };
             };
             break;
        default :
          cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
               << Nom_type_deformation(type_de_deformation);
          cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3 \n";
          Sortie(1);
       };
    }
   else
    // dans le cas où on travaille dans le repère gi en 3D on dispose directement des infos de déformation
    {//récup de eps_mecaBB
     Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
     Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
     // on sauve ce qu'il y a dans la direction 11 pour éviter les erreurs d'arrondis
     double eps11= (*ptintmeca.EpsBB())(1,1);
     double Deps11= (*ptintmeca.DepsBB())(1,1);
     double delta_eps11= (*ptintmeca.DeltaEpsBB())(1,1);
     // les manips entre tenseurs
     (*ptintmeca.EpsBB()) = eps_mecaBB;
     (*ptintmeca.DeltaEpsBB()) = eps_mecaBB - eps_mecaBB_t;
     (*ptintmeca.DepsBB()) = (*ptintmeca.DeltaEpsBB()) * unSurDeltat;
     // récup des valeurs sauvegardées (normalement ne devrait pas changer...)
     (*ptintmeca.EpsBB()).Coor(1,1) = eps11;
     (*ptintmeca.DeltaEpsBB()).Coor(1,1) = delta_eps11;
     (*ptintmeca.DepsBB()).Coor(1,1) = Deps11;
    };

   if(d_jacobien_tdt != NULL)
     { Vecteur& d_hsurh0 = save_resul.d_hsurh0;
       Vecteur& d_bsurb0 = save_resul.d_bsurb0;
       // jacobien_tdt_3D = hsurh0 * jacobien_1D
       // --> d_jacobien_tdt_3D = d_hsurh0 * jacobien_1D + hsurh0 * d_jacobien_1D
       d_jacobien_tdt_3D = d_hsurh0 * save_resul.bsurb0 * jacobien
                            + save_resul.hsurh0 * d_bsurb0 * jacobien
                            + save_resul.hsurh0 * save_resul.bsurb0 * (*d_jacobien_tdt);
     };
  
//---debug
//cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3, eps_BB_3D(3,3) = "
//     << eps_BB_3D(3,3) << " jacobien_tdt_3D= " <<  jacobien_tdt_3D;
// --- fin debug
  
      
 };


// mise à jour des informations liées à la déformation  de 1 vers 3
// uniquement méthode 1
void LoiContraintesPlanesDouble::Mise_a_jour_deformations_et_Jacobien_en_3D
								()

 { // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   // -- calcul du jacobien 3D
   // première phase on se contente de recopier, ce qui suppose une épaisseur = 1. et une largeur = 1.
   jacobien_tdt_3D = sauve_jacobien1D_tdt * save_resul.hsurh0 * save_resul.bsurb0;
   // -- 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);
             // idem pour la direction 2
             eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
             delta_eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - save_resul.b_tsurb0 * save_resul.b_tsurb0);
          };
          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);
             // idem pour la direction 2
             eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0);
             delta_eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0) - log(save_resul.b_tsurb0);
          };
          break;
     default :
       cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
            << Nom_type_deformation(type_de_deformation);
       cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_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(2,2) = delta_eps_BB_3D(2,2) * unSurDeltat;
   Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3) * unSurDeltat;
  
 };

   
// calcul de la déformation d'épaisseur et de largeur
// cas de la méthode 1
void LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume1(Tenseur3BH& sigBH_t,double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
                                        ,Tenseur3BH& sigBH,double& jaco_1D_t)
 { // récup des infos spécifiques
   SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   // on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
   // avec trace(sigma) en 1D = idem en 3D car on est en contraintes doublement planes
   // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_0*h_0*b_0) = delta_V
   // on considère de plus que l'on est isotrope donc h/h_0 = b/b_0
   // en appelant B= trace(sigma)/3 /K_moy on a:
   // (h/h_0)^2 = jaco1D_0 / jaco1D_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 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
       // on test les mini-maxi
       if (hSurh_0 < mini_hsurh0)
        {if (Permet_affichage() > 3)
           {cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
                 << " on le ramene a la limite "
                 << "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume(... " << endl;
           };
         hSurh_0 = mini_hsurh0;
        };
       if (hSurh_0 > maxi_hsurh0)
        {if (Permet_affichage() > 3)
           {cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
                 << " on le ramene a la limite "
                 << "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume(... " << endl;
           };
         hSurh_0 = maxi_hsurh0;
        };
       // on met à jour épaisseur et largeur d'une même valeur
       save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
       save_resul.bsurb0 = hSurh_0; // élongation de largeur

////--- debug
//       save_resul.hsurh0 = save_resul.h_tsurh0; // élongation d'épaisseur
//       save_resul.bsurb0 = save_resul.b_tsurb0; // élongation d'épaisseur
////-----fin debug


      }
   else
     { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
  
// ---- la suite est un essai, mais ne semble pas non plus fonctionner avec la loi critère
// car le jacobien à t n'est pas calculé ?? ce qui n'est pas normal
// ticket #120
  
//   // a priori cela ne marche pas bien: on va tenter avec le delta jacobien
//   // on utilise la relation: (V(tdt)-V(t))/V(tdt) = trace(delta_sigma)/3 /K_moy
//   // avec trace(delta_sigma) en 1D = idem en 3D car on est en contraintes doublement planes
//   // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_t*h_t*b_t) = delta_V
//   // on considère de plus que l'on est isotrope donc delta_h/h = delta_b/b d'où h/h_t = b/b_t
//   // en appelant B= trace(delta_sigma)/3 /K / (h_t*b_t) on a:
//   // (h/h_t)^2 = jaco1D_t / (jaco1D_tdt * (1.-B))
//   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 B = (sigBH(1,1)-sigBH_t(1,1)) / (3. * module_compressibilite );
//       double hSurh_t = sqrt((jaco_1D_t / jaco_1D_tdt) / (1.-B)) ;
//      
//       save_resul.hsurh0 = hSurh_t * save_resul.h_tsurh0;
//       save_resul.bsurb0 =  hSurh_t * save_resul.b_tsurb0; // élongation d'épaisseur
//      }
//   else
//     { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
  
  
 };

// calcul de la déformation d'épaisseur et de largeur
// cas de la méthode 2
void LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume2(Tenseur3BH& sigBH_t,double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
                                        ,Tenseur3BH& sigBH,double& jaco_1D_t)
 { // récup des infos spécifiques
   SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   // on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
   // avec trace(sigma) en 1D = idem en 3D car on est en contraintes doublement planes
   // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_0*h_0*b_0) = delta_V
   // on considère de plus que l'on est isotrope donc h/h_0 = b/b_0
   // en appelant B= trace(sigma)/3 /K_moy on a:
   // (h/h_0)^2 = jaco1D_0 / jaco1D_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 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
       // on test les mini-maxi
       if (hSurh_0 < mini_hsurh0)
        {if (Permet_affichage() > 3)
           {cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
                 << " on le ramene a la limite "
                 << "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume2(... " << endl;
           };
         hSurh_0 = mini_hsurh0;
        };
       if (hSurh_0 > maxi_hsurh0)
        {if (Permet_affichage() > 3)
           {cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
                 << " on le ramene a la limite "
                 << "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume1(... " << endl;
           };
         hSurh_0 = maxi_hsurh0;
        };
       // on met à jour épaisseur et largeur d'une même valeur
       save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
       save_resul.bsurb0 = hSurh_0; // élongation de largeur

////--- debug
//       save_resul.hsurh0 = save_resul.h_tsurh0; // élongation d'épaisseur
//       save_resul.bsurb0 = save_resul.b_tsurb0; // élongation d'épaisseur
////-----fin debug


      }
   else
     { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
  
// ---- la suite est un essai, mais ne semble pas non plus fonctionner avec la loi critère
// car le jacobien à t n'est pas calculé ?? ce qui n'est pas normal
// ticket #120
  
//   // a priori cela ne marche pas bien: on va tenter avec le delta jacobien
//   // on utilise la relation: (V(tdt)-V(t))/V(tdt) = trace(delta_sigma)/3 /K_moy
//   // avec trace(delta_sigma) en 1D = idem en 3D car on est en contraintes doublement planes
//   // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_t*h_t*b_t) = delta_V
//   // on considère de plus que l'on est isotrope donc delta_h/h = delta_b/b d'où h/h_t = b/b_t
//   // en appelant B= trace(delta_sigma)/3 /K / (h_t*b_t) on a:
//   // (h/h_t)^2 = jaco1D_t / (jaco1D_tdt * (1.-B))
//   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 B = (sigBH(1,1)-sigBH_t(1,1)) / (3. * module_compressibilite );
//       double hSurh_t = sqrt((jaco_1D_t / jaco_1D_tdt) / (1.-B)) ;
//
//       save_resul.hsurh0 = hSurh_t * save_resul.h_tsurh0;
//       save_resul.bsurb0 =  hSurh_t * save_resul.b_tsurb0; // élongation d'épaisseur
//      }
//   else
//     { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
  
  
 };

// calcul des déformations d'épaisseur et de largeur et des variations
// méthode 1
void LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
                    ,TenseurHH& sigHH_,Vecteur& d_jaco_1D
                    ,Tableau <TenseurHH *>& d_sigHH,Tableau <TenseurBB *>& d_gijBB_tdt,TenseurBB & gijBB_)
  { // récup des infos spécifiques
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
    // on utilise la relation: log(V/V0) = trace(sigma)/3 /K_moy
    // avec trace(sigma) en 1D = idem en 3D car on est en contraintes planes double
    // jacobien 1D = la longueur, donc le volume: V = jacobien 1D * h * b -> (jaco1D_tdt * h * b)/(jaco1D_0*h_0*b_0) = V/V0
    // en appelant log_VsurV0= trace(sigma)/3 /K_moy on a:
    // h/h_0 * b/b_0= jaco1D_0 / jaco1D_tdt * exp(log_VsurV0)
    // ici sans info sup, on considère une loi isotrope d'où h/h_0 = b/b_0
    // -> h/h_0 = (jaco1D_0 / jaco1D_tdt * exp(log_VsurV0))^(0.5)
    if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
      { const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_); //   pour simplifier
        const Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_); //   pour simplifier
        Tenseur1BH sigBH = gijBB * sigHH;
        double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
        double exp_log_VsurV0 = exp(log_VsurV0);
        double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
        // on test les mini-maxi
        if (hSurh_0 < mini_hsurh0)
         {if (Permet_affichage() > 3)
            {cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
                  << " on le ramene a la limite "
                  << "\n LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(... " << endl;
            };
          hSurh_0 = mini_hsurh0;
         };
        if (hSurh_0 > maxi_hsurh0)
         {if (Permet_affichage() > 3)
            {cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
                  << " on le ramene a la limite "
                  << "\n LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(... " << endl;
            };
          hSurh_0 = maxi_hsurh0;
         };
        // on met à jour épaisseur et largeur d'une même valeur
        save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
        save_resul.bsurb0 = hSurh_0; // élongation de largeur
       
        // 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
        //2.* h/h_0 * d h/h0 = - djaco_1D * (jaco_1D_0 / ((jaco_1D_tdt)^2)) * exp(log_VsurV0)
        //          + d_(sigBH.Trace()) * (jaco_1D_0 / jaco_1D_tdt) * exp(log_VsurV0)  / (3. * module_compressibilite )
        int taille = d_hsurh0.Taille();
        double fact1 = - jaco_1D_0  / (jaco_1D_tdt*jaco_1D_tdt) * exp_log_VsurV0 / (2. * hSurh_0);
        double fact2 = (jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0 / (3. * module_compressibilite / (2. * hSurh_0));
////---debug
//cout << "\n Calcul_d_eps33_parVarVolume()  module_compressibilite= " << module_compressibilite
//     << " log_VsurV0= "<< log_VsurV0
//     << " sigBH.Trace()= " <<  sigBH.Trace()
//     << "  hSurh_0 =" << hSurh_0
//     << "\n  fact1 =" << fact1 << "  fact2 =" << fact2;
//// --- fin debug
       
        for (int i=1;i<=taille;i++)
         {
          Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i)));  // passage en dim 1
          const Tenseur1BB & dgijBB = *((Tenseur1BB*)(d_gijBB_tdt(i))) ; // pour simplifier l'ecriture
          Tenseur1BH d_sigBH = gijBB * dsigHH + dgijBB * sigHH;
          d_hsurh0(i) = fact1 * d_jaco_1D(i) + fact2 * (d_sigBH.Trace());
//---debug
//cout << "\n Calcul_d_eps33_parVarVolume()   "
//     << " d_hsurh0("<< i<<")= "<< d_hsurh0(i);
// --- fin debug
         };
        // en isotrope les axes h et b sont équivalents
        save_resul.d_bsurb0 = save_resul.d_hsurh0;
       }
    else
      { save_resul.hsurh0 = save_resul.bsurb0 = 1.;
        save_resul.d_hsurh0.Zero();save_resul.d_bsurb0.Zero();
      };
   
  };

// calcul de la variation des déformations d'hauteur et de largeur, par rapport à la def_11, à cause
// des contraintes planes doubles: première méthode
void LoiContraintesPlanesDouble::Calcul_d_eps_eg_11()
{ // on prépare la résolution
//    ,d_sig_ef_gh(2,2),d_eps_ef_11(2)
  // la matrice
  double a = d_sig_ef_gh(1,1) = d_sigma_deps_3D(2,2,2,2);
  double b = d_sig_ef_gh(1,2) = d_sigma_deps_3D(2,2,3,3);
  double c = d_sig_ef_gh(2,1) = d_sigma_deps_3D(3,3,2,2);
  double d = d_sig_ef_gh(2,2) = d_sigma_deps_3D(3,3,3,3);
  double determinant = a*d - b*c;
//  //----- debug
//  cout << "\n debug LoiContraintesPlanesDouble::Calcul_d_eps_eg_11() "<< "\n d_sig_ef_gh= ";
//  d_sig_ef_gh.Affiche();
//  //---fin debug
  if (Dabs(determinant) > ConstMath::petit)
   { // le second membre
     d_eps_ef_11(1)=-d_sigma_deps_3D(2,2,1,1);
     d_eps_ef_11(2)=-d_sigma_deps_3D(3,3,1,1);
//     //----- debug
//     cout << "\n  second membre= "<< d_eps_ef_11(1) << ", " << d_eps_ef_11(2);
//     //---fin debug
     // résolution
     d_sig_ef_gh.Resol_systID(d_eps_ef_11);
//     //----- debug
//     cout << "\n  solution=  "<< d_eps_ef_11(1) << ", " << d_eps_ef_11(2);
//       // vérif
//     a = d_sig_ef_gh(1,1) = d_sigma_deps_3D(2,2,2,2);
//     b = d_sig_ef_gh(1,2) = d_sigma_deps_3D(2,2,3,3);
//     c = d_sig_ef_gh(2,1) = d_sigma_deps_3D(3,3,2,2);
//     d = d_sig_ef_gh(2,2) = d_sigma_deps_3D(3,3,3,3);
//     Vecteur truc(2);
//     truc = d_sig_ef_gh * d_eps_ef_11;
//     cout << "\n  second membre correspondant  "<< truc(1) << ", " << truc(2)<< flush;
//     //---fin debug
   }
  else // sinon on met à 0
   {d_eps_ef_11.Inita(0.);};
};


// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D: valable pour les 2 méthodes
void LoiContraintesPlanesDouble::Calcul_invariants_et_def_cumul()
  { // récup des infos spécifiques
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) 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: sig^ef(h/h0,b/b0) = 0
// avec "e,f" = 2 et 3
// h/h0, bb0 sont les inconnues du problème et sont l'élément d'entrée: <x(1),x(2)>
// 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& LoiContraintesPlanesDouble::Residu_constitutif (const double & alpha,const Vecteur & x, int& test)
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) 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 bsurb0 = x(1);double hsurh0 = x(2);
    // vérif des grandeurs
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    if (Limitation_h_b(hsurh0,bsurb0))
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanesDouble::Residu_constitutif(... " << endl;
       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
    save_resul.bsurb0 = bsurb0; // sauvegarde de l'élongation de largeur
    CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
    CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
    giB_tdt_3D.CoordoB(3) = giBtdt_3 * hsurh0; giH_tdt_3D.CoordoH(3) = (giBtdt_3 / hsurh0).Bas_haut();
    giB_tdt_3D.CoordoB(2) = giBtdt_2 * bsurb0; giH_tdt_3D.CoordoH(2) = (giBtdt_2 / bsurb0).Bas_haut();
    //  - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs
    // et en largeurs sont normées initialement, seule l'intensité
    //    des normales est à mettre à jour compte tenue de la variation d'épaisseurs et de largeur
    double hsurh0_2 = hsurh0 * hsurh0;
    gijBB_tdt_3D.Coor(3,3) = hsurh0_2; gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
    double bsurb0_2 = bsurb0 * bsurb0;
    gijBB_tdt_3D.Coor(2,2) = bsurb0_2; gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_2;
   
    // mise à jour des informations liées à la déformation  de 1 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(2,2);
    residu(2) = sig_HH_3D(3,3);
    return residu;

    };

// calcul de la matrice tangente de la résolution de l'équation constitutive: sig^ef(h/h0,b/b0) = 0
// avec "e,f" = 2 et 3
// h/h0, bb0 sont les inconnues du problème et sont l'élément d'entrée: <x(1),x(2)>
// 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& LoiContraintesPlanesDouble::Mat_tangente_constitutif
                              (const double & alphap,const Vecteur & x, Vecteur& residu, int& test)
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) 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 et de largeur
    // pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = b0 = 1
    // donc à tdt il faut multiplier par h/h0 à tdt et idem pour b/b0
    double bsurb0 = x(1);double hsurh0 = x(2);
    // vérif des grandeurs
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    if (Limitation_h_b(hsurh0,bsurb0))
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif(... " << endl;
       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
    save_resul.bsurb0 = bsurb0; // sauvegarde de l'élongation de largeur
    CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
    CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
    giB_tdt_3D.CoordoB(3) = giBtdt_3 * hsurh0;
    giH_tdt_3D.CoordoH(3) = (giBtdt_3 / hsurh0).Bas_haut();
    giB_tdt_3D.CoordoB(2) = giBtdt_2 * bsurb0;
    giH_tdt_3D.CoordoH(2) = (giBtdt_2 / bsurb0).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;
    double bsurb0_2 = bsurb0 * bsurb0;
    gijBB_tdt_3D.Coor(2,2) = bsurb0_2;
    gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_2;

    // mise à jour des informations liées à la déformation  de 1 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(2,2);
    residu(2) = sig_HH_3D(3,3);
    // la matrice tangente dépend du type de mesure de déformation, avec la relation
    // en appelant hb/hb0 les deux ddl: h/h0 et b/b0, i,j,k et l variant de 2 à 3
    // d_sig^ij(h/h0,b/b0) / d_ddl = d_sig^ij(h/h0,b/b0)/d_eps_kl * d_eps_kl/d_ddl
    // comme les repères restent orthogonales à tous moments il n'y a pas de couplage entre les directions 2 et 3
    Vecteur d_eps_kl_sur_d_ddl(2,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.)
              // idem pour la direction 22
              d_eps_kl_sur_d_ddl(1) = bsurb0;d_eps_kl_sur_d_ddl(2) = 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_eps_kl_sur_d_ddl(1) = log(bsurb0);d_eps_kl_sur_d_ddl(2) = 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 LoiContraintesPlanesDouble::Mat_tangente_constitutif \n";
        Sortie(1);
     };
   
    // d'où la dérivée du résidu, en considérant qu'il n'y a pas de couplage entre les directions 2 et 3
    // c-a-d d_eps33_sur_d_bSurb0 = 0 et d_eps22_sur_d_hSurh0 = 0
    // non c'est faux, du à la partie sphérique des contraintes par exemple, il faut tout considérer
    derResidu(1,1) = d_sigma_deps_3D(2,2,2,2) * d_eps_kl_sur_d_ddl(1);
    derResidu(1,2) = d_sigma_deps_3D(2,2,3,3) * d_eps_kl_sur_d_ddl(2);
    derResidu(2,1) = d_sigma_deps_3D(3,3,2,2) * d_eps_kl_sur_d_ddl(1);
    derResidu(2,2) = d_sigma_deps_3D(3,3,3,3) * d_eps_kl_sur_d_ddl(2);
    return derResidu;
   };

            //--- cas de la résolution de l'équation sigij(eps_meca)*V_j{.l}=0, l=2 et 3
  // calcul de la fonction résidu de la résolution de l'équation constitutive
  // 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& LoiContraintesPlanesDouble::Residu_constitutif_3D(const double & alpha,const Vecteur & x, int& test)
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
    test = 1; // init par défaut
   
    //récup de eps_mecaBB_t
    Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
   
    // les inconnues sont les déformations mécaniques {epsilon'}_kl;
    // (k,l) =  (22)  (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
    // c-a-d les def qui vont servir dans l'appel de la loi 3D
    eps_P_BB.Coor(1,2) = x(3);eps_P_BB.Coor(2,2) = x(1);eps_P_BB.Coor(3,3) = x(2);
    eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
    if (Permet_affichage() > 4)
      {cout << "\n eps_P_BB (dans Va) avant limitation : " ;eps_P_BB.Ecriture(cout);}

    // vérif des grandeurs et calcul des variations d'épaisseur et de largeur
    double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
       // la méthode va mettre à jour les variations bsurb0 et hsurh0 et le jacobien
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    int info_limit = Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D);
    if (info_limit)
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(... " << endl;
       erreur_sortie = true;test = info_limit - 2;
      };
    if (erreur_sortie && (Permet_affichage() > 4))
     {cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
      cout << "\n eps_BB_3D: ";eps_BB_3D.Ecriture(cout);
      cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
      cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) << flush;
     };
   
    // dans la direction V_1 c'est la déformation issue de la cinématique
    // elle a déjà été calculée, et les autres composantes demeurent nulle
    // donc eps_P_BB est au complet,
    eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
    eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
    // on s'occupe maintenant de l'incrément de déformation
    delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
    // au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
    delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
    // recup de l'incrément de temps
    double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
    double unSurDeltat;
    if (Abs(deltat) >= ConstMath::trespetit)
     {unSurDeltat = 1./deltat;}
    else
     // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
     { // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
       #ifdef MISE_AU_POINT
       if (unSurDeltat < 0)
        { cout << "\n le pas de temps est négatif !! "
               << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(...";
        };
       #endif
       unSurDeltat = ConstMath::tresgrand;
     };
    // ici on considère que la déformation est de type Almansi
    Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
   
    // 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();
   
/*    // --- on s'occupe des jacobiens
    // à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
    // ici il s'agit du jacobien relatif au comportement mécanique
    // |J| = sqrt(det(2 * eps_ij + g_ij(0)))
    // calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
    gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
    gij_meca_BB.Matrice_composante(mat_inter);
    // calcul du jacobien
    jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
*/
    // 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 en orthonormee
                   // car on travaille dans la base naturelle
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
     {cout << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(..."
           << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
      cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
      cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
      cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) ;
      // en base orthonormee
      Tenseur3BB epsBB_3D_inter;
      eps_BB_3D.BaseAbsolue(epsBB_3D_inter,*(umat_cont_3D->giH_tdt));
      cout << "\n eps_BB_3D en absolue 3D (dans Va): "<<epsBB_3D_inter;
      cout << "\n jacobien_tdt_3D= "<< jacobien_tdt_3D
           << "\n jacobien_0= " << jacobien_0_3D
           << flush;
     };
    #endif
    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);
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
     {cout << "\n sig_HH_3D: " << sig_HH_3D;
      cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
     };
    #endif

    //-------- calcul du residu -----------
    // residu(mi) = V_{mj} * \sigma^{ij}(\epsilon'_{kl}) : m=2 et 3
   
    CoordonneeH V_sig_H(sig_HH_3D * (*ViB_3D)(2));
    // on considère que les densités d'effort interne dans le plan, la composante 23 devant
    // être nulle compte tenu du tenseur des déformations
//    residu_3D(1) = V_sig_H(1);residu_3D(3) = V_sig_H(2);
    residu_3D(1) = V_sig_H * (*ViB_3D)(2); // composante 22
    residu_3D(3) = (*ViB_3D)(1) * V_sig_H ; // composante 12
// essai avec le mixte
//    residu_3D(3) = 0.5 * ((*ViB_3D)(1) * V_sig_H + V_sig_H * (*ViB_3D)(1)); // composante 12
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
      {cout << "\n en g_i: sig_HH * ViB(2)= " << V_sig_H
            << " df_ds_22= " << residu_3D(1) << " df_ds_21= " << residu_3D(3)
            << flush;
      };
    #endif

    V_sig_H = sig_HH_3D * (*ViB_3D)(3);
//    residu_3D(2) = V_sig_H(3); // seule la composante dans la direction 3 est utilisée
    residu_3D(2) = V_sig_H * (*ViB_3D)(3); // composante 33
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
      {cout << "\n en g_i: sig_HH * ViB(3)= " << V_sig_H << " df_ds_33= " << residu_3D(2);
       // pour info on calcul la projection symétrique
       double f_12 = (*ViB_3D)(1)* sig_HH_3D * (*ViB_3D)(2);
       cout << " df_ds_12= " << f_12 << flush;
      };
    #endif
   
    return residu_3D;

    };
            //--- cas de la résolution de l'équation sigij(eps_meca)*V_j{.l}=0, l=2 et 3
  // calcul de la matrice tangente de la résolution de l'équation constitutive
  // 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& LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D
                              (const double & alphap,const Vecteur & x, Vecteur& residu_local, int& test)
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
    test = 1; // init par défaut

    //récup de eps_mecaBB_t
    Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
   
    // les inconnues sont les déformations mécaniques {epsilon'}_kl;
    // (k,l) =  (22)  (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
    // c-a-d les def qui vont servir dans l'appel de la loi 3D
    eps_P_BB.Coor(1,2) = x(3);eps_P_BB.Coor(2,2) = x(1);eps_P_BB.Coor(3,3) = x(2);
    eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;

    // vérif des grandeurs et calcul des variations d'épaisseur et de largeur
    double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
       // la méthode va mettre à jour les valeurs de bsurb0 et hsurh0, et le jacobien
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    if (Permet_affichage() > 4)
      {cout << "\n eps_P_BB (dans Va) avant limitation : " ;eps_P_BB.Ecriture(cout);};
    int info_limit = Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D);
    if (info_limit)
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(... " << endl;
       erreur_sortie = true;test = info_limit - 2;
      };
    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);
      cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) << flush;
     };

    // dans la direction V_1 c'est la déformation issue de la cinématique
    // elle a déjà été calculée, et les autres composantes demeurent nulle
    // donc eps_P_BB est au complet,
    eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
    eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
    // on s'occupe maintenant de l'incrément de déformation
    delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
    // au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
    delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
    // recup de l'incrément de temps
    double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
    double unSurDeltat;
    if (Abs(deltat) >= ConstMath::trespetit)
     {unSurDeltat = 1./deltat;}
    else
     // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
     { // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
       #ifdef MISE_AU_POINT
       if (unSurDeltat < 0)
        { cout << "\n le pas de temps est négatif !! "
               << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(...";
        };
       #endif
       unSurDeltat = ConstMath::tresgrand;
     };
    // ici on considère que la déformation est de type Almansi
    Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;

    // 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();
   
/*    // --- on s'occupe des jacobiens
    // à t=0 on considère que c'est le même jacobien que transmis
    jacobien_0_3D = *(umat_cont_3D->jacobien_0);
    // à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
    // ici il s'agit du jacobien relatif au comportement mécanique
    // |J| = sqrt(det(2 * eps_ij + g_ij(0)))
    // calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
    gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
    gij_meca_BB.Matrice_composante(mat_inter);
    // calcul du jacobien
    jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
*/
    // 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
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
     {cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(..."
           << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
      cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
      cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
      cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) ;
      // en base orthonormee
      Tenseur3BB epsBB_3D_inter;
      eps_BB_3D.BaseAbsolue(epsBB_3D_inter,*(umat_cont_3D->giH_tdt));
      cout << "\n eps_BB_3D en absolue 3D (dans Va): "<<epsBB_3D_inter;
      cout << "\n jacobien_tdt_3D= "<< jacobien_tdt_3D
           << "\n jacobien_0= " << jacobien_0_3D
           << flush;
     };
    #endif
    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);
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
     {cout << "\n sig_HH_3D: " << sig_HH_3D;
      cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
     };
    #endif

    //-------- calcul du residu -----------
    // residu(mi) = V_{mj} * \sigma^{ij}(\epsilon'_{kl}) : m=2 et 3
   
    //residu(1) = V_2 . \sigma . V_2 = V_{2i} ~\sigma^{ij} ~V_{2j} = ( V_2 \otimes  V_2) :  \sigma
    //residu(2) = V_3 . \sigma . V_3 = V_{3i} ~\sigma^{ij} ~V_{3j} = ( V_3 \otimes  V_3) :  \sigma
    //residu(3) = V_1 . \sigma . V_2 = V_{1i} ~\sigma^{ij} ~V_{2j} = ( V_1 \otimes  V_2) :  \sigma

   
    CoordonneeH V_sig_H(sig_HH_3D * (*ViB_3D)(2)); // = \sigma . V_2
    // ici on utilise le vecteur passé en argument
    residu_local(1) = V_sig_H * (*ViB_3D)(2); // composante 22 = V_2 . \sigma . V_2
//    residu(3) = V_sig_H * (*ViB_3D)(1); // composante 21
    residu_local(3) = (*ViB_3D)(1) * V_sig_H ; // composante 12 = V_1 . \sigma . V_2
// essai avec le mixte
//    residu_local(3) = 0.5 * ((*ViB_3D)(1) * V_sig_H + V_sig_H * (*ViB_3D)(1)); // composante 12
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
      {cout << "\n en g_i: sig_HH * ViB(2)= " << V_sig_H
            << " df_ds_22= " << residu_3D(1) << " df_ds_21= " << residu_3D(3)
            << flush;
      };
    #endif

    V_sig_H = sig_HH_3D * (*ViB_3D)(3); // = \sigma . V_3
    residu_local(2) = V_sig_H * (*ViB_3D)(3); // composante 33 = V_3 . \sigma . V_3
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
      {cout << "\n en g_i: sig_HH * ViB(3)= " << V_sig_H << " df_ds_33= " << residu_3D(2);
       // pour info on calcul la projection symétrique
       double f_12 = (*ViB_3D)(1)* sig_HH_3D * (*ViB_3D)(2);
       cout << " df_ds_12= " << f_12 << flush;
      };
    #endif

    // calcul de la matrice tangente
    // formule générique
    //\frac{\partial \left( \sigma^{ij}(\epsilon'_{kl}) ~ V_{mj} \right )}{\partial \epsilon'_{rs}}
    // = V_{mj}~\frac{\partial \sigma^{ij}(\epsilon_{kl}) }{\partial \epsilon_{op}}~V^r_{~.o}~ V^s_{~.p}
    // ici on doit utiliser uniquement:
    //\frac{\partial \left( \sigma^{3j}(\epsilon'_{kl}) ~ V_{3j} \right ) }{\partial \epsilon'_{rs}}
    //\frac{\partial \left( \sigma^{\alpha j}(\epsilon'_{kl}) ~ V_{2j} \right ) }{\partial \epsilon'_{rs}}
    //  avec  \alpha  =1,2  ; (k,l) et (r,s) } = (2,2) ; (3,3) ; (1,2)
   
    // d  (V_{mi} ~\sigma^{ij} ~V_{nj} \right ) / d \epsilon'_{gh}
    // =  V_{mi}~V_{nj}~ d \sigma^{ij} / d \epsilon_{op} ~V^r_{~.o}~ V^s_{~.p}
    //  pour (m,n) et (g,h)  = (2,2); (3,3); et (1,2)


    derResidu_3D.Initialise(0.); // init
//    for (int j=1;j<4;j++) // j représente la composante de l'axe de projection
//     { for (int i=1;i<4;i++) // i représente le classement de (\epsilon'_{kl})
//               // i= 1 pour (k,l) = (2,2), i=2 pour (k,l) = (3,3), i=3 pour (k,l) = (1,2)
//         {// d (sig(2,j) V(2)(j)) / d (eps')(k,l)
//          derResidu_3D(1,i) +=  (*ViB_3D)(2)(j) *
//                                 (d_sigma_deps_3D(2,j,1,1) * (*ViH_3D)(indx(i))(1)  * (*ViH_3D)(indy(i))(1)
//                                  + d_sigma_deps_3D(2,j,1,2) * (*ViH_3D)(indx(i))(1)  * (*ViH_3D)(indy(i))(2)
//                                  + d_sigma_deps_3D(2,j,2,1) * (*ViH_3D)(indx(i))(2)  * (*ViH_3D)(indy(i))(1)
//                                  + d_sigma_deps_3D(2,j,3,3) * (*ViH_3D)(indx(i))(3)  * (*ViH_3D)(indy(i))(3)
//                                 );
//          // d (sig(1,j) V(2)(j)) / d (eps')(k,l)
//          derResidu_3D(3,i) +=  (*ViB_3D)(2)(j) *
//                                 (d_sigma_deps_3D(1,j,1,1) * (*ViH_3D)(indx(i))(1)  * (*ViH_3D)(indy(i))(1)
//                                  + d_sigma_deps_3D(1,j,1,2) * (*ViH_3D)(indx(i))(1)  * (*ViH_3D)(indy(i))(2)
//                                  + d_sigma_deps_3D(1,j,2,1) * (*ViH_3D)(indx(i))(2)  * (*ViH_3D)(indy(i))(1)
//                                  + d_sigma_deps_3D(1,j,3,3) * (*ViH_3D)(indx(i))(3)  * (*ViH_3D)(indy(i))(3)
//                                 );
//          // d (sig(3,j) V(2)(j)) / d (eps')(k,l)
//          derResidu_3D(2,i) +=  (*ViB_3D)(3)(j) *
//                                 (d_sigma_deps_3D(3,j,1,1) * (*ViH_3D)(indx(i))(1)  * (*ViH_3D)(indy(i))(1)
//                                  + d_sigma_deps_3D(3,j,1,2) * (*ViH_3D)(indx(i))(1)  * (*ViH_3D)(indy(i))(2)
//                                  + d_sigma_deps_3D(3,j,2,1) * (*ViH_3D)(indx(i))(2)  * (*ViH_3D)(indy(i))(1)
//                                  + d_sigma_deps_3D(3,j,3,3) * (*ViH_3D)(indx(i))(3)  * (*ViH_3D)(indy(i))(3)
//                                 );
//
//        };
//     };

// pour (m,n) et (g,h)  = (2,2); (3,3); et (1,2)
// d (V_{mi} ~\sigma^{ij} ~V_{nj}  ) / d\epsilon'_{gh} =
//  (V_{mi} V_{nj}) (d \sigma^{ij} / d \epsilon_{op}) (V^g_{~.o}  V^h_{~.p})
// pour 2,2 et 3,3 pas de pb de symétrie
// pour 1,2
   
   
    for (int mn=1;mn<4;mn++) // mn représente la composante de l'axe de projection
      {// d'abord les termes diagonaux
       for (int gh=1;gh<3;gh++) // gh représente le classement de (\epsilon'_{gh})
          derResidu_3D(mn,gh) = VhVg_BB(mn) && (d_sigma_deps_3D && VhVg_BB(gh));
       // maintenant les termes non diagonaux à doubler
       derResidu_3D(mn,3) = 2.* VhVg_BB(mn) && (d_sigma_deps_3D && VhVg_BB(3));
      };

    return derResidu_3D;
   };


// limitation des variations d'épaisseurs et largeurs
// ramène true s'il y a eu une modif: sert pour la méthode 1
bool LoiContraintesPlanesDouble::Limitation_h_b(double& hsurh0, double& bsurb0)
 {  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;
      };
    // il ne faut pas que la nouvelle largeur soit inf au mini
    if (bsurb0 < mini_bsurb0)
      {if (Permet_affichage() > 3)
         {cout << "\n *** attention bsurb0 " << bsurb0 << ", est inferieur a la limite fixee "<<mini_bsurb0
               << " on le ramene a la limite ";
         };
       bsurb0 = mini_bsurb0;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;
     };
    if (bsurb0 > maxi_bsurb0)
     {if (Permet_affichage() > 3)
        {cout << "\n *** attention bsurb0 " << bsurb0 << ", est superieur a la limite fixee "<<maxi_bsurb0
              << " on le ramene a la limite ";
        };
      bsurb0 = maxi_bsurb0;retour = true;
     };
    // retour
    return retour;
 };


  // calcul de l'opérateur tangent  (2ième méthode)
  // d sigma^{ij} / d epsilon_{kl}
void LoiContraintesPlanesDouble::Calcul2_dsigma_deps(TenseurHHHH& d_sigma_deps)
  {

    #ifdef MISE_AU_POINT
    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 calcule deux grandeurs intermédiaires
    // a) d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}
    // on récupère que le pointeur, ce qui évite d'utiliser l'opérateur de copie du tenseur
    // !!!!**** a vérifier -> oui c'est ok
    //     const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3

    Tenseur3HH dSigDeps11P_HH(d_sigma_deps_3D && V1V1_BB);
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 7)
       { cout << "\n (*ViB_3D)(1)= "<<(*ViB_3D)(1);
         cout << "\n V1V1_BB: "<<  V1V1_BB
              << "\n dSigDeps11P_HH: "<< dSigDeps11P_HH
              << flush;
       };
    #endif

    // b) d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^g_{~.o'}~ V^h_{~.p'}
    Tableau <Tenseur3HH* > dsigdeps_gh_P(3);
    for (int gh=1;gh<4;gh++)
       {dsigdeps_gh_P(gh) = ((Tenseur3HH*) &( d_sigma_deps_3D  && VhVg_BB(gh)));
        #ifdef MISE_AU_POINT
        if (Permet_affichage() > 7)
          {  cout << "\n VhVg_BB("<<gh<<")= "<< VhVg_BB(gh);
             cout << "\n dsigdeps_gh_P(gh): "<< (*dsigdeps_gh_P(gh)) << flush;
          }
        #endif
       };

    // c) mat1(mn,gh) = [ V_{mi}~(d sigma^{ij} / d epsilon_{op}(3D)~~V^g_{~.o}~ V^h_{~.p}) ~~V_{nj}]
    //    SM(mn) =  ( V_{mi}~~ d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}~~V_{nj}  )
    Mat_pleine mat1(3,3);
    Vecteur SM(3);
    for (int mn=1;mn<4;mn++)
     {SM(mn) = VhVg_BB(mn) && dSigDeps11P_HH;
      for (int gh=1;gh<4;gh++)
        mat1(mn,gh) = VhVg_BB(mn) && (*dsigdeps_gh_P(gh));
     };
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 7)
      {  cout << "\n mat1 "<< mat1;
         cout << "\n SM: "<< SM << flush;
      }
    #endif

    // d)  calcul des variations des déformations transversales
    //  en fonction de la déformation longitudinale
    //  d epsilon'_gh / d epsilon'_11 =
    //  - [ V_{mi}~(d sigma^{ij} / d epsilon_{op}(3D)~~V^g_{~.o}~ V^h_{~.p}) ~~V_{nj}]^{-1}~
    //  * ( V_{mi}~~ d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}~~V_{nj}  )
    Vecteur  depsP_gh_depsP_11_BB(3);
    depsP_gh_depsP_11_BB = - mat1.Inverse() * SM;
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 7)
      {  cout << "\n depsP_gh_depsP_11_BB "<< depsP_gh_depsP_11_BB << flush;
      }
    #endif

    // e) d sigma'^{11} / d epsilon'_11 (1D~CP) =
    // = V_{1i}~V_{1j}~d sigma^{ij} / d epsilon_{op}(3D)~~V^1_{~.o}~ V^1_{~.p}
    // + V_{1i}~V_{1j}~d sigma^{ij} / d epsilon_{op}(3D)~~~V^g_{~.o}~ V^h_{~.p}
    //    * d epsilon'_gh / d epsilon'_11
    // ici on se sert du fait que l'on a déjà calculé:
    // dSigDeps11P_HH = d sigma^{ij} / d epsilon_{op}(3D)~~V^1_{~.o}~ V^1_{~.p}
   
    Tenseur1HHHH d_sigP11_depsP11_HHHH
        (  (V1V1_BB && dSigDeps11P_HH)
          + (V1V1_BB &&
                           (  (*dsigdeps_gh_P(1))*depsP_gh_depsP_11_BB(1)
                            + (*dsigdeps_gh_P(2))*depsP_gh_depsP_11_BB(2)
                            + (*dsigdeps_gh_P(3))*depsP_gh_depsP_11_BB(3)
                           )
            )
        );
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 7)
        cout << "\n d_sigP11_depsP11_HHHH "<< d_sigP11_depsP11_HHHH << flush;
    #endif
    // f)  d_sigma^{ij} / d_epsilon_{kl} \approx
    // V_1^{~.i}~V_1^{~.j}~d sigma'^{11} / d epsilon'_11 (1D~CP)~\gamma^{k}_{~.1} ~\gamma^{l}_{~.1}
        //fonctions static définissant le produit tensoriel de deux vecteurs
    // si les vecteurs sont égaux le tenseur est symétrique sinon il est non symétrique
    TenseurHH & inter1_HH = Tenseur3HH::Prod_tensoriel((*ViH_3D)(1),(*ViH_3D)(1));
    Coordonnee3H u_H(gamma3D.Colonne(1)); // qui correspond aux \gamma^{k}_{~.1}
    TenseurHH & inter2_HH = Tenseur3HH::Prod_tensoriel(u_H,u_H);
    d_sigma_deps = Tenseur3HHHH::Prod_tensoriel(inter1_HH,inter2_HH);
    d_sigma_deps *= d_sigP11_depsP11_HHHH(1,1,1,1);
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
        { Tenseur3HHHH  inter_HHHH;
          cout << "\n dans le repere locale d_sigma_deps=  \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(i,j,k,l) << " ; ";
               if (e>6) {cout << "\n"; e=1;}
             };
          Tenseur3HHHH  inter1_HHHH(d_sigma_deps);
          inter1_HHHH.ChangeBase(inter_HHHH,giB_tdt_3D);
          cout << "\n dans le repere orthonormee d_sigma_deps=  \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

  };


  // vérif des grandeurs et calcul des valeurs de l'épaisseur et de la largeur
  // ainsi que calcul du jacobien final
  //   = 0 pas de  modif
  //   = 1 modif due aux maxi permis de h et b, = 2 si def d'entrée incorrecte
  // cas de la 2ième méthode
int LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b
      (double& hsurh0,Tenseur3BB& eps_P_BB_3D, double& bsurb0,double& jacobien_tdt_3DD, const double & jacobien_0_3DD)
 {  int retour = 0; // init par défaut
    bool modif_hsurh0 = false; bool modif_bsurb0 = false;  // init

    // à partir de la valeur des déformations en cours on regarde les variations transversales correspondantes
    // ici les déformations eps_BB_3D sont exprimées dans le repère V_a donc orthonormée
    switch (type_de_deformation)
      {case  DEFORMATION_STANDART :  case DEFORMATION_POUTRE_PLAQUE_STANDART :
            // cas d'une déformation d'Almansi
            { // on est dans le repère V_a, V_3 est direction principale du tenseur de def -> formule classique
              // epsBB33 = 1/2 * (1. -(h_0/h_tdt)**2) avec un maxi de o.5
              if (eps_P_BB_3D(3,3) > 0.5)
                 {if (Permet_affichage() > 0)
                    cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
                       << "\n **erreur: la deformation d'epaisseur dans V_a "<<eps_P_BB_3D(3,3)
                       << ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
                       << " la nouvelle epaisseur, on garde donc l'ancienne " << flush;
                  retour = 2; // là le pb est sérieux
                 }
               else // sinon cas normal
                {hsurh0 = 1./sqrt(1.-2.*eps_P_BB_3D(3,3));
                 // 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; modif_hsurh0 = true;
                    retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
                   };
                 // 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;modif_hsurh0 = true;
                   retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
                  };
                 // si on a modifié via les extrèmes il faut mettre à jour la def
                 if (modif_hsurh0)
                   eps_BB_3D.Coor(3,3) = 0.5 * (1.-hsurh0 * hsurh0);
                };
               // -- suite seulement si retour est différent de 2 ---
               //    c-a-d que l'on n'a pas eu de pb sérieux sur h
               if (retour < 2)
                {// pour les deux autres def on n'est pas en directions principales mais c'est un plan principal
                 // on calcule les def principales via la technique de Mohr (possible car on est en repère ortho)
                 double epsP_c = 0.5 * (eps_P_BB_3D(1,1) + eps_P_BB_3D(2,2)); // le centre
                 double R = sqrt(Sqr(eps_P_BB_3D(2,2) - epsP_c) + Sqr(eps_P_BB_3D(1,2)));
                 double epsP_I = epsP_c + R; // la plus grande
                 double epsP_II = epsP_c - R; // la deuxième mais qui en toute rigueur pourrait être intercalé avec eps_P_BB_3D(3,3)
                 // d'où la variation de jacobien à condition de ne pas avoir de pb
                 bool pb = false;
   //              if ((epsP_I > 0.5) || (epsP_II > 0.5))
                 if (epsP_I > 0.5)
                   {if (Permet_affichage() > 0)
                      cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
                           << "\n **erreur: la deformation eps'_I  "<<epsP_I
                           << ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
                           << " la variation de volume " << flush;
                    epsP_I = 0.5-ConstMath::trespetit; // on le met au maxi
                    retour = 2; // là le pb est sérieux
                    pb = true;
                   };
                 if (epsP_II > 0.5)
                   {if (Permet_affichage() > 0)
                      cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
                           << "\n **erreur: la deformation eps'_II  "<<epsP_I
                           << ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
                           << " la variation de volume " << flush;
                    epsP_II = 0.5-ConstMath::trespetit; // on le met au maxi
                    retour = 2; // là le pb est sérieux
                    pb = true;
                   };
                 // on considère que la direction tranversale est celle correspondante à celle différente de
                 // eps_P_BB_3D(1,1) donc on en déduit une variation de b mais qui ne sera celle de V_2 que
                 // dans un cas isotrope, sinon c'est celle de la direction principale de def dans le plan !
                 if (retour < 2) // là encore on ne modifie que si retour est < 2
                         // car si retour == 2 on va obtenir des bsurb0 infinie
                  
                  // eps_P_BB_3D(1,1), c'est la déformation cinématique imposée selon V1
                  // si elle est > eps_P_BB_3D(2,2), cela veut dire que epsP_I est de son côté
                  // par contre si c'est < , cela veut dire que c'est epsP_II qui est de son côté
                  // d'où le choix qui suit
                  { if (eps_P_BB_3D(1,1) > eps_P_BB_3D(2,2))
                     {bsurb0 = 1./sqrt(1.-2.*epsP_II);}
                    else // cas d'une compression
                     {bsurb0 = 1./sqrt(1.-2.*epsP_I);};
                  };
               
                 // il ne faut pas que la nouvelle largeur soit inf au mini
                 if (bsurb0 < mini_bsurb0)
                   {if (Permet_affichage() > 3)
                      {cout << "\n *** attention bsurb0 " << bsurb0 << ", est inferieur a la limite fixee "<<mini_bsurb0
                            << " on le ramene a la limite ";
                      };
                    bsurb0 = mini_bsurb0; modif_bsurb0 = true;
                    retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
                   };
                 if (bsurb0 > maxi_bsurb0)
                  {if (Permet_affichage() > 3)
                     {cout << "\n *** attention bsurb0 " << bsurb0 << ", est superieur a la limite fixee "<<maxi_bsurb0
                           << " on le ramene a la limite ";
                     };
                   bsurb0 = maxi_bsurb0;modif_bsurb0 = true;
                   retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
                  };
                 // si on a modifié via les extrèmes il faut mettre à jour la def
                 // ici il s'agit uniquement de limiter eps_BB_3D(2,2) à des valeurs non absurdes donc on
                 // fait l'approximation (fausse en ortho mais cela ne devrait pas être un gros pb)
                 if (modif_bsurb0)
                     {eps_BB_3D.Coor(2,2) = 0.5 * (1.-bsurb0 * bsurb0);
                      // on est obligé de refaire une passe sur les def principales qui vont ^tre utilisées
                      // pour la variation de jacobien
                      epsP_c = 0.5 * (eps_P_BB_3D(1,1) + eps_P_BB_3D(2,2)); // le centre
                      R = sqrt(Sqr(eps_P_BB_3D(2,2) - epsP_c) + Sqr(eps_P_BB_3D(1,2)));
                      epsP_I = epsP_c + R; // la plus grande
                      epsP_II = epsP_c - R; // la deuxième mais qui en toute rigueur pourrait être intercalé avec eps_P_BB_3D(3,3)
                      // d'où la variation de jacobien à condition de ne pas avoir de pb
                      if (epsP_I > 0.5)
                        {if (Permet_affichage() > 0)
                           cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
                                << "\n **erreur: la deformation eps'_I (deuxieme passe!!) "<<epsP_I
                                << ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
                                << " la variation de volume " << flush;
                         epsP_I = 0.5-ConstMath::trespetit; // on le met au maxi
                         retour = 2; // là le pb est sérieux
                         pb = true;
                        };
                      if (epsP_II > 0.5)
                        {if (Permet_affichage() > 0)
                           cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
                                << "\n **erreur: la deformation eps'_II  (deuxieme passe!!)  "<<epsP_I
                                << ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
                                << " la variation de volume " << flush;
                         epsP_II = 0.5-ConstMath::trespetit; // on le met au maxi
                         retour = 2; // là le pb est sérieux
                         pb = true;
                        };
                     };

                 if (!pb) // cas normal
                  {double var_vol = 1./sqrt((1-2.*eps_P_BB_3D(3,3))*(1-2.*epsP_I)*(1-2.*epsP_II));
                   //  jacobien_0_3DD contient uniquement le jacobien_0_1D et on part par défaut
                   // avec ici h0 = b0 = 1, c-a-d ici dans la loi de comp: on ne calcule en fait que h/h0 et b/b0
                   // la véritable épaisseur et largeur est obtenue au niveau de l'élément fini
                   jacobien_tdt_3DD = jacobien_0_3DD * var_vol;
                  };
                 // sinon on laise en l'état
               }; // fin du cas retour < 2
            };
            break;
// -- pour l'instant je laisse le log en attente
//       case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
//           // cas d'une def logarithmique ou une approximation
//            { hsurh0 = exp(eps_BB_3D(3,3));
//              bsurb0 = exp(eps_BB_3D(2,2));
//            };
//            break;
       default :
         cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
              << Nom_type_deformation(type_de_deformation);
         cout << "\n LoiContraintesPlanesDouble::Limitation2_h_b(.. \n";
         Sortie(1);
      };
//---- debug
//jacobien_tdt_3DD = jacobien_0_3DD; // pour voir si c'est un pb de jacob ien
    // retour
    return retour;
 };


  // calcul de l'état final, dans le cas de la méthode de perturbation: version 2
  void LoiContraintesPlanesDouble::Cal_avec_perturbation2()
  {
    // récup du conteneur spécifique
    SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   
    //récup de eps_mecaBB_t
    // 1) dans le repère gi
    Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
    // 2) dans le dernier repère V_a: a utiliser en sortie post, pas pour le calcul !!
    Tenseur3BB& eps_P_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_P_mecaBB_t));
    Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_resul.eps_P_mecaBB));

    // les inconnues sont les déformations mécaniques {epsilon'}_kl;
    // (k,l) =  (22)  (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
    // c-a-d les def qui vont servir dans l'appel de la loi 3D
    // contrairement au cas de la méthode de Newton, ici on utilise les grandeurs stockées à t
    // sauf pour la cinématique
    // --- initialisation des déformations mécaniques dans le repère V_a
    eps_P_BB = eps_mecaBB_t; // là on est dans le repère g^i
    // on passe dans le repère V_a qui est donc le nouveau repère, pour ce passage
    eps_P_BB.ChBase(betaP_3D);
   
    // on calcul la déformation cinématique dans la direction V_1
    // \epsilon'_{11}(t+\Delta t) = \epsilon_{kl}(t+\Delta t) ~\gamma^k_{~.1} ~\gamma^l_{~.1}
    // idem pour l'accroissement de déformation et la vitesse de déformation
    eps_P_BB.Coor(1,1) = 0.; Deps_P_BB.Coor(1,1)  = 0.; delta_eps_P_BB.Coor(1,1) = 0.;
    eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
    for (int k=1;k<4;k++) // calcul de la matrice gamma
     for (int l=1;l<4;l++)
       {eps_P_BB.Coor(1,1) += (*epsBB_tdt_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
        Deps_P_BB.Coor(1,1) += (*DepsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
        delta_eps_P_BB.Coor(1,1) += (*delta_epsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
       };
  
    // vérif des grandeurs et calcul de la nouvelle épaisseur et  largeur
    double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
       // la méthode va mettre à jour les valeurs bsurb0 et hsurh0 et du jacobien
    bool erreur_sortie = false; // pour gestion  de sortie d'infos supplémentaires
    if (Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D))
      {if (Permet_affichage() > 3)
         cout << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(... " << endl;
       erreur_sortie = true;
      };
    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);
     };
   
    // dans la direction V_1 c'est la déformation issue de la cinématique
    // elle a déjà été calculée, et les autres composantes demeurent nulle
    // donc eps_P_BB est au complet,
    eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
    eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
    // on s'occupe maintenant de l'incrément de déformation
    delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
    // au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
    delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
    // recup de l'incrément de temps
    double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
    double unSurDeltat;
    if (Abs(deltat) >= ConstMath::trespetit)
     {unSurDeltat = 1./deltat;}
    else
     // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
     { // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
       #ifdef MISE_AU_POINT
       if (unSurDeltat < 0)
        { cout << "\n le pas de temps est négatif !! "
               << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(...";
        };
       #endif
       unSurDeltat = ConstMath::tresgrand;
     };
    // ici on considère que la déformation est de type Almansi
    Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
   
    // 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();
   
    // --- on s'occupe des jacobiens
    // à t=0 on considère que c'est le même jacobien que transmis
    jacobien_0_3D = *(umat_cont_3D->jacobien_0);
    // à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
    // ici il s'agit du jacobien relatif au comportement mécanique
    // |J| = sqrt(det(2 * eps_ij + g_ij(0)))
    // calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
    gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
    gij_meca_BB.Matrice_composante(mat_inter);
    // calcul du jacobien
    jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
   
    // 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 en orthonormee
                   // car on travaille dans la base naturelle
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
     {cout << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(..."
           << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
      cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
      cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
     };
    #endif
    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);
    #ifdef MISE_AU_POINT
    if (Permet_affichage() > 8)
     {cout << "\n sig_HH_3D: " << sig_HH_3D;
      cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
     };
    #endif
    // sauvegarde des résultats
    save_resul.def_P(1) = eps_P_BB(2,2); // récup de la solution dans le repère V_a
    save_resul.def_P(2) = eps_P_BB(3,3); //    "
    save_resul.def_P(3) = eps_P_BB(1,2); //    "
    // on doit sauvegarder les def méca dans le repère g^i, ce sont directement
    // les dernières composantes utilisées pour l'appel de la loi 3D
    Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
    eps_mecaBB = eps_BB_3D; // en g^i
    eps_P_mecaBB = eps_P_BB; // dans V_a pour le post traitement
   };


  // passage entre le calcul classique en 1D de Calcul_dsigma_deps et le calcul 3D
  // méthode 2
  void LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_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
            )
  {
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 7)
     cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_deps";
   #endif
   // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   module_compressibilite=module_cisaillement=0.; // init
   energ.Inita(0.); // initialisation des énergies mises en jeux
      
   // --- pour les contraintes passage en 3D
   bool plusZero = true;
   // on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
   // sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
   sig_HH_t_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_t),plusZero);
   // --- pour la cinématique
   // passage des informations spécifique à la loi le_SaveResul
   // initialisation du tenseurs contrainte
   sig_HH_3D.Inita(0.);
   sig_HH_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_tdt),plusZero);
   // passage des métriques de l'ordre 1 vers 3
   Passage_metrique_ordre1_vers_3(ex);
   // passage des informations liées à la déformation  de 1 vers 3, et variation de volume
   Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
   Vecteur* titi = NULL;                // "    "    "
   Passage_deformation_volume_ordre1_vers_3(DepsBB,epsBB_tdt,toto,delta_epsBB
                                          ,jacobien_0,jacobien,titi,*(ex.jacobien_t));
   
   //  -- construction de la base ortho associée à la direction de traction/compression
   // ici il s'agit pour le premier vecteur de la direction 1,
   // et compte tenu de Passage_metrique_ordre1_vers_3, le repère gi est orthogonale
   // on va donc le normer
   BaseB Vi_B(3); BaseH Vi_H(3);
   BaseB * ViB = &Vi_B;BaseH * ViH = &Vi_H;
   // tout d'abord on récupère la base gi actuelle et on la norme:
   // on se sert de Vi_B comme stockage  transitoire
//   cout << "\n base giB(1) "<< (* ex.giB_tdt)(1)
//          << "\n base giB(2) "<< (* ex.giB_tdt)(2);
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 8)
    { cout << "\n base giB_tdt_3D(1) "<< giB_tdt_3D(1)
           << "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
           << "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
    };
   #endif

   
   Vi_B.CoordoB(1) = giB_tdt_3D(1);Vi_B.CoordoB(1).Normer();
   Vi_B.CoordoB(2) = giB_tdt_3D(2);Vi_B.CoordoB(2).Normer();
   Vi_B.CoordoB(3) = giB_tdt_3D(3);Vi_B.CoordoB(3).Normer();
   // calcul des coordonnées locales de V_a en contravariants:
   // V_a^{~i} = \vec V_a . \vec g^j
   for (int a=1;a<4;a++) for (int i=1;i<4;i++)
     Vi_H.CoordoH(a)(i) = Vi_B.CoordoB(a) * giH_tdt_3D(i);
   // calcul des coordonnées contra
   // Vi_B(a) = g_ij * Vi_H(a)
   for (int a=1;a<4;a++)
    Vi_B.CoordoB(a) = gijBB_t_3D * Vi_H(a);
   
   // appel de la résolution 3D
     // calcul des contraintes et ses variations  par rapport aux déformations a t+dt
   // ceci dans le cas où l'axe de traction simple est quelconque
   // ici on suppose que:
   //   - tous les calculs sont en dim 3
   //   - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
   //
   // ViB et ViH : correspond à la base de traction telle que ViB(1) = ViH(1) = l'axe de traction
   //              La base Vi est supposée orthonormée
   // en_base_orthonormee:  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
   Calcul_dsigma_deps (en_base_orthonormee,ViB,ViH,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 1
   ((Tenseur1HH*) &sigHH_tdt)->Affectation_trans_dimension(sig_HH_3D,false);
   // maintenant on renseigne le tenseur de sortie
   bool pluszero = true;
   d_sigma_deps.Affectation_trans_dimension(d_sigma_deps_3D, pluszero);

  };

  // passage entre le calcul classique en 1D de Calcul_DsigmaHH_tdt et le calcul 3D
  // méthode 2
  void LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_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)
  {
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 7)
     cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_td";
   #endif
   // récup du conteneur spécifique
   SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
   module_compressibilite=module_cisaillement=0.; // init
   energ.Inita(0.); // initialisation des énergies mises en jeux
     
   // --- pour les contraintes passage en 3D
   bool plusZero = true;
   // on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
   // sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
   sig_HH_t_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_t),plusZero);
   // --- pour la cinématique
   // passage des informations spécifique à la loi le_SaveResul
   // initialisation du tenseurs contrainte
   sig_HH_3D.Inita(0.);
   sig_HH_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_tdt),plusZero);
   // passage des métriques de l'ordre 1 vers 3
   Passage_metrique_ordre1_vers_3(ex);
   // passage des informations liées à la déformation  de 1 vers 3, et variation de volume
   Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
   Vecteur* titi = NULL;                // "    "    "
   Passage_deformation_volume_ordre1_vers_3(DepsBB,epsBB_tdt,toto,delta_epsBB
                                         ,jacobien_0,jacobien,titi,*(ex.jacobien_t));
   
   //  -- construction de la base ortho associée à la direction de traction/compression
   // ici il s'agit pour le premier vecteur de la direction 1,
   // et compte tenu de Passage_metrique_ordre1_vers_3, le repère gi est orthogonale
   // on va donc le normer
   BaseB Vi_B(3); BaseH Vi_H(3);
   BaseB * ViB = &Vi_B;BaseH * ViH = &Vi_H;
   // tout d'abord on récupère la base gi actuelle et on la norme:
   // on se sert de Vi_B comme stockage  transitoire
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 8)
    { cout << "\n base giB_tdt_3D(1) "<< giB_tdt_3D(1)
           << "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
           << "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
    };
   #endif

   Vi_B.CoordoB(1) = giB_tdt_3D(1);Vi_B.CoordoB(1).Normer();
   Vi_B.CoordoB(2) = giB_tdt_3D(2);Vi_B.CoordoB(2).Normer();
   Vi_B.CoordoB(3) = giB_tdt_3D(3);Vi_B.CoordoB(3).Normer();
   // calcul des coordonnées locales de V_a en contravariants:
   // V_a^{~i} = \vec V_a . \vec g^j
   for (int a=1;a<4;a++) for (int i=1;i<4;i++)
     Vi_H.CoordoH(a)(i) = Vi_B.CoordoB(a) * giH_tdt_3D(i);
   // calcul des coordonnées contra
   // Vi_B(a) = g_ij * Vi_H(a)
   for (int a=1;a<4;a++)
    Vi_B.CoordoB(a) = gijBB_t_3D * Vi_H(a);
   
   #ifdef MISE_AU_POINT
   if (Permet_affichage() > 4)
     {cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_td( ";
      cout << "\n base Vi_B(1) "<< Vi_B(1) << "\n base Vi_B(2) "<< Vi_B(2)
           << "\n base Vi_B(3) "<< Vi_B(3)
           << "\n base Vi_H(1) "<< Vi_H(1) << "\n base Vi_H(2) "<< Vi_H(2)
           << "\n base Vi_H(3) "<< Vi_H(3) << flush;
     };
   #endif

   
   // appel de la résolution 3D
     // calcul des contraintes et ses variations  par rapport aux déformations a t+dt
   // ceci dans le cas où l'axe de traction simple est quelconque
   // ici on suppose que:
   //   - tous les calculs sont en dim 3
   //   - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
   //
   // ViB et ViH : correspond à la base de traction telle que ViB(1) = ViH(1) = l'axe de traction
   //              La base Vi est supposée orthonormée
   // en_base_orthonormee:  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
   bool en_base_orthonormee = false;
   Calcul_dsigma_deps (en_base_orthonormee,ViB,ViH,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 1
   ((Tenseur1HH*) &sigHH_tdt)->Affectation_trans_dimension(sig_HH_3D,false);
   // maintenant on renseigne le tenseur de sortie
   bool pluszero = true;
   
   // d_sigma / d_ddl = d_sigma / d_eps * d_eps / d_ddl
   // ici l'opérateur tangent correspond à la sensibilité de d sigma'^11 / d eps'_11
   // et la direction V_1 c'est la direction de g_1 (à la norme près donc
   // d_sigma / d_eps normalement c'est == à d sigma'^11 / d eps'_11
   // du coup ici on n'a pas à itérer sur toutes les déformations
   int nb_ddlPlus1 = d_epsBB.Taille()+1;
   for (int i_ddl=1;i_ddl<nb_ddlPlus1;i_ddl++)
     d_sigHH(i_ddl)->Coor(1,1) = (d_sigma_deps_3D(1,1,1,1) * (*d_epsBB(i_ddl)))(1,1);

  };