// FICHIER : Poly_hyper3D.cc
// CLASSE : Poly_hyper3D


// 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 <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "CharUtil.h"

#include "Poly_hyper3D.h"
#include "MathUtil.h"

Poly_hyper3D::Poly_hyper3D () :  // Constructeur par defaut
  Hyper_W_gene_3D(POLY_HYPER3D,CAT_MECANIQUE,3)
  ,Cij(),K(-ConstMath::trespetit),Cij_nD(),K_nD(NULL)
  ,K_use(0.),Cij_use()
  ,Cij_temperature(NULL),K_temperature(NULL),type_pot_vol(1)
  ,avec_courbure(false),a_courbure(0.),r_courbure(0)  ,W_d(0.),W_v(0.)
  ,W_d_J1(0.),W_d_J2(0.),W_v_J3(0.),W_v_J3J3(0.)
  ,W_d_J1_2(0.),W_d_J1_J2(0.),W_d_J2_2(0.)
   {  };
	

// Constructeur de copie
Poly_hyper3D::Poly_hyper3D (const Poly_hyper3D& loi) :
	Hyper_W_gene_3D(loi)
   ,Cij(loi.Cij),K(loi.K),Cij_nD(loi.Cij_nD),K_nD(loi.K_nD)
   ,K_use(loi.K),Cij_use(loi.Cij)
   ,Cij_temperature(loi.Cij_temperature)
   ,K_temperature(loi.K_temperature),type_pot_vol(loi.type_pot_vol)
   ,W_d(loi.W_d),W_v(loi.W_v)
   ,W_d_J1(loi.W_d_J1),W_d_J2(loi.W_d_J2),W_v_J3(loi.W_v_J3),W_v_J3J3(loi.W_v_J3J3)
   ,W_d_J1_2(loi.W_d_J1_2),W_d_J1_J2(loi.W_d_J1_J2),W_d_J2_2(loi.W_d_J2_2)
   ,avec_courbure(loi.avec_courbure),a_courbure(loi.a_courbure),r_courbure(loi.r_courbure)
	 {   // on regarde s'il s'agit de courbes locales ou de courbes globales
      int tai = Cij.Taille();
      for (int i=1;i<= tai; i++)
  	   { int taj = Cij(i).Taille();
  	     for (int j=1;j<=taj;j++)
           if (Cij_temperature(i)(j) != NULL)
              if (Cij_temperature(i)(j)->NomCourbe() == "_") 
                Cij_temperature(i)(j) = Courbe1D::New_Courbe1D(*(loi.Cij_temperature(i)(j)));
  	   };
      if (K_temperature != NULL)   
         if (K_temperature->NomCourbe() == "_") 
            K_temperature = Courbe1D::New_Courbe1D(*(loi.K_temperature));;
      if(avec_courbure)
       { if (a_temperature != NULL)
           if (a_temperature->NomCourbe() == "_")
              a_temperature = Courbe1D::New_Courbe1D(*(loi.a_temperature));
         if (r_temperature != NULL)
            if (r_temperature->NomCourbe() == "_")
               r_temperature = Courbe1D::New_Courbe1D(*(loi.r_temperature));
       };
          
      // idem pour les fonctions nD
      if (K_nD != NULL)
        if (K_nD->NomFonction() == "_")
            {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
             string non_courbe("_");
             K_nD = Fonction_nD::New_Fonction_nD(*loi.K_nD);
            };
      int bai = Cij_nD.Taille();
      for (int i=1;i<= bai; i++)
        { int baj = Cij_nD(i).Taille();
          for (int j=1;j<=baj;j++)
            {if (Cij_nD(i)(j)->NomFonction() == "_")
              {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
               string non_courbe("_");
               Cij_nD(i)(j) = Fonction_nD::New_Fonction_nD(*loi.Cij_nD(i)(j));
              };
            };
        };
	  };

Poly_hyper3D::~Poly_hyper3D ()
// Destructeur 
{ int tai = Cij.Taille();
  for (int i=1;i<= tai; i++)
  	{ int taj = Cij(i).Taille();
  	  for (int j=1;j<=taj;j++)
        if (Cij_temperature(i)(j) != NULL)
           if (Cij_temperature(i)(j)->NomCourbe() == "_") delete Cij_temperature(i)(j);
  	};
  if (K_temperature != NULL)
      if (K_temperature->NomCourbe() == "_") delete K_temperature;
  if (a_temperature != NULL)
      if (a_temperature->NomCourbe() == "_") delete a_temperature;
  if (r_temperature != NULL)
      if (r_temperature->NomCourbe() == "_") delete r_temperature;

  // idem pour les fonctions nD
  if (K_nD != NULL) if (K_nD->NomFonction() == "_") delete K_nD;
  int bai = Cij_nD.Taille();
  for (int i=1;i<= bai; i++)
    { int baj = Cij_nD(i).Taille();
      for (int j=1;j<=baj;j++)
         if (Cij_nD(i)(j) != NULL) if (Cij_nD(i)(j)->NomFonction() == "_") delete Cij_nD(i)(j);
    };
  };
  
// initialise les donnees particulieres a l'elements
// de matiere traite ( c-a-dire au pt calcule)
// Il y a creation d'une instance de SaveResul particuliere
// a la loi concernee
// la SaveResul classe est remplie par les instances heritantes
// le pointeur de SaveResul est sauvegarde au niveau de l'element
// c'a-d que les info particulieres au point considere sont stocke
// au niveau de l'element et non de la loi.
Loi_comp_abstraite::SaveResul * Poly_hyper3D::New_et_Initialise()
  { // on sortie_post
    int para=0;
    if (sortie_post)
      {para++; // premier element K
       int tai = Cij.Taille();
       for (int i=1;i<= tai; i++)
          para += Cij(i).Taille();
      };

    // stockage en conséquence
    SaveResulHyper_W_gene_3D * pt = new SaveResulHyper_W_gene_3D(para);
    // insertion éventuelle de conteneurs de grandeurs quelconque
    this->Insertion_conteneur_dans_save_result(pt);
    return pt;
  };

