// 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 "IsoHyper3DFavier3.h"
//#include "ComLoi_comp_abstraite.h"

# include <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "ConstMath.h"
#include "CharUtil.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"


//================== initialisation des variables static ======================
// indicateur utilisé par Verif_Potentiel_et_var  
int IsoHyper3DFavier3::indic_Verif_PoGrenoble_et_var = 0;
double IsoHyper3DFavier3::limite_co2=700.; // limite à partir de laquelle on considère que cosh(co2) = infini
                                           // ici cosh(700.) = 5.0712e+303 !!

//================== fin d'initialisation des variables static ================
//================== stockage particulier SaveResul ================

IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::SaveResulIsoHyper3DFavier3(bool avec_para): // constructeur par défaut :
   SaveResulHyperD(),para_loi(NULL),para_loi_t(NULL)
      {if (avec_para)
        {para_loi = new Vecteur(4,-ConstMath::trespetit);
         para_loi_t = new Vecteur(4,-ConstMath::trespetit);
        };
      };
IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::SaveResulIsoHyper3DFavier3(const SaveResulIsoHyper3DFavier3& sav): // de copie
   SaveResulHyperD(sav),para_loi(NULL),para_loi_t(NULL)
     {if (sav.para_loi != NULL)
       {para_loi = new Vecteur(*sav.para_loi);
        para_loi_t = new Vecteur(*sav.para_loi_t);
       };
     };

IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::~SaveResulIsoHyper3DFavier3()  // destructeur
  { if (para_loi != NULL) delete para_loi;
    if (para_loi_t != NULL) delete para_loi_t;
  };
      // affectation
Loi_comp_abstraite::SaveResul & IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3
     ::operator = ( const Loi_comp_abstraite::SaveResul & a)
   { SaveResulHyperD::operator=(a);
     SaveResulIsoHyper3DFavier3& sav = *((SaveResulIsoHyper3DFavier3*) &a);
     if (sav.para_loi != NULL)
      {if (para_loi == NULL)
        {para_loi = new Vecteur(4,-ConstMath::trespetit);
         para_loi_t = new Vecteur(4,-ConstMath::trespetit);
        };
       *para_loi = *(sav.para_loi);
       *para_loi_t = *(sav.para_loi_t);
      };
     map_type_quelconque = sav.map_type_quelconque;
     return *this;
   };
      //============= lecture écriture dans base info ==========
         // cas donne le niveau de la récupération
         // = 1 : on récupère tout
         // = 2 : on récupère uniquement les données variables (supposées comme telles)
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::Lecture_base_info (ifstream& ent,const int cas)
   {SaveResulHyperD::Lecture_base_info (ent,cas);
    string toto;
    ent >> toto;
    if (toto == "sans_para")
     {if (para_loi != NULL)
        {delete para_loi;para_loi=NULL;
         delete para_loi_t;para_loi_t=NULL;
        };
     }
    else
     {if (para_loi == NULL)
       {para_loi = new Vecteur(4);
        para_loi_t = new Vecteur(4);
       };
      ent >> *para_loi;
      *para_loi_t=*para_loi;
     };
   };
         // cas donne le niveau de sauvegarde
         // = 1 : on sauvegarde tout
         // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::Ecriture_base_info(ofstream& sort,const int cas)
   {SaveResulHyperD::Ecriture_base_info (sort,cas);
    if (para_loi == NULL)
     {sort << " sans_para ";}
    else
     {sort << " avec_para ";
      sort << *para_loi;
     };
   };
 
      // mise à jour des informations transitoires
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::TdtversT()
   {SaveResulHyperD::TdtversT();
    if (para_loi != NULL)
      {*para_loi_t = *para_loi;};
    // mise à jour de la liste des grandeurs quelconques internes
    Mise_a_jour_map_type_quelconque();
   };
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::TversTdt()
   {SaveResulHyperD::TversTdt();
    if (para_loi != NULL)
      {*para_loi = *para_loi_t;};
    // mise à jour de la liste des grandeurs quelconques internes
    Mise_a_jour_map_type_quelconque();
   };

      // affichage à l'écran des infos
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::Affiche() const
   {SaveResulHyperD::Affiche();
    if (para_loi == NULL)
     {cout << " sans_para ";}
    else
     {cout << " avec_para ";
      cout << *para_loi;
     };
   };


//================== fin stockage particulier SaveResul ================

IsoHyper3DFavier3::IsoHyper3DFavier3 ()  : // Constructeur par defaut
 Hyper3D(ISOHYPER3DFAVIER3,CAT_MECANIQUE,false),K(ConstMath::trespetit),Qor(ConstMath::trespetit)
 ,mur(ConstMath::trespetit),mu_inf(ConstMath::trespetit)
 ,nQor(ConstMath::trespetit),gammaQor(ConstMath::trespetit)
 ,n_mu_inf(ConstMath::trespetit),gamma_mu_inf(ConstMath::trespetit)
 // dépendance nD éventuelle
 ,K_nD(NULL),Qor_nD(NULL),mur_nD(NULL),mu_inf_nD(NULL)
    {};
// Constructeur de copie
IsoHyper3DFavier3::IsoHyper3DFavier3 (const IsoHyper3DFavier3& loi) :
  Hyper3D (loi),K(loi.K),Qor(loi.Qor),mur(loi.mur),mu_inf(loi.mu_inf)
  ,nQor(loi.nQor),gammaQor(loi.gammaQor)
  ,n_mu_inf(loi.n_mu_inf),gamma_mu_inf(loi.gamma_mu_inf)
 // dépendance nD éventuelle
 ,K_nD(loi.K_nD),Qor_nD(loi.Qor_nD),mur_nD(loi.mur_nD),mu_inf_nD(loi.mu_inf_nD)

 {// on regarde s'il s'agit d'une fonction locale ou d'une fonction globale
  if (K_nD != NULL) // 1)
   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);
       };
  if (Qor_nD != NULL) // 2)
   if (Qor_nD->NomFonction() == "_")
       {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
        string non_courbe("_");
        Qor_nD = Fonction_nD::New_Fonction_nD(*loi.Qor_nD);
       };
  if (mur_nD != NULL) // 3)
   if (mur_nD->NomFonction() == "_")
       {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
        string non_courbe("_");
        mur_nD = Fonction_nD::New_Fonction_nD(*loi.mur_nD);
       };
  if (mu_inf_nD != NULL) // 4)
   if (mu_inf_nD->NomFonction() == "_")
       {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
        string non_courbe("_");
        mu_inf_nD = Fonction_nD::New_Fonction_nD(*loi.mu_inf_nD);
       };

  // --- on remplit avec les grandeurs succeptible d'être utilisées
//  // acces à listdeTouslesQuelc_dispo_localement
//  list <EnumTypeQuelconque >& list_tousQuelc =  ListdeTouslesQuelc_dispo_localement();
//  list_tousQuelc.push_back(E_YOUNG);
//  list_tousQuelc.push_back(NU_YOUNG);
//  // on supprime les doublons localement
//  list_tousQuelc.sort(); // on ordonne la liste
//  list_tousQuelc.unique(); // suppression des doublons
 };

  

