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

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


//================== initialisation des variables static ======================
                       

// ---------- classe Hyper3D ---------

Hyper3D::Hyper3D ()  : // Constructeur par defaut
 HyperD(RIEN_COMP,RIEN_CATEGORIE_LOI_COMP,3,false)
 ,indic_Verif_Potentiel_et_var(0),limite_inf_Qeps(sqrt(2)*6.e-5)
 ,limite_inf_bIIb(36.e-10)
    {};
// Constructeur utile si l'identificateur du nom de la loi
// de comportement et le paramètre phase sont connus
Hyper3D::Hyper3D (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp,bool avec_ph) :
	 HyperD(id_compor,categorie_comp,3,avec_ph)
   ,indic_Verif_Potentiel_et_var(0),limite_inf_Qeps(sqrt(2)*6.e-5)
   ,limite_inf_bIIb(36.e-10)
	  {} ;
// Constructeur utile si l'identificateur du nom de la loi
// de comportement et le paramètre phase sont connus			
Hyper3D::Hyper3D (char* nom,Enum_categorie_loi_comp categorie_comp,bool avec_ph) :
 HyperD(nom,categorie_comp,3,avec_ph)
   ,indic_Verif_Potentiel_et_var(0),limite_inf_Qeps(sqrt(2)*6.e-5)
   ,limite_inf_bIIb(36.e-10)
	  {} ;
// Constructeur de copie
Hyper3D::Hyper3D (const Hyper3D& loi) :
  HyperD (loi)
  ,indic_Verif_Potentiel_et_var(loi.indic_Verif_Potentiel_et_var)
  ,limite_inf_Qeps(loi.limite_inf_Qeps)
  ,limite_inf_bIIb(loi.limite_inf_bIIb)
	{};



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

 // Calcul des invariants, de epsBH, 
  // retour de IdGBH qui pointe sur le bon tenseur
//++ TensBH * Invariants (TenseurBB&  epsBB_t,TenseurBB& gijBB_t,
TenseurBH * Hyper3D::Invariants (const TenseurBB&  epsBB_t,const TenseurBB& ,
           const TenseurHH & gijHH_t,const double& jacobien_0,const double& jacobien_t,
//++           Invariant & invariant,TensBH & epsBH);
           Invariant & inv,TenseurBH & epsBH_t)
  {
    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_t); // passage en dim 3
    Tenseur3BH & epsBH = *((Tenseur3BH*) &epsBH_t); // passage en dim 3
    const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_t); // passage en dim 3

    epsBH = epsBB * gijHH;  // deformation en mixte
    inv.Ieps = epsBH.Trace(); // trace de la déformation
    Tenseur3BH  eps_barreBH = epsBH - (inv.Ieps/3.) * IdBH3;
    inv.V = (jacobien_t/jacobien_0);
    // dans le cas de la prise en compte de la dilatation il faut changer V
    if (dilatation)
     // on considère que la variation de volume thermique = la trace du tenseur de def thermique
     {double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace();
      inv.V -= var_vol_thermique;
     };
    double inv_V = inv.V; // init pour être éventuellement modifié
    // s'il n'y a pas de régularisation : fact_regularisation = 0.
    if (Dabs(inv_V) < Dabs(fact_regularisation))
       inv_V += fact_regularisation;
    inv.bIIb = eps_barreBH.II() / 2.; // second invariant
//  la première méthode de calcul est la méthode traditionnelle, mais elle nécessite beaucoup
//  plus de calcul que la seconde (qui donne le même résultat à 10^{-18} près environ
//  une fois que V est calculé    
//    inv.bIIIb = eps_barreBH.III() / 3.; // troisième invariant 
    inv.bIIIb =  -0.125/(inv_V * inv_V ) +0.125  - 0.25*inv.Ieps
                 + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
                 - 0.5*inv.bIIb  +1./3. * inv.bIIb  * inv.Ieps;    // troisième invariant 
//    if (avec_regularisation)
//     {// inv.V += fact_regularisation;
//       inv.bIIb += fact_regularisation;
//       inv.bIIIb += fact_regularisation;
//     };
    return & IdBH3;
   };
   
   
// calcul des invariants et de leurs variations, de epsBH,
 // et de sa variation, puis  retour de IdGBH qui pointe sur le bon tenseur 
//++ TensBH *  Invariants_et_var (TenseurBB&  epsBB_tdt,
TenseurBH *  Hyper3D::Invariants_et_var (const TenseurBB&  epsBB_tdt,
           const TenseurBB& ,const Tableau <TenseurBB *>& d_gijBB_tdt,
           const TenseurHH & gijHH_tdt,const Tableau <TenseurHH *>& d_gijHH_tdt,
           const double& jacobien_0,const double& jacobien_tdt,const Vecteur& d_jacobien_tdt,
           InvariantVarDdl& inv,
//++           TensBH & epsBH_tdt,Tableau<TensBH> & depsBH_tdt); 
           TenseurBH & epsBH_tdt,Tableau<TenseurBH*> & depsBH_tdt) 
  {
    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
    Tenseur3BH & epsBH = *((Tenseur3BH*) &epsBH_tdt); // passage en dim 3
    const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // passage en dim 3

    epsBH = epsBB * gijHH;  // deformation en mixte
    inv.Ieps = epsBH.Trace(); // trace de la déformation
    Tenseur3BH  eps_barreBH = epsBH - (inv.Ieps/3.) * IdBH3;
    inv.V = (jacobien_tdt/jacobien_0);
    // dans le cas de la prise en compte de la dilatation il faut changer V
    if (dilatation)
     // on considère que la variation de volume thermique = la trace du tenseur de def thermique
     {double var_vol_thermique = (gijHH_tdt * (*epsBB_therm)).Trace();
      inv.V -= var_vol_thermique;
     };
    double inv_V = inv.V; // init pour être éventuellement modifié
    // s'il n'y a pas de régularisation : fact_regularisation = 0.
    if (Dabs(inv_V) < Dabs(fact_regularisation))
       inv_V += fact_regularisation;

    inv.bIIb = eps_barreBH.II() / 2.; // second invariant

//  la première méthode de calcul est la méthode traditionnelle, mais elle nécessite beaucoup
//  plus de calcul que la seconde (qui donne le même résultat à 10^{-18} près environ
//  une fois que V est calculé    
//    inv.bIIIb = eps_barreBH.III() / 3.; // troisième invariant
//    cout << "\n inv.bIIIb= "<< inv.bIIIb;
    inv.bIIIb =  -0.125/(inv_V * inv_V ) +0.125  - 0.25*inv.Ieps
                 + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
                 - 0.5*inv.bIIb  +1./3. * inv.bIIb  * inv.Ieps;    // troisième invariant 
//    cout << "  ,   "<< inv.bIIIb << " debug Hyper3D::Invariants_et_var " << endl;
//    if (avec_regularisation)
//     { //inv.V += fact_regularisation;
//       inv.bIIb += fact_regularisation;
//       inv.bIIIb += fact_regularisation;
//     };
   
// --- debug
///   pour info la bonne expression de bIIIb est donnée par tutu !!
// mais le résultat est identique (au erreur de calcul près) avec titi (c'est d'après la formule de denis
// thèse chapitre, II.2.3, formule 2.19
// et idem avec toto1 , en fait il n'y a que toto2 qui est complètement faux 
//
//   double toto2 = inv.bIIIb - (-1./(8.*inv.V * inv.V ) +1./8. - 0.25*inv.Ieps
//                               + 1./3. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps*11/54.
//                               -0.5*inv.bIIb  + inv.bIIb  * inv.Ieps);
//                               
//   double je = 1.-2./3.* inv.Ieps;
//   double je3=je*je*je;
//   double titi = 1./(inv.V*inv.V) - (je3 - 4.*inv.bIIb*je-8.*inv.bIIIb); 
//   double tutu = 1./(8.*inv.V*inv.V) - 1./8.+1./4.*inv.Ieps-1./6.*inv.Ieps*inv.Ieps
//                                    +1./27.*inv.Ieps*inv.Ieps*inv.Ieps 
//                                    + 0.5*inv.bIIb*1.-1./3.*inv.bIIb* inv.Ieps+inv.bIIIb;                        
//   double toto1 = inv.bIIIb -( -1./(8.*inv.V * inv.V ) +1./8. - 0.25*inv.Ieps
//                               + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
//                               - 0.5*inv.bIIb  +1./3. * inv.bIIb  * inv.Ieps);
//
//   if ((Abs(toto1)>ConstMath::petit) || (Abs(toto2)>ConstMath::petit)|| (Abs(titi)>ConstMath::petit)) 
//      cout << "\n ** toto1= " << toto1 << " toto2= " << toto2 << " titi= " << titi<< " tutu= " << tutu;
//
//// -- fin debug   
//   

    // cas le la variation 
    int nbddl = d_gijBB_tdt.Taille();
    Tenseur3BH  deps_barreBH; // grandeurs intermédiaires
    for (int i = 1; i<= nbddl; i++)
    {const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ;  // pour simplifier l'ecriture
     const Tenseur3BB & dgijBB = *((Tenseur3BB*)(d_gijBB_tdt(i)));  // passage en dim 3
     Tenseur3BH & depsBH = *((Tenseur3BH*)(depsBH_tdt(i)));  // passage en dim 3
     depsBH = 0.5 * dgijBB * gijHH_tdt + epsBB * dgijHH;
     double dIeps = (depsBH).Trace();inv.dIeps(i) = dIeps;
     deps_barreBH = depsBH - (dIeps/3.) * IdBH3; // en mixte les coordonnées de IdBH3 sont fixes
     double dbIIb = inv.dbIIb(i) = deps_barreBH && eps_barreBH;
     double dV = inv.dV(i) = (d_jacobien_tdt(i)/jacobien_0) ; 
     
//   Pour minimiser le temps de calcul on calcul la dérivée à l'aide de la formule scalaire 
//     inv.dbIIIb(i) = (deps_barreBH * eps_barreBH) && eps_barreBH;
//  -- formule scalaire de base 
//     inv.bIIIb =  -0.125/(inv.V * inv.V ) +0.125  - 0.25*inv.Ieps
//                 + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
//                 - 0.5*inv.bIIb  +1./3. * inv.bIIb  * inv.Ieps;    // troisième invariant 
//     d'où la variation
     inv.dbIIIb(i) =  0.25*dV/(inv_V * inv_V * inv_V)   - 0.25*dIeps
                     + 1./3. * dIeps * inv.Ieps - dIeps*inv.Ieps*inv.Ieps/9.
                     - 0.5*dbIIb  +1./3. * ( dbIIb  * inv.Ieps + inv.bIIb  * dIeps ) ;    // troisième invariant 

    };

// --- debug
//  affichage des invariants pour comparer
//   cout << "\n Hyper3D::Invariants_et_var ( V:" << inv.V <<" Ieps:" << inv.Ieps 
//	     << " bIIb:" << inv.bIIb << " bIIIb:" << inv.bIIIb 
//		  << "\n dV:" << inv.dV <<"\n dIeps:" << inv.dIeps 
//	     << "\n dbIIb:" << inv.dbIIb << "\n dbIIIb:" << inv.dbIIIb
//		  << endl;
//
//// -- fin debug   
       
    return & IdBH3;
   };
   
 // calcul des invariants  et de leurs variations par rapport aux déformations
 // et de sa variation, puis  retour de IdGBH qui pointe sur le bon tenseur 