// Lecture des donnees de la classe sur fichier
void Poly_hyper3D::LectureDonneesParticulieres 
                (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
  { string toto,nom; 
    // lecture du degré maxi du polynôme
    *(entreePrinc->entree) >> nom;   
    if(nom != "degre_maxi=")
         { cout << "\n erreur en lecture du parametre degre_maxi, on attendait la chaine degre_maxi= et on a lu: " << nom;
           cout << "\n Poly_hyper3D::LectureDonneesParticulieres(...  " << endl ;
           entreePrinc->MessageBuffer("erreur 1   ");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);          
          };
    int tai=0; *(entreePrinc->entree) >> tai;
    if (tai < 1)
         { cout << "\n erreur en lecture du parametre degre_maxi,  on a lu: " << tai
                << " qui est plus petit que 1 !! , le minimum c'est 1. ";
           cout << "\n Poly_hyper3D::LectureDonneesParticulieres(...  " << endl ;
           entreePrinc->MessageBuffer("erreur 2   ");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);          
          };
    Cij.Change_taille(tai); // on dimensionne correctement la première dimension de Cij
    Cij_use.Change_taille(tai);
    Cij_temperature.Change_taille(tai); // idem pour la dépendance éventuelle à la température
    // puis on dimensionne la deuxième dimension
    for (int i=1;i<= tai; i++)
      {Cij(i).Change_taille(i+1,(-ConstMath::trespetit)); // init à (-ConstMath::trespetit)
       Cij_use(i).Change_taille(i+1,(-ConstMath::trespetit));
       Cij_temperature(i).Change_taille(i+1,NULL); // idem init à NULL
      };
      
    // maintenant on lit les coefficients, ils doivent être tous présents
    for (int i=1;i<= tai; i++)
      for (int j=0;j<= i; j++)
      	{ // constitution de l'indicateur
      	  string indicateur("C");   indicateur += ChangeEntierSTring(j);indicateur +=ChangeEntierSTring(i-j);
      	  indicateur =indicateur + "=";
      	  // lecture du coefficient Cij 
          *(entreePrinc->entree) >> nom;   
          if(nom != indicateur)
            { cout << "\n erreur en lecture de l'indicateur "<<indicateur<< ",  on a lu: " << nom;
              cout << "\n Poly_hyper3D::LectureDonneesParticulieres(...  " << endl ;
              entreePrinc->MessageBuffer("erreur 3   ");
              throw (UtilLecture::ErrNouvelleDonnee(-1));
              Sortie(1);          
            };
          // on regarde si le coefficient  est thermo dépendante
          indicateur="C";   indicateur += ChangeEntierSTring(j);indicateur +=ChangeEntierSTring(i-j);
          indicateur +="_thermo_dependant_";
          if(strstr(entreePrinc->tablcar,indicateur.c_str())!=0)
           { thermo_dependant=true;
             *(entreePrinc->entree) >> nom;
             if (nom != indicateur) 
              { cout << "\n erreur en lecture de la thermodependance , on aurait du lire le mot cle " << indicateur 
                     << " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
                entreePrinc->MessageBuffer("**erreur4 Poly_hyper3D::LectureDonneesParticulieres (...**");
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);     
               };
             // lecture de la loi d'évolution en fonction de la température
             *(entreePrinc->entree) >>  nom;    
             // on regarde si la courbe existe, si oui on récupère la référence 
             if (lesCourbes1D.Existe(nom)) 
              { Cij_temperature(i)(j+1) = lesCourbes1D.Trouve(nom);}
             else
              { // sinon il faut la lire maintenant
                string non_courbe("_");   
                Cij_temperature(i)(j+1) = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
                // lecture de la courbe
                Cij_temperature(i)(j+1)->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
               }; 
             // prepa du flot de lecture
             entreePrinc->NouvelleDonnee(); 
           }// fin du cas où il y a une thermo dépendance
          else  // cas sans thermodépendance
            // lecture du coef Cij
           { *(entreePrinc->entree) >> Cij(i)(j+1) ;};
      	}; //-- fin de la boucle sur j
              
    // ---- lecture du coefficient K (dernier coefficient obligatoire) -----
    *(entreePrinc->entree) >> nom;   
    if(nom != "K=")
         { cout << "\n erreur en lecture du parametre K, on attendait la chaine K= et on a lu: " << nom;
           cout << "\n Poly_hyper3D::LectureDonneesParticulieres  "
                << "(UtilLecture * entreePrinc) " << endl ;
           entreePrinc->MessageBuffer("erreur 5   ");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);          
          }
    // on regarde si le coefficient  est thermo dépendante
    if(strstr(entreePrinc->tablcar,"K_thermo_dependant_")!=0)
     { thermo_dependant=true;
       *(entreePrinc->entree) >> nom;
       if (nom != "K_thermo_dependant_") 
         { cout << "\n erreur en lecture de la thermodependance de K, on aurait du lire le mot cle K_thermo_dependant_"
                << " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
           entreePrinc->MessageBuffer("**erreur6 Poly_hyper3D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
          }
       // lecture de la loi d'évolution en fonction de la température
       *(entreePrinc->entree) >>  nom;    
       // on regarde si la courbe existe, si oui on récupère la référence 
       if (lesCourbes1D.Existe(nom)) 
         { K_temperature = lesCourbes1D.Trouve(nom);
          }
       else
        { // sinon il faut la lire maintenant
          string non_courbe("_");   
          K_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
          // lecture de la courbe
          K_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
         } 
       // prepa du flot de lecture
       entreePrinc->NouvelleDonnee(); 
      }
    else
     { // lecture de K 
       *(entreePrinc->entree) >> K ;
       // s'il n'y a plus rien n'a lire, il faut passer un enregistrement
       if((strstr(entreePrinc->tablcar,"type_potvol_")==0)
          && (strstr(entreePrinc->tablcar,"avec_courbure_")==0)
          && (strstr(entreePrinc->tablcar,"avec_nD_")==0)
          && (strstr(entreePrinc->tablcar,"sortie_post_")==0)
         )
          entreePrinc->NouvelleDonnee();
      };
    // lecture éventuelle du type de potentiel  
    if(strstr(entreePrinc->tablcar,"type_potvol_")!=0)
     { *(entreePrinc->entree) >> nom;
       if (nom != "type_potvol_") 
         { cout << "\n erreur en lecture du type de variation volumique, on aurait du lire le mot cle type_potvol_"
                << " suivi d'un nombre entier ";
           entreePrinc->MessageBuffer("**erreur7 Poly_hyper3D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
          }
       // lecture du type de variation volumique
       *(entreePrinc->entree) >>  type_pot_vol;    
       // on regarde si ce type est possible
       switch (type_pot_vol) 
       	{case 1: case 2: case 3: case 4 : break; // ok
       	 default:
         { cout << "\n erreur en lecture du type de variation volumique: valeur lue: " << type_pot_vol 
                << " , actuellement seule les types 1, 2 ,3 et 4 sont implantes ";
           entreePrinc->MessageBuffer("**erreur8 Poly_hyper3D::LectureDonneesParticulieres (...**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
          };
       	};
       // s'il n'y a plus rien n'a lire, il faut passer un enregistrement
       if((strstr(entreePrinc->tablcar,"avec_nD_")==0)
         && (strstr(entreePrinc->tablcar,"avec_courbure_")==0)
         && (strstr(entreePrinc->tablcar,"sortie_post_")==0)
        )
           entreePrinc->NouvelleDonnee();
     };
                
    
     // ---- lecture éventuelle d'un terme de courbure
     if(strstr(entreePrinc->tablcar,"avec_courbure_")!=0)
      { *(entreePrinc->entree) >> nom;
        if (nom != "avec_courbure_")
          { cout << "\n erreur en lecture du mot cle, on aurait du lire le mot cle avec_courbure_"
                 << " et on a lue: " << nom ;
            entreePrinc->MessageBuffer("**erreur9 Poly_hyper3D::LectureDonneesParticulieres (...**");
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);
          };
        *(entreePrinc->entree) >> nom;
        if (nom != "a_courbure=")
          { cout << "\n erreur en lecture du parametre a, on aurait du lire le mot cle a_courbure="
                 << " suivi d'un nombre reel et on a lue: " << nom ;
            entreePrinc->MessageBuffer("**erreur10 Poly_hyper3D::LectureDonneesParticulieres (...**");
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);
          };
        avec_courbure=true;
        // on regarde si le coefficient  est thermo dépendante
        if(strstr(entreePrinc->tablcar,"a_thermo_dependant_")!=0)
         { thermo_dependant=true;
           *(entreePrinc->entree) >> nom;
           if (nom != "a_thermo_dependant_")
             { cout << "\n erreur en lecture de la thermodependance de a, on aurait du lire le mot cle a_thermo_dependant_"
                    << " suivi du nom d'une courbe de charge ou de la courbe elle meme, et on a lue: " << nom ;
               entreePrinc->MessageBuffer("**erreur10 Poly_hyper3D::LectureDonneesParticulieres (...**");
               throw (UtilLecture::ErrNouvelleDonnee(-1));
               Sortie(1);
             };
           // lecture de la loi d'évolution en fonction de la température
           *(entreePrinc->entree) >>  nom;
           // on regarde si la courbe existe, si oui on récupère la référence
           if (lesCourbes1D.Existe(nom))
             { a_temperature = lesCourbes1D.Trouve(nom);}
           else
             { // sinon il faut la lire maintenant
               string non_courbe("_");
               a_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
               // lecture de la courbe
               a_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
              };
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
         }
        else
         { // lecture de a
           *(entreePrinc->entree) >> a_courbure ;
         };
        //-- on lit maintenant le paramètre r
        *(entreePrinc->entree) >> nom;
        if (nom != "r_courbure=")
          { cout << "\n erreur en lecture du parametre r, on aurait du lire le mot cle r_courbure="
                 << " suivi d'un nombre reel et on a lue: " << nom ;
            entreePrinc->MessageBuffer("**erreur11 Poly_hyper3D::LectureDonneesParticulieres (...**");
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);
          };
        // on regarde si le coefficient  est thermo dépendante
        if(strstr(entreePrinc->tablcar,"r_thermo_dependant_")!=0)
         { thermo_dependant=true;
           *(entreePrinc->entree) >> nom;
           if (nom != "r_thermo_dependant_")
             { cout << "\n erreur en lecture de la thermodependance de r, on aurait du lire le mot cle r_thermo_dependant_"
                    << " suivi du nom d'une courbe de charge ou de la courbe elle meme, et on a lue: " << nom ;
               entreePrinc->MessageBuffer("**erreur12 Poly_hyper3D::LectureDonneesParticulieres (...**");
               throw (UtilLecture::ErrNouvelleDonnee(-1));
               Sortie(1);
             };
           // lecture de la loi d'évolution en fonction de la température
           *(entreePrinc->entree) >>  nom;
           // on regarde si la courbe existe, si oui on récupère la référence
           if (lesCourbes1D.Existe(nom))
             { r_temperature = lesCourbes1D.Trouve(nom);}
           else
             { // sinon il faut la lire maintenant
               string non_courbe("_");
               r_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
               // lecture de la courbe
               r_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
              }
           // prepa du flot de lecture
           entreePrinc->NouvelleDonnee();
         }
        else
         { // lecture de a
           *(entreePrinc->entree) >> r_courbure ;
         };
        // prepa du flot de lecture
        entreePrinc->NouvelleDonnee();
      } //-- fin de la lecture éventuelle d'un terme de courbure
     else
      {avec_courbure=false;};
                
    // -------- on regarde s'il y a une dépendance à une fonction nD
    if (strstr(entreePrinc->tablcar,"avec_nD_")!=NULL)
    { string nom1,nom2;
      string nom_class_methode = "Poly_hyper3D::LectureDonneesParticulieres";
      // lecture des  fonctions nD,
      entreePrinc->NouvelleDonnee();  // passage d'une ligne
      
      // on commence par dimensionner le tableau des pointeurs de fonctions nD
      int tai = Cij.Taille(); Cij_nD.Change_taille(tai);
      for (int i=1;i<= tai; i++)
        Cij_nD(i).Change_taille( Cij(i).Taille(),NULL);
      
      // on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_avec_nD_"
      // ou un nouveau mot clé global auquel cas il y a pb !!
      MotCle motCle; // ref aux mots cle
      while (strstr(entreePrinc->tablcar,"fin_parametres_avec_nD_")==0)
      {
        // si on a  un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur
        if ( motCle.SimotCle(entreePrinc->tablcar))
         { cout << "\n erreur de lecture des parametre de dependance a une fonction nD: on n'a pas trouve le mot cle "
                << " fin_parametres_avec_nD_ et par contre la ligne courante contient un mot cle global  ";
           entreePrinc->MessageBuffer
                ("** erreur51 des parametres de reglage de la loi de comportement Poly_hyper3D**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
        
        // lecture d'un mot clé
        *(entreePrinc->entree) >> nom1;
        
        if ((entreePrinc->entree)->rdstate() == 0)
          {} // lecture normale
        #ifdef ENLINUX
        else  if ((entreePrinc->entree)->fail())
           // on a atteind la fin de la ligne et on appelle un nouvel enregistrement
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
              *(entreePrinc->entree) >>nom1;
           }
        #else
        else  if ((entreePrinc->entree)->eof())
          // la lecture est bonne mais on a atteind la fin de la ligne
          { if(nom1 != "fin_parametres_avec_nD_")
              {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom1;};
          }
        #endif
        else // cas d'une erreur de lecture
         { cout << "\n erreur de lecture inconnue  ";
           entreePrinc->MessageBuffer
             ("** erreur61 des parametres de reglage de la loi de comportement Maheo_hyper**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
       
        // on récupère les indices, on doit avoir la forme suivante C_IJ_nD=
        // ou alors il s'agit de K_nD
        if (   (nom1 != "K_nD=")
             &&
               (   (nom1[0] != 'C') || (nom1[1] != '_' ) || (nom1[4] != '_' )
                || (nom1[5] != 'n' ) || (nom1[6] != 'D' ) || (nom1[7] != '=' )
                || isNumeric(&nom1[2],10) || isNumeric(&nom1[3],10)
               )
            )
          { cout << "\n erreur en lecture, l'indicateur lue de la fonction nD = " << nom1
                 << " n'est pas correct, il devrait etre de la forme  C_IJ_nD= avec I et J un entier  "
                 << " entre 0 et 9 "
                 << " il n'est donc pas utilisable !! ";
            string message("\n**erreur011** \n"+nom_class_methode+"(...");
            entreePrinc->MessageBuffer(message);
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);
          };
          
        // deux cas suivant qu'il s'agit d'un coef ou de K
        if(nom1 == "K_nD=")
         {string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {K_nD = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             K_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             K_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (K_nD->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << K_nD->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = K_nD->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         }
        else // sinon il s'agit d'un coef
         {int i = ChangeEntier(&nom1[2]);
          int j = ChangeEntier(&nom1[3]);
          string nom_fonct;
          // lecture du nom de la fonction
          *(entreePrinc->entree) >>  nom_fonct;
          
          // maintenant on définit la fonction
          if (lesFonctionsnD.Existe(nom_fonct))
            {Cij_nD(i)(j) = lesFonctionsnD.Trouve(nom_fonct);
            }
          else
            {// sinon il faut la lire maintenant
             string non("_");
             Cij_nD(i)(j) = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
             // lecture de la fonction
             Cij_nD(i)(j)->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
             // maintenant on vérifie que la fonction est utilisable
             if (Cij_nD(i)(j)->NbComposante() != 1 )
              { cout << "\n erreur en lecture, la fonction " << nom_fonct
                     << " est une fonction vectorielle a   " << Cij_nD(i)(j)->NbComposante()
                     << " composante alors qu'elle devrait etre scalaire ! "
                     << " elle n'est donc pas utilisable !! ";
                string message("\n**erreur031** \n"+nom_class_methode+"(...");
                entreePrinc->MessageBuffer(message);
                throw (UtilLecture::ErrNouvelleDonnee(-1));
                Sortie(1);
              };
            };
          // on regarde si la fonction nD intègre la température
          const Tableau <Ddl_enum_etendu>& tab_enu = Cij_nD(i)(j)->Tab_enu_etendu();
          if (tab_enu.Contient(TEMP))
            thermo_dependant=true;
          // prepa du flot de lecture
          entreePrinc->NouvelleDonnee();
         };
       }; //-- fin du while
    };
    
    // on définit le nombre de paramètre total de la loi
    tai = Cij.Taille();
    nb_para_loi = 1; // le K
    for (int i=1;i<= tai; i++)
       { int baj = Cij(i).Taille();
         for (int j=1;j<=baj;j++)
           nb_para_loi++;
        };
        
    // lecture de l'indication éventuelle du post traitement
    string nom_class_methode = "Poly_hyper3D::LectureDonneesParticulieres";string le_mot_cle = "sortie_post_";
    if (entreePrinc->Lecture_un_parametre_int(0,nom_class_methode,0,1,le_mot_cle,sortie_post))
       entreePrinc->NouvelleDonnee(); // passage éventuel du mot clé sortie_post_;

    // --- appel au niveau de la classe mère
    // ici  il n'y a pas de type de déformation associé
    // mais on prend la def standart d'almansi, pour les fonctions associées éventuelles
    Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
                 (*entreePrinc,lesFonctionsnD,false);
   };

 // affichage de la loi
void Poly_hyper3D::Affiche()  const 
  { cout << " \n loi de comportement hyper elastique 3D de type polynomiale en invariants "; 
    // affichage du degré maxi
    int tai = Cij.Taille();
    cout << " degre_maxi_polynome= " <<  tai << " ";
    // on balaie les coefficients
    for (int i=1;i<= tai; i++)
      for (int j=1;j<= i+1; j++)
      	{ // constitution de l'indicateur
      	  string indicateur("C"); indicateur += ChangeEntierSTring(j);indicateur +=ChangeEntierSTring(i-j);
      	   indicateur +="=";
      	  cout << indicateur << " ";
          if ( Cij_temperature(i)(j) != NULL) { cout << " thermo dependant " 
                                         << " courbe f(T): " << Cij_temperature(i)(j)->NomCourbe() <<" ";}
          else                          { cout <<  Cij(i)(j) ;};
      	};
    if ( K_temperature != NULL) { cout << " K thermo dependant " 
                                         << " courbe K=f(T): " << K_temperature->NomCourbe() <<" ";}
    else                          { cout << " K= " << K << " ";}; 
    cout << " type de potentiel= " << type_pot_vol << " ";
    if (avec_courbure)
      { if ( a_temperature != NULL) { cout << " a thermo dependant " 
                                             << " courbe a=f(T): " << a_temperature->NomCourbe() <<" ";}
          else                          { cout << " a= " << a_courbure ;} 
        if ( r_temperature != NULL) { cout << " r thermo dependant " 
                                             << " courbe r=f(T): " << r_temperature->NomCourbe() <<" ";}
          else                          { cout << " r= " << r_courbure ;} 
      };

    // les dépendances à une fonction nD
    if (K_nD != NULL)
     {cout << " fonction nD K_nD: ";
      if (K_nD->NomFonction() != "_")
           cout << K_nD->NomFonction();
      else
           K_nD->Affiche();
     };
    int bai = Cij_nD.Taille();
    for (int i=1;i<= bai; i++)
       { int baj = Cij_nD(i).Taille();
         for (int j=1;j<=baj;j++)
            if (Cij_nD(i)(j) != NULL)
              {cout << " fonction nD C_"<<i<<"_"<<j<<" ";
               if (Cij_nD(i)(j)->NomFonction() != "_")
                    cout << Cij_nD(i)(j)->NomFonction();
               else
                    Cij_nD(i)(j)->Affiche();
              };
       };
    // appel de la classe mère
    Loi_comp_abstraite::Affiche_don_classe_abstraite();
  };
            
// affichage et definition interactive des commandes particulières à chaques lois
void Poly_hyper3D::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
 {  ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
	   cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? ";
	   string rep = "_";
    // procédure de lecture avec prise en charge d'un retour chariot
    rep = lect_return_defaut(true,"o");

    sort << "\n# ----------------------------------------------------------------------------------"
         << "\n# |... loi de comportement hyper elastique 3D de type polynomiale en invariants .. |"
         << "\n# |                   .. coefficients de type Cij  ..                              |"
         << "\n# ----------------------------------------------------------------------------------"
         << "\n\n#           exemple de definition de loi"
         << "\n   degre_maxi= 2 C01= 0.0167 C10= 0.145 C02= 0.0167 C11= 0.145 C20= 0.0167  K= 100   "
         << "\n#   .. fin de la definition de la loi polynomiale \n" << endl;
	   if ((rep != "o") && (rep != "O" ) && (rep != "0") )
      { sort << "\n#   " 
         << "\n#  on note que l'on met tout d'abord le groupe des coeff de degre 1, puis le groupe de degre 2  " 
         << "\n#  tous les coefficients doivent etre presents " 
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  il est possible de definir des parametres thermo-dependants (1 ou 2 ou 3 parametres)"
         << "\n#  par exemple pour les trois parametres on ecrit: "
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  C01= C01_thermo_dependant_ courbe1 "
         << "\n#  C10= C10_thermo_dependant_ courbe2 "
         << "\n#  C02= C02_thermo_dependant_ courbe2 "
         << "\n#  C11= C10_thermo_dependant_ courbe2 "
         << "\n#  C20= 0.0167  K=  K_thermo_dependant_ courbe4 "
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  noter qu'apres la definition de chaque courbe, on change de ligne, a l'inverse "
         << "\n#  si la valeur du parametre est fixe, on poursuit sur la meme ligne. "
         << "\n# "
         << "\n#---------------------------------- type_potvol_ ----------------------------------------"
         << "\n#  un  parametre facultatif permet d'indiquer le type de variation volumique "
         << "\n#  que l'on desire: par defaut il s'agit de : K(1-(1+log(V))/V) qui correspond au type 1"
         << "\n#   mais on peut choisir:                     K/2(V-1)          qui correspond au type 2 "
         << "\n#                     ou:                     K/2(log(V))^2     qui correspond au type 3 "
         << "\n#                     ou:                     K/2(V-1)^2          qui correspond au type 4 "
         << "\n#   en indiquant (en fin de ligne) le mot cle: type_potvol_ suivi du type"
         << "\n#  par defaut type_potvol_ a la valeur 1 "
         << "\n# "
         << "\n#---------------------------------- avec_courbure_ -------------------------------------"
         << "\n#  apres le type de variation volumique on peut indiquer facultativement l'ajout au potentiel  "
         << "\n#  d'un terme permettant de raidir le comportement a partir d'un certain niveau de chargement "
         << "\n#  pour cela on indique le mot cle: avec_courbure_ puis on change de ligne "
         << "\n#  le terme additionnelle depend de deux parametres: a et r, a positionne la valeur de J1 "
         << "\n#  a partir de laquelle il y a durcissement, r controle la courbure du changement de regime "
         << "\n#  exemple de declaration :  "
         << "\n#   ....  type_potvol_ 2  avec_courbure_ "
         << "\n#      a_courbure= 94 r_courbure= 1.24 "
         << "\n#  ces deux parametres peuvent etre l'un et/ou l'autre dependant de la temperature  "
         << "\n#  dans ce cas la declaration de dependance suit les regles habituelles "
         << "\n#  exemple de declaration dans le cas d'une dependance a la temperature :  "
         << "\n#   ...  type_potvol_ 2  avec_courbure_ "
         << "\n#      a_courbure= a_thermo_dependant_  "
         << "\n#      r_courbure= r_thermo_dependant_  "
         << "\n# "
         << "\n#--------------------- dependance explicite a une fonction nD -----------------------------------"
         << "\n# il est egalement possible de definir une dependance  des parametres a une fonction nD "
         << "\n#  dans ce cas, a la suite de la dependance explicite a la temperature via une courbr 1D "
         << "\n#  ou a la suite du dernier parametre fixe, on met le mot cle: avec_nD_ "
         << "\n#  puis sur les lignes qui suivent on met des fonctions nD multiplicatives des parametres "
         << "\n#  initiaux, puis un mot cle signalant la fin des fonctions: fin_parametres_avec_nD_ "
         << "\n# suit un exemple de declaration (il est possible d'ommettre certaine fonction ! "
         << "\n# la presence de fct nD est facultative), de meme il est possible soit d'utiliser une fonction "
         << "\n# deja defini soit de donner directement la fonction nD (comme pour les autres lois) "
         << "\n# Enfin, apres chaque definition de fonction  il faut passer a une nouvelle ligne "
         << "\n# "
         << "\n#           exemple de definition de loi"
         << "\n# degre_maxi= 2 C01= 0.0167 C10= 0.145 C02= 0.0167 C11= 0.145 C20= 0.0167  K= 100   avec_nD_ "
         << "\n#   K_nD= fct_nD_1 "
         << "\n#   C_01_nD= fct_nD_2 "
         << "\n#   C_11_nD= fct_nD_3 "
         << "\n# fin_parametres_avec_nD_ "
         << "\n# "
         << "\n# Remarques: "
         << "\n# - les fonction nD sont des fonctions multiplicatives  "
         << "\n# des parametres initiaux (ou de ceux obtenus apres dependance explicite a la temperature  "
         << "\n#  ( bien noter que la loi obtenue peut-etre quelconque, en particulier plus du tout "
         << "\n#  hyperelastique, tout depend des choix des fcts nD !) "
         << "\n# "
         << "\n# "
         << "\n#----------------------------------- sortie_post_ -------------------------------------"
         << "\n#  il est possible de recuperer differentes grandeurs de travail par exemple "
         << "\n#  l'intensite du potentiel, comme ces grandeurs sont calculees au moment de la resolution "
         << "\n#  et ne sont pas stockees. On peut egalement recuperer la valeur des coefficients de la loi "
         << "\n#  (utile s'ils varient) sous forme d'un tableau de n uplets dont le premier est K suivi des CIJ "
         << "\n#  dont l'ordre est celui d'entree des parametres CIJ a la lecture "
         << "\n#  Il faut indiquer le mot cle sortie_post_ suivi de 1 (par defaut = 0) "
         << "\n#  ensuite au moment de la constitution du .CVisu on aura acces aux grandeurs de travail "
         << "\n#  ex: "
         << "\n#    sortie_post_  1 "
         << "\n#  ce mot cle est le dernier des parametres specifiques de la loi il doit se situe "
         << "\n#  a la fin de la derniere ligne de donnees "
         << "\n#"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#------------------------------------------------------------------------------------"
         << endl;
       };
    // appel de la classe Hyper_W_gene_3D
    Hyper_W_gene_3D::Info_commande_LoisDeComp_hyper3D(entreePrinc);
    // appel de la classe mère
    Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);     
  };  		  	  

// test si la loi est complete
int Poly_hyper3D::TestComplet()
    { int ret = LoiAbstraiteGeneral::TestComplet();
      int tai = Cij.Taille();
      string indicateur;
      for (int i=1;i<= tai; i++)
        for (int j=0;j<= i; j++)
          { indicateur="C"; indicateur += ChangeEntierSTring(j);indicateur +=ChangeEntierSTring(i-j); 
            if ((Cij_temperature(i)(j+1) == NULL) && (Cij(i)(j+1) == (-ConstMath::trespetit)))
             { cout << " \n le coefficient " << indicateur << " de la loi Poly_hyper3D n'est pas defini pour la loi " 
                    << Nom_comp(id_comp)
                    << '\n' << endl;
               ret = 0;
             };                
          };
      if (avec_courbure)
      { if ((a_temperature == NULL) && (a_courbure == (-ConstMath::trespetit)))
              { cout << " \n le coefficient de courbure a (loi de Hart Smith + raidissement) n'est pas defini pour la loi " 
                     << Nom_comp(id_comp)
                     << '\n' << endl;
        ret = 0;
       };
        if ((r_temperature == NULL) && (r_courbure == (-ConstMath::trespetit)))
              { cout << " \n le coefficient de courbure r (loi de Hart Smith + raidissement) n'est pas defini pour la loi " 
                     << Nom_comp(id_comp)
                     << '\n' << endl;
        ret = 0;
               };                
      };               
       //
       if (ret == 0)
        {this-> Affiche();
         ret = 0;
        };
       return ret; 
    }; 

	   //----- lecture écriture de restart -----
	   // cas donne le niveau de la récupération
       // = 1 : on récupère tout
       // = 2 : on récupère uniquement les données variables (supposées comme telles)
void Poly_hyper3D::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
                                             ,LesFonctions_nD& lesFonctionsnD)
  { string nom; bool test;
    if (cas == 1)
     { // le degré maxi
       int tai=0;
       ent >> nom >> tai;
       Cij.Change_taille(tai); // on dimensionne correctement la première dimension de Cij
       Cij_use.Change_taille(tai);
       // pour les températures on commence par vider  le tableau
       int tait = Cij_temperature.Taille();
       for (int i=1;i<=tait;i++)
        {int taj = Cij_temperature(i).Taille();
         for (int j=1;j<=taj;j++)
           if (Cij_temperature(i)(j) != NULL)
             if (Cij_temperature(i)(j)->NomCourbe() == "_") delete Cij_temperature(i)(j);
        };
       // maintenant on peut adapter la dimension
       Cij_temperature.Change_taille(tai); // idem pour la dépendance éventuelle à la température
       // puis on dimensionne la deuxième dimension
       for (int i=1;i<= tai; i++)
        {Cij(i).Change_taille(i+1,(-ConstMath::trespetit)); // init à (-ConstMath::trespetit)
         Cij_use(i).Change_taille(i+1,(-ConstMath::trespetit));
         Cij_temperature(i).Change_taille(i+1,NULL); // idem init à NULL
         // on lit les informations
         for (int j=1;j<= i+1; j++)
          {ent >> nom >> test;
           if (!test)
             { ent >> Cij(i)(j);}
           else
             { ent >> nom;
               Cij_temperature(i)(j) = lesCourbes1D.Lecture_pour_base_info(ent,cas,Cij_temperature(i)(j));
             };
          };
        }; //-- fin de la boucle i
     
       // K
       ent >> nom >> test;
       // vérification
       #ifdef MISE_AU_POINT
       if (nom != "K")
         { cout << "\n erreur en lecture du parametres K de la loi de mooney rivlin 3D"
                << " on devait lire K= avant le second parametre "
                << " et on a lue: " << nom << " "
                << "\n Poly_hyper3D::Lecture_base_info_loi(..."
                << endl;
           Sortie(1);
         };
       #endif
       if (!test)
         { ent >> K;
           if (K_temperature != NULL) {if (K_temperature->NomCourbe() == "_")
              delete K_temperature; K_temperature = NULL;};
         }
       else
         { ent >> nom; K_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,K_temperature); };
       // le type de potentiel
       ent >>  nom  >> type_pot_vol ;
       // --- potentiel de gestion de courbure, éventuellement
       ent >>  nom   ;
       #ifdef MISE_AU_POINT
       if (nom != "avec_courbure=")
         { cout << "\n erreur en lecture de l'indicateur de courbure de la loi de Poly_hyper3D "
                << " on devait lire avec_courbure= et on a lue: " << nom << " "
                << "\n Poly_hyper3D::Lecture_base_info_loi(..."
                << endl;
           Sortie(1);
         };
       #endif
       ent >> avec_courbure;
       if (avec_courbure)
        {// a_courbure
         ent >> nom >> test;
         // vérification
         #ifdef MISE_AU_POINT
         if (nom != "a=")
          { cout << "\n erreur en lecture du parametres a de la gestion eventuelle de courbure de la loi de Poly_hyper3D  3D"
                 << " on devait lire a= avant le premier parametre "
                 << " et on a lue: " << nom << " "
                 << "\n Poly_hyper3D::Lecture_base_info_loi(..."
                 << endl;
            Sortie(1);
          };
         #endif
         if (!test)
          { ent >> a_courbure;
            if (a_temperature != NULL) {if (a_temperature->NomCourbe() == "_")
               delete a_temperature; a_temperature = NULL;};
          }
         else
          { ent >> nom; a_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,a_temperature); };
         // r_courbure
         ent >> nom >> test;
         // vérification
         #ifdef MISE_AU_POINT
         if (nom != "r=")
          { cout << "\n erreur en lecture du parametres r de la gestion eventuelle de courbure de la loi de Poly_hyper3D "
                 << " on devait lire r= avant le premier parametre "
                 << " et on a lue: " << nom << " "
                 << "\n Poly_hyper3D::Lecture_base_info_loi(..."
                 << endl;
            Sortie(1);
          };
         #endif
         if (!test)
           {ent >> r_courbure;
            if (r_temperature != NULL) {if (r_temperature->NomCourbe() == "_")
               delete r_temperature; r_temperature = NULL;};
           }
         else
          { ent >> nom; r_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,r_temperature); };
        };
       // --- fonctions nD éventuelles
       ent >> nom;
       if (nom == "avec_fct_nD")
        { ent >> nom;
          // d'abord K
          if (nom == "K_nD:")
            K_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,K_nD);
          // puis les coefficients
          int bai = Cij_nD.Taille();
          // comme il s'agit du cas==1, on considère que ce ne sera pas exécuté souvent
          // on ne cherche pas a optimiser, du coup on commence par supprimer les fct éventuellement existantes
          // ceci pour faciliter la lecture
          for (int i=1;i<= bai; i++)
            { int baj = Cij_nD(i).Taille();
              for (int j=1;j<=baj;j++)
                 if (Cij_nD(i)(j) != NULL) if (Cij_nD(i)(j)->NomFonction() == "_")
                    {delete Cij_nD(i)(j);Cij_nD(i)(j)=NULL;};
            };
          // on dimensionne si nécessaire
          if (!bai)
            // cas où le tableau est vide
            {int tai = Cij.Taille(); Cij_nD.Change_taille(tai);
             for (int i=1;i<= tai; i++)
               Cij_nD(i).Change_taille( Cij(i).Taille(),NULL);
            };
          // on lit et on rempli
          for (int i=1;i<= bai; i++)
            { int baj = Cij_nD(i).Taille();
              for (int j=1;j<=baj;j++)
               { int a,b;
                 ent >> nom >> a >> b >> nom;
                 Cij_nD(i)(j) = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Cij_nD(i)(j));
               };

           };
        };
        
        // on définit le nombre de paramètre total de la loi
        tai = Cij.Taille();
        nb_para_loi = 1; // le K
        for (int i=1;i<= tai; i++)
           { int baj = Cij(i).Taille();
             for (int j=1;j<=baj;j++)
               nb_para_loi++;
            };

     }
	   // appel de la classe mère
	   Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
 };
       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Poly_hyper3D::Ecriture_base_info_loi(ofstream& sort,const int cas)
  { if (cas == 1) 
     { // le degré maxi 
       int tai = Cij.Taille();
       sort << "\n degre_maxi_polynome= " <<  tai << " ";
     
       // on balaie les coefficients
       sort << "\n ";
       for (int i=1;i<= tai; i++)
        for (int j=0;j<= i; j++)
      	  { // constitution de l'indicateur
      	    string indicateur("C"); indicateur += ChangeEntierSTring(j);indicateur +=ChangeEntierSTring(i-j);
           indicateur +="=";
      	    sort << indicateur << " ";
           if ( Cij_temperature(i)(j+1) != NULL)
                { sort << true << "  courbe_f(T): ";
                  LesCourbes1D::Ecriture_pour_base_info(sort,cas,Cij_temperature(i)(j+1));
                }
           else { sort << false <<  Cij(i)(j+1) ;};
      	  };
       // la partie volumique
	      sort << "  K= ";
       if (K_temperature == NULL)
         { sort << false << " " << K << " ";}
       else 
         { sort << true << " fonction_K_temperature ";
           LesCourbes1D::Ecriture_pour_base_info(sort,cas,K_temperature);
         };
       // le potentiel   
       sort << " typ_pot= " << type_pot_vol << " "; // le type de potentiel   
       // --- potentiel de gestion de courbure, éventuellement
       sort << " avec_courbure= " << avec_courbure;
       if (avec_courbure)
        { sort << "  a= ";             
          if (a_temperature == NULL)
            { sort << false << " " << a_courbure << " ";}
          else
            { sort << true << " fonction_a_temperature ";
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,a_temperature);
            };
          sort << "\n  r= ";             
          if (r_temperature == NULL)
            { sort << false << " " << r_courbure << " ";}
          else
            { sort << true << " fonction_r_temperature ";
              LesCourbes1D::Ecriture_pour_base_info(sort,cas,r_temperature);
            };
        };
       // --- cas de fonction nD éventuelles
       int bai = Cij_nD.Taille();
       if (!bai)
        { sort << " sans_fctnD ";}
       else
        { sort << " avec_fct_nD ";
          // d'abord K
          if (K_nD != NULL)
             {sort << " K_nD: " << " ";
              LesFonctions_nD::Ecriture_pour_base_info(sort, cas,K_nD);
             }
          else {sort << " non_K_nD "; };
          // puis les coefficients
          for (int i=1;i<= bai; i++)
            { int baj = Cij_nD(i).Taille();
              for (int j=1;j<=baj;j++)
               if (Cij_nD(i)(j) != NULL)
                 { sort << " coeff " << i <<" "<<j<<" fctnD ";
                   LesFonctions_nD::Ecriture_pour_base_info(sort, cas,Cij_nD(i)(j));
                 }
               else
                 { sort << " coeff " << i <<" "<<j<<" sans_fctnD "; };
            };
        };
     }
	   // appel de la classe mère
    Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);  
  };
     

 // ========== codage des METHODES VIRTUELLES  protegees:================
        // calcul des contraintes a t+dt