// Lecture des donnees de la classe sur fichier
void IsoHyper3DFavier3::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
             ,LesFonctions_nD& lesFonctionsnD)
  { // --- lecture des quatres coefficients de la loi
    string nom_class_methode("IsoHyper3DFavier3::LectureDonneesParticulieres");
    // on regarde si K dépend d'une fonction nD
    if(strstr(entreePrinc->tablcar,"K_fonction_nD:")!=0)
     { string nom;
       string mot_cle2="K_fonction_nD:";

       // on lit le nom de la fonction
       string nom_fonct;
       bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
       if (!lec )
          { entreePrinc->MessageBuffer("**erreur01 en lecture**  "+mot_cle2);
            throw (UtilLecture::ErrNouvelleDonnee(-1));
            Sortie(1);
          };
       // 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 courbe
          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**erreur02** \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;
     }
    //  sinon si le module K n'a pas été lue, on le récupère en version scalaire
    else
     { // lecture du module
       *(entreePrinc->entree) >> K ;
     };
   // idem  si Qor dépend d'une fonction nD
   if(strstr(entreePrinc->tablcar,"Qor_fonction_nD:")!=0)
    { string nom;
      string mot_cle2="Qor_fonction_nD:";

      // on lit le nom de la fonction
      string nom_fonct;
      bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
      if (!lec )
         { entreePrinc->MessageBuffer("**erreur03 en lecture**  "+mot_cle2);
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
      // maintenant on définit la fonction
      if (lesFonctionsnD.Existe(nom_fonct))
        {Qor_nD = lesFonctionsnD.Trouve(nom_fonct);
        }
      else
        {// sinon il faut la lire maintenant
         string non("_");
         Qor_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
         // lecture de la courbe
         Qor_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
         // maintenant on vérifie que la fonction est utilisable
         if (Qor_nD->NbComposante() != 1 )
          { cout << "\n erreur en lecture, la fonction " << nom_fonct
                 << " est une fonction vectorielle a   " << Qor_nD->NbComposante()
                 << " composante alors qu'elle devrait etre scalaire ! "
                 << " elle n'est donc pas utilisable !! ";
            string message("\n**erreur04** \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 = Qor_nD->Tab_enu_etendu();
      if (tab_enu.Contient(TEMP))
        thermo_dependant=true;
    }
   //  sinon si le module Qor n'a pas été lue, on le récupère en version scalaire
   else
    { // lecture du module
      *(entreePrinc->entree) >> Qor ;
    };
   // idem si mur dépend d'une fonction nD
   if(strstr(entreePrinc->tablcar,"mur_fonction_nD:")!=0)
    { string nom;
      string mot_cle2="mur_fonction_nD:";

      // on lit le nom de la fonction
      string nom_fonct;
      bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
      if (!lec )
         { entreePrinc->MessageBuffer("**erreur05 en lecture**  "+mot_cle2);
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
      // maintenant on définit la fonction
      if (lesFonctionsnD.Existe(nom_fonct))
        {mur_nD = lesFonctionsnD.Trouve(nom_fonct);
        }
      else
        {// sinon il faut la lire maintenant
         string non("_");
         mur_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
         // lecture de la courbe
         mur_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
         // maintenant on vérifie que la fonction est utilisable
         if (mur_nD->NbComposante() != 1 )
          { cout << "\n erreur en lecture, la fonction " << nom_fonct
                 << " est une fonction vectorielle a   " << mur_nD->NbComposante()
                 << " composante alors qu'elle devrait etre scalaire ! "
                 << " elle n'est donc pas utilisable !! ";
            string message("\n**erreur06** \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 = mur_nD->Tab_enu_etendu();
      if (tab_enu.Contient(TEMP))
        thermo_dependant=true;
    }
   //  sinon si le module mur n'a pas été lue, on le récupère en version scalaire
   else
    { // lecture du module
      *(entreePrinc->entree) >> mur ;
    };
   // idem si mu_inf dépend d'une fonction nD
   if(strstr(entreePrinc->tablcar,"mu_inf_fonction_nD:")!=0)
    { string nom;
      string mot_cle2="mu_inf_fonction_nD:";

      // on lit le nom de la fonction
      string nom_fonct;
      bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
      if (!lec )
         { entreePrinc->MessageBuffer("**erreur07 en lecture**  "+mot_cle2);
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
      // maintenant on définit la fonction
      if (lesFonctionsnD.Existe(nom_fonct))
        {mu_inf_nD = lesFonctionsnD.Trouve(nom_fonct);
        }
      else
        {// sinon il faut la lire maintenant
         string non("_");
         mu_inf_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
         // lecture de la courbe
         mu_inf_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
         // maintenant on vérifie que la fonction est utilisable
         if (mu_inf_nD->NbComposante() != 1 )
          { cout << "\n erreur en lecture, la fonction " << nom_fonct
                 << " est une fonction vectorielle a   " << mu_inf_nD->NbComposante()
                 << " composante alors qu'elle devrait etre scalaire ! "
                 << " elle n'est donc pas utilisable !! ";
            string message("\n**erreur08** \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 = mu_inf_nD->Tab_enu_etendu();
      if (tab_enu.Contient(TEMP))
        thermo_dependant=true;
    }
   //  sinon si le module mu_inf n'a pas été lue, on le récupère en version scalaire
   else
    { // lecture du module
      *(entreePrinc->entree) >> mu_inf ;
    };
//    *(entreePrinc->entree) >> K >> Qor >> mur >> mu_inf ;
    // on regarde ensuite s'il y a la phase
    if (strstr(entreePrinc->tablcar,"avec_phase")!=NULL)
     { // lecture des 4 paramètres de phase
       entreePrinc->NouvelleDonnee();  // passage d'une ligne
       *(entreePrinc->entree) >> nQor >> gammaQor >> n_mu_inf >> gamma_mu_inf ;
       avec_phase=true;
     };
      
    // Lecture des paramètre specifique sur fichier (exe: regularisation sortie_post...
    Hyper3D::LectParaSpecifiques(entreePrinc,lesCourbes1D,lesFonctionsnD);
    
    // appel au niveau de la classe mère
    Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
                 (*entreePrinc,lesFonctionsnD);

    };
// affichage de la loi
void IsoHyper3DFavier3::Affiche() const 
  { cout << " \n loi de comportement 3D hyperelastique isotrope favier3 : " << Nom_comp(id_comp)
         << " paramètres : \n";
    if (K_nD != NULL)
       {cout << "\n K fonction nD:  ";
        cout << "K_fonction_nD:" << " ";
        if (K_nD->NomFonction() != "_")
             cout << K_nD->NomFonction();
        else
             K_nD->Affiche();
       }
    else
       { cout << " K= " << K ;};
    if (Qor_nD != NULL)
       {cout << "\n Qor fonction nD:  ";
        cout << "Qor_fonction_nD:" << " ";
        if (Qor_nD->NomFonction() != "_")
             cout << Qor_nD->NomFonction();
        else
             Qor_nD->Affiche();
       }
    else
       { cout << " Qor= " << Qor ;};
    if (mur_nD != NULL)
       {cout << "\n mur fonction nD:  ";
        cout << "mur_fonction_nD:" << " ";
        if (mur_nD->NomFonction() != "_")
             cout << mur_nD->NomFonction();
        else
             mur_nD->Affiche();
       }
    else
       { cout << " mur_nD= " << mur_nD ;};
    if (mu_inf_nD != NULL)
       {cout << "\n mu_inf fonction nD:  ";
        cout << "mu_inf_fonction_nD:" << " ";
        if (mu_inf_nD->NomFonction() != "_")
             cout << mu_inf_nD->NomFonction();
        else
             mu_inf_nD->Affiche();
       }
    else
       { cout << " mu_inf= " << mu_inf ;};
//    cout << " K= " << K << " ; Qor = " << Qor << " ; mur = " << mur << " ; mu_inf = " << mu_inf ;
    if (avec_phase)
     { cout << "\n cas de la loi avec phase, parametre de la phase: "
            << " nQor= " << nQor << " gammaQor= " << gammaQor << " n_mu_inf= " 
            << n_mu_inf << " gamma_mu_inf= " << gamma_mu_inf;
      }
    cout << endl;
    // appel de la classe mère
    Loi_comp_abstraite::Affiche_don_classe_abstraite();
  };
            
// affichage et definition interactive des commandes particulières à chaques lois
void IsoHyper3DFavier3::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# .......  loi de comportement 3D hyperelastique isotrope favier3  ........";
    if ((rep != "o") && (rep != "O" ) && (rep != "0") )
      { sort
         << "\n# **** exemple de cas AVEC la phase :"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#     K         |     Qor     |       mur     |   mu_inf  | avec_phase(falculatif) |"
         << "\n#------------------------------------------------------------------------------------"
         << "\n     160000           150             17000       220           avec_phase"
         << "\n#-------------------------------------------------"
         << "\n#    n_Q    |    gamma_Q  |  n_mu  |   gamma_mu  |"
         << "\n#-------------------------------------------------"
         << "\n      0.25          0.4        0.25        0.4  "
         << "\n# **** exemple de cas AVEC utilisation de fonction nD :"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#     K         |     Qor     |       mur     |   mu_inf  | avec_phase(falculatif) |"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#   16000    Qor_fonction_nD: fct1 17000  mu_inf_fonction_nD: fct2"
         << "\n#"
         << "\n#  une fonction nD est possible pour chacun des 4 coefficients, via les mots cle "
         << "\n# K_fonction_nD:  "
         << "\n# Qor_fonction_nD: "
         << "\n# mur_fonction_nD: "
         << "\n# mu_inf_fonction_nD: "
         << "\n# on peut utiliser 1 ou 2 ou 3 ou 4 fct nD, au choix "
         << "\n# ( bien noter que la loi obtenue peut-etre quelconque, en particulier plus du tout "
         << "\n# hyperelastique, tout depend des choix de fct nD !) "
         << "\n#"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  il est possible d'indiquer 2 parametres specifiques de limite inf "
         << "\n#  limite_inf_Qeps_ : si Qeps < limite_inf_Qeps on considere que l'angle de phase et ses derivees sont nulles "
         << "\n#  limite_inf_bIIb_ : si Dabs(inv.bIIb) < limite_inf_bIIb   "
         << "\n#    les derivees du potentiel sont calculees pas difference fini, sauf ceux relatif a la phase qui sont mises a 0 "
         << "\n#  les  mots cle limite_inf_Qeps_ suivi du facteur voulu et limite_inf_bIIb_ suivi du facteur "
         << "\n#  doivent etre dans cet ordre "
         << "\n#  ex: "
         << "\n#    limite_inf_Qeps_  8.5e-5  limite_inf_bIIb_ 36.e-10 "
         << "\n#  ces mots cle doivent  se situer avant le mot cle avec_regularisation_ "
         << "\n#------------------------------------------------------------------------------------"
         << "\n#  il est possible d'indiquer un facteur de regularisation qui permet d'eviter "
         << "\n#  de potentiels  problemes de NaN, de type division par 0 par exemple "
         << "\n#  1/a est remplace par 1/(a+fact_regularisation), par defaut fact_regularisation = 1.e-12 "
         << "\n#  pour indiquer un facteur de regulation non nul on indique en dernier parametre "
         << "\n#  le mot cle avec_regularisation_ suivi du facteur voulu "
         << "\n#  ex: "
         << "\n#    avec_regularisation_  1.e-12 "
         << "\n#  ce mot cle doit se situer avant le mot cle sortie_post_ "
         << "\n#------------------------------------------------------------------------------------"
         << "\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, 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#------------------------------------------------------------------------------------";         
      };
    sort << "\n# **** exemple de cas SANS la phase :"
         << "\n#------------------------------------------------------------------------------------"
         << "\n#     K         |     Qor     |       mur     |   mu_inf  | avec_phase(falculatif) |"
         << "\n#------------------------------------------------------------------------------------"
         << "\n    160000           150             17000       220           "
         << endl;
    // appel de la classe Hyper3D
    Hyper3D::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 IsoHyper3DFavier3::TestComplet()
    { int ret = LoiAbstraiteGeneral::TestComplet();
      if ((K == ConstMath::trespetit) && (K_nD == NULL))
       { cout << " \n Le parametre K n'est pas defini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         Affiche();
         ret = 0;
       }
      if ((Qor == ConstMath::trespetit) && (Qor_nD == NULL))
       { cout << " \n Le parametre Qor n'est pas defini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         Affiche();
         ret = 0;
       }
      if ((mur == ConstMath::trespetit) && (mur_nD == NULL))
       { cout << " \n Le parametre mur n'est pas defini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         Affiche();
         ret = 0;
       }
      if ((mu_inf == ConstMath::trespetit) && (mu_inf_nD == NULL))
       { cout << " \n Le parametre mu_inf n'est pas defini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         Affiche();
         ret = 0;
       }
      if (avec_phase)
        if ((nQor == ConstMath::trespetit) || (gammaQor == ConstMath::trespetit)  
        ||  (n_mu_inf == ConstMath::trespetit)  ||  (gamma_mu_inf == ConstMath::trespetit))
         { cout << " \n Les paramètres de la phase ne sont pas défini pour  la loi " << Nom_comp(id_comp)
              << '\n';
         Affiche();     
         ret = 0;
       }
                
       return ret; 
    }; 

  // récupération des grandeurs particulière (hors ddl )
  // correspondant à liTQ
  // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
  void IsoHyper3DFavier3::Grandeur_particuliere
                (bool absolue,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list<int>& decal) const
  {  // tout d'abord on appelle la class mère
     HyperD::Grandeur_particuliere(absolue,liTQ,saveDon,decal);
     //  --- puis la  loi elle-même
     // on passe en revue la liste
     List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end();
     list<int>::iterator idecal=decal.begin();
     for (itq=liTQ.begin();itq!=itqfin;itq++,idecal++)
       {TypeQuelconque& tipParticu = (*itq); // pour simplifier
        if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
         switch (tipParticu.EnuTypeQuelconque().EnumTQ())
           {case PARA_LOI_FAVIER:
             //  ----- cas des paramètres de la loi
              { SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveDon);
                Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
                if (save_resul.para_loi != NULL)
                   {tyTQ(1+(*idecal))=  *(save_resul.para_loi);}
                else {tyTQ(1+(*idecal)).Zero();};
                (*idecal)++; break;
              }
            default: ;// on ne fait rien
           };
       };
  };

  // récupération et création de la liste de tous les grandeurs particulières
  // ces grandeurs sont ajoutées à la liste passées en paramètres
  // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
  void IsoHyper3DFavier3::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
   { // tout d'abord on appelle la class mère
     HyperD::ListeGrandeurs_particulieres(absolue,liTQ);
     //  --- puis la  loi elle-même
     // $$$ cas des paramètres de la loi de comportement
      {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
       for (itq=liTQ.begin();itq!=itqfin;itq++)
         if ((*itq).EnuTypeQuelconque() == PARA_LOI_FAVIER)
          {Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
           int taille = tyTQ.Taille()+1;
           tyTQ.Change_taille(taille); nexistePas = false;
          };
       if (nexistePas)
         {Vecteur  tens(4); // les 4 paramètres
          Tab_Grandeur_Vecteur gr_vec_para(tens,1);
          // def d'un type quelconque représentatif
          TypeQuelconque typQ(PARA_LOI_FAVIER,SIG11,gr_vec_para);
          liTQ.push_back(typQ);
        };
      };

   };

	   //----- 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 IsoHyper3DFavier3::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef
                              ,LesCourbes1D& lesCourbe1D,LesFonctions_nD& lesFonctionsnD)
  { string toto,nom;
    if (cas == 1)
     { ent >> toto >> thermo_dependant >> toto;
       // tout d'abord la lecture de K
       int type =0;
       ent >> type;
       switch (type)
        { case 1 :
            {ent >> nom;
             if (nom != " K_fonction_nD: ")
                 { cout << "\n erreur en lecture de la fonction nD, on attendait "
                        << " K_fonction_nD: et on a lue " << nom
                        << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                   Sortie(1);
                 };
               K_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,K_nD);
               // on regarde si la fonction nD intègre la température
               // car l'utilisateur peut éventuellement lors d'un restart changer
               // de  paramètres pour la fonction nD
               const Tableau <Ddl_enum_etendu>& tab_enu = K_nD->Tab_enu_etendu();
               if (tab_enu.Contient(TEMP))
                 thermo_dependant=true;
             break;
            };
          case 2 :
            {ent >> K;
             break;
            };
          default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
                        << " reconnus sont uniquement: 1 c-a-d K fonction nD,  "
                        << " 2 K une valeur fixe "
                        << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                   Sortie(1);
                   break;
        };
       // puis lecture de Qor
       type =0;
       ent >> toto >> type;
       switch (type)
        { case 1 :
            {ent >> nom;
             if (nom != " Qor_fonction_nD: ")
                 { cout << "\n erreur en lecture de la fonction nD, on attendait "
                        << " Qor_fonction_nD: et on a lue " << nom
                        << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                   Sortie(1);
                 };
               Qor_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Qor_nD);
               // on regarde si la fonction nD intègre la température
               // car l'utilisateur peut éventuellement lors d'un restart changer
               // de  paramètres pour la fonction nD
               const Tableau <Ddl_enum_etendu>& tab_enu = Qor_nD->Tab_enu_etendu();
               if (tab_enu.Contient(TEMP))
                 thermo_dependant=true;
             break;
            };
          case 2 :
            {ent >> Qor;
             break;
            };
          default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
                        << " reconnus sont uniquement: 1 c-a-d Qor fonction nD,  "
                        << " 2 Qor une valeur fixe "
                        << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                   Sortie(1);
                   break;
        };

      // puis lecture de mur
      type =0;
      ent >> toto >> type;
      switch (type)
       { case 1 :
           {ent >> nom;
            if (nom != " mur_fonction_nD: ")
                { cout << "\n erreur en lecture de la fonction nD, on attendait "
                       << " mur_fonction_nD: et on a lue " << nom
                       << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                  Sortie(1);
                };
              mur_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mur_nD);
              // on regarde si la fonction nD intègre la température
              // car l'utilisateur peut éventuellement lors d'un restart changer
              // de  paramètres pour la fonction nD
              const Tableau <Ddl_enum_etendu>& tab_enu = mur_nD->Tab_enu_etendu();
              if (tab_enu.Contient(TEMP))
                thermo_dependant=true;
            break;
           };
         case 2 :
           {ent >> mur;
            break;
           };
         default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
                       << " reconnus sont uniquement: 1 c-a-d mur fonction nD,  "
                       << " 2 mur une valeur fixe "
                       << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                  Sortie(1);
                  break;
       };

      // puis lecture de mu_inf
      type =0;
      ent >> toto >> type;
      switch (type)
       { case 1 :
           {ent >> nom;
            if (nom != " mu_inf_fonction_nD: ")
                { cout << "\n erreur en lecture de la fonction nD, on attendait "
                       << " mu_inf_fonction_nD: et on a lue " << nom
                       << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                  Sortie(1);
                };
              mu_inf_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu_inf_nD);
              // on regarde si la fonction nD intègre la température
              // car l'utilisateur peut éventuellement lors d'un restart changer
              // de  paramètres pour la fonction nD
              const Tableau <Ddl_enum_etendu>& tab_enu = mu_inf_nD->Tab_enu_etendu();
              if (tab_enu.Contient(TEMP))
                thermo_dependant=true;
            break;
           };
         case 2 :
           {ent >> mu_inf;
            break;
           };
         default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
                       << " reconnus sont uniquement: 1 c-a-d mu_inf fonction nD,  "
                       << " 2 mu_inf une valeur fixe "
                       << "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
                  Sortie(1);
                  break;
       };
     
   //    ent >> toto >> K >> toto >> Qor >> toto >> mur >> toto >> mu_inf >> avec_phase;
       // avec phase éventuelle
       ent  >> avec_phase;
       if (avec_phase)
        { ent >> toto >> nQor >> toto >> gammaQor >> toto
              >> n_mu_inf >> toto >> gamma_mu_inf;
        };
       // prise en compte éventuelle de paramètres particulier: ex: régularisation
       // géré par HyperD et Hyper3D
       Hyper3D::Lecture_para_specifiques(ent,cas);
       // gestion du post-traitement
       ent >> toto >> sortie_post ;
     };
    // appel class mère
    Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbe1D,lesFonctionsnD);
   };
       // cas donne le niveau de sauvegarde
       // = 1 : on sauvegarde tout
       // = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void IsoHyper3DFavier3::Ecriture_base_info_loi(ofstream& sort,const int cas)
  { if (cas == 1) 
     { sort << " ISOHYPER3DFAVIER3,thermodependance= " << thermo_dependant ;
       sort << "  module_K ";
       if (K_nD != NULL)
         {sort << " 1 K_fonction_nD: " << " ";
          LesFonctions_nD::Ecriture_pour_base_info(sort, cas,K_nD);
         }
       else
         { sort << " 2 " << K ; };
       sort << " seuil ";
       if (Qor_nD != NULL)
         {sort << " 1 Qor_fonction_nD: " << " ";
          LesFonctions_nD::Ecriture_pour_base_info(sort, cas,Qor_nD);
         }
       else
         { sort << " 2 " << Qor ; };
       sort << "  pente_origine ";
       if (mur_nD != NULL)
         {sort << " 1 mur_fonction_nD: " << " ";
          LesFonctions_nD::Ecriture_pour_base_info(sort, cas,mur_nD);
         }
       else
         { sort << " 2 " << mur ; };
       sort << " pente_infini ";
       if (mu_inf_nD != NULL)
         {sort << " 1 mu_inf_fonction_nD: " << " ";
          LesFonctions_nD::Ecriture_pour_base_info(sort, cas,mu_inf_nD);
         }
       else
         { sort << " 2 " << mu_inf ; };
//       sort << "  module_dilatation " << K << " seuil " << Qor
//            << "  pente_origine " << mur << " pente_infini " << mu_inf << " avec_phase " << avec_phase;
       sort << " avec_phase " << avec_phase;
       if (avec_phase)
        { sort << "\n nQor= " << nQor << " gammaQor= " << gammaQor << " n_mu_inf= " 
               << n_mu_inf << " gamma_mu_inf= " << gamma_mu_inf;
        };
       // prise en compte éventuelle de paramètres particulier: ex: régularisation
       // géré par HyperD et Hyper3D
       Hyper3D::Ecriture_para_specifiques(sort,cas);
       // gestion du post-traitement
       sort << " sortie_post= "<< sortie_post << " ";
     };
    // appel de la classe mère
    Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);  
  };

         
// calcul d'un module d'young équivalent à la loi, ceci pour un
// chargement nul  si rien n'est calculé, sinon c'est
// une valeur correspondante au dernier calcul
double IsoHyper3DFavier3::Module_young_equivalent(Enum_dure temps,const Deformation & def,SaveResul * saveDon)
  { // compte tenu du fait que l'on ne connait pas la métrique etc... on ramène le module en cours
    
    // on récupère le conteneur
    SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);

    double K_sortie=K; // init
    double mur_sortie=mur; // init
    double mu_inf_sortie=mu_inf; // init
    
    // traitement du cas particulier de l'existence de fct  nD
    if (save_resul.para_loi != NULL)
     {bool recup = false;
      switch (temps)
       {
        case TEMPS_t : // on utilise les grandeurs stockées à t
           if ((*save_resul.para_loi_t)(1) != (-ConstMath::trespetit))
            {K_sortie= (*save_resul.para_loi_t)(1);recup = true;}
           break;
        case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
           if ((*save_resul.para_loi)(1) != (-ConstMath::trespetit))
             {K_sortie= (*save_resul.para_loi)(1);recup = true;}
           break;
        case TEMPS_0 : // rien n'a été calculé
          recup = false;
       };
      // dans tous les autres cas on utilise soit les fonctions
      // si elles existent sinon on concerve les valeurs par défaut
      if (!recup)
        {//on essaie un calcul
         // cas des courbes d'évolution
         if (K_nD != NULL)
          // là il faut calcul la fonction nD
          {
            // 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 = K_nD->Li_enu_etendu_scalaire();
            List_io <TypeQuelconque >& li_quelc = K_nD->Li_equi_Quel_evolue();

            // on initialise les grandeurs du tableau pour les valeurs numériques
            Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
            // init par défaut des types quelconques
            List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
            for (il = li_quelc.begin();il != ilfin; il++)
               (*il).Grandeur_pointee()->InitParDefaut();

            // calcul de la valeur et retour dans tab_ret
            Tableau <double> & tab_val
                    = K_nD->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 relative au module d'Young "
                      << " doit calculer un scalaire or le tableau de retour est de taille "
                      << tab_val.Taille() << " ce n'est pas normal !\n";
                 cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
                 Sortie(1);
               };
            #endif
            // on récupère le premier élément du tableau uniquement
            K_sortie = tab_val(1);
          };
        };
        
      // -- idem pour mur
      recup = false;
      switch (temps)
       {
        case TEMPS_t : // on utilise les grandeurs stockées à t
           if ((*save_resul.para_loi_t)(3) != (-ConstMath::trespetit))
            {mur_sortie= (*save_resul.para_loi_t)(3);recup = true;}
           break;
        case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
           if ((*save_resul.para_loi)(3) != (-ConstMath::trespetit))
             {mur_sortie= (*save_resul.para_loi)(3);recup = true;}
           break;
             
        case TEMPS_0 : // rien n'a été calculé
          recup = false;
       };
      // dans tous les autres cas on utilise soit les fonctions
      // si elles existent sinon on concerve les valeurs par défaut
      if (!recup)
        {//on essaie un calcul
         // cas des courbes d'évolution
         if (mur_nD != NULL)
          // là il faut calcul la fonction nD
          {
            // 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 = mur_nD->Li_enu_etendu_scalaire();
            List_io <TypeQuelconque >& li_quelc = mur_nD->Li_equi_Quel_evolue();

            // on initialise les grandeurs du tableau pour les valeurs numériques
            Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
            // init par défaut des types quelconques
            List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
            for (il = li_quelc.begin();il != ilfin; il++)
               (*il).Grandeur_pointee()->InitParDefaut();

            // calcul de la valeur et retour dans tab_ret
            Tableau <double> & tab_val
                    = mur_nD->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 relative au module d'Young "
                      << " doit calculer un scalaire or le tableau de retour est de taille "
                      << tab_val.Taille() << " ce n'est pas normal !\n";
                 cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
                 Sortie(1);
               };
            #endif
            // on récupère le premier élément du tableau uniquement
            mur_sortie = tab_val(1);
          };
        };
      
      // -- idem pour mu_inf
      recup = false;
      switch (temps)
       {
        case TEMPS_t : // on utilise les grandeurs stockées à t
           if ((*save_resul.para_loi_t)(4) != (-ConstMath::trespetit))
            {mu_inf_sortie= (*save_resul.para_loi_t)(4);recup = true;}
           break;
        case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
           if ((*save_resul.para_loi)(4) != (-ConstMath::trespetit))
             {mu_inf_sortie= (*save_resul.para_loi)(4);recup = true;}
           break;
             
        case TEMPS_0 : // rien n'a été calculé
          recup = false;
       };
      // dans tous les autres cas on utilise soit les fonctions
      // si elles existent sinon on concerve les valeurs par défaut
      if (!recup)
        {//on essaie un calcul
         // cas des courbes d'évolution
         if (mu_inf_nD != NULL)
          // là il faut calcul la fonction nD
          {
            // 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 = mu_inf_nD->Li_enu_etendu_scalaire();
            List_io <TypeQuelconque >& li_quelc = mu_inf_nD->Li_equi_Quel_evolue();

            // on initialise les grandeurs du tableau pour les valeurs numériques
            Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
            // init par défaut des types quelconques
            List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
            for (il = li_quelc.begin();il != ilfin; il++)
               (*il).Grandeur_pointee()->InitParDefaut();

            // calcul de la valeur et retour dans tab_ret
            Tableau <double> & tab_val
                    = mu_inf_nD->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 relative au module d'Young "
                      << " doit calculer un scalaire or le tableau de retour est de taille "
                      << tab_val.Taille() << " ce n'est pas normal !\n";
                 cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
                 Sortie(1);
               };
            #endif
            // on récupère le premier élément du tableau uniquement
            mu_inf_sortie = tab_val(1);
          };
        };
     };
      
    // retour du module
    double module = (9.*K_sortie*(mur_sortie+mu_inf_sortie)/((mur_sortie+mu_inf_sortie)+3.*K_sortie));
    return module;
  };