TenseurBH *  Hyper3D::Invariants_et_varEps (const TenseurBB&  epsBB_tdt,
           const TenseurBB& ,const TenseurHH & gijHH_tdt,
           const double& jacobien_0,const double& jacobien_tdt,
           InvariantVarEps& inv,TenseurBH & epsBH_tdt) 
  {
    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
    Tenseur3BH & epsBH = *((Tenseur3BH*) &epsBH_tdt); // passage en dim 3
    const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // passage en dim 3

    epsBH = epsBB * gijHH;  // deformation en mixte
    inv.Ieps = epsBH.Trace(); // trace de la déformation
    Tenseur3BH  eps_barreBH = epsBH - (inv.Ieps/3.) * IdBH3;
    inv.V = (jacobien_tdt/jacobien_0);
    // dans le cas de la prise en compte de la dilatation il faut changer V
    if (dilatation)
     // on considère que la variation de volume thermique = la trace du tenseur de def thermique
     {double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace();
      inv.V -= var_vol_thermique;
     };
    double inv_V = inv.V; // init pour être éventuellement modifié
    // s'il n'y a pas de régularisation : fact_regularisation = 0.
    if (Dabs(inv_V) < Dabs(fact_regularisation))
       inv_V += fact_regularisation;

    inv.bIIb = eps_barreBH.II() / 2.; // second invariant

//  la première méthode de calcul est la méthode traditionnelle, mais elle nécessite beaucoup
//  plus de calcul que la seconde (qui donne le même résultat à 10^{-18} près environ
//  une fois que V est calculé    
    inv.bIIIb = eps_barreBH.III() / 3.; // troisième invariant
// cout << "\n inv.bIIIb= "<< inv.bIIIb;
    inv.bIIIb =  -0.125/(inv_V * inv_V ) +0.125  - 0.25*inv.Ieps
                 + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
                 - 0.5*inv.bIIb  +1./3. * inv.bIIb  * inv.Ieps;    // troisième invariant 
// cout << "  ,   "<< inv.bIIIb << " debug Hyper3D::Invariants_et_varEps " << endl;
//    if (avec_regularisation)
//     { //inv.V += fact_regularisation;
//       inv.bIIb += fact_regularisation;
//       inv.bIIIb += fact_regularisation;
//     };
    // cas le la variation
    Tenseur3HH epsHH = gijHH_tdt * epsBH;
    inv.dIeps_deps_HH = gijHH_tdt - 2. * epsHH;
    // on a : bIIb = bII - 1/6 Ieps^2 ; et dbII = epsHH - 2. * epsHH * epsBH
    Tenseur3HH d_bII_HH =  epsHH - 2. * epsHH * epsBH ;
    inv.dbIIb_deps_HH = d_bII_HH - (inv.Ieps/3.) * inv.dbIIb_deps_HH;
    // on a = bIIIb = bIII - 2/3 (I * bII - I^3/9) : thèse denis, page II.2.1
    Tenseur3HH d_bIII = epsHH * epsBH - 2. * epsHH * epsBH * epsBH;
    inv.dbIIIb_deps_HH = d_bIII + ((inv.Ieps * inv.Ieps / 3.)- (2.*inv.bIIb/3.) ) // mise en facteur du coeff 
                                                 * inv.dIeps_deps_HH            // scalaire
                                          - (2.*inv.Ieps/3.) * inv.dbIIb_deps_HH;
    // formule identique au cas de Mooney Rivlin , voir document de synthèse sur les lois hyper-élastiques chap 9                                      
    inv.dV_deps_HH = inv.V * gijHH_tdt;
    // on libère les tenseurs 
    LibereTenseur();LibereTenseurQ();
    return & IdBH3;
   };
             
    // calcul du potentiel et de ses dérivées non compris la phase
Hyper3D::PotenSansPhaseSansVar Hyper3D::Potentiel(const Invariant & inv,const double& jacobien0)
  { // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
    // mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
    // directement une technique de différence finie pour obtenir une dérivée numérique
    // def du retour
    PotenSansPhaseSansVar potret;
    // récup du conteneur spécifique du point, pour sauvegarde éventuelle
    SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
//    if ((Dabs(inv.bIIb) > limite_inf_bIIb)||(avec_regularisation))
    if (Dabs(inv.bIIb) > limite_inf_bIIb)
      { // calcul du nouvel  invariant sous la nouvelle forme
        InvariantQeps inv_spes = InvariantSpecifSansPhase(inv);
        // calcul du potentiel sous la forme polaire
        PoGrenobleSansPhaseSansVar pot = PoGrenoble(inv_spes,inv);
        // calcul de potret
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V *  pot.EV  + pot.E);
        potret.EbIIb = jacobien0 * inv.V * pot.EQ * inv_spes.dQepsdbIIb;
        // stockage éventuel des invariants et du potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= inv_spes.Qeps;
           save_resulHyperD.invP->cos3phi= 0.;
           save_resulHyperD.invP->potentiel= pot.E;
         };
       }
    else
      // cas d'un calcul par différence finie
      { // calcul des invariants sous la nouvelle forme
        double Qeps = Invariant0Qeps(inv);
        // calcul du potentiel sous la forme polaire au point actuel
        // ainsi que la variation suivant V, la phase ne peut pas être prise en compte
        // car Qeps est nulle
        PoGrenoble_V pot = PoGrenoble_et_V(Qeps,inv);
        // calcul du potentiel pour une petite variation suivant bIIb
        Invariant inv_n = inv; inv_n.bIIb += limite_inf_bIIb;
        double Qeps_n = Invariant0Qeps(inv_n);
        double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
        // stockage éventuel des invariants et potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv_n.V;
           save_resulHyperD.invP->Qeps= Qeps_n;
           save_resulHyperD.invP->cos3phi= 0.;
           save_resulHyperD.invP->potentiel= pot.E;
         };
        double E_n = PoGrenoble(Qeps_n,inv_n);
        // calcul de potret
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V * pot.EV  + pot.E);
        {double inter = (E_n - pot.E  );
         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
         // le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
         if (inter < ConstMath::trespetit)
            { potret.EbIIb = 0.; }
         else
            {potret.EbIIb = jacobien0 * inv.V * inter * unsurdelta;};
        };

       }
    return potret;             
   }; 

    // calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PotenAvecPhaseSansVar Hyper3D::PotentielPhase(const Invariant& inv,const double& jacobien0)
  { // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
    // mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
    // directement une technique de différence finie pour obtenir une dérivée numérique
        // def du retour
    PotenAvecPhaseSansVar potret;
    // récup du conteneur spécifique du point, pour sauvegarde éventuelle
    SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
//    if ((Dabs(inv.bIIb) > limite_inf_bIIb)||(avec_regularisation))
    if (Dabs(inv.bIIb) > limite_inf_bIIb)
      { // calcul des invariants sous la nouvelle forme
        InvariantQepsCosphi inv_spes = InvariantSpecif(inv);
        // calcul du potentiel sous la forme polaire
        PoGrenobleAvecPhaseSansVar pot = PoGrenoblePhase(inv_spes,inv);
        // stockage éventuel des invariants et potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= inv_spes.Qeps;
           save_resulHyperD.invP->cos3phi= inv_spes.cos3phi;
           save_resulHyperD.invP->potentiel= pot.E;
         };
        // calcul de potret
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V * ( pot.EV + pot.Ecos3phi * inv_spes.dcos3phidV) + pot.E);
        potret.EIeps = jacobien0 * inv.V * (pot.Ecos3phi * inv_spes.dcos3phidIeps);
        potret.EbIIb = jacobien0 * inv.V * (pot.EQ * inv_spes.dQepsdbIIb 
                       + pot.Ecos3phi * inv_spes.dcos3phidbIIb);
       }
    else
      // cas d'un calcul par différence finie, on considère que toute les dérivées suivant la phase
      // sont nulles dans ce cas
      { // calcul des invariants sous la nouvelle forme
        double Qeps = Invariant0Qeps(inv);
        // calcul du potentiel sous la forme polaire au point actuel
        // ainsi que la variation suivant V
        PoGrenoble_V pot = PoGrenoble_et_V(Qeps,inv);
        // calcul du potentiel pour une petite variation suivant bIIb
        Invariant inv_n = inv; inv_n.bIIb += limite_inf_bIIb;
        Invariant0QepsCosphi inv_spes_n = Invariant0Specif(inv_n);
        double E_n = PoGrenoble(inv_spes_n,inv_n);
        double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
        // stockage éventuel des invariants et du potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= inv_spes_n.Qeps;
           save_resulHyperD.invP->cos3phi= inv_spes_n.cos3phi;
           save_resulHyperD.invP->potentiel= pot.E;
         };
        // calcul de potret
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V * pot.EV  + pot.E);
        potret.EIeps = 0.;