void Poly_hyper3D::Calcul_SigmaHH(TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl,
             TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB & epsBB_,
             TenseurBB & ,
             TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,
             double& jacobien_0,double& jacobien,TenseurHH & sigHH_
		  	,EnergieMeca & energ,const EnergieMeca & ,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Expli_t_tdt& ex )
 {   
    #ifdef MISE_AU_POINT
    if (epsBB_.Dimension() != 3)
       { cout << "\nErreur : la dimension devrait etre 3 !\n";
         cout << " Poly_hyper3D::Calcul_SigmaHH\n";
         Sortie(1);
       };
    if (tab_ddl.NbDdl() != d_gijBB_.Taille())
       { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_  !\n";
         cout << " Poly_hyper3D::Calcul_SigmaHH\n";
         Sortie(1);
       };
    #endif
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_);  //   passage explicite en tenseur dim 3
    
    // init des coefficients utilisés
    Cij_use = Cij;
    K_use = K;
    
    // informations éventuelles
    #ifdef MISE_AU_POINT
    if  (Permet_affichage() > 4)
      { cout << "\nPoly_hyper3D::Calcul_SigmaHH(... ";
      };
    #endif
    // cas de la thermo dépendance, on calcul les grandeurs
    int tai = Cij.Taille();

    // on balaie les coefficients
    for (int i=1;i<= tai; i++)
      for (int j=1;j<= i+1; j++)
      	 if (Cij_temperature(i)(j) != NULL) 
      	     {Cij_use(i)(j) = Cij_temperature(i)(j)->Valeur(*temperature);
             #ifdef MISE_AU_POINT
             if  (Permet_affichage() > 4)
                cout << "\n C_"<<i<<"_"<<j<<"(T)= " << Cij_use(i)(j) << " ";
             #endif
            };
    if (K_temperature != NULL)
      {K_use = K_temperature->Valeur(*temperature);
       #ifdef MISE_AU_POINT
       if  (Permet_affichage() > 4)
          cout << "\n K_(T)="<<K_use << " ";
       #endif
      };
    if (avec_courbure)
     { if (a_temperature != NULL)
         {a_courbure = a_temperature->Valeur(*temperature);
          #ifdef MISE_AU_POINT
          if  (Permet_affichage() > 4)
             cout << "\n a_courbure(T)="<<a_courbure << " ";
          #endif
         };
       if (r_temperature != NULL)
         {r_courbure = r_temperature->Valeur(*temperature);
          #ifdef MISE_AU_POINT
          if  (Permet_affichage() > 4)
             cout << "\n r_courbure(T)="<<r_courbure << " ";
          #endif
         };
     };
  
    // récup du conteneur spécifique du point, pour sauvegarde éventuelle
    SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
    // dans le cas de l'utilisation de fct nD
    // --- cas de fonction nD éventuelles
    int bai = Cij_nD.Taille();
    if ((!bai) || (K_nD != NULL)) // si diff de 0
      { // opération de transmission de la métrique
        const Met_abstraite::Impli* ex_impli = NULL;
        const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &ex;
        const Met_abstraite::Umat_cont* ex_expli = NULL;
        // tout d'abord K
        if (K_nD != NULL)
          // là il faut calculer la fonction nD
          { // on utilise la méthode générique de loi abstraite
            list <SaveResul*> list_save; // inter pour l'appel de la fonction
            list_save.push_back(saveResul);
            Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                (K_nD,1 // une seule valeur attendue en retour
                 ,ex_impli,ex_expli_tdt,ex_expli
                 ,NULL
                 ,NULL
                 ,&list_save
                );
            // on récupère le premier élément du tableau uniquement
            K_use *= tab_val(1);
          };

        // puis les coef: on boucle sur le tableau
        for (int i=1;i<= bai; i++)
          { int baj = Cij_nD(i).Taille();
            for (int j=1;j<=baj;j++)
             if (Cij_nD(i)(j) != NULL)
                // là il faut calculer la fonction nD
              { // on utilise la méthode générique de loi abstraite
                list <SaveResul*> list_save; // inter pour l'appel de la fonction
                list_save.push_back(saveResul);
                Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                    (Cij_nD(i)(j),1 // une seule valeur attendue en retour
                     ,ex_impli,ex_expli_tdt,ex_expli
                     ,NULL
                     ,NULL
                     ,&list_save
                    );
                // on récupère le premier élément du tableau uniquement
                Cij_use(i)(j) *= tab_val(1);
              };
          };
      };
    if (sortie_post)
     {// sauvegarde des paramètres matériau
      int ia = 1; // init
      // d'abord K
      (*save_resulHyper_W.para_loi)(ia) = K_use;
      // puis les coef
      int tai = Cij.Taille();
      for (int i=1;i<= tai; i++)
        { int taj = Cij(i).Taille();
          for (int j=1;j<=taj;j++,ia++)
           (*save_resulHyper_W.para_loi)(ia) = Cij_use(i)(j);
        };
     };

    // calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D)
    Invariants_et_var1(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_,gijHH_,jacobien_0,jacobien);
    // calcul du potentiel et de ses dérivées premières / aux invariants J_r
    Potentiel_et_var(module_compressibilite);
    // stockage éventuel pour du post-traitement
    if (sortie_post)
     { // récup du conteneur spécifique du point, pour sauvegarde éventuelle
       SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
       save_resulHyper_W.invP->potentiel= W_d + W_v + W_c;
     };
  
    // calcul du tenseur des contraintes
    sigHH = ((W_d_J1+W_c_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2) 
               + ((W_v_J3+W_c_J3)/V) * d_J_r_epsBB_HH(3);
    #ifdef MISE_AU_POINT
    if  (Permet_affichage() > 4) cout << "\n sigHH " << sigHH;
    #endif

//    // calcul du module de compressibilité
//    module_compressibilite = 2. * V * (W_v_J3+W_c_J3) / (MaX(ConstMath::petit, log(V) ));
    // pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
    module_cisaillement = 0.;

    // traitement des énergies
    energ.Inita(0.); 
    energ.ChangeEnergieElastique((W_d+W_v)/V);
    
    LibereTenseur(); 
 };
        
        // calcul des contraintes a t+dt et de ses variations  