// récupération d'un module de compressibilité équivalent à la loi, ceci pour un chargement nul
// si rien n'est calculé, sinon c'est
// une valeur correspondante au dernier calcul
double IsoHyper3DFavier3::Module_compressibilite_equivalent(Enum_dure temps,const Deformation & def,SaveResul * saveDon)
   { // compte tenu du fait que l'on ne connait pas la métrique etc... on ramène le module en cours
     
     // on récupère le conteneur
     SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);

     double K_sortie=K; // init
     // traitement du cas particulier de l'existence de fct  nD
     if (save_resul.para_loi != NULL)
      {bool recup = false;
       switch (temps)
        {
         case TEMPS_t : // on utilise les grandeurs stockées à t
            if ((*save_resul.para_loi_t)(1) != (-ConstMath::trespetit))
             {K_sortie= (*save_resul.para_loi_t)(1);recup = true;}
            break;
         case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
            if ((*save_resul.para_loi)(1) != (-ConstMath::trespetit))
              {K_sortie= (*save_resul.para_loi)(1);recup = true;}
            break;
         case TEMPS_0 : // rien n'a été calculé
           recup = false;
        };
       // dans tous les autres cas on utilise soit les fonctions
       // si elles existent sinon on concerve les valeurs par défaut
       if (!recup)
         {//on essaie un calcul
          // cas des courbes d'évolution
          if (K_nD != NULL)
           // là il faut calcul la fonction nD
           {
             // 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 = K_nD->Li_enu_etendu_scalaire();
             List_io <TypeQuelconque >& li_quelc = K_nD->Li_equi_Quel_evolue();

             // on initialise les grandeurs du tableau pour les valeurs numériques
             Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
             // init par défaut des types quelconques
             List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
             for (il = li_quelc.begin();il != ilfin; il++)
                (*il).Grandeur_pointee()->InitParDefaut();

             // calcul de la valeur et retour dans tab_ret
             Tableau <double> & tab_val
                     = K_nD->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 relative au module d'Young "
                       << " doit calculer un scalaire or le tableau de retour est de taille "
                       << tab_val.Taille() << " ce n'est pas normal !\n";
                  cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
                  Sortie(1);
                };
             #endif
             // on récupère le premier élément du tableau uniquement
             K_sortie = tab_val(1);
           };
         };
      };

     // retour du module
     double K= K_sortie/3.;
     return K;
   };