//        potret.EbIIb = jacobien0 * inv.V * (-pot.E + E_n)/limite_inf_bIIb;
        {double inter = (E_n-pot.E );
         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
         // le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
         if (inter < ConstMath::trespetit)
            { potret.EbIIb = 0.; }
         else
            {potret.EbIIb = jacobien0 * inv.V * inter *unsurdelta;};
        };
       }
                  
    return potret;             
   }; 
    // calcul  du potentiel sans phase et dérivées avec  ses variations par rapport aux invariants
Hyper3D::PotenSansPhaseAvecVar Hyper3D::Potentiel_et_var(const Invariant& inv,const double& jacobien0 )
  { // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
    // mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
    // directement une technique de différence finie pour obtenir une dérivée numérique
    // def du retour
    PotenSansPhaseAvecVar potret;
    // récup du conteneur spécifique du point, pour sauvegarde éventuelle
    SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
//    if ((Dabs(inv.bIIb) > ConstMath::petit)||(avec_regularisation))
    if (Dabs(inv.bIIb) > limite_inf_bIIb)
      { // calcul des invariants sous la nouvelle forme
        Invariant2Qeps inv_spes = Invariant2SpecifSansPhase(inv);
        // calcul du potentiel sous la forme polaire
        PoGrenobleSansPhaseAvecVar pot = PoGrenoble_et_var(inv_spes,inv);
//debug
//cout << "\n debug: Hyper3D::"
//     << " E="<<pot.E << "  pot.EV=" << pot.EV << " pot.EVV="<<pot.EVV<<endl;
//if (pot.E < 0)
//   cout << "\n pb ";
//fin debug
        // calcul de potret
        //dérivées premières
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V *  pot.EV  + pot.E);
        potret.EbIIb = jacobien0 * inv.V * pot.EQ * inv_spes.dQepsdbIIb; 
        // dérivées secondes
        potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
        potret.EbIIb2 = jacobien0 * inv.V * ( pot.EQQ * inv_spes.dQepsdbIIb * inv_spes.dQepsdbIIb 
                        + pot.EQ * inv_spes.dQepsdbIIb2 ); 
        potret.EbIIbV = jacobien0 * ((pot.EQ + inv.V * pot.EQV) * inv_spes.dQepsdbIIb);
        // stockage éventuel des invariants et du potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= inv_spes.Qeps;
           save_resulHyperD.invP->cos3phi= 0.;
           save_resulHyperD.invP->potentiel= pot.E;
         };
        // appel d'une vérification éventuelle
//+++verif***        Verif_Potentiel_et_var(inv,jacobien0,potret);
       }
    else
      // cas d'un calcul par différence finie, on considère que toute les dérivées suivant la phase
      // sont nulles dans ce cas
      { // calcul des invariants sous la nouvelle forme
        Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb);
        double Qeps = Invariant0Qeps(inv_n0);
        // calcul du potentiel sous la forme polaire au point actuel
        // ainsi que la variation première et seconde suivant V
        PoGrenoble_VV pot = PoGrenoble_et_VV(Qeps,inv_n0);
        // calcul du potentiel pour une petite variation suivant bIIb
        Invariant inv_n = inv_n0;
        inv_n.bIIb += limite_inf_bIIb;
        double delta = limite_inf_bIIb; double unsurdelta = 1./delta;

        double Qeps_n = Invariant0Qeps(inv_n);
        // stockage éventuel des invariants et du potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= Qeps_n;
           save_resulHyperD.invP->cos3phi= 0.;
           save_resulHyperD.invP->potentiel= pot.E;
         };
        PoGrenoble_V pot_V = PoGrenoble_et_V(Qeps_n,inv_n);
        // et pour la dérivée seconde de nouveau une petite variation
        inv_n.bIIb += limite_inf_bIIb;
        Qeps_n = Invariant0Qeps(inv_n);
        double E_n = PoGrenoble(Qeps_n,inv_n);
        // calcul de potret
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V * pot.EV  + pot.E);
        {double inter = (pot_V.E - pot.E);
         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
         // le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
         if (inter < ConstMath::trespetit)
            { potret.EbIIb = 0.; }
         else
            {potret.EbIIb = jacobien0 * inv.V * inter *unsurdelta;};
        };
//        potret.EbIIb = jacobien0 * inv.V * (pot_V.E - pot.E)/limite_inf_bIIb;
        // dérivées secondes
        potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
        {double inter = ((E_n + pot.E) - 2.*pot_V.E );
         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
         // la dérivée seconde du potentiel doit-être elle aussi > 0 (au moins à l'origine)
         if (inter < ConstMath::trespetit)
            { potret.EbIIb2 = 0.; }
         else
            {potret.EbIIb2 = jacobien0 * inv.V * inter
                             * unsurdelta * unsurdelta;};
 //cout << "\n debug ((E_n + pot.E) - 2.*pot_V.E )= " << ((E_n + pot.E) - 2.*pot_V.E ) << endl;
        };
        // pour la variation de la dérivée / V, a priori on peut avoir les deux signes
        // donc on ne fait pas de filtrage
//        {double inter = (pot_V.EV - pot.EV);
//         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
//         if (inter < ConstMath::trespetit)
//            { potret.EbIIbV = potret.EbIIb / inv.V; }
//         else
//            {potret.EbIIbV = potret.EbIIb / inv.V
//                          + jacobien0 * inv.V * inter * unsurdelta;};
//        };
        potret.EbIIbV = potret.EbIIb / inv.V
              + jacobien0 * inv.V * (pot_V.EV - pot.EV) * unsurdelta;
//         potret.EbIIbV =  potret.EbIIb / inv.V
//                          + jacobien0 * inv.V * (pot_V.EV - pot.EV) / limite_inf_bIIb;

        // appel d'une vérification éventuelle
//+++verif***        Verif_Potentiel_et_var(inv,jacobien0,potret);
       }
    // retour                                 
    return potret;             
   }; 
    // calcul  du potentiel avec phase et dérivées avec  ses variations par rapport aux invariants