void Poly_hyper3D::Calcul_DsigmaHH_tdt (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl
            ,BaseB& ,TenseurBB & ,TenseurHH &
            ,BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> &
            ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,TenseurBB &
		          ,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt_
		          ,Tableau <TenseurBB *>& d_gijBB_tdt
		  	       ,Tableau <TenseurHH *>& ,double& jacobien_0,double& jacobien
		  	       ,Vecteur& ,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 (epsBB_tdt.Dimension() != 3)
       { cout << "\nErreur : la dimension devrait etre 3 !\n";
      cout << " Poly_hyper3D::Calcul_DsigmaHH_tdt\n";
      Sortie(1);
    };
    if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille())
       { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt  !\n";
      cout << " Poly_hyper3D::Calcul_SDsigmaHH_tdt\n";
      Sortie(1);
    };
    #endif
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt);  //  passage en dim 3 explicite
    Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_);  //  passage en dim 3 explicite
    
    // init des coefficients utilisés
    Cij_use = Cij;
    K_use = K;
    
    // informations éventuelles
    #ifdef MISE_AU_POINT
    if  (Permet_affichage() > 4)
      { cout << "\nPoly_hyper3D::Calcul_DsigmaHH_tdt(... ";
      };
    #endif
    // cas de la thermo dépendance, on calcul les grandeurs
    int tai = Cij.Taille();
  
    // on balaie les coefficients
    for (int i=1;i<= tai; i++)
      for (int j=1;j<= i+1; j++)
        if (Cij_temperature(i)(j) != NULL)
            {Cij_use(i)(j) = Cij_temperature(i)(j)->Valeur(*temperature);
             #ifdef MISE_AU_POINT
             if  (Permet_affichage() > 4)
                cout << "\n C_"<<i<<"_"<<j<<"(T)= " << Cij_use(i)(j) << " ";
             #endif
            };
    if (K_temperature != NULL)
      {K_use = K_temperature->Valeur(*temperature);
       #ifdef MISE_AU_POINT
       if  (Permet_affichage() > 4)
          cout << "\n K_(T)="<<K_use << " ";
       #endif
      };
    if (avec_courbure)
     { if (a_temperature != NULL)
         {a_courbure = a_temperature->Valeur(*temperature);
          #ifdef MISE_AU_POINT
          if  (Permet_affichage() > 4)
             cout << "\n a_courbure(T)="<<a_courbure << " ";
          #endif
         };
       if (r_temperature != NULL)
         {r_courbure = r_temperature->Valeur(*temperature);
          #ifdef MISE_AU_POINT
          if  (Permet_affichage() > 4)
             cout << "\n r_courbure(T)="<<r_courbure << " ";
          #endif
         };
     };
  
    // récup du conteneur spécifique du point, pour sauvegarde éventuelle
    SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
    // dans le cas de l'utilisation de fct nD
    int bai = Cij_nD.Taille();
    if ((!bai) || (K_nD != NULL)) // si diff de 0
      { // opération de transmission de la métrique
        const Met_abstraite::Impli* ex_impli = &ex;
        const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
        const Met_abstraite::Umat_cont* ex_expli = NULL;
        if (K_nD != NULL)
          // là il faut calculer la fonction nD
          { // on utilise la méthode générique de loi abstraite
            list <SaveResul*> list_save; // inter pour l'appel de la fonction
            list_save.push_back(saveResul);
            Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                (K_nD,1 // une seule valeur attendue en retour
                 ,ex_impli,ex_expli_tdt,ex_expli
                 ,NULL
                 ,NULL
                 ,&list_save
                );
            // on récupère le premier élément du tableau uniquement
            K_use *= tab_val(1);
          };
          
        // puis les coef: on boucle sur le tableau
        for (int i=1;i<= bai; i++)
          { int baj = Cij_nD(i).Taille();
            for (int j=1;j<=baj;j++)
             if (Cij_nD(i)(j) != NULL)
                // là il faut calculer la fonction nD
              { // on utilise la méthode générique de loi abstraite
                list <SaveResul*> list_save; // inter pour l'appel de la fonction
                list_save.push_back(saveResul);
                Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                    (Cij_nD(i)(j),1 // une seule valeur attendue en retour
                     ,ex_impli,ex_expli_tdt,ex_expli
                     ,NULL
                     ,NULL
                     ,&list_save
                    );
                // on récupère le premier élément du tableau uniquement
                Cij_use(i)(j) *= tab_val(1);
              };
          };

      };
    if (sortie_post)
     {// sauvegarde des paramètres matériau
      int ia = 1; // init
      // d'abord K
      (*save_resulHyper_W.para_loi)(ia) = K_use;
      // puis les coef
      int tai = Cij.Taille();
      for (int i=1;i<= tai; i++)
        { int taj = Cij(i).Taille();
          for (int j=1;j<=taj;j++,ia++)
           (*save_resulHyper_W.para_loi)(ia) = Cij_use(i)(j);
        };
     };

    // calcul des invariants et de leurs variations premières et secondes
    Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
    // calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
    Potentiel_et_var2(module_compressibilite);
    // stockage éventuel pour du post-traitement
    if (sortie_post)
     { // récup du conteneur spécifique du point, pour sauvegarde éventuelle
       SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
       save_resulHyper_W.invP->potentiel= W_d + W_v + W_c;
     };
    // calcul du tenseur des contraintes
    double unSurV=1./V;
    sigHH = ((W_d_J1+W_c_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2) 
                + ((W_v_J3+W_c_J3)/V) * d_J_r_epsBB_HH(3);
    #ifdef MISE_AU_POINT
    if  (Permet_affichage() > 4) cout << "\n sigHH " << sigHH;
    #endif

    // calcul de la variation seconde du potentiel par rapport à epsij epskl
    Tenseur3HHHH d2W_d2epsHHHH 
        =  // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
      (W_d_J1_2 + W_c_J1_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
         + W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
         + W_d_J2_2  * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
           // puis les dérivées premières du potentiel déviatoire + courbure éventuellement
       + (W_d_J1 + W_c_J1) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
           // enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
       + (W_v_J3 + W_c_J3) * d_J_3_eps2BB_HHHH
       + (W_v_J3J3 + W_c_J3_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
    if (avec_courbure)
        d2W_d2epsHHHH += W_c_J1_J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(3));

    // calcul de la variation du tenseur des contraintes par rapports aux déformations
    // on tient compte du fait que V*sigHH = d W/ d epsij
    Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV);
    
    //	Tenseur3HHHH dSigdepsHHHH(1,interHH,d_J_r_epsBB_HH(3));
    Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt);
    //	dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH;
    
    Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH);  // cas des tenseurs généraux
    dSigdepsHHHH += interHHHH;                        // cas des tenseurs généraux
    
    //---------------------------------------------------------------------------------
    // vérif numérique de l'opérateur tangent
    //    Cal_dsigma_deps_num (*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien,dSigdepsHHHH);
    //---------------------------------------------------------------------------------

	   // calcul des variations / aux ddl
    int nbddl = d_gijBB_tdt.Taille();
    for (int i = 1; i<= nbddl; i++)
     { // on fait  uniquement une égalité d'adresse et de ne pas utiliser
       // le constructeur d'ou la profusion d'* et de ()
       Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i)));  // passage en dim 3
       const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); //    "
       dsigHH = dSigdepsHHHH && depsBB;
       };  
  
    // pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
    module_cisaillement = 0.;

    // traitement des énergies
    energ.Inita(0.); 
    energ.ChangeEnergieElastique((W_d+W_v)/V);
    
    LibereTenseur(); 
    LibereTenseurQ();
     
 };