// ===========  METHODES Protégées dérivant de virtuelles : ==============

//  METHODES internes spécifiques à l'hyperélasticité isotrope découlant de 
//  méthodes virtuelles de Hyper3D

    // calcul du potentiel tout seul sans la phase car Qeps est  nul
    // ou très proche de 0
double IsoHyper3DFavier3::PoGrenoble
             (const double & Qeps,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   // des variables intermédiaires
   double a = Qor/2./mur;
   double co1 = Qor*a;
   double co2 = Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder
   {// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées
   double E = K/6. * (logV)*(logV) 
       + co1*log_A + mu_inf * Qeps * Qeps;
   // retour
   return E;
   };

    // calcul du potentiel tout seul avec la phase donc dans le cas où Qeps est non nul
double IsoHyper3DFavier3::PoGrenoble
             (const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };

   // dans le cas de l'existence de la phase, on modifie Qr et muinf
   double bmuinf,bQr,muinfF,Qr;
   if (avec_phase)
     {bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
      bQr = (1.+gammaQor*inv_spec.cos3phi);
      muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
      Qr = Qor/pow(bQr,nQor);
      }
   else
     {bmuinf = 1.;
      bQr = 1.;
      muinfF = mu_inf;
      Qr = Qor;
      }
   // des variables intermédiaires
   double a = Qr/2./mur;
   double co1 = Qr*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées
   double E = K/6. * (logV)*(logV) 
       + co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
   // retour
   return E;
   };

    // calcul du potentiel tout seul sans la phase car Qeps est  nul
    // ou très proche de 0, et de sa variation suivant V uniquement