Hyper3D::PotenAvecPhaseAvecVar Hyper3D::PotentielPhase_et_var(const Invariant& inv,const double& jacobien0)
  { // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
    // mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
    // directement une technique de différence finie pour obtenir une dérivée numérique
    // def du retour
    PotenAvecPhaseAvecVar potret; // mise à 0 de tous les termes
    // récup du conteneur spécifique du point, pour sauvegarde éventuelle
    SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
// cout << "\n " << potret.E; // pour avoir a cces en debug ??
//    if ((Dabs(inv.bIIb) > limite_inf_bIIb)||avec_regularisation)
    if (Dabs(inv.bIIb) > limite_inf_bIIb)
//    if (Dabs(inv.bIIb) > limite_inf_Qeps*limite_inf_Qeps)
      { // calcul des invariants sous la nouvelle forme
        Invariant2QepsCosphi inv_spes = Invariant2Specif(inv);
        // calcul du potentiel sous la forme polaire
        PoGrenobleAvecPhaseAvecVar pot = PoGrenoblePhase_et_var(inv_spes,inv);
        // stockage éventuel des invariants et du potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= inv_spes.Qeps;
           save_resulHyperD.invP->cos3phi= inv_spes.cos3phi;
           save_resulHyperD.invP->potentiel= pot.E;
         };
//cout << "\n calcul normal";
//debug
//cout << "\n debug: Hyper3D::"
//     << " E="<<pot.E << "  pot.EV=" << pot.EV << " pot.EVV="<<pot.EVV<<endl;
//if (pot.E < 0)
//   cout << "\n pb ";
//fin debug
        // calcul de potret
        //dérivées premières
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        potret.EV = jacobien0 * (inv.V * ( pot.EV + pot.Ecos3phi * inv_spes.dcos3phidV) + pot.E);
        potret.EIeps = jacobien0 * inv.V * (pot.Ecos3phi * inv_spes.dcos3phidIeps);
        potret.EbIIb = jacobien0 * inv.V * (pot.EQ * inv_spes.dQepsdbIIb 
                      + pot.Ecos3phi * inv_spes.dcos3phidbIIb);
        // dérivées secondes
        double dcos3phidV_2 = inv_spes.dcos3phidV * inv_spes.dcos3phidV; // simplification
        double dcos3phidIeps_2 = inv_spes.dcos3phidIeps * inv_spes.dcos3phidIeps; //   "             
        double dcos3phidbIIb_2 = inv_spes.dcos3phidbIIb * inv_spes.dcos3phidbIIb; //   "             
        potret.EVV = jacobien0 * (2. * ( pot.Ecos3phi + inv.V * pot.EVcos3phi) * inv_spes.dcos3phidV
                     + 2. * pot.EV + inv.V * pot.EVV + inv.V * pot.Ecos3phi2 * dcos3phidV_2 
                     + inv.V * pot.Ecos3phi * inv_spes.dcos3phidV2);
        potret.EIeps2 = jacobien0 * inv.V * ( pot.Ecos3phi2 * dcos3phidIeps_2 + pot.Ecos3phi * inv_spes.dcos3phidIeps2);
        potret.EbIIb2 = jacobien0 * inv.V * ( pot.EQQ * inv_spes.dQepsdbIIb * inv_spes.dQepsdbIIb  
                        + 2. * pot.EQcos3phi * inv_spes.dQepsdbIIb * inv_spes.dcos3phidbIIb
                        + pot.Ecos3phi2 * dcos3phidbIIb_2
                        + pot.EQ * inv_spes.dQepsdbIIb2 +  pot.Ecos3phi * inv_spes.dcos3phidbIIb2); 
        potret.EIepsV = jacobien0 * (( pot.Ecos3phi + inv.V * pot.EVcos3phi) * inv_spes.dcos3phidIeps
                        + inv.V * pot.Ecos3phi2 * inv_spes.dcos3phidV * inv_spes.dcos3phidIeps);
        potret.EbIIbV = jacobien0 * ((pot.EQ + inv.V * pot.EQV) * inv_spes.dQepsdbIIb 
                        + (pot.Ecos3phi + inv.V * pot.EVcos3phi ) * inv_spes.dcos3phidbIIb
                        + inv.V * pot.EQcos3phi * inv_spes.dcos3phidV * inv_spes.dQepsdbIIb
                        + inv.V * pot.Ecos3phi2 * inv_spes.dcos3phidV * inv_spes.dcos3phidbIIb
                        + inv.V * pot.Ecos3phi * inv_spes.dcos3phidVdbIIb);
        potret.EbIIbIeps = jacobien0 * inv.V * ( pot.EQcos3phi * inv_spes.dcos3phidIeps * inv_spes.dQepsdbIIb
                        + pot.Ecos3phi2 * inv_spes.dcos3phidIeps * inv_spes.dcos3phidbIIb 
                        + pot.Ecos3phi * inv_spes.dcos3phidIepsdbIIb);

//   PotenAvecPhaseAvecVar toto = potret; // pour le debug
//   cout << "\n " <<  toto.E << " " << toto.EV; 
//  cout << "\n cas analytique ";
// --- pour le débug   
//   switch(fpclassify(toto.EV))
//    { case FP_NAN: case FP_INFINITE:
//        { // cas d'un nombre infini 
//          cout << "\n erreur valeur nan pour EV! ";
//          cout << " poret= " << potret.E; // pour que potret existe pour le débug
//        };
//        break;
//    };
    
// --- fin débug
        // appel d'une vérification éventuelle
//      Verif_Potentiel_et_var(inv,jacobien0,potret);
      }
      
// ---------------------------------------------------------------
// non, le calcul par différence finies vraiment ne marche pas !!! c'est la merde
// car le calcul normale est extrèmement sensible aux limites initiales limite_inf_Qeps et limite_inf_bIIb 
// la ça marche pour le cu mais peut-être pour un autre matériau ça ne marchera pas !!! 
// peut-être mettre en paramètre, mais c'est vraiment difficile à utiliser 
// bref à suivre ...
/*       
    else if (Dabs(inv.bIIb) > limite_inf_bIIb /10.)
      { // calcul des dérivées par différences finies
  cout << "\n cas difference fini avec phase ";
    
    Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl 
    Invariant0QepsCosphi is_sp0 = Invariant0Specif(inv_n0);
    double E_n = PoGrenoble(is_sp0,inv_n0);  // la valeur du potentiel de référence
    // ici on est avec phase donc trois invariants indépendant V et Qeps et cos3phi,
    // 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/100.;
    Invariant inv_et_dbIIb = inv_n0;inv_et_dbIIb.bIIb += delta;
    Invariant0QepsCosphi is_et_dbIIb = InvariantDe_V_bIIb_Ieps(inv_et_dbIIb);
    double E_et_dbIIb = PoGrenoble(is_et_dbIIb,inv_et_dbIIb);
    Invariant inv_et_dbIIb2 = inv_n0;inv_et_dbIIb2.bIIb += 2. * delta;
    Invariant0QepsCosphi is_et_dbIIb2 = InvariantDe_V_bIIb_Ieps(inv_et_dbIIb2);
    double E_et_dbIIb2 = PoGrenoble(is_et_dbIIb2,inv_et_dbIIb2);
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    Invariant0QepsCosphi is_et_dV = InvariantDe_V_bIIb_Ieps(inv_et_dV);
    double E_et_dV = PoGrenoble(is_et_dV,inv_et_dV);
    Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
    Invariant0QepsCosphi is_et_dV2 = InvariantDe_V_bIIb_Ieps(inv_et_dV2);
    double E_et_dV2 = PoGrenoble(is_et_dV2,inv_et_dV2);
    Invariant inv_et_dVdbIIb = inv_n0;
    inv_et_dVdbIIb.V += delta;
    inv_et_dVdbIIb.bIIb += delta;
    Invariant0QepsCosphi is_et_dVdbIIb = InvariantDe_V_bIIb_Ieps(inv_et_dVdbIIb);
    double E_et_dVdbIIb = PoGrenoble(is_et_dVdbIIb,inv_et_dVdbIIb);
    // cas avec Ieps
    
    Invariant inv_et_dIeps = inv_n0;inv_et_dIeps.Ieps += delta;
    Invariant0QepsCosphi is_et_dIeps = InvariantDe_V_bIIb_Ieps(inv_et_dIeps);
    double E_et_dIeps = PoGrenoble(is_et_dIeps,inv_et_dIeps);
    
    Invariant inv_et_dIeps2 = inv_n0;inv_et_dIeps2.Ieps += 2. * delta;
    Invariant0QepsCosphi is_et_dIeps2 = InvariantDe_V_bIIb_Ieps(inv_et_dIeps2);
    double E_et_dIeps2 = PoGrenoble(is_et_dIeps2,inv_et_dIeps2);
    
    Invariant inv_et_dIepsdbIIb = inv_n0;
    inv_et_dIepsdbIIb.Ieps += delta;
    inv_et_dIepsdbIIb.bIIb += delta;
    Invariant0QepsCosphi is_et_dIepsdbIIb = InvariantDe_V_bIIb_Ieps(inv_et_dIepsdbIIb);
    double E_et_dIepsdbIIb = PoGrenoble(is_et_dIepsdbIIb,inv_et_dIepsdbIIb);
    
    Invariant inv_et_dIepsdV = inv_n0;
    inv_et_dIepsdV.Ieps += delta;
    inv_et_dIepsdV.V += delta;
    Invariant0QepsCosphi is_et_dIepsdV = InvariantDe_V_bIIb_Ieps(inv_et_dIepsdV);
    double E_et_dIepsdV = PoGrenoble(is_et_dIepsdV,inv_et_dIepsdV);
    
    
    // calcul des dérivées premières
    potret.EV = jacobien0 *(E_et_dV * inv_et_dV.V - E_n * inv.V)/delta;
    potret.EbIIb = (jacobien0 * inv.V)*(E_et_dbIIb - E_n)/delta;
    potret.EIeps = (jacobien0 * inv.V)*(E_et_dIeps  - E_n )/delta;
    // calcul des dérivées secondes hors phases
    potret.EVV = jacobien0 *(E_et_dV2 * inv_et_dV2.V - 2.*E_et_dV * inv_et_dV.V + E_n * inv.V)
                   /(delta * delta);
    potret.EbIIb2 = (jacobien0 * inv.V)*(E_et_dbIIb2 - 2.*E_et_dbIIb + E_n) /(delta * delta);

    potret.EIeps2 = (jacobien0 * inv.V) *(E_et_dIeps2  - 2.*E_et_dIeps  + E_n )/(delta * delta);

    double E_V_a_dbIIb = jacobien0 * (E_et_dVdbIIb * inv_et_dVdbIIb.V - E_et_dbIIb * inv_et_dbIIb.V) /delta;
    potret.EbIIbV = ( E_V_a_dbIIb - potret.EV)/delta;
    
    double E_Ieps_a_dbIIb = (jacobien0 * inv.V) * (E_et_dIepsdbIIb  - E_et_dbIIb) /delta;
    potret.EbIIbIeps = ( E_Ieps_a_dbIIb - potret.EIeps)/delta;
    
    double E_Ieps_a_dV = jacobien0 * (E_et_dIepsdV * inv_et_dIepsdV.V - E_et_dV * inv.V) /delta;
    potret.EIepsV = ( E_Ieps_a_dV - potret.EIeps)/delta;
    
    		
    	}*/