// calcul des contraintes et ses variations  par rapport aux déformations a t+dt
// 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
void Poly_hyper3D::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & ,TenseurBB&
            ,TenseurBB & epsBB_tdt,TenseurBB &, 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 (epsBB_tdt.Dimension() != 3)
       { cout << "\nErreur : la dimension devrait etre 3 !\n";
      cout << " Poly_hyper3D::Calcul_dsigma_deps\n";
      Sortie(1);
    };
    #endif
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt);  // // passage en dim 3 explicite
    Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) ex.gijHH_tdt);  //  passage en dim 3 explicite
    
    // init des coefficients utilisés
    Cij_use = Cij;
    K_use = K;
    
    // informations éventuelles
    #ifdef MISE_AU_POINT
    if  (Permet_affichage() > 4)
      { cout << "\nPoly_hyper3D::Calcul_dsigma_deps(... ";
      };
    #endif
    // cas de la thermo dépendance, on calcul les grandeurs
    int tai = Cij.Taille();
  
    // on balaie les coefficients
    for (int i=1;i<= tai; i++)
      for (int j=1;j<= i+1; j++)
        if (Cij_temperature(i)(j) != NULL)
            {Cij_use(i)(j) = Cij_temperature(i)(j)->Valeur(*temperature);
             #ifdef MISE_AU_POINT
             if  (Permet_affichage() > 4)
                cout << "\n C_"<<i<<"_"<<j<<"(T)= " << Cij_use(i)(j) << " ";
             #endif
            };
    if (K_temperature != NULL)
      {K_use = K_temperature->Valeur(*temperature);
       #ifdef MISE_AU_POINT
       if  (Permet_affichage() > 4)
          cout << "\n K_(T)="<<K_use << " ";
       #endif
      };
    if (avec_courbure)
     { if (a_temperature != NULL)
         {a_courbure = a_temperature->Valeur(*temperature);
          #ifdef MISE_AU_POINT
          if  (Permet_affichage() > 4)
             cout << "\n a_courbure(T)="<<a_courbure << " ";
          #endif
         };
       if (r_temperature != NULL)
         {r_courbure = r_temperature->Valeur(*temperature);
          #ifdef MISE_AU_POINT
          if  (Permet_affichage() > 4)
             cout << "\n r_courbure(T)="<<r_courbure << " ";
          #endif
         };
     };
  
  // récup du conteneur spécifique du point, pour sauvegarde éventuelle
  SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
  // dans le cas de l'utilisation de fct nD
  int bai = Cij_nD.Taille();
  if ((!bai) || (K_nD != NULL)) // si diff de 0
    { // opération de transmission de la métrique
      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 (K_nD != NULL)
        // là il faut calculer la fonction nD
        { // on utilise la méthode générique de loi abstraite
          list <SaveResul*> list_save; // inter pour l'appel de la fonction
          list_save.push_back(saveResul);
          Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
              (K_nD,1 // une seule valeur attendue en retour
               ,ex_impli,ex_expli_tdt,ex_expli
               ,NULL
               ,NULL
               ,&list_save
              );
          // on récupère le premier élément du tableau uniquement
          K_use *= tab_val(1);
        };

      // puis les coef: on boucle sur le tableau
      for (int i=1;i<= bai; i++)
        { int baj = Cij_nD(i).Taille();
          for (int j=1;j<=baj;j++)
           if (Cij_nD(i)(j) != NULL)
              // là il faut calculer la fonction nD
            { // on utilise la méthode générique de loi abstraite
              list <SaveResul*> list_save; // inter pour l'appel de la fonction
              list_save.push_back(saveResul);
              Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
                  (Cij_nD(i)(j),1 // une seule valeur attendue en retour
                   ,ex_impli,ex_expli_tdt,ex_expli
                   ,NULL
                   ,NULL
                   ,&list_save
                  );
              // on récupère le premier élément du tableau uniquement
              Cij_use(i)(j) *= tab_val(1);
            };
        };
    };
    if (sortie_post)
     {// sauvegarde des paramètres matériau
      int ia = 1; // init
      // d'abord K
      (*save_resulHyper_W.para_loi)(ia) = K_use;
      // puis les coef
      int tai = Cij.Taille();
      for (int i=1;i<= tai; i++)
        { int taj = Cij(i).Taille();
          for (int j=1;j<=taj;j++,ia++)
           (*save_resulHyper_W.para_loi)(ia) = Cij_use(i)(j);
        };
     };


    // calcul des invariants et de leurs variations premières et secondes
    Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),*(ex.gijBB_tdt),gijHH_tdt,jacobien_0,jacobien);
    // calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
    Potentiel_et_var2(module_compressibilite);
    // stockage éventuel pour du post-traitement
    if (sortie_post)
     { // récup du conteneur spécifique du point, pour sauvegarde éventuelle
       SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
       save_resulHyper_W.invP->potentiel= W_d + W_v + W_c;
     };
    // calcul du tenseur des contraintes, on travaille ici dans le repère matériel finale correspondant
    // aux coordonnées initiales X_(0)^a, on obtient donc un tenseur dans la base naturelle finale
    double unSurV=1./V;
    Tenseur3HH  sig_localeHH  = ((W_d_J1+W_c_J1)*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2)
	        + ((W_v_J3+W_c_J3)*unSurV) * d_J_r_epsBB_HH(3);
	   // passage éventuelle dans la base I_a
    if (en_base_orthonormee)
       {sig_localeHH.BaseAbsolue(sigHH,*(ex.giB_tdt));}
    else {sigHH = sig_localeHH;}; // sinon la base locale est la bonne

    #ifdef MISE_AU_POINT
    if  (Permet_affichage() > 4)
      {cout << "\n sigHH en local" << sigHH;
       Tenseur3HH sig_3D_inter;
       sigHH.BaseAbsolue(sig_3D_inter,*(ex.giB_tdt));
       cout << "\n sigHH en absolue 3D : "<<sig_3D_inter;
      };
    #endif

    // --- modif 9 fevrier 2015 --- la formule n'était pas complète
    // calcul de la variation seconde du potentiel par rapport à epsij epskl