Hyper3D::PoGrenoble_V IsoHyper3DFavier3::PoGrenoble_et_V
             (const double & Qeps,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   PoGrenoble_V ret;
   // des variables intermédiaires
   double a = Qor/2./mur;
   double co1 = Qor*a;
   double co2 = Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + mu_inf * Qeps * Qeps;
   ret.EV = K/3.*logV/inv.V;
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

    // calcul du potentiel et de sa variation suivant V uniquement
Hyper3D::PoGrenoble_V IsoHyper3DFavier3::PoGrenoble_et_V
             (const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   Hyper3D::PoGrenoble_V ret;
   // dans le cas de l'existence de la phase, on modifie Qr et muinf
   double bmuinf,bQr,muinfF,Qr;
   if (avec_phase)
     {bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
      bQr = (1.+gammaQor*inv_spec.cos3phi);
      muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
      Qr = Qor/pow(bQr,nQor);
      }
   else
     {bmuinf = 1.;
      bQr = 1.;
      muinfF = mu_inf;
      Qr = Qor;
      }
   //  des variables intermédiaires
   double a = Qr/2./mur;
   double co1 = Qr*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
   ret.EV = K/3.*logV/inv.V;
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

    // calcul du potentiel tout seul sans la phase car Qeps est  nul
    // ou très proche de 0, et de ses variations première et seconde suivant V uniquement
Hyper3D::PoGrenoble_VV IsoHyper3DFavier3::PoGrenoble_et_VV
             (const double & Qeps,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   PoGrenoble_VV ret;
   // des variables intermédiaires
   double a = Qor/2./mur;
   double co1 = Qor*a;
   double co2 = Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + mu_inf * Qeps * Qeps;
   ret.EV = K/3.*logV/inv.V;
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

    // calcul du potentiel et de sa variation première et seconde suivant V uniquement
Hyper3D::PoGrenoble_VV IsoHyper3DFavier3::PoGrenoble_et_VV
             (const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   Hyper3D::PoGrenoble_VV ret;
   // dans le cas de l'existence de la phase, on modifie Qr et muinf
   double bmuinf,bQr,muinfF,Qr;
   if (avec_phase)
     {bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
      bQr = (1.+gammaQor*inv_spec.cos3phi);
      muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
      Qr = Qor/pow(bQr,nQor);
      }
   else
     {bmuinf = 1.;
      bQr = 1.;
      muinfF = mu_inf;
      Qr = Qor;
      }
   // des variables intermédiaires
   double a = Qr/2./mur;
   double co1 = Qr*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
   ret.EV = K/3.*logV/inv.V;
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };

 // calcul du potentiel et de ses dérivées non compris la phase
Hyper3D::PoGrenobleSansPhaseSansVar IsoHyper3DFavier3::PoGrenoble
            (const InvariantQeps & inv_spec,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   PoGrenobleSansPhaseSansVar ret;
   // des variables intermédiaires
   double a = Qor/2./mur;
   double co1 = Qor*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + mu_inf * inv_spec.Qeps * inv_spec.Qeps;
   ret.EV = K/3.*logV/inv.V;
//   // dans le cas où l'on est très près de l'origine (voir nul)
//   // il faut un traitement particulier
//   if (co2 > ConstMath::unpeupetit)
    // cas normal
     ret.EQ = Qor * tanh(co2) + 2.* mu_inf * inv_spec.Qeps;
/*   else
    // cas ou  Qeps est très proche de zéro, on utilise un développement
    // limité  
     ret.EQ = (Qor/a + mu_inf)* inv_spec.Qeps 
              - 1./3. * Qor/(a*a*a) 
                      * inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps; */ 
       
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
   
   // calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PoGrenobleAvecPhaseSansVar IsoHyper3DFavier3::PoGrenoblePhase
             (const InvariantQepsCosphi& inv_spec,const Invariant & inv) 
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   PoGrenobleAvecPhaseSansVar ret;
   // dans le cas de l'existence de la phase, on modifie Qr et muinf
   double bmuinf,bQr,muinfF,Qr;
   if (avec_phase)
     {bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
      bQr = (1.+gammaQor*inv_spec.cos3phi);
      muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
      Qr = Qor/pow(bQr,nQor);
      }
   else
     {bmuinf = 1.;
      bQr = 1.;
      muinfF = mu_inf;
      Qr = Qor;
      }
   // des variables intermédiaires
   double a = Qr/2./mur;
   double co1 = Qr*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
   	 // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	}
   else
   	{ log_A = log(A);};
   double logV = log(inv.V);
   // le potentiel et ses dérivées
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
   ret.EV = K/3.*logV/inv.V;
//   // dans le cas où l'on est très près de l'origine (voir nul)
//   // il faut un traitement particulier
//   if (co2 > ConstMath::unpeupetit)
    // cas normal
     ret.EQ = Qr * tanh(co2) + 2.* muinfF * inv_spec.Qeps;
/*   else
    // cas ou  Qeps est très proche de zéro, on utilise un développement
    // limité  
     ret.EQ = (Qor/a + mu_inf)* inv_spec.Qeps 
              - 1./3. * Qor/(a*a*a) 
                      * inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps; */ 
       
   // cas de la variation faisant intervenir  la phase
   double dEdmuinf= inv_spec.Qeps * inv_spec.Qeps;
   double dmuinf_dcos3phi = (- n_mu_inf * gamma_mu_inf * muinfF/bmuinf);
   double dQdcos3phi = (- nQor * gammaQor * Qr / bQr);
   ret.Ecos3phi = dEdmuinf * dmuinf_dcos3phi + ret.EQ * dQdcos3phi;
   ret.Ks = K / 3. *(logV/2. + 1.);
   
   // retour
   return ret;
   };
      

// calcul  du potentiel sans phase et dérivées avec  ses variations par rapport aux invariants
Hyper3D::PoGrenobleSansPhaseAvecVar IsoHyper3DFavier3::PoGrenoble_et_var
             (const Invariant2Qeps& inv_spec,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   PoGrenobleSansPhaseAvecVar ret;
   // des variables intermédiaires
   double a = Qor/2./mur;
   double co1 = Qor*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2); 
   double log_A=0.;
   double tanhco2 = 1.; // init par défaut
   if ((Abs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est //   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
     // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	 tanhco2 = 1.; // et dans ce cas on est à la limite : th(A) =(environ)= exp(co2)/exp(co2) = 1
   	}
   else
   	{ log_A = log(A);
      tanhco2 = tanh(co2);
    };
   double logV = log(inv.V);
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) 
       + co1*log_A + mu_inf * inv_spec.Qeps * inv_spec.Qeps;
   ret.EV = K/3.*logV/inv.V;
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
   ret.EQV = 0.;
   // dans le cas où l'on est très près de l'origine (voir nul)
   // il faut un traitement particulier
//   if (co2 > ConstMath::unpeupetit)
    // cas normal
    { ret.EQ = Qor * tanhco2 + 2.* mu_inf * inv_spec.Qeps;
      // dérivée seconde
      ret.EQQ = 2.*mur*(1. - tanhco2 * tanhco2) + 2.* mu_inf;
     }
/*   else
    // cas ou  Qeps est très proche de zéro, on utilise un développement
    // limité  
    { ret.EQ = (Qor/a + mu_inf)* inv_spec.Qeps 
              - 1./3. * Qor/(a*a*a) 
                      * inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;  
      // dérivée seconde
      ret.EQQ = (Qor/a + mu_inf) 
                - Qor/(a*a*a) 
                      * inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;  
     }*/
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
// ---      Verif_PoGrenoble_et_var(inv_spec.Qeps,inv,ret );

//--- débug: autre vérification à décommenter si on veut vérifier un nan événtuel !! --------------------
//   if ((Dabs(ret.E) > ConstMath::tresgrand) || (Dabs(ret.EV) > ConstMath::tresgrand)
//       || (Dabs(ret.EVV) > ConstMath::tresgrand) || (Dabs(ret.EQV) > ConstMath::tresgrand)
//       || (Dabs(ret.EQQ) > ConstMath::tresgrand) || (Dabs(ret.EQ) > ConstMath::tresgrand) )
//   	{ cout << "\n attention *** on a detecter un comportement bizarre sur le potentiel de la loi Favier3D "
//   	      << " potentiel= " << ret.E << " var/V= " << ret.EV << " var/Q= " << ret.EQ
//   	      << " var/VV= " << ret.EVV << " var/QQ= " << ret.EQQ  << " var/QV= " << ret.EQV << endl;
//   	};
//---- fin débug -----------------------------------------------------------------------------------------   	
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
    
    // calcul  du potentiel avec phase et dérivées avec  ses variations par rapport aux invariants
Hyper3D::PoGrenobleAvecPhaseAvecVar IsoHyper3DFavier3::PoGrenoblePhase_et_var
             (const Invariant2QepsCosphi& inv_spec,const Invariant & inv)
 { // calcul éventuelle des paramètres de la loi
   if (K_nD != NULL)
    {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (Qor_nD != NULL)
    {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mur_nD != NULL)
    {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   if (mu_inf_nD != NULL)
    {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
     (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
    };
   // sauvegarde éventuelle des paramètres matériau
   SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
   if (save_resul.para_loi != NULL)
     {(*save_resul.para_loi)(1) = K;
      (*save_resul.para_loi)(2) = Qor;
      (*save_resul.para_loi)(3) = mur;
      (*save_resul.para_loi)(4) = mu_inf;
     };
   Hyper3D::PoGrenobleAvecPhaseAvecVar ret;
   // dans le cas de l'existence de la phase, on modifie Qr et muinf
   double bmuinf,bQr,muinfF,Qr;
   if (avec_phase)
     {bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
      bQr = (1.+gammaQor*inv_spec.cos3phi);
      muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
      Qr = Qor/pow(bQr,nQor);
      }
   else
     {bmuinf = 1.;
      bQr = 1.;
      muinfF = mu_inf;
      Qr = Qor;
      }
//debug
//cout << "\n debug: IsoHyper3DFavier3::PoGrenoblePhase_et_var "
//     << " Qr="<<Qr << "  Qor=" << Qor << " muinfF="<<muinfF <<"  mu_inf="<<mu_inf<<endl;
//fin debug
   // des variables intermédiaires
   double a = Qr/2./mur;
   double co1 = Qr*a;
   double co2 = inv_spec.Qeps/a ;
   double A = cosh(co2);
   double log_A=0.;
   double tanhco2 = 1.; // init par défaut
   if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est infini
//   if ((Dabs(A)>ConstMath::pasmalgrand)||(!isfinite(A)))
//   if (Dabs(A)>ConstMath::tresgrand)
   	{// si co2 est  grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
   	 // ou exp(-co2) = inf
   	 // mais en fait on n'utilise que le log(A) donc on peut faire une approximation
   	 // soit co2 est grand 
     // on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
   	 // d'où log(A) =(environ)= co2-log(2)
   	 log_A = MaX(0.,Dabs(co2)-log(2.));       
   	 tanhco2 = 1.; // et dans ce cas on est à la limite : th(A) =(environ)= exp(co2)/exp(co2) = 1
   	}
   else
   	{ log_A = log(A);
      tanhco2 = tanh(co2);
    };
//debug
//cout << "\n debug: 3DFavier3::Po..ePhase_et_var "
//     << " A="<< A << " co2="<< co2<< " log_A="<<log_A << "  log(A)=" << log(A) << " tanhco2="<<tanhco2 <<"  tanh(co2)="<<tanh(co2) <<endl;
//fin debug
   double logV = log(inv.V);
   double logCosh = log_A;
   double unsurmur=1./mur;
   // le potentiel et ses dérivées premières
   ret.E = K/6. * (logV)*(logV) 
       + co1*logCosh + muinfF * inv_spec.Qeps * inv_spec.Qeps;
   ret.EV = K/3.*logV/inv.V;
   // dérivées secondes
   ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
   ret.EQV = 0.;
   // dans le cas où l'on est très près de l'origine (voir nul)
   // il faut un traitement particulier
//   if (co2 > ConstMath::unpeupetit)
    // cas normal
    { ret.EQ = Qr * tanhco2 + 2.* muinfF * inv_spec.Qeps;
      // dérivée seconde
      ret.EQQ = 2.*mur*(1. - tanhco2 * tanhco2) + 2.* muinfF;
     }
/*   else
    // cas ou  Qeps est très proche de zéro, on utilise un développement
    // limité  
    { ret.EQ = (Qr/a + muinfF)* inv_spec.Qeps 
              - 1./3. * Qr/(a*a*a) 
                      * inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;  
      // dérivée seconde
      ret.EQQ = (Qr/a + muinfF) 
                - Qr/(a*a*a) 
                      * inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;  
     }*/

   // cas de la variation première par rapport à la phase
   double dEdQr= Qr*unsurmur * logCosh - inv_spec.Qeps * tanhco2;
   double QepsSurQr = inv_spec.Qeps/Qr;
   double tanhco2_2=tanhco2*tanhco2;
   double dEdmuinf= inv_spec.Qeps * inv_spec.Qeps;
   double dmuinf_dcos3phi = (- n_mu_inf * gamma_mu_inf * muinfF/bmuinf);
   double dQdcos3phi = (- nQor * gammaQor * Qr / bQr);
   ret.Ecos3phi = dEdmuinf * dmuinf_dcos3phi + dEdQr * dQdcos3phi;
   // cas des variations seconde faisant intervenir la phase la phase
   double d2EdQr2=unsurmur*logCosh-2.* QepsSurQr * tanhco2 + co2*QepsSurQr*(1.-tanhco2_2);
   double d2muinf_dcos3phi_2 = n_mu_inf*(1.+n_mu_inf)* muinfF*gamma_mu_inf*gamma_mu_inf/bmuinf/bmuinf;
   double d2Qr_dcos3phi_2  = nQor*(1.+nQor)*Qr*gammaQor*gammaQor/bQr/bQr;
   double d2E_dQ_dmuinf = 2. * inv_spec.Qeps ;
   double d2E_dQ_dQr = tanhco2 - co2 * (1.-tanhco2_2);
   
   ret.EVcos3phi= 0.;
   ret.EQcos3phi= d2E_dQ_dmuinf * dmuinf_dcos3phi + d2E_dQ_dQr * dQdcos3phi;
   ret.Ecos3phi2= d2EdQr2 * ret.Ecos3phi * ret.Ecos3phi + dEdmuinf * d2muinf_dcos3phi_2 
                  + dEdQr * d2Qr_dcos3phi_2;
/*   ret.EVcos3phi= 0.;
   ret.EQcos3phi= 0.;
   ret.Ecos3phi2= 0.;
   ret.Ecos3phi = 0.;*/

   // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
// ---   Verif_PoGrenoble_et_var(inv_spec.Qeps,inv,inv_spec.cos3phi,ret );
// ---   Invariant2Qeps toto(inv_spec.Qeps,inv_spec.dQepsdbIIb,inv_spec.dQepsdbIIb2);
// ---   Hyper3D::PoGrenobleSansPhaseAvecVar truc = PoGrenoble_et_var(toto,inv);
   ret.Ks = K / 3. *(logV/2. + 1.);
   // retour
   return ret;
   };
 
              
///=============== fonctions pour la vérification et la mise au point ===============
    
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DFavier3::Verif_PoGrenoble_et_var
            (const double & Qeps,const Invariant & inv
               ,const PoGrenobleSansPhaseAvecVar& potret )
  { // calcul éventuelle des paramètres de la loi
    if (K_nD != NULL)
     {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    if (Qor_nD != NULL)
     {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    if (mur_nD != NULL)
     {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    if (mu_inf_nD != NULL)
     {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    // sauvegarde éventuelle des paramètres matériau
    SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
    if (save_resul.para_loi != NULL)
      {(*save_resul.para_loi)(1) = K;
       (*save_resul.para_loi)(2) = Qor;
       (*save_resul.para_loi)(3) = mur;
       (*save_resul.para_loi)(4) = mu_inf;
      };
    // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenoble_et_var == 0)
     { cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var \n";
      } 
    indic_Verif_PoGrenoble_et_var++;
    // potret contiend le potentiel et ses variations première et seconde
    // on va calculer ces mêmes variations par différences finies et comparer les deux résultats
    // calcul des invariants sous la nouvelle forme
    
    double toto = potret.E; // pour que le débugger affiche potret
    
    Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl 
    double E_n = PoGrenoble(Qeps,inv_n0);  // la valeur du potentiel de référence
    // ici on est sans phase donc deux invariants indépendant V et Qeps,
    // pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
    // pour les dérivées premières et secondes
    double delta = ConstMath::unpeupetit;
    double Qeps_et_dQeps = Qeps+delta;
    double E_et_dQeps = PoGrenoble(Qeps_et_dQeps,inv_n0);
    
    double Qeps_et_dQeps2 = Qeps + 2.*delta;
    double E_et_dQeps2 = PoGrenoble(Qeps_et_dQeps2,inv_n0);
    
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    double E_et_dV = PoGrenoble(Qeps,inv_et_dV);
    
    Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
    double E_et_dV2 = PoGrenoble(Qeps,inv_et_dV2);
    
    Invariant inv_et_dVdQeps = inv_n0;
    inv_et_dVdQeps.V += delta;
    double Qeps_et_dVdQeps = Qeps+delta;
    double E_et_dVdQeps = PoGrenoble(Qeps_et_dVdQeps,inv_et_dVdQeps);
    // calcul des dérivées premières
    double E_V = (E_et_dV  - E_n)/delta;
    double E_Qeps = (E_et_dQeps - E_n)/delta;
    // calcul des dérivées secondes
    double E_V2 = (E_et_dV2 - 2.*E_et_dV  + E_n ) /(delta * delta);
    double E_Qeps2 = (E_et_dQeps2 - 2.*E_et_dQeps + E_n)/(delta * delta);
    double E_V_a_dQeps = (E_et_dVdQeps  - E_et_dQeps )/delta;
    double E_VQeps = ( E_V_a_dQeps - E_V)/delta;
    // comparaison avec les valeurs de dérivées analytiques
    bool erreur = false;
    if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
     {if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
        {if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = true;}
      else 
        erreur = true;}  
    if (diffpourcent(potret.EQ,E_Qeps,MaX(Dabs(E_Qeps),Dabs(potret.EQ)),0.05)) 
     {if (MiN(Dabs(E_Qeps),Dabs(potret.EQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps),Dabs(potret.EQ)) > 50.*delta) erreur = true;}
      else 
        erreur = true;}  
    if (diffpourcent(potret.EQQ,E_Qeps2,MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)),0.05)) 
     {if (MiN(Dabs(E_Qeps2),Dabs(potret.EQQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)) > 50.*delta) erreur = true;}
      else 
        erreur = true;}  
    if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05)) 
     {if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
        {if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = true;}
      else 
        erreur = true;}
    if (diffpourcent(potret.EQV,E_VQeps,MaX(Dabs(E_VQeps),Dabs(potret.EQV)),0.05)) 
     {if (MiN(Dabs(E_VQeps),Dabs(potret.EQV)) == 0.)
      {if (Dabs(potret.EQV) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps)) erreur = true;
          }
       else 
         if ( MaX(Dabs(E_VQeps),Dabs(potret.EQV)) > 150.*delta) erreur = true;
       }  
      else 
        erreur = true;
      }  
    if (erreur)
     { cout << "\n erreur dans le calcul analytique des derivees du potentiel";
       cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var(.."
            << " , numero d'increment = " << indic_Verif_PoGrenoble_et_var;
       Sortie(1);
      } 
 };
    // vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DFavier3::Verif_PoGrenoble_et_var
            (const double & Qeps,const Invariant & inv,const double& cos3phi
               ,const PoGrenobleAvecPhaseAvecVar& potret )
  { // calcul éventuelle des paramètres de la loi
    if (K_nD != NULL)
     {K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    if (Qor_nD != NULL)
     {Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    if (mur_nD != NULL)
     {mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    if (mu_inf_nD != NULL)
     {mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
      (mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
     };
    // sauvegarde éventuelle des paramètres matériau
    SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
    if (save_resul.para_loi != NULL)
      {(*save_resul.para_loi)(1) = K;
       (*save_resul.para_loi)(2) = Qor;
       (*save_resul.para_loi)(3) = mur;
       (*save_resul.para_loi)(4) = mu_inf;
      };
    // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_PoGrenoble_et_var == 0)
     { cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
       cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var \n";
      } 
    indic_Verif_PoGrenoble_et_var++;
    // potret contiend le potentiel et ses variations première et seconde
    // on va calculer ces mêmes variations par différences finies et comparer les deux résultats
    // calcul des invariants sous la nouvelle forme
    
    double toto = potret.E; // pour que le débugger affiche potret
    
    Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl 
    Invariant0QepsCosphi invcos3phi(Qeps,cos3phi);
    double E_n = PoGrenoble(invcos3phi,inv_n0);  // la valeur du potentiel de référence
    // ici on est sans phase donc deux invariants indépendant V et Qeps,
    // pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
    // pour les dérivées premières et secondes
    double delta = 10.*ConstMath::unpeupetit;
    double Qeps_et_dQeps = Qeps+delta;
    Invariant0QepsCosphi invcos3phi1(Qeps_et_dQeps,cos3phi);
    double E_et_dQeps = PoGrenoble(invcos3phi1,inv_n0);
    
    double Qeps_moins_dQeps = Qeps - delta;
    Invariant0QepsCosphi invcos3phi2(Qeps_moins_dQeps,cos3phi);
    double E_moins_dQeps = PoGrenoble(invcos3phi2,inv_n0);
    
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    double E_et_dV = PoGrenoble(invcos3phi,inv_et_dV);
    
    Invariant inv_moins_dV = inv_n0;inv_moins_dV.V -=  delta;
    double E_moins_dV = PoGrenoble(invcos3phi,inv_moins_dV);
    
    Invariant inv_et_dVdQeps = inv_n0;
    inv_et_dVdQeps.V += delta;
    double Qeps_et_dVdQeps = Qeps+delta;
    Invariant0QepsCosphi invcos3phi3(Qeps_et_dVdQeps,cos3phi);
    double E_et_dVdQeps = PoGrenoble(invcos3phi3,inv_et_dVdQeps);
    
    // cas des variations avec cos3phi    
    double cos3phi_et_dcos3phi = cos3phi+delta;
    Invariant0QepsCosphi invcos3phi_et_dcos3phi(Qeps,cos3phi_et_dcos3phi);
    double E_et_dcos3phi = PoGrenoble(invcos3phi_et_dcos3phi,inv_n0);
    
    double cos3phi_moins_dcos3phi = cos3phi-delta;
    Invariant0QepsCosphi invcos3phi_moins_dcos3phi(Qeps,cos3phi_moins_dcos3phi);
    double E_moins_dcos3phi = PoGrenoble(invcos3phi_moins_dcos3phi,inv_n0);
    
    Invariant0QepsCosphi invcos3phi_et_dcos3phi_a_dQeps(Qeps_et_dQeps,cos3phi_et_dcos3phi);
    double E_et_dcos3phi_a_dQeps = PoGrenoble(invcos3phi_et_dcos3phi_a_dQeps,inv_n0);
    
    double E_et_dcos3phi_a_dV = PoGrenoble(invcos3phi_et_dcos3phi,inv_et_dV);
    
    
    // calcul des dérivées premières
    double E_V = (E_et_dV  - E_n)/delta;
    double E_Qeps = (E_et_dQeps - E_n)/delta;
    double E_cos3phi = (E_et_dcos3phi -E_n)/delta;
    // calcul des dérivées secondes
    double E_V2 = (E_et_dV - 2.*E_n  + E_moins_dV ) /(delta * delta);
    double E_Qeps2 = (E_et_dQeps - 2.*E_n + E_moins_dQeps)/(delta * delta);
    double E_V_a_dQeps = (E_et_dVdQeps  - E_et_dQeps )/delta;
    double E_VQeps = ( E_V_a_dQeps - E_V)/delta;
    double E_cos3phi_2 = (E_et_dcos3phi - 2.*E_n  + E_moins_dcos3phi ) /(delta * delta);
    double E_dQeps_a_dcos3phi = (E_et_dcos3phi_a_dQeps  - E_et_dcos3phi )/delta;
    double E_dcos3phi_dQeps = (E_dQeps_a_dcos3phi  - E_Qeps )/delta;
    double E_dV_a_dcos3phi = (E_et_dcos3phi_a_dV  - E_et_dcos3phi )/delta;
    double E_dcos3phi_dV = (E_dV_a_dcos3phi  - E_V )/delta;

 // dans les dérivées secondes par rapport à cos3phi on a des erreurs, sans doute à cause des petites valeurs ??
  
    // comparaison avec les valeurs de dérivées analytiques
    int erreur = 0;
    if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
     {if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
        {if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = 1;}
      else 
       erreur = 2; } 
    if (diffpourcent(potret.EQ,E_Qeps,MaX(Dabs(E_Qeps),Dabs(potret.EQ)),0.05)) 
     {if (MiN(Dabs(E_Qeps),Dabs(potret.EQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps),Dabs(potret.EQ)) > 50.*delta) erreur = 3;}
      else 
        erreur = 4;}        
    if (diffpourcent(potret.Ecos3phi,E_cos3phi,MaX(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)),0.05)) 
     {if (MiN(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)) == 0.)
        {if ( MaX(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)) > 50.*delta) erreur = 31;}
      else 
        erreur = 41;}
    if (diffpourcent(potret.EQQ,E_Qeps2,MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)),0.05)) 
     {if (MiN(Dabs(E_Qeps2),Dabs(potret.EQQ)) == 0.)
        {if ( MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)) > 50.*delta) erreur = 5;}
      else 
        erreur = 6;}  
    if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05)) 
     {if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
        {if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = 7;}
      else 
        erreur = 8;}        
    if (diffpourcent(potret.Ecos3phi2,E_cos3phi_2,MaX(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)),0.05)) 
     {if (MiN(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)) == 0.)
        {if ( MaX(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)) > 50.*delta) erreur = 71;}
      else 
        erreur = 81;}
    if (diffpourcent(potret.EQV,E_VQeps,MaX(Dabs(E_VQeps),Dabs(potret.EQV)),0.05)) 
     {if (MiN(Dabs(E_VQeps),Dabs(potret.EQV)) == 0.)
      {if (Dabs(potret.EQV) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps)) erreur = 9;
          }
       else 
         if ( MaX(Dabs(E_VQeps),Dabs(potret.EQV)) > 150.*delta) erreur = 10;
        } 
      else 
        erreur = 11;
      }
      
    if (diffpourcent(potret.EVcos3phi,E_dcos3phi_dV,MaX(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)),0.05)) 
     {if (MiN(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)) == 0.)
      {if (Dabs(potret.EVcos3phi) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           /*Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps))*/ erreur = 91;
          }
       else 
         if ( MaX(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)) > 150.*delta) erreur = 101;
        } 
      else 
        erreur = 111;
      }
      
    if (diffpourcent(potret.EQcos3phi,E_dcos3phi_dQeps,MaX(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)),0.05)) 
     {if (MiN(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)) == 0.)
      {if (Dabs(potret.EQcos3phi) == 0.)
         // si la valeur de la dérivée numérique est nulle cela signifie peut-être
         // que l'on est à un extréma, la seule chose que le l'on peut faire est de
         // vérifier que l'on tend numériquement vers cette extréma ici un minima
         {
           /*Invariant inv_et_ddVddQeps = inv_n0;
           inv_et_ddVddQeps.V += delta/10.;
           double Qeps_et_ddVddQeps = Qeps+delta/10.;
           double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
           double Qeps_et_ddQeps = Qeps+delta/10.;
           double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
           double E_V_a_ddQeps = (E_et_ddVddQeps  - E_et_ddQeps )/(delta*10.);
           if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps))*/ erreur = 92;
          }
       else 
         if ( MaX(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)) > 150.*delta) erreur = 102;
        } 
      else 
        erreur = 112;
      }
      
        
    if (erreur > 0)
     { cout << "\n erreur dans le calcul analytique des derivees du potentiel";
       cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var(.."
            << " , numero d'increment = " << indic_Verif_PoGrenoble_et_var
            << " , numero d'erreur : " << erreur ;
 //      Sortie(1);
      } 
 };