/*    else if (Dabs(inv.bIIb) > ConstMath::petit)
      { // ici on calcul comme si on n'avait pas de phase
        // calcul des invariants sous la nouvelle forme
        Invariant2Qeps inv_spes = Invariant2SpecifSansPhase(inv);
        // calcul du potentiel sous la forme polaire
        PoGrenobleSansPhaseAvecVar pot = PoGrenoble_et_var(inv_spes,inv);
//  cout << "\n cas analytique sans phase ";
        // calcul de potret
        //dérivées premières
        potret.E = pot.E * inv.V * jacobien0;
        potret.EV = jacobien0 * (inv.V *  pot.EV  + pot.E);
        potret.EbIIb = jacobien0 * inv.V * pot.EQ * inv_spes.dQepsdbIIb; 
        // dérivées secondes
        potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
        potret.EbIIb2 = jacobien0 * inv.V * ( pot.EQQ * inv_spes.dQepsdbIIb * inv_spes.dQepsdbIIb 
                        + pot.EQ * inv_spes.dQepsdbIIb2 ); 
        potret.EbIIbV = jacobien0 * ((pot.EQ + inv.V * pot.EQV) * inv_spes.dQepsdbIIb);
        // appel d'une vérification éventuelle
//+++verif***        Verif_Potentiel_et_var(inv,jacobien0,potret);
       }*/
    else  // dernier cas avec un inv.bIIb vraiment petit
      // cas d'un calcul par différence finie, on considère que toute les dérivées suivant la phase
      // sont nulles dans ce cas, car de toutes manière on ne peut pas calculer la phase
      { // calcul des invariants sous la nouvelle forme
//  cout << "\n cas difference fini sans phase ";
        Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb);
        double Qeps = Invariant0Qeps(inv_n0);
        // calcul du potentiel sous la forme polaire au point actuel
        // ainsi que la variation première et seconde suivant V
        PoGrenoble_VV pot = PoGrenoble_et_VV(Qeps,inv_n0);
        // calcul du potentiel pour une petite variation suivant bIIb
        Invariant inv_n = inv_n0;
        double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
        inv_n.bIIb += delta;
        Invariant0QepsCosphi inv_spes_n = Invariant0Specif(inv_n);
        // stockage éventuel des invariants et du potentiel
        if (sortie_post)
         { save_resulHyperD.invP->V= inv.V;
           save_resulHyperD.invP->Qeps= inv_spes_n.Qeps;
           save_resulHyperD.invP->cos3phi= inv_spes_n.cos3phi;
           save_resulHyperD.invP->potentiel= pot.E;
         };
        PoGrenoble_V pot_V = PoGrenoble_et_V(inv_spes_n,inv_n);
        // et pour la dérivée seconde de nouveau une petite variation
        inv_n.bIIb += delta;
        inv_spes_n = Invariant0Specif(inv_n);
        double E_n = PoGrenoble(inv_spes_n,inv_n);
        // calcul de potret
        potret.E = pot.E * inv.V * jacobien0;
        potret.Ks = pot.Ks;
        // dérivées premières
        potret.EV = jacobien0 * (inv.V * pot.EV  + pot.E);
        potret.EbIIb = jacobien0 * inv.V * ( pot_V.E - pot.E)*unsurdelta;
        {double inter = (pot_V.E - pot.E);
         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
         // le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
         if (inter < ConstMath::trespetit)
            { potret.EbIIb = 0.; }
         else
            {potret.EbIIb = jacobien0 * inv.V * inter * unsurdelta;};
        };

        // dérivées secondes
        potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
        {double inter = ((E_n + pot.E) - 2.*pot_V.E ); //(E_n - 2.*pot_V.E + pot.E)
         // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
         // la dérivée seconde du potentiel doit-être elle aussi > 0 (au moins à l'origine)
         if (inter < ConstMath::trespetit)
            { potret.EbIIb2 = 0.; }
         else
            {potret.EbIIb2 = jacobien0 * inv.V * inter
                             * unsurdelta * unsurdelta;};
       //cout << "\n debug ((E_n + pot.E) - 2.*pot_V.E )= " << ((E_n + pot.E) - 2.*pot_V.E ) << endl;
        }
        // pour la variation de la dérivée / V, a priori on peut avoir les deux signes
        // donc on ne fait pas de filtrage
        potret.EbIIbV =  potret.EbIIb / inv.V
                         + jacobien0 * inv.V * (pot_V.EV - pot.EV) *unsurdelta;
        // pour les dérivées selon l'angle tous les termes sont considéré nul
        potret.EIeps = 0.;
        potret.EIeps2 = 0.; potret.EIepsV = 0.; potret.EbIIbIeps = 0.;         
        // appel d'une vérification éventuelle
//        Verif_Potentiel_et_var(inv,jacobien0,potret);
      };
    // retour 
    
// --- debug --------
//cout << "\n E= " << potret.E << " EV= " << potret.EV << " EbIIb= " << potret.EbIIb << " EIeps= " << potret.EIeps
//     << "\n EVV= " << potret.EVV << " EbIIb2= " << potret.EbIIb2 << " EIeps2= " << potret.EIeps2
//     << "\n EbIIbIeps = " << potret.EbIIbIeps << " EIepsV= " << potret.EIepsV  << " EbIIbV= " << potret.EbIIbV << endl;

// ---- fin débug -------    
     
    return potret;             
   }; 
   
    // calcul de l'invariant Qeps seul valable même si Qeps est nul
double Hyper3D::Invariant0Qeps(const Invariant & inv)
 { // les nouveaux invariants
   double Qeps =  sqrt(2. * inv.bIIb);
//   if (avec_regularisation)
//    {Qeps += fact_regularisation;}
   return Qeps;
  };
  
    // calcul des invariants Qeps,cos3phi seuls