//    Tenseur3HHHH d2W_d2epsHHHH 
//        =  // tout d'abord les dérivées secondes du potentiel déviatoire
//           W_d_J1_2  * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
//         + W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
//         + W_d_J2_2  * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
//           // puis les dérivées premières du potentiel déviatoire
//          + W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
//           // enfin les dérivées seconde et première du potentiel sphérique
//          + W_v_J3 * d_J_3_eps2BB_HHHH
//          + W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));

    // calcul de la variation seconde du potentiel par rapport à epsij epskl
    // calcul de la variation seconde du potentiel par rapport à epsij epskl
    Tenseur3HHHH d2W_d2epsHHHH 
        =  // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
           (W_d_J1_2 + W_c_J1_2)  * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
         + W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
         + W_d_J2_2  * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
           // puis les dérivées premières du potentiel déviatoire + courbure éventuellement
       + (W_d_J1 + W_c_J1) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
           // enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
       + (W_v_J3 + W_c_J3) * d_J_3_eps2BB_HHHH
       + (W_v_J3J3 + W_c_J3_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
    //	// calcul de la variation seconde du potentiel par rapport à epsij epskl
    //	Tenseur3HHHH d2W_d2epsHHHH 
    //	    = W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
    //	      + W_v_J3 * d_J_3_eps2BB_HHHH
    //	      + W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
    // --- fin modif ---
  
    // calcul de la variation du tenseur des contraintes par rapports aux déformations
    // on tient compte du fait que V*sigHH = d W/ d epsij
    Tenseur3HH interHH = -sig_localeHH ; //* (-0.5*unSurV*unSurV);
    
    // calcul de la variation du tenseur des contraintes par rapports aux déformations
    // on tient compte du fait que V*sigHH = d W/ d epsij
    Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt);
    //	dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH;
    
    Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH);  // cas des tenseurs généraux
    dSigdepsHHHH += interHHHH;                        // cas des tenseurs généraux

      // transfert des informations: on pas d'un tenseur de 81 composantes à 36         
      // avec des symétries par rapport aux deux premiers indices et par rapport aux deux derniers
    ///    Tenseur3HHHH   d_sigma_depsHHHH; d_sigma_depsHHHH.TransfertDunTenseurGeneral(dSigdepsHHHH.Symetrise1et2_3et4());
    // calcul de la première partie de l'opérateur tangent (correspond au changement de repère
    // gi_tdt -> Ia de l'opérateur calculer précédemment
    Tenseur3HHHH & d_sigma_depsFinHHHH =  *((Tenseur3HHHH*) &d_sigma_deps_); // pour accés directe
	   // passage éventuelle dans la base I_a
    if (en_base_orthonormee)
       {dSigdepsHHHH.ChangeBase(d_sigma_depsFinHHHH,*(ex.giB_tdt));}
    else
       {d_sigma_depsFinHHHH = dSigdepsHHHH;};
  
//    // calcul du module de compressibilité
//    module_compressibilite = 2. * V * (W_v_J3+W_c_J3) / (MaX(ConstMath::petit, log(V) ));
    // pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
    module_cisaillement = 0.;

    // traitement des énergies
    energ.Inita(0.); 
    energ.ChangeEnergieElastique((W_d+W_v)/V);
    
    LibereTenseur(); 
    LibereTenseurQ();
     
 };

//---------------------- méthodes privées -------------------------------
	   // calcul du potentiel et de ses dérivées premières / aux invariants J_r
void Poly_hyper3D::Potentiel_et_var(double & module_compressibilite)
	{ // calcul de grandeurs intermédiaires
	  double unSurV=1./V; // le cas V=0 a normalement été traité dans Hyper_W_gene_3D
	  double unSurV2=unSurV*unSurV;
	  // calcul du potentiel
	  double logV = log(V);
	  
   W_d=0.;       // init W_d
	  W_d_J1 = 0.;  // init variation / J1
	  W_d_J2 = 0.;  // init variation / J2
	  int tai=Cij_use.Taille();
   for (int n=1;n<= tai; n++) // n c'est le degré
      for (int k=1;k<=n+1;k++)
      	  { W_d += Cij_use(n)(k) * PuissPoten((J_r(1)-3.),k-1) * PuissPoten((J_r(2)-3.),n-k+1);
      		   W_d_J1 += Cij_use(n)(k) * (k-1) * PuissPoten((J_r(1)-3.),k-2) * PuissPoten((J_r(2)-3.),n-k+1);
      		   W_d_J2 += Cij_use(n)(k) * PuissPoten((J_r(1)-3.),k-1) * (n-k+1) * PuissPoten((J_r(2)-3.),n-k);
      		 };
   #ifdef MISE_AU_POINT
   if  (Permet_affichage() > 4)
     { cout << "\n Poly_hyper3D::Potentiel_et_var((..."
            << "\n W_d= " << W_d << ", W_d_J1= "<< W_d_J1 <<", W_d_J2= "<< W_d_J2 << " ";
     };
   #endif
   
   // puis partie volumique
   double VmoinsUn = V-1.;
   switch (type_pot_vol)
    {case 1: {W_v = K_use * (1- (1+logV)*unSurV);         // potentiel 1
              W_v_J3 = 0.5 * K_use * logV * unSurV * unSurV2;    // variation / J3
              module_compressibilite = K_use  *  unSurV2;
              break;
              }
     case 2: {W_v = 0.5 * K_use * VmoinsUn;         // potentiel 2
              W_v_J3 = 0.25 * K_use * unSurV;    // variation / J3
              if (Dabs(VmoinsUn) > ConstMath::petit )
                {module_compressibilite = 0.5 * K_use  / logV;}
              else // on fait un développement limité du log pour V proche de 1.
                {module_compressibilite = 0.5 * K_use  / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
              break;
              }
     case 3: {W_v = K_use * 0.5 * logV * logV; // potentiel 3
              W_v_J3 = K_use * 0.5 * unSurV2 * logV; // variation / J3
              module_compressibilite =  K_use / V ;
              break;
              }
     case 4: {W_v = 0.5 * K_use * VmoinsUn * VmoinsUn; // potentiel
              W_v_J3 = 0.5 * K_use * unSurV * VmoinsUn;    // variation / J3
              if (Dabs(VmoinsUn) > ConstMath::petit )
                {module_compressibilite = VmoinsUn * K_use /logV ;}
              else // on fait un développement limité du log pour V proche de 1.
                {module_compressibilite = VmoinsUn * K_use  / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
              break;
              }
    };
   #ifdef MISE_AU_POINT
   if  (Permet_affichage() > 4)
     {switch (type_pot_vol)
       {case 1: case 2: case 3: case 4:
          cout << "\n W_v= " << W_v << ", W_v_J3= "<< W_v_J3 << ", K_s= " << module_compressibilite << " ";
          break;
       };
     };
   #endif
   
	  //--- partie courbure éventuelle
	  if (avec_courbure)
	   {double unSurV3=unSurV2*unSurV;
     double aP2r = pow(a_courbure,2.*r_courbure);
     double J_rMoins3P2r = pow((J_r(1)-3.),(2.*r_courbure));
     W_c = unSurV * ( J_rMoins3P2r * (J_r(1)-3.) / aP2r / (2.*r_courbure+1) );
     W_c_J1 = unSurV * ( J_rMoins3P2r/aP2r );
     W_c_J3 = -0.5*unSurV3 * ( J_rMoins3P2r * (J_r(1)-3.) /aP2r/(2.*r_courbure+1) );
     module_compressibilite += 2. * V * W_c_J3 / (MaX(ConstMath::petit, log(V) ));
     #ifdef MISE_AU_POINT
     if  (Permet_affichage() > 4)
       {cout << "\n W_c= " << W_c << ", W_c_J1= "<< W_c_J1 << ", W_c_J3= " << W_c_J3 ;
       };
     #endif
    }
	  else
	   {W_c = W_c_J1 = W_c_J3 = 0.;};	
	};

	   // calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
void Poly_hyper3D::Potentiel_et_var2(double & module_compressibilite)
	{ // calcul de grandeurs intermédiaires
	  double unSurV=1./V; // le cas V=0 a normalement été traité dans Hyper_W_gene_3D
	  // calcul du potentiel
	  double logV = log(V);
	  
   W_d=0.;       // init W_d
	  W_d_J1 = 0.;  // init variation / J1
	  W_d_J2 = 0.;  // init variation / J2
	  W_d_J1_2 = W_d_J1_J2 = W_d_J2_2 = 0.; // init des dérivées secondes
	  int tai=Cij_use.Taille();
   for (int n=1;n<= tai; n++) // n c'est le degré
     for (int k=1;k<=n+1;k++)
       { W_d += Cij_use(n)(k) * PuissPoten((J_r(1)-3.),k-1) * PuissPoten((J_r(2)-3.),n-k+1);
      		 W_d_J1 += Cij_use(n)(k) * (k-1) * PuissPoten((J_r(1)-3.),k-2) * PuissPoten((J_r(2)-3.),n-k+1);
      		 W_d_J2 += Cij_use(n)(k) * PuissPoten((J_r(1)-3.),k-1) * (n-k+1) * PuissPoten((J_r(2)-3.),n-k);
      		 // dérivées secondes
      		 W_d_J1_2 += Cij_use(n)(k) * (k-1)*(k-2) * PuissPoten((J_r(1)-3.),k-3) * PuissPoten((J_r(2)-3.),n-k+1);
      		 W_d_J1_J2 += Cij_use(n)(k) * (k-1) * PuissPoten((J_r(1)-3.),k-2) * (n-k+1) * PuissPoten((J_r(2)-3.),n-k);
      		 W_d_J2_2 += Cij_use(n)(k) * PuissPoten((J_r(1)-3.),k-1) * (n-k+1)*(n-k) * PuissPoten((J_r(2)-3.),n-k-1);
       };
   #ifdef MISE_AU_POINT
   if  (Permet_affichage() > 4)
     { cout << "\n Poly_hyper3D::Potentiel_et_var2(..."
            << "\n W_d= " << W_d << ", W_d_J1= "<< W_d_J1 <<", W_d_J2= "<< W_d_J2
            << ", W_d_J1_2= "<< W_d_J1_2 <<", W_d_J1_J2= "<< W_d_J1_J2
            <<", W_d_J2_2= "<< W_d_J2_2 << " ";
     };
   #endif
   
   // puis partie volumique
   double unSurV2=unSurV*unSurV;
   double VmoinsUn = V-1.;
   switch (type_pot_vol)
    {case 1: {W_v = K_use * (1- (1+logV)*unSurV);         // potentiel 1
              W_v_J3 = 0.5 * K_use * logV * unSurV * unSurV2;    // variation / J3
              // calcul des variations secondes / (J3)^2
              W_v_J3J3 = 0.25 * K_use * unSurV2 * unSurV2 * unSurV *(1.-3.*logV);
              module_compressibilite = K_use  *  unSurV2;
              break;
              }
     case 2: {W_v = 0.5 * K_use * VmoinsUn;         // potentiel 2
              W_v_J3 = 0.25 * K_use * unSurV;    // variation / J3
              // calcul des variations secondes / (J3)^2
              W_v_J3J3 = -K_use * 0.125 * unSurV2 * unSurV;
              if (Dabs(VmoinsUn) > ConstMath::petit )
                {module_compressibilite = 0.5 * K_use  / logV;}
              else // on fait un développement limité du log pour V proche de 1.
                {module_compressibilite = 0.5 * K_use  / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
              break;
              }
     case 3: {W_v = K_use * 0.5 * logV * logV; // potentiel 3
              W_v_J3 = K_use * 0.5 * unSurV2 * logV; // variation / J3
              // calcul des variations secondes / (J3)^2
              W_v_J3J3 = K_use * 0.25 * unSurV2 * unSurV2 * (1.-2.*logV);
              module_compressibilite =  K_use / V ;
              break;
              }
     case 4: {W_v = 0.5 * K_use * VmoinsUn * VmoinsUn; // potentiel
              W_v_J3 = 0.5 * K_use * unSurV * VmoinsUn;    // variation / J3
              // calcul des variations secondes / (J3)^2
              W_v_J3J3 = K_use * 0.25 * unSurV2 * unSurV;
              if (Dabs(VmoinsUn) > ConstMath::petit )
                {module_compressibilite = VmoinsUn * K_use /logV ;}
              else // on fait un développement limité du log pour V proche de 1.
                {module_compressibilite = VmoinsUn * K_use  / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
              break;
              }
    };
   #ifdef MISE_AU_POINT
   if  (Permet_affichage() > 4)
     {switch (type_pot_vol)
       {case 1: case 2: case 3: case 4:
          cout << "\n W_v= " << W_v << ", W_v_J3= "<< W_v_J3 << ", W_v_J3J3= "<< W_v_J3J3
               << ", K_s= " << module_compressibilite << " ";
          break;
       };
     };
   #endif

	  // -- partie courbure éventuelle ---
	  if (avec_courbure)
	   {double unSurV3=unSurV2*unSurV;
     double aP2r = pow(a_courbure,2.*r_courbure);
     double J_rMoins3P2rMoins1 = pow((J_r(1)-3.),(2.*r_courbure-1));
     double J_rMoins3P2r = J_rMoins3P2rMoins1 * (J_r(1)-3.);
     W_c = unSurV * ( J_rMoins3P2r * (J_r(1)-3.) / aP2r / (2.*r_courbure+1) );
     W_c_J1 = unSurV * ( J_rMoins3P2r/aP2r );
     W_c_J3 = -0.5*unSurV3 * ( J_rMoins3P2r * (J_r(1)-3.) /aP2r/(2.*r_courbure+1) );
     W_c_J1_2 = unSurV * r_courbure * J_rMoins3P2rMoins1 / aP2r;
     W_c_J3_2 = 3./4. * unSurV3 * unSurV2 * ( J_rMoins3P2r * (J_r(1)-3.) /aP2r/(2.*r_courbure+1) ); 
     W_c_J1_J3 = -0.5*unSurV3 * ( J_rMoins3P2r/aP2r );
     module_compressibilite += 2. * V * W_c_J3 / (MaX(ConstMath::petit, log(V) ));
     #ifdef MISE_AU_POINT
     if  (Permet_affichage() > 4)
       {cout << "\n W_c= " << W_c << ", W_c_J1= "<< W_c_J1 << ", W_c_J3= " << W_c_J3
             << ", W_c_J1_2= "<< W_c_J1_2 << ", W_c_J3_2= " << W_c_J3_2
             << ", W_c_J1_J3= " << W_c_J1_J3;
       };
     #endif
    }
	  else
	   {W_c = W_c_J1 = W_c_J3 = W_c_J1_2 = W_c_J3_2 = W_c_J1_J3 =0.;};		
  
//	  cout << "\n Jr= " << J_r(1) << " " << J_r(2) << " " << J_r(3) << " pot " << W_d << " " << W_v ;
	};


	   // calcul de la dérivée numérique de la contrainte
void Poly_hyper3D::Cal_dsigma_deps_num (const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
		                       ,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
		                       ,const double& jacobien_0,const double& jacobien
		                       ,Tenseur3HHHH& dSigdepsHHHH)
		  	
		  	
	{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
   const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); //  "
   const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); //  "
   const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); //  "
	
   Tenseur3BB  gijBBtdt_N; // tenseur  modifié
   Tenseur3HH  gijHHtdt_N; // idem_0
	  double delta = ConstMath::unpeupetit*10.;
	  double unSurDelta = 1./delta;
   
   // init des coefficients utilisés
   Cij_use = Cij;
   K_use = K;

   // cas de la thermo dépendance, on calcul les grandeurs
   int tai = Cij.Taille();
   // on balaie les coefficients
   for (int i=1;i<= tai; i++)
     for (int j=1;j<= i+1; j++)
     	 if (Cij_temperature(i)(j) != NULL) 
     	     Cij_use(i)(j) = Cij_temperature(i)(j)->Valeur(*temperature);
   if (K_temperature != NULL) K_use = K_temperature->Valeur(*temperature);
	  if (avec_courbure)
    { if (a_temperature != NULL) a_courbure = a_temperature->Valeur(*temperature);
      if (r_temperature != NULL) r_courbure = r_temperature->Valeur(*temperature);
		  };

	  // cas des contraintes et de ses variations analytiques
//	  Tenseur3HHHH dSigdepsHHHH; // le tenseur contenant les dérivées analytiques
	  Tenseur3HH SigmaHH_deb;
	  Cal_sigmaEtDer_pour_num(gijBB_0_,gijHH_0_,gijBB_tdt_,gijHH_tdt_
		                       ,jacobien_0,jacobien,SigmaHH_deb,dSigdepsHHHH);      
	  // dimensionnement pour la matrice numérique
	  Tenseur3HHHH dSigdepsHHHH_num;
	  
	  // on va boucler sur les composantes de gijBB
	  for (int i=1;i<=3;i++) 
	    for (int j=1;j<=3;j++)
	     { gijBBtdt_N = gijBB_tdt;
	       gijBBtdt_N.Coor(i,j) += delta;
	       // en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i)
	       // car le tenseur est symétrique
	       // on a donc en variation numérique la somme des deux dérivées
	       // on définit un coeff multiplicatif qui vaut 1 ou 0.5
	       double coef=1.; if (i != j) coef = 0.5;
	       gijHHtdt_N = gijBBtdt_N.Inverse();
	       double jacobien_N=sqrt(gijBBtdt_N.Det());
	         
	       // cas des contraintes 
	       Tenseur3HH SigmaHH_N;
	       Cal_sigma_pour_num(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
		                       ,jacobien_0,jacobien_N,SigmaHH_N);
		                       	  		  	
	       // calcul des dérivées numériques et comparaisons
	       for (int k=1;k<=3;k++)
	         for (int l=1;l<=3;l++)
	           { // 
	             double derSigNum = coef * 2.*(SigmaHH_N(k,l) - SigmaHH_deb(k,l) )*unSurDelta;
	             dSigdepsHHHH_num.Change(k,l,i,j,derSigNum);
	             double derSigAna = dSigdepsHHHH(k,l,i,j);//0.5*(dSigdepsHHHH(k,l,i,j) + dSigdepsHHHH(k,l,j,i));
              bool erreur = false;
              if (diffpourcent(derSigNum,derSigAna,MaX(Dabs(derSigNum),Dabs(derSigAna)),0.1))
              if (MaX(Dabs(derSigNum),Dabs(derSigAna)) > 200.)
               {if (MiN(Dabs(derSigNum),Dabs(derSigAna)) == 0.)
                 {if ( MaX(Dabs(derSigNum),Dabs(derSigAna)) > 50.*delta) erreur = true;}
                else erreur = true; 
               };
              erreur = false; // a virer
              if (erreur)
               { 
                 // calcul des dérivées d'éléments intermédiaires pour voir
                 
                 //
                 cout << "\n erreur dans le calcul analytique de l'operateur tangent " 
                      << "\n derSigNum= " << derSigNum << " derSigAna= " << derSigAna 
                      << " klij= " << k << " " << l << " " << i << " " << j
                      << " SigmaHH_N(k,l)= " << SigmaHH_N(k,l);
                 cout << "\n Poly_hyper3D::Calcul_derivee_numerique(..";
                 cout << "\n un caractere ";
                 // --- pour le débug ----
                 // calcul des invariants et de leurs variations premières en numérique
	                Invariants_et_var1_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
		                               ,jacobien_0,jacobien_N);
		               // calcul des invariants et de leurs variations premières et secondes
		               Invariants_et_var2_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
		                               ,jacobien_0,jacobien_N);
                 string toto;
//                 toto= lect_chaine();
               };
  
	           };
	     };
	  // passage des dérivées numériques aux dérivées finales
	  dSigdepsHHHH= dSigdepsHHHH_num;  
	};

	   // calcul de la contrainte avec le minimum de variable de passage, utilisé pour le numérique