Hyper3D::Invariant0QepsCosphi Hyper3D::Invariant0Specif(const Invariant & inv)
 { Hyper3D::Invariant0QepsCosphi inv_ret; // def du retour
   // les nouveaux invariants
   inv_ret.Qeps = sqrt(2. * inv.bIIb);
   double inv_ret_Qeps = inv_ret.Qeps; // init pour être éventuellement modifié
     // ensuite avec la régularisation
   if (avec_regularisation)
    {inv_ret_Qeps += fact_regularisation;}
   double unsurQeps = 1./inv_ret_Qeps;  // pour simplifier
   // --- angle de phase
   double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
   double troisracine6 = 3. * sqrt(6.);
   inv_ret.cos3phi = troisracine6 * inv.bIIIb * unsurQeps3;
   if (inv_ret.cos3phi > 1.) inv_ret.cos3phi=1.;
   if (inv_ret.cos3phi < -1.) inv_ret.cos3phi=-1.;
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::Invariant0Specif(" << endl;
//fin debug
   // retour
   return inv_ret;                     
  };

// calcul des invariants Qeps,cos3phi en fonction de V, bIIb Ieps
// on utilise pas l'invariant bIIIb, 
Hyper3D::Invariant0QepsCosphi Hyper3D::InvariantDe_V_bIIb_Ieps(const Invariant & inv)
 { Hyper3D::Invariant0QepsCosphi inv_ret; // def du retour
///******** en fait ça merdouille, le bIIIb calculé donne n'importe quoi quand c'est petit donc
/// a éviter !!! 
/// a défaut d'autre chose pour l'instant on utilise  Invariant0Specif
///
   // les nouveaux invariants
   inv_ret.Qeps = sqrt(2. * inv.bIIb); 
   double inv_ret_Qeps = inv_ret.Qeps; // init pour être éventuellement modifié
   double inv_V = inv.V; // init pour être éventuellement modifié
     // ensuite avec la régularisation
   if (avec_regularisation)
    {inv_ret_Qeps += fact_regularisation;
     inv_V += fact_regularisation;
    };
   double unsurQeps = 1./inv_ret_Qeps;  // pour simplifier
   // on recalcul bIIIb
   double bIIIb=-1./(8.*inv_V*inv_V) + 1./8.-1./4.*inv.Ieps+1./6.*inv.Ieps*inv.Ieps
                                    -1./27.*inv.Ieps*inv.Ieps*inv.Ieps 
                                    - 0.5*inv.bIIb+1./3.*inv.bIIb* inv.Ieps;   
   // --- angle de phase avec la valeur locale de bIIIb
   double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
   double troisracine6 = 3. * sqrt(6.);
   inv_ret.cos3phi = troisracine6 * bIIIb * unsurQeps3;
   // on limite les valeurs au normal
   if (inv_ret.cos3phi > 1.) { inv_ret.cos3phi = 1.;}
   else if (inv_ret.cos3phi < -1.) { inv_ret.cos3phi = -1.;};
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::InvariantDe_V_bIIb_Ieps(" << endl;
//fin debug
   // retour
   return inv_ret;                     
  };
    

    // calcul de l'invariants Qeps ainsi que sa variation par rapport
    // à  bIIb (Qeps doit être non nul)
Hyper3D::InvariantQeps Hyper3D::InvariantSpecifSansPhase(const Invariant & inv)
 { Hyper3D::InvariantQeps inv_ret; // def du retour
   // tout d'abord le nouvel invariant
   inv_ret.Qeps = sqrt(2. * inv.bIIb); 
   // maintenant la variation
   if (avec_regularisation)
    {inv_ret.dQepsdbIIb = 1./(inv_ret.Qeps+fact_regularisation);}
   else
    {inv_ret.dQepsdbIIb = 1./inv_ret.Qeps;};
   // retour
   return inv_ret;                     
  };

    // calcul des invariants Qeps,cos3phi, ainsi que leur variation par rapport
    // aux trois invariants Ieps, V et bIIb
Hyper3D::InvariantQepsCosphi Hyper3D::InvariantSpecif(const Invariant & inv)
 { Hyper3D::InvariantQepsCosphi inv_ret; // def du retour, init à 0
   // tout d'abord les nouveaux invariants
   inv_ret.Qeps = sqrt(2. * inv.bIIb); 
   double inv_V = inv.V; // init pour être éventuellement modifié
   double unsurQeps; // init = 1./inv_ret.Qeps;  // pour simplifier
   if (avec_regularisation)
        {unsurQeps = 1./(inv_ret.Qeps+fact_regularisation);
         inv_V += fact_regularisation;
        }
   else {unsurQeps = 1./inv_ret.Qeps;};
   // maintenant les variations / pour Qeps
   double untiers = 1./3.;
   inv_ret.dQepsdbIIb = unsurQeps;
   // --- angle de phase
   // dans le cas ou Qeps est trop petit on laisse les 0
//   if ((inv_ret.Qeps >= limite_inf_Qeps)||(avec_regularisation))
   if (inv_ret.Qeps >= limite_inf_Qeps)
    { double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
      double unsurQeps5 = unsurQeps3 * unsurQeps * unsurQeps;
      double troisracine6 = 3. * sqrt(6.);
      double neufracine6 = 3. * troisracine6;
      inv_ret.cos3phi = troisracine6 * inv.bIIIb * unsurQeps3;
      if (inv_ret.cos3phi > 1.) inv_ret.cos3phi=1.;
      if (inv_ret.cos3phi < -1.) inv_ret.cos3phi=-1.;
      // puis les variations
      inv_ret.dcos3phidV = troisracine6 / 4. * unsurQeps3 /(inv_V * inv_V * inv_V);
      inv_ret.dcos3phidIeps = troisracine6  * unsurQeps3 * (-0.25 + inv.bIIb*untiers 
                              + inv.Ieps*untiers - inv.Ieps * inv.Ieps/9.);
      inv_ret.dcos3phidbIIb = - neufracine6  * unsurQeps5 * inv.bIIIb
                             + troisracine6  * unsurQeps3 * (inv.Ieps*untiers - 0.5); 
    };
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::InvariantSpecif(" << endl;
//fin debug
   // retour
   return inv_ret;                                       
  };

    // calcul de l'invariants Qeps ainsi que ses variations première et seconde
    // par rapport à bIIb (Qeps doit être non nul)
Hyper3D::Invariant2Qeps Hyper3D::Invariant2SpecifSansPhase(const Invariant & inv)
 { Hyper3D::Invariant2Qeps inv_ret; // def du retour et init à 0
   // tout d'abord les nouveaux invariants
   inv_ret.Qeps = sqrt(2. * inv.bIIb); 
   double unsurQeps; // init = 1./inv_ret.Qeps;  // pour simplifier
   if (avec_regularisation)
        {unsurQeps = 1./(inv_ret.Qeps+fact_regularisation);}
   else {unsurQeps = 1./inv_ret.Qeps;};
//   double unsurQeps = 1./inv_ret.Qeps;  // pour simplifier
   // maintenant les variations
   inv_ret.dQepsdbIIb = unsurQeps;
   // et les variations secondes
   inv_ret.dQepsdbIIb2 = -1. * unsurQeps * unsurQeps * unsurQeps; 
   // retour
   return inv_ret;                     
  };


    // calcul des invariants Qeps,cos3phi, ainsi que leur variation première et seconde
    // par rapport aux trois invariants Ieps, V et bIIb
Hyper3D::Invariant2QepsCosphi Hyper3D::Invariant2Specif(const Invariant & inv)
 { Hyper3D::Invariant2QepsCosphi inv_ret; // def du retour, init à 0
   // tout d'abord les nouveaux invariants
   inv_ret.Qeps = sqrt(2. * inv.bIIb); 
   double unsurQeps; // init = 1./inv_ret.Qeps;  // pour simplifier
   double inv_V = inv.V; // init pour être éventuellement modifié
   if (avec_regularisation)
        {unsurQeps = 1./(inv_ret.Qeps+fact_regularisation);
         inv_V += fact_regularisation;
        }
   else {unsurQeps = 1./inv_ret.Qeps;};
//   double unsurQeps = 1./inv_ret.Qeps;  // pour simplifier
   // et les variations / Qeps
   double untiers = 1./3.;
   inv_ret.dQepsdbIIb = unsurQeps;
   double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
   inv_ret.dQepsdbIIb2 = -1. * unsurQeps3; 
   // --- angle de phase
   // dans le cas ou Qeps est trop petit on laisse les 0
//   if ((inv_ret.Qeps >= limite_inf_Qeps)||(avec_regularisation))
   if (inv_ret.Qeps >= limite_inf_Qeps)
    { double unsurQeps5 = unsurQeps3 * unsurQeps * unsurQeps;
      double unsurQeps7 = unsurQeps5 * unsurQeps * unsurQeps;
      double racine6=sqrt(6.);
      double troisracine6 = 3. * racine6;
      double neufracine6 = 3. * troisracine6;
      inv_ret.cos3phi = troisracine6 * inv.bIIIb * unsurQeps3;
      if (inv_ret.cos3phi > 1.) inv_ret.cos3phi=1.;
      if (inv_ret.cos3phi < -1.) inv_ret.cos3phi=-1.;
     
      double unsurV3 = 1. / (inv_V * inv_V * inv_V);
      
      // maintenant les variations
      inv_ret.dcos3phidV = troisracine6 * unsurQeps3 * unsurV3 / 4. ;
      inv_ret.dcos3phidIeps = troisracine6  * unsurQeps3 * (-0.25 + inv.bIIb*untiers
                              + inv.Ieps*untiers - inv.Ieps * inv.Ieps/9.);
      inv_ret.dcos3phidbIIb = - neufracine6  * unsurQeps5 * inv.bIIIb
                         + troisracine6  * unsurQeps3 * (inv.Ieps*untiers - 0.5);
      // les variations secondes
      inv_ret.dcos3phidV2 = - neufracine6 * unsurQeps3 * unsurV3 / (4. * inv_V);
      inv_ret.dcos3phidIeps2 = troisracine6 * unsurQeps3 * (1./3. - 2./9. * inv.Ieps);
      inv_ret.dcos3phidbIIb2 = 15. * troisracine6 * inv.bIIIb * unsurQeps7
                              - 2.*neufracine6  * unsurQeps5 * (untiers * inv.Ieps - 0.5);
      inv_ret.dcos3phidVdbIIb = - neufracine6 * unsurQeps5 * unsurV3 * 0.25;
      inv_ret.dcos3phidIepsdbIIb = racine6 * unsurQeps3
                                   - neufracine6 * unsurQeps5
                               * (-0.25 + inv.bIIb*untiers + untiers * inv.Ieps
                                  - 1./9. * inv.Ieps * inv.Ieps );
    };
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::Invariant2Specif(" << endl;
//fin debug
   // retour
   return inv_ret;                     
  };

// Lecture des paramètre specifique sur fichier
void Hyper3D::LectParaSpecifiques(UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D,LesFonctions_nD& lesFonctionsnD)
  { string nom_class_methode("Hyper3D::LectParaSpecifiques");
    
    // cas limite_inf_Qeps (variable stockée dans Hyper3D.h)
    string nom1;
    if (strstr(entreePrinc->tablcar,"limite_inf_Qeps_")!=NULL)
     { *(entreePrinc->entree) >> nom1 ;
       if (nom1 != "limite_inf_Qeps_")
         { cout << "\n erreur en lecture du drapeau de limite_inf_Qeps, on attendait le mot cle "
             << " limite_inf_Qeps_ on a lue: " << nom1 ;
           entreePrinc->MessageBuffer("**Hyper3D::LectParaSpecifiques((.....**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
       *(entreePrinc->entree) >> limite_inf_Qeps;
     };
   // cas limite_inf_Qeps (variable stockée dans Hyper3D.h)
   if (strstr(entreePrinc->tablcar,"limite_inf_bIIb_")!=NULL)
    { *(entreePrinc->entree) >> nom1 ;
      if (nom1 != "limite_inf_bIIb_")
        { cout << "\n erreur en lecture du drapeau de limite_inf_bIIb, on attendait le mot cle "
            << " limite_inf_bIIb_ on a lue: " << nom1 ;
          entreePrinc->MessageBuffer("**Hyper3D::LectParaSpecifiques((.....**");
          throw (UtilLecture::ErrNouvelleDonnee(-1));
          Sortie(1);
        };
      *(entreePrinc->entree) >> limite_inf_bIIb;
    };
    
    // cas avec régularisation (variable stockée dans Hyper3D.h)
    avec_regularisation = false; // par défaut
    if (strstr(entreePrinc->tablcar,"avec_regularisation_")!=NULL)
     { *(entreePrinc->entree) >> nom1 ;
       if (nom1 != "avec_regularisation_")
         { cout << "\n erreur en lecture du drapeau de regularisation, on attendait le mot cle "
             << " avec_regularisation_ on a lue: " << nom1 ;
           entreePrinc->MessageBuffer("**Hyper3D::LectParaSpecifiques((.....**");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);
         };
       *(entreePrinc->entree) >> fact_regularisation;
       avec_regularisation=true;
     };
    // lecture de l'indication éventuelle du post traitement
    string le_mot_cle = "sortie_post_";
    entreePrinc->Lecture_un_parametre_int(0,nom_class_methode,0,1,le_mot_cle,sortie_post);
    
  };

// I/O des paramètres spécifiques: ex: régularisation, limite inf ...
void Hyper3D::Ecriture_para_specifiques(ofstream& sort,const int cas)
   { if (cas == 1)
      {// la régularisation éventuelle
       if (avec_regularisation)
         sort << " avec_regularisation " << fact_regularisation;
       else
         sort << " sans_regularisation ";
       // limites
       sort << " limite_inf_Qeps= " << limite_inf_Qeps
            << " limite_inf_bIIb= " << limite_inf_bIIb << " " ;
      };
   };
void Hyper3D::Lecture_para_specifiques(ifstream& ent,const int cas)
   { if (cas == 1)
      {// la régularisation éventuelle
       string nom;
       ent >> nom;
       if (nom == "avec_regularisation")
         {ent >> fact_regularisation;avec_regularisation=true;}
       else if (nom == "sans_regularisation")
         avec_regularisation=false;
       else
        {cout << "\n erreur en lecture du facteur de regularisation (loi hyperelastique) "
              << " on a lue: " << nom << " au lieu du mot cle: avec_regularisation ou  sans_regularisation";
         Sortie(1);
        };
       // limites
       ent >> nom;
       if (nom == "limite_inf_Qeps=")
         ent >> limite_inf_Qeps;
       else
        {cout << "\n erreur en lecture de limite_inf_Qeps (loi hyperelastique) "
              << " on a lue: " << nom << " au lieu du mot cle: limite_inf_Qeps= ";
         Sortie(1);
        };
       ent >> nom;
       if (nom == "limite_inf_bIIb=")
         ent >> limite_inf_bIIb;
       else
        {cout << "\n erreur en lecture de limite_inf_bIIb (loi hyperelastique) "
              << " on a lue: " << nom << " au lieu du mot cle: limite_inf_bIIb= ";
         Sortie(1);
        };
      };
   };

///=============== 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 Hyper3D::Verif_Potentiel_et_var(const InvariantVarDdl& inv,const double& jacobien0
                                     ,const PotenSansPhaseAvecVar& potret )
  { // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_Potentiel_et_var == 0)
     { cout << "\n ****verification des derivees du potentiels par rapport aux invariants****";
       cout << "\n Hyper3D::Verif_Potentiel_et_var \n";
      } 
    indic_Verif_Potentiel_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 Qeps = Invariant0Qeps(inv_n0);  
    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 = limite_inf_bIIb;
    double mini_val = ConstMath::pasmalpetit;
    Invariant inv_et_dbIIb = inv_n0;inv_et_dbIIb.bIIb += delta;
    double Qeps_et_dbIIb = Invariant0Qeps(inv_et_dbIIb);
    double E_et_dbIIb = PoGrenoble(Qeps_et_dbIIb,inv_et_dbIIb);
    Invariant inv_et_dbIIb2 = inv_n0;inv_et_dbIIb2.bIIb += 2. * delta;
    double Qeps_et_dbIIb2 = Invariant0Qeps(inv_et_dbIIb2);
    double E_et_dbIIb2 = PoGrenoble(Qeps_et_dbIIb2,inv_et_dbIIb2);
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    double Qeps_et_dV = Invariant0Qeps(inv_et_dV);
    double E_et_dV = PoGrenoble(Qeps_et_dV,inv_et_dV);
    Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
    double Qeps_et_dV2 = Invariant0Qeps(inv_et_dV2);
    double E_et_dV2 = PoGrenoble(Qeps_et_dV2,inv_et_dV2);
    Invariant inv_et_dVdbIIb = inv_n0;
    inv_et_dVdbIIb.V += delta;
    inv_et_dVdbIIb.bIIb += delta;
    double Qeps_et_dVdbIIb = Invariant0Qeps(inv_et_dVdbIIb);
    double E_et_dVdbIIb = PoGrenoble(Qeps_et_dVdbIIb,inv_et_dVdbIIb);
    // calcul des dérivées premières
    double E_V = jacobien0 *(E_et_dV * inv_et_dV.V - E_n * inv.V)/delta;
    double E_bIIb = (jacobien0 * inv.V)*(E_et_dbIIb - E_n)/delta;
    // calcul des dérivées secondes
    double E_V2 = jacobien0 *(E_et_dV2 * inv_et_dV2.V - 2.*E_et_dV * inv_et_dV.V + E_n * inv.V)
                   /(delta * delta);
    double E_bIIb2 = (jacobien0 * inv.V)*(E_et_dbIIb2 - 2.*E_et_dbIIb + E_n)
                     /(delta * delta);
    double E_V_a_dbIIb = jacobien0 * (E_et_dVdbIIb * inv_et_dVdbIIb.V - E_et_dbIIb * inv.V)
                         /delta;
    double E_VbIIb = ( E_V_a_dbIIb - E_V)/delta;
    // comparaison avec les valeurs de dérivées analytiques
    int erreur = 0;
    
    if (diffpourcent(potret.EV,E_V,E_V,0.05)) 
     { if (MiN(Dabs(potret.EV),Dabs(E_V)) <= mini_val)
         {if( MaX(Dabs(potret.EV),Dabs(E_V)) > 50.*delta)
           if (Dabs(potret.EV) == 0.)
           // on regarde si la dérivée numérique tend vers zéro lorque l'incrément tend vers 0
            { 
             Invariant inv_et_ddV = inv_n0;inv_et_ddV.V += delta/10.;
             double Qeps_et_ddV = Invariant0Qeps(inv_et_ddV);
             double E_et_ddV = PoGrenoble(Qeps_et_ddV,inv_et_ddV);
             double E_ddV = jacobien0 *(E_et_ddV * inv_et_ddV.V - E_n * inv.V)/(delta*10.);
             if (10.* Dabs(E_ddV) > Dabs(E_V)) erreur = 5;
             }
            else {erreur = 5;cout << " \n err5: potret.EV= " << potret.EV << " E_V= " << E_V;}
          else {erreur = 52;cout << " \n err52: potret.EV= " << potret.EV << " E_V= " << E_V;};
         } 
      else {erreur = 53;cout << " \n err53: potret.EV= " << potret.EV << " E_V= " << E_V;};
     };
      
// cout << " \n info: potret.EV= " << potret.EV << " E_V= " << E_V;      
      
      
    if (diffpourcent(potret.EbIIb,E_bIIb,E_bIIb,0.05)) 
     {if (MiN(Dabs(potret.EbIIb),Dabs(E_bIIb)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIb),Dabs(E_bIIb)) > 50.*delta) 
         	{erreur = 1;cout << " \n err 1: potret.EbIIb= " << potret.EbIIb << " E_bIIb= " << E_bIIb;}
         }
      else {erreur = 11;cout << " \n err11: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
     };
    if (diffpourcent(potret.EbIIb2,E_bIIb2,E_bIIb2,0.05)) 
     {if (MiN(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) > 50.*delta) 
         	{erreur = 2;cout << " \n err 2: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;}
         }
      else {erreur = 21;cout << " \n err21: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;};
     };
    if (diffpourcent(potret.EVV,E_V2,E_V2,0.05)) 
     {if (MiN(Dabs(potret.EVV),Dabs(E_V2)) <= mini_val)
         {if( MaX(Dabs(potret.EVV),Dabs(E_V2)) > 50.*delta) 
         	{erreur = 3;cout << " \n err 3: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;}
         }
      else {erreur = 31;cout << " \n err31: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;};
     };
   if (diffpourcent(potret.EbIIbV,E_VbIIb,E_VbIIb,0.05)) 
     {if (MiN(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) > 50.*delta) 
         	{erreur = 4;cout << " \n err 4: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;}
         }
      else {erreur = 41;cout << " \n err41: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
     };

//    if (erreur)
//     { cout << "\n erreur dans le calcul analytique des dérivées du potentiel";
//       cout << "\n Hyper3D::Verif_Potentiel_et_var, numéro d'appel = " 
//            << indic_Verif_Potentiel_et_var;
//       Sortie(1);
//      } 
 };

// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
// cas avec la phase
void Hyper3D::Verif_Potentiel_et_var(const InvariantVarDdl& inv,const double& jacobien0
                                 ,const PotenAvecPhaseAvecVar& potret )
  { // dans le cas du premier passage on indique qu'il y a vérification
    if (indic_Verif_Potentiel_et_var == 0)
     { cout << "\n ****verification des derivees du potentiels par rapport aux invariants****";
       cout << "\n Hyper3D::Verif_Potentiel_et_var \n";
      } 
    indic_Verif_Potentiel_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 is_sp0 = Invariant0Specif(inv_n0);
    double E_n = PoGrenoble(is_sp0,inv_n0);  // la valeur du potentiel de référence
    // ici on est avec phase donc trois invariants indépendant V et Qeps et cos3phi,
    // 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 = limite_inf_bIIb;
    double mini_val = ConstMath::pasmalpetit;
    Invariant inv_et_dbIIb = inv_n0;inv_et_dbIIb.bIIb += delta;
    Invariant0QepsCosphi is_et_dbIIb = Invariant0Specif(inv_et_dbIIb);
    double E_et_dbIIb = PoGrenoble(is_et_dbIIb,inv_et_dbIIb);
    Invariant inv_et_dbIIb2 = inv_n0;inv_et_dbIIb2.bIIb += 2. * delta;
    Invariant0QepsCosphi is_et_dbIIb2 = Invariant0Specif(inv_et_dbIIb2);
    double E_et_dbIIb2 = PoGrenoble(is_et_dbIIb2,inv_et_dbIIb2);
    Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
    Invariant0QepsCosphi is_et_dV = Invariant0Specif(inv_et_dV);
    double E_et_dV = PoGrenoble(is_et_dV,inv_et_dV);
    Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
    Invariant0QepsCosphi is_et_dV2 = Invariant0Specif(inv_et_dV2);
    double E_et_dV2 = PoGrenoble(is_et_dV2,inv_et_dV2);
    Invariant inv_et_dVdbIIb = inv_n0;
    inv_et_dVdbIIb.V += delta;
    inv_et_dVdbIIb.bIIb += delta;
    Invariant0QepsCosphi is_et_dVdbIIb = Invariant0Specif(inv_et_dVdbIIb);
    double E_et_dVdbIIb = PoGrenoble(is_et_dVdbIIb,inv_et_dVdbIIb);
    // cas avec Ieps
    
    Invariant inv_et_dIeps = inv_n0;inv_et_dIeps.Ieps += delta;
    Invariant0QepsCosphi is_et_dIeps = Invariant0Specif(inv_et_dIeps);
    double E_et_dIeps = PoGrenoble(is_et_dIeps,inv_et_dIeps);
    
    Invariant inv_et_dIeps2 = inv_n0;inv_et_dIeps2.Ieps += 2. * delta;
    Invariant0QepsCosphi is_et_dIeps2 = Invariant0Specif(inv_et_dIeps2);
    double E_et_dIeps2 = PoGrenoble(is_et_dIeps2,inv_et_dIeps2);
    
    Invariant inv_et_dIepsdbIIb = inv_n0;
    inv_et_dIepsdbIIb.Ieps += delta;
    inv_et_dIepsdbIIb.bIIb += delta;
    Invariant0QepsCosphi is_et_dIepsdbIIb = Invariant0Specif(inv_et_dIepsdbIIb);
    double E_et_dIepsdbIIb = PoGrenoble(is_et_dIepsdbIIb,inv_et_dIepsdbIIb);
    
    Invariant inv_et_dIepsdV = inv_n0;
    inv_et_dIepsdV.Ieps += delta;
    inv_et_dIepsdV.V += delta;
    Invariant0QepsCosphi is_et_dIepsdV = Invariant0Specif(inv_et_dIepsdV);
    double E_et_dIepsdV = PoGrenoble(is_et_dIepsdV,inv_et_dIepsdV);
    
    
    // calcul des dérivées premières

    // calcul des dérivées premières
    double E_V = jacobien0 *(E_et_dV * inv_et_dV.V - E_n * inv.V)/delta;
    double E_bIIb = (jacobien0 * inv.V)*(E_et_dbIIb - E_n)/delta;
    double E_Ieps = (jacobien0 * inv.V)*(E_et_dIeps  - E_n )/delta;
    // calcul des dérivées secondes hors phases
    double E_V2 = jacobien0 *(E_et_dV2 * inv_et_dV2.V - 2.*E_et_dV * inv_et_dV.V + E_n * inv.V)
                   /(delta * delta);
    double E_bIIb2 = (jacobien0 * inv.V)*(E_et_dbIIb2 - 2.*E_et_dbIIb + E_n) /(delta * delta);

    double E_Ieps2 = (jacobien0 * inv.V) *(E_et_dIeps2  - 2.*E_et_dIeps  + E_n )/(delta * delta);

    double E_V_a_dbIIb = jacobien0 * (E_et_dVdbIIb * inv_et_dVdbIIb.V - E_et_dbIIb * inv_et_dbIIb.V) /delta;
    double E_VbIIb = ( E_V_a_dbIIb - potret.EV)/delta;
    
    double E_Ieps_a_dbIIb = (jacobien0 * inv.V) * (E_et_dIepsdbIIb  - E_et_dbIIb) /delta;
    double E_bIIbIeps = ( E_Ieps_a_dbIIb - potret.EIeps)/delta;
    
    double E_Ieps_a_dV = jacobien0 * (E_et_dIepsdV * inv_et_dIepsdV.V - E_et_dV * inv.V) /delta;
    double E_IepsV = ( E_Ieps_a_dV - potret.EIeps)/delta;
    
    // comparaison avec les valeurs de dérivées analytiques
    int erreur = 0;
    
    if (diffpourcent(potret.EV,E_V,E_V,0.05)) 
     {if (MiN(Dabs(potret.EV),Dabs(E_V)) <= mini_val)
         {if( MaX(Dabs(potret.EV),Dabs(E_V)) > 50.*delta)
           if (Dabs(potret.EV) == 0.)
           // on regarde si la dérivée numérique tend vers zéro lorque l'incrément tend vers 0
            { 
             Invariant inv_et_ddV = inv_n0;inv_et_ddV.V += delta/10.;
             double Qeps_et_ddV = Invariant0Qeps(inv_et_ddV);
             double E_et_ddV = PoGrenoble(Qeps_et_ddV,inv_et_ddV);
             double E_ddV = jacobien0 *(E_et_ddV * inv_et_ddV.V - E_n * inv.V)/(delta*10.);
             if (10.* Dabs(E_ddV) > Dabs(E_V)) erreur = 5;
             }
            else {erreur = 5;cout << " \n err5: potret.EV= " << potret.EV << " E_V= " << E_V;}
          else {erreur = 52;cout << " \n err52: potret.EV= " << potret.EV << " E_V= " << E_V;};
         } 
      else 
       {erreur = 53;cout << " \n err53: potret.EV= " << potret.EV << " E_V= " << E_V;};
     };
    if (diffpourcent(potret.EbIIb,E_bIIb,E_bIIb,0.05))
     {if (MiN(Dabs(potret.EbIIb),Dabs(E_bIIb)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIb),Dabs(E_bIIb)) > 50.*delta) 
         	{erreur = 1;cout << " \n err 1: potret.EbIIb= " << potret.EbIIb << " E_bIIb= " << E_bIIb;}
         }
      else {erreur = 11;cout << " \n err 11: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
     };
    if (diffpourcent(potret.EbIIb2,E_bIIb2,E_bIIb2,0.05))
     {if (MiN(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) > 50.*delta) 
         	{erreur = 2;cout << " \n err 2: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;}
         }
      else {erreur = 21;cout << " \n err 21: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;};
     };
    if (diffpourcent(potret.EVV,E_V2,E_V2,0.05))
     {if (MiN(Dabs(potret.EVV),Dabs(E_V2)) <= mini_val)
         {if( MaX(Dabs(potret.EVV),Dabs(E_V2)) > 50.*delta) 
         	{erreur = 3;cout << " \n err 3: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;}
         }
      else {erreur = 31;cout << " \n err 31: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;};
     };
    if (diffpourcent(potret.EbIIbV,E_VbIIb,E_VbIIb,0.05))
     {if (MiN(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) > 50.*delta) 
         	{erreur = 4;cout << " \n err 4: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;}
         }
      else {erreur = 41;cout << " \n err 41: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
     };

    if (diffpourcent(potret.EIeps,E_Ieps,E_Ieps,0.05)) 
     {if (MiN(Dabs(potret.EIeps),Dabs(E_Ieps)) <= mini_val)
         {if( MaX(Dabs(potret.EIeps),Dabs(E_Ieps)) > 50.*delta) 
         	{erreur = 6;cout << " \n err 6: potret.EIeps= " << potret.EIeps << " E_Ieps= " << E_Ieps;}
         }
      else {erreur = 61;cout << " \n err 61: potret.EIeps= " << potret.EIeps << " E_Ieps= " << E_Ieps;};
     };

    if (diffpourcent(potret.EIeps2,E_Ieps2,E_Ieps2,0.05)) 
     {if (MiN(Dabs(potret.EIeps2),Dabs(E_Ieps2)) <= mini_val)
         {if( MaX(Dabs(potret.EIeps2),Dabs(E_Ieps2)) > 50.*delta) 
         	{erreur = 7;cout << " \n err 7: potret.EIeps2= " << potret.EIeps2 << " E_Ieps2= " << E_Ieps2;}
         }
      else {erreur = 71;cout << " \n err 71: potret.EIeps2= " << potret.EIeps2 << " E_Ieps2= " << E_Ieps2;};
     };
   
    if (diffpourcent(potret.EbIIbIeps,E_bIIbIeps,E_bIIbIeps,0.05)) 
     {if (MiN(Dabs(potret.EbIIbIeps),Dabs(E_bIIbIeps)) <= mini_val)
         {if( MaX(Dabs(potret.EbIIbIeps),Dabs(E_bIIbIeps)) > 50.*delta) 
         	{erreur = 8;cout << " \n err 8: potret.EbIIbIeps= " << potret.EbIIbIeps << " E_bIIbIeps= " << E_bIIbIeps;}
         }
      else {erreur = 81;cout << " \n err 81: potret.EbIIbIeps= " << potret.EbIIbIeps << " E_bIIbIeps= " << E_bIIbIeps;};
     };

    if (diffpourcent(potret.EIepsV,E_IepsV,E_IepsV,0.05)) 
     {if (MiN(Dabs(potret.EIepsV),Dabs(E_IepsV)) <= mini_val)
         {if( MaX(Dabs(potret.EIepsV),Dabs(E_IepsV)) > 50.*delta) 
         	{erreur = 9;cout << " \n err 9: potret.EIepsV= " << potret.EIepsV << " E_IepsV= " << E_IepsV;}
         }
      else {erreur = 91;cout << " \n err 91: potret.EIepsV= " << potret.EIepsV << " E_IepsV= " << E_IepsV;};
     };
/*    if (erreur)
     { cout << "\n erreur dans le calcul analytique des derivees du potentiel";
       cout << "\n Hyper3D::Verif_Potentiel_et_var, numero d'appel = " 
            << indic_Verif_Potentiel_et_var << " erreur = " << erreur << endl;
       Sortie(1);
      }*/ 
 };