void Poly_hyper3D::Cal_sigma_pour_num(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0
		                       ,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt
		                       ,const double& jacobien_0,const double& jacobien,TenseurHH & sigHH_)
	{    
   Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_);  //  passage en dim 3 explicite
	  // calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D)
//    Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// pour vérif on appelle var2, mais c'est à virer
   Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
   // calcul du potentiel et de ses dérivées premières / aux invariants J_r
   double module_compressibilite; // ne sert pas ici
   Potentiel_et_var(module_compressibilite);
	
//	// calcul du tenseur des contraintes
//	sigHH = (W_d_J1/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2) 
//	        + (W_v_J3/V) * d_J_r_epsBB_HH(3);
		
   // calcul du tenseur des contraintes
   sigHH = ((W_d_J1+W_c_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2) 
           + ((W_v_J3+W_c_J3)/V) * d_J_r_epsBB_HH(3);
	};

	   // idem avec la variation 
void Poly_hyper3D::Cal_sigmaEtDer_pour_num(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0
		                       ,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt
		                       ,const double& jacobien_0,const double& jacobien
		                       ,TenseurHH & sigHH_,Tenseur3HHHH& dSigdepsHHHH)	  		  	
	{    
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_);  //  passage en dim 3 explicite
	// calcul des invariants et de leurs variations premières et seconde (méthode de Hyper_W_gene_3D)
    Invariants_et_var2(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
	   // calcul du potentiel et de ses dérivées premières / aux invariants J_r
    double module_compressibilite; // ne sert pas ici
	   Potentiel_et_var2(module_compressibilite);
	
//	// calcul du tenseur des contraintes
//	double unSurV=1./V;
//	sigHH = (W_d_J1*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2) 
//	        + (W_v_J3*unSurV) * d_J_r_epsBB_HH(3);
   double unSurV=1./V;
   sigHH = ((W_d_J1+W_c_J1)*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2) 
           + ((W_v_J3+W_c_J3)*unSurV) * d_J_r_epsBB_HH(3);

	// calcul de la variation seconde du potentiel par rapport à epsij epskl
	Tenseur3HHHH d2W_d2epsHHHH 
	    = W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
	      + W_v_J3 * d_J_3_eps2BB_HHHH
	      + W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
//	cout << "\n d2W_d2epsHHHH(1,1,1,1)= " << d2W_d2epsHHHH(1,1,1,1);      
	
	// calcul de la variation du tenseur des contraintes par rapports aux déformations
	// on tient compte du fait que V*sigHH = d W/ d epsij
	Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV);
	
	Tenseur3HHHH d_igdepsHHHH(1,interHH,gijHH_tdt);
	d_igdepsHHHH += (unSurV) * d2W_d2epsHHHH;
	dSigdepsHHHH = d_igdepsHHHH;
	};