// FICHIER : Hyper_W_gene_3D.cc
// CLASSE : Hyper_W_gene_3D


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

//#include "Debug.h"

# include <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "ParaGlob.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"
#include "CharUtil.h"
  
// CONSTRUCTEURS :
// ---------- classe de stockage des grandeurs spécifiques pour la loi ---------

// constructeur par défaut
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::SaveResulHyper_W_gene_3D() :
           invP(NULL),invP_t(NULL),map_type_quelconque()
           ,para_loi(NULL),para_loi_t(NULL)
{ // mise à jour de la liste des grandeurs quelconques internes
  Mise_a_jour_map_type_quelconque();
};
// avec_para : le nombre indique si on sauvegarde des paramètres matériaux qui varient
//             et également des invariants.
// si == 0 : pas de stockage (joue en fait le role de sortie_post)
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::SaveResulHyper_W_gene_3D(int avec_para) :
   invP(NULL),invP_t(NULL)
   ,para_loi(NULL),para_loi_t(NULL)
   ,map_type_quelconque()
{ if (avec_para)
    {invP = new Invariantpost3D();
     invP_t = new Invariantpost3D();
    };
  if (avec_para)
    {para_loi = new Vecteur(avec_para,-ConstMath::trespetit);
     para_loi_t = new Vecteur(avec_para,-ConstMath::trespetit);
    };
  // mise à jour de la liste des grandeurs quelconques internes
  Mise_a_jour_map_type_quelconque();
};

// constructeur de copie	           
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::SaveResulHyper_W_gene_3D(const SaveResulHyper_W_gene_3D& sav):
     invP(NULL),invP_t(NULL),para_loi(NULL),para_loi_t(NULL)
     ,map_type_quelconque(sav.map_type_quelconque)
 { if (sav.invP != NULL)
     {invP = new Invariantpost3D(*(sav.invP));
      invP_t = new Invariantpost3D(*(sav.invP_t));
     };
   if (sav.para_loi != NULL)
     {para_loi = new Vecteur(*sav.para_loi);
      para_loi_t = new Vecteur(*sav.para_loi_t);
     };
   // mise à jour de la liste des grandeurs quelconques internes
   Mise_a_jour_map_type_quelconque();

 };

//  destructeur
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::~SaveResulHyper_W_gene_3D()
  { if (invP != NULL)
     { delete invP;
       delete invP_t;
     };
    if (para_loi != NULL) delete para_loi;
    if (para_loi_t != NULL) delete para_loi_t;
  };
  
// affectation
Loi_comp_abstraite::SaveResul & Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::operator = ( const Loi_comp_abstraite::SaveResul & a)
  { Hyper_W_gene_3D::SaveResulHyper_W_gene_3D& sav = *((Hyper_W_gene_3D::SaveResulHyper_W_gene_3D*) &a);
    if (sav.invP != NULL)
     { if (invP == NULL)
         {invP = new Invariantpost3D(*(sav.invP));
          invP_t = new Invariantpost3D(*(sav.invP_t));
         }
       else
        { *invP = (*(sav.invP)); *invP_t = (*(sav.invP_t));
        };
     }
    else
     // sinon cela veut dire que invP et invP_t ne doivent pas être affecté
     { if (invP != NULL)
         {delete invP;invP=NULL;
          delete invP_t;invP_t=NULL;
         };
     };
    map_type_quelconque = sav.map_type_quelconque;
    if (sav.para_loi != NULL)
     {if (para_loi == NULL)
       {para_loi = new Vecteur(sav.para_loi->Taille(),-ConstMath::trespetit);
        para_loi_t = new Vecteur(sav.para_loi->Taille(),-ConstMath::trespetit);
       };
      *para_loi = *(sav.para_loi);
      *para_loi_t = *(sav.para_loi_t);
     };

    return *this;
  };

//------- lecture écriture dans base info -------
   // cas donne le niveau de la récupération
   // = 1 : on récupère tout
   // = 2 : on récupère uniquement les données variables (supposées comme telles)
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Lecture_base_info (ifstream& ent,const int )
  { string toto; int presence=0;
    ent >> toto >> presence;
    #ifdef MISE_AU_POINT
      if (toto != "dat_sp_hyperelas:")
         { cout << "\n erreur dans la lecture des donnees specifiques hyperelastiques, on ne trouve pas le mot cle: dat_sp_hyperelas: "
                << "\n  Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Lecture_base_info (... ";
           Sortie(1);
          }       
    #endif
    if (presence == 1)
     {if (invP == NULL)
        {invP = new Invariantpost3D();
         invP_t = new Invariantpost3D();
        };
      // lecture
      ent >> (*invP); (*invP_t)= (*invP); //toto >>  invP->V >> toto >>  invP->Qeps >> toto >>  invP->cos3phi;
     }
    else
     { if (invP != NULL)
         { delete invP;invP=NULL;
           delete invP_t;invP_t=NULL;
         };
     };
    // cas des paramètres de lois éventuelles
    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 (toto == "avec_para_taille=")
     {int taille;
      ent >> taille;
      if (para_loi == NULL)
       {para_loi = new Vecteur(taille);
        para_loi_t = new Vecteur(taille);
       };
      ent >> *para_loi;
      *para_loi_t=*para_loi;
     };
    // mise à jour de la liste des grandeurs quelconques internes
    Mise_a_jour_map_type_quelconque();

   };
   
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables 
//(supposées comme telles)
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Ecriture_base_info(ofstream& sort,const int )
  { sort << "\n dat_sp_hyperelas: ";
    if (invP != NULL)
      { sort << " 1 " << (*invP); }
    else { sort << " 0  ";};
    // cas des paramètres de lois éventuelles
    if (para_loi == NULL)
     {sort << " sans_para ";}
    else
     {sort << " avec_para_taille= "<<para_loi->Taille() << " ";
      sort << *para_loi;
     };

  };
 
      // mise à jour des informations transitoires
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::TdtversT()
   {if (invP != NULL)
      (*invP_t) = (*invP);
    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 Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::TversTdt()
   {if (invP != NULL)
      (*invP) = (*invP_t);
    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 Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Affiche() const
   { cout << "\n data specifiques hyperelastiques : " ;
     if (invP != NULL)
       { cout << (*invP); };
     if (para_loi == NULL)
      {cout << " sans_para ";}
     else
      {cout << " avec_para ";
       cout << *para_loi;
      };
   };


// mise à jour de la liste des grandeurs quelconques internes
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Mise_a_jour_map_type_quelconque()
  { // on parcours la map
    map <  EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >::iterator
      il,ilfin=map_type_quelconque.end();
    for (il=map_type_quelconque.begin();il != ilfin;il++)
    
     { switch ((*il).first)
       {  case POTENTIEL:
           { Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) (*il).second.Grandeur_pointee()); // pour simplifier
             (*tyTQ.ConteneurDouble()) = invP->potentiel;
             break;
           }
          // -----cas des paramètres de la loi
          case PARA_LOI_TYPE_MOONEY:
           { Grandeur_Vecteur& gr= *((Grandeur_Vecteur*) map_type_quelconque[PARA_LOI_TYPE_MOONEY].Grandeur_pointee()); // pour simplifier
            if (para_loi != NULL)
               {*(gr.ConteneurVecteur()) =  *(para_loi);}
            else {gr.ConteneurVecteur()->Zero();};
              break;
            }

          default:
            cout << "\n *** erreur on demande l'acces  a : "
                 << NomTypeQuelconque((*il).first)
                 << " or celui-ci n'est pas dispo pour la loi ";
            this->Affiche();
            cout << " revoir la mise en donnees ! " << endl;
            Sortie(1);
       };

     };

};

// ---------- classe Hyper_W_gene_3D ---------

		// Constructeur par defaut
Hyper_W_gene_3D::Hyper_W_gene_3D ():
     Loi_comp_abstraite()
     ,I_B(0.),I_BB(0.),II_B(0.),III_B(0.)
     ,sortie_post(0),nb_para_loi(0)
    // I_B(0.),I_BB(0.),I_BBB(0.),II_B(0.),III_B(0.)
     ,J_r(3)
     ,d_I_B_epsBB_HH(),d_II_B_epsBB_HH(),d_III_B_epsBB_HH()
     ,d_J_r_epsBB_HH(3)
    // ,d_I_B_eps2BB_HHHH(),d_II_B_eps2BB_HHHH(),d_III_B_eps2BB_HHHH()
     ,d_J_1_eps2BB_HHHH(),d_J_2_eps2BB_HHHH(),d_J_3_eps2BB_HHHH()
     ,V(0.),BB_HH(),BB_HB()
     ,IxI_HHHH(),IxbarreI_HHHH(),IxB_HHHH(),BxI_HHHH()
	{
		
	};
		
 
		// Constructeur utile si l'identificateur du nom de la loi
		// de comportement et la dimension sont connus
		// vit_def indique si oui ou non la loi utilise la vitesse de déformation
Hyper_W_gene_3D::Hyper_W_gene_3D (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp
		                    ,int dimension,bool vit_def):
     Loi_comp_abstraite(id_compor,categorie_comp,dimension,vit_def)
     ,I_B(0.),I_BB(0.),II_B(0.),III_B(0.)
     ,sortie_post(0),nb_para_loi(0)
    // I_B(0.),I_BB(0.),I_BBB(0.),II_B(0.),III_B(0.)
     ,J_r(3)
     ,d_I_B_epsBB_HH(),d_II_B_epsBB_HH(),d_III_B_epsBB_HH()
     ,d_J_r_epsBB_HH(3)
    // ,d_I_B_eps2BB_HHHH(),d_II_B_eps2BB_HHHH(),d_III_B_eps2BB_HHHH()
     ,d_J_1_eps2BB_HHHH(),d_J_2_eps2BB_HHHH(),d_J_3_eps2BB_HHHH()
     ,V(0.),BB_HH(),BB_HB()
     ,IxI_HHHH(),IxbarreI_HHHH(),IxB_HHHH(),BxI_HHHH()
	{
		
	};

		// Constructeur de copie
Hyper_W_gene_3D::Hyper_W_gene_3D (const Hyper_W_gene_3D& loi):
     Loi_comp_abstraite(loi)
     ,nb_para_loi(loi.nb_para_loi)
     ,I_B(loi.I_B),I_BB(loi.I_BB),II_B(II_B),III_B(III_B)
    // I_B(loi.I_B),I_BB(loi.I_BB),I_BBB(loi.I_BBB),II_B(II_B),III_B(III_B)
     ,J_r(loi.J_r)
     ,d_I_B_epsBB_HH(loi.d_I_B_epsBB_HH),d_II_B_epsBB_HH(loi.d_II_B_epsBB_HH)
    // ,d_III_B_epsBB_HH(loi.d_III_B_epsBB_HH)
     ,d_J_r_epsBB_HH(loi.d_J_r_epsBB_HH)
    // ,d_I_B_eps2BB_HHHH(loi.d_I_B_eps2BB_HHHH)
     ,d_J_1_eps2BB_HHHH(loi.d_J_1_eps2BB_HHHH)
     ,d_J_2_eps2BB_HHHH(loi.d_J_2_eps2BB_HHHH)
     ,d_J_3_eps2BB_HHHH(loi.d_J_3_eps2BB_HHHH)
    // ,d_II_B_eps2BB_HHHH(loi.d_II_B_eps2BB_HHHH)
    // ,d_III_B_eps2BB_HHHH(loi.d_III_B_eps2BB_HHHH)
     ,V(loi.V),BB_HH(loi.BB_HH),BB_HB(loi.BB_HB)
     ,IxI_HHHH(),IxbarreI_HHHH(),IxB_HHHH(),BxI_HHHH()
     ,sortie_post(loi.sortie_post)
	{
		
	};
		
		// DESTRUCTEUR :
		
Hyper_W_gene_3D::~Hyper_W_gene_3D ()
	{
		
	};

// 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 Hyper_W_gene_3D::Grandeur_particuliere
              (bool ,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list<int>& decal) const
{  // ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas
   // 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 POTENTIEL:
           { SaveResulHyper_W_gene_3D & save_resul = *((SaveResulHyper_W_gene_3D*) saveDon);
             Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
             if (sortie_post)
                   tyTQ(1+(*idecal)) = (save_resul.invP)->potentiel;
             else
                   tyTQ(1+(*idecal)) = 0.;
             (*idecal)++; break;
           }
          case PARA_LOI_TYPE_MOONEY:
           //  ----- cas des coefficients de loi
            { SaveResulHyper_W_gene_3D & save_resul = *((SaveResulHyper_W_gene_3D*) 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 Hyper_W_gene_3D::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
 {// ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas 
  Tableau <double> tab_1(1);
  Tab_Grandeur_scalaire_double grand_courant(tab_1);
  // on ne propose des grandeurs que si elles ont été stockée
  // non car le choix est un choix global, par contre les valeurs ne seront différentes de 0 que si sortie_post
  // est activé !!
//  if (sortie_post)
   { // def d'un type quelconque représentatif à chaque grandeur
     // a priori ces grandeurs sont défini aux points d'intégration identique à la contrainte par exemple
     // enu_ddl_type_pt est définit dans la loi Abtraite générale
     // on n'ajoute que si sortie_post est vraie, sinon aucune grandeur n'est sauvegardé, donc on ne peut
     // plus y accèder
     //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
     // cas du potentiel
     {List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
      for (itq=liTQ.begin();itq!=itqfin;itq++)
        if ((*itq).EnuTypeQuelconque() == POTENTIEL)
         { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
           int taille = tyTQ.Taille()+1;
           tyTQ.Change_taille(taille); nexistePas = false;
         };
      if (nexistePas)
        {TypeQuelconque typQ1(POTENTIEL,enu_ddl_type_pt,grand_courant);
         liTQ.push_back(typQ1);
        };
     };
   };
  // $$$ 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_TYPE_MOONEY)
      {Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
       int taille = tyTQ.Taille()+1;
       tyTQ.Change_taille(taille); nexistePas = false;
       tyTQ(taille).Change_taille(nb_para_loi);
      };
   if (nexistePas)
     {Vecteur  tens(nb_para_loi); //
      Tab_Grandeur_Vecteur gr_vec_para(tens,1);
      // def d'un type quelconque représentatif
      TypeQuelconque typQ(PARA_LOI_TYPE_MOONEY,EPS11,gr_vec_para);
      liTQ.push_back(typQ);
    };
  };

 };
 
// activation du stockage de grandeurs quelconques qui pourront ensuite être récupéré
// via le conteneur SaveResul, si la grandeur n'existe pas ici, aucune action
void Hyper_W_gene_3D::Activation_stockage_grandeurs_quelconques(list <EnumTypeQuelconque >& listEnuQuelc)
 {  // récup de la liste de stockage
    list <EnumTypeQuelconque >& listlocale = ListQuelc_mis_en_acces_localement();
  
    // on parcours la liste des grandeurs à activer
    // et on remplit la liste locale
    list <EnumTypeQuelconque >::iterator il, ilfin = listEnuQuelc.end();
    for (il = listEnuQuelc.begin();il != ilfin; il++)
//    for (EnumTypeQuelconque enu : listEnuQuelc)
     // on ne remplit que s'il s'agit d'une grandeur qui peut-être accessible
     {switch (*il)
       {case POTENTIEL: case PARA_LOI_TYPE_MOONEY:
          listlocale.push_back(*(il));
          break;
        default: ;  // pour les autres cas on ne fait rien
       };
     };
  
    // on supprime les doublons localement
    listlocale.sort(); // on ordonne la liste
    listlocale.unique(); // suppression des doublons
  
 };
 
// insertion des conteneurs ad hoc concernant le stockage de grandeurs quelconques
// passée en paramètre, dans le save result: ces conteneurs doivent être valides
// c-a-d faire partie de listdeTouslesQuelc_dispo_localement
void Hyper_W_gene_3D::Insertion_conteneur_dans_save_result(SaveResul * a)
 {
    // -- autre stockage éventuel en fonction des grandeurs quelconques demandées par d'autres lois
    // on va regarder s'il faut activer ou non invP et invP_t, pour le cas
    // de la récupération de grandeur quelconque au moment de l'utilisation de la loi
    // récup de la liste de stockage
    list <EnumTypeQuelconque >& listlocale = ListQuelc_mis_en_acces_localement();
    // def
    Hyper_W_gene_3D::SaveResulHyper_W_gene_3D& sav = *((Hyper_W_gene_3D::SaveResulHyper_W_gene_3D*) a);
  
    List_io <EnumTypeQuelconque>::iterator jk,jkfin = listlocale.end();
    for (jk=listlocale.begin();jk != jkfin;jk++)
//    for (EnumTypeQuelconque enu : listQuelc_dispo_localement)
//     { switch (enu)
      {EnumTypeQuelconque enu = *jk;
       switch (enu)
        { case POTENTIEL: // cas on des grandeurs sont demandées
            { if (sav.invP == NULL) // def si ils n'existent pas
               {sav.invP = new Invariantpost3D();
                sav.invP_t = new Invariantpost3D();
                if ((Permet_affichage() > 0 ) && (sortie_post == 0))
                 {cout << "\n *** warning: on demande pour une loi hyper elastique, l'utilisation de la grandeur "
                       << NomTypeQuelconque(enu) << " or on n'a pas demande sortie_post_ different de 0 "
                       << " on met arbitrairement sortie_post_ = 1 ";
                 };
               if (sortie_post == 0)
                 sortie_post = 1;
               };
              // on crée le conteneur ad hoc pour le passage d'info
              // def d'un conteneur de grandeurs quelconques, initialisée à 0
              Grandeur_scalaire_double grand_courant(0.);
              TypeQuelconque typQ1(enu,EPS11,grand_courant);
              sav.map_type_quelconque[enu]=(typQ1);
              break;
            }
          case PARA_LOI_TYPE_MOONEY:  // cas on des grandeurs sont demandées
            { if (sav.para_loi == NULL) // def si ils n'existent pas
               {sav.para_loi = new Vecteur; // la dimension sera ensuite précisée
                sav.para_loi_t = new Vecteur; // dans la loi dérivée
               };
              // on crée le conteneur ad hoc pour le passage d'info
              // def d'un conteneur de grandeurs quelconques, initialisée à 0
              Grandeur_Vecteur grand_courant(3);grand_courant.ConteneurVecteur()->Zero();
              TypeQuelconque typQ1(enu,EPS11,grand_courant);
              sav.map_type_quelconque[enu]=(typQ1);
              break;
            }

          default:
            cout << "\n *** erreur on demande l'acces  a : "
                 << NomTypeQuelconque(enu)
                 << " or celui-ci n'est pas dispo pour la loi ";
            this->Affiche();
            cout << " revoir la mise en donnees ! " << endl;
            Sortie(1);
        };
     };
 
 };


// --------------- méthodes internes -------------
// affichage et definition interactive des commandes particulières à chaques lois
void Hyper_W_gene_3D::Info_commande_LoisDeComp_hyper3D(UtilLecture& entreePrinc)
 {  ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
    sort << "\n# il est egalement possible de stocker et donc de recuperer en post-traitement, differentes grandeurs "
         << "\n#  specifiquement utilises : ici la valeur du potentiel ,"
         << "\n#  pour cela sur la derniere ligne des donnees de la loi il "
         << "\n#  faut la presence du mot cle: sortie_post_ suivi de 1, par defaut sa valeur est 0 "
         << "\n# ex: sortie_post_ 1 ";
 };


		// calcul des invariants et de leurs variations premières
void Hyper_W_gene_3D::InvariantsEtVar1(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
		                                ,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
		                                ,const double& jacobien_0,const double& jacobien)
	{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
   const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); //  "
   const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); //  "
   const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); //  "
   // tout d'abord on définit deux tenseurs qui nous seront utiles par la suite
   const Tenseur3HH & B_HH = gijHH_0; // le tenseur B
   B_HB = B_HH * gijBB_tdt; // le tenseur B  en mixte
   BB_HB = B_HB * B_HB;
   BB_HH = BB_HB * gijHH_tdt; // le tenseur B . B en deux fois contra
// BB_HH = B_HB * B_HH ;
   // calcul des trois invariants traces primaires de B
   I_B= B_HB.Trace(); //gijBB_tdt && B_HH ; //gijHH_0;
   I_BB= BB_HB.Trace(); //gijBB_tdt && BB_HH; //((gijBB_tdt * gijHH_0) * gijBB_tdt) && gijHH_0;
//	 I_BBB= gijBB_tdt && ((BB_HH * gijBB_tdt) * B_HH); // ((((gijBB_tdt * gijHH_0) * gijBB_tdt) * gijHH_0) * gijBB_tdt) && gijHH_0;
   // calcul des trois invariants en I de B: I_B, II_B, III_B
   II_B=0.5*(I_B*I_B-I_BB);
   const double untiers=1./3.;
//	 III_B=untiers*I_BBB - 0.5*I_B*I_BB + I_B*I_B*I_B/6;
   V= jacobien/jacobien_0;  // variation relative de volume
   III_B=V*V;               // calcul plus rapide que le précédent
   // calcul des trois invariants en J
   if (III_B <= ConstMath::petit)
     { if (Permet_affichage() >= 1)
         cout << "\n **** erreur, le volume final sur volume initial est quasi-nulle ! III_B=" << (III_B)
              << " on limite (V**(2)) a : " << ConstMath::petit;
       if (Permet_affichage() >= 5)
         cout << "\n Hyper_W_gene_3D::Invariants_et_var1( ";
       III_B = ConstMath::petit;     
     };
   // variables intermédiaires 	
   J3_puiss_untiers=pow(III_B,untiers);
   unSurJ3_puissuntiers=1./J3_puiss_untiers;
   unSurJ3_puissuntiers2=unSurJ3_puissuntiers*unSurJ3_puissuntiers;
   unSurJ3_puissuntiers4=unSurJ3_puissuntiers2*unSurJ3_puissuntiers2;
   // les J
   J_r(1)= I_B*unSurJ3_puissuntiers;
   J_r(2)=II_B*unSurJ3_puissuntiers2;
   J_r(3)=III_B;
//	 V = sqrt(III_B); // autre calcul de la variation relative de volume
	 // calcul des variations premières des I_r par rapports aux composantes epsBB d'Almansi
// ( en fait on ne les calcul pas pour optimiser le temps)	 
//	 d_I_B_epsBB_HH=2.* B_HH ; //gijHH_0;
//	 d_II_B_epsBB_HH= 2.* (I_B * B_HH - BB_HH);      
//	                // - 2. * ((gijHH_0 * gijBB_tdt) * gijHH_0);
//	 d_III_B_epsBB_HH=(2.*III_B)*gijHH_tdt;
	 // calcul des variation premières des J_r par rapports aux composantes epsBB d'Almansi

   d_J_r_epsBB_HH(1)=(2.*unSurJ3_puissuntiers)*(B_HH - (untiers*I_B)*gijHH_tdt);
   d_J_r_epsBB_HH(2)=(2.*unSurJ3_puissuntiers2)*(I_B * B_HH - BB_HH - (2.*untiers*II_B)*gijHH_tdt);
   d_J_r_epsBB_HH(3)=(2.*III_B)*gijHH_tdt;
   #ifdef MISE_AU_POINT
   if  (Permet_affichage() > 4)
     { cout << "\nHyper_W_gene_3D::InvariantsEtVar1(..."
            << "\n I_B= " << I_B << ", II_B= "<< II_B <<", III_B= "<< III_B
            << "\n J_r(1)= " << J_r(1) <<", J_r(2)= " << J_r(2) << ", J_r(3)= " << J_r(3)
            << "\n d_J_r_epsBB_HH(1)= " << d_J_r_epsBB_HH(1)
            << ", d_J_r_epsBB_HH(2)= " << d_J_r_epsBB_HH(2)
            << ", d_J_r_epsBB_HH(3)= " << d_J_r_epsBB_HH(3);
     };
   #endif
	};

		// calcul des invariants et de leurs variations premières et secondes
void Hyper_W_gene_3D::InvariantsEtVar2(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
		                       ,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
		                       ,const  double& jacobien_0,const double& jacobien)
	{ // calcul des invariants et de leurs variations premières
	  Invariants_et_var1(gijBB_0_,gijHH_0_,gijBB_tdt_,gijHH_tdt_,jacobien_0,jacobien);

	  const double untiers=1./3.;
	  // calcul des variations secondes
	  const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
   const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); //  "
   const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); //  "
   const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); //  "
   const Tenseur3HH & B_HH = gijHH_0; // récup du tenseur B
   // le tenseur BB est déjà calculé
   // calcul de variables intermédiaires
   IxI_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,gijHH_tdt);
   IxbarreI_HHHH=Tenseur3HHHH::Prod_tensoriel_barre(gijHH_tdt,gijHH_tdt);
   IxB_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,B_HH);
   BxI_HHHH=Tenseur3HHHH::Prod_tensoriel(B_HH,gijHH_tdt);
   
   // variations secondes des I_r 
   //(finalement on ne les calcule pas pour optimiser)
   ////  d_I_B_eps2BB_HHHH.Inita(0.);
   ////  d_II_B_eps2BB_HHHH = 4.*(Tenseur3HHHH::Prod_tensoriel(B_HH,B_HH)
   ////                        - Tenseur3HHHH::Prod_tensoriel_barre(B_HH,B_HH));
   ////  d_III_B_eps2BB_HHHH = (4. * III_B) * 
   ////           ( Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,gijHH_tdt)
   ////            - Tenseur3HHHH::Prod_tensoriel_barre(gijHH_tdt,gijHH_tdt));

   // variations secondes des J_r
//       d_J_1_eps2BB_HHHH = (-2.*untiers*unSurJ3_puissuntiers4)*
//                         Tenseur3HHHH::Prod_tensoriel(d_I_B_epsBB_HH,d_III_B_epsBB_HH)
//                             - (untiers * I_B * unSurJ3_puissuntiers4) * d_III_B_eps2BB_HHHH;
                             
   d_J_1_eps2BB_HHHH = (-4.*untiers*unSurJ3_puissuntiers)*
                     ( BxI_HHHH + IxB_HHHH - I_B*(untiers*IxI_HHHH + IxbarreI_HHHH) );

   d_J_2_eps2BB_HHHH = (-8. * untiers * unSurJ3_puissuntiers2 )*
                          (I_B* BxI_HHHH
                           - Tenseur3HHHH::Prod_tensoriel(BB_HH,gijHH_tdt)
                           - (2.*untiers)*II_B * IxI_HHHH
                           )
                       +(4. * unSurJ3_puissuntiers2 )*    
                         (Tenseur3HHHH::Prod_tensoriel(B_HH,B_HH)
                           - Tenseur3HHHH::Prod_tensoriel_barre(B_HH,B_HH)
                         -(2.*untiers )*(I_B * IxB_HHHH
                           - Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,BB_HH))
                         +(2.*untiers * II_B)*IxbarreI_HHHH
                         );
                           
   d_J_3_eps2BB_HHHH = (4. * III_B) *  ( IxI_HHHH - IxbarreI_HHHH);                   
		
	};

//-------------- vérification de la dérivée avec une dérivée numérique -------------


// calcul des dérivées numériques premières
void   Hyper_W_gene_3D::Calcul_derivee_numerique(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
		                       ,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
		                       ,const  double& jacobien_0,const double& jacobien) 
	{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
   const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); //  "
   const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); //  "
   const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); //  "
	
   Tenseur3BB  gijBBtdt_N; // tenseur  modifié
   Tenseur3HH  gijHHtdt_N; // idem_0 
   const Tenseur3HH & B_HH = gijHH_0;	  
	  // on choisit un incrément
	  double delta = ConstMath::unpeupetit;///10.;
	  double unSurDelta = 1./delta;
	  // sauvegarde des invariants actuelles et de leur variations
	  Vecteur J_r_sauve(J_r) ; 
	  double I_B_sauve=I_B; double II_B_sauve=II_B; double III_B_sauve=III_B;
	  Tableau <Tenseur3HH> d_J_r_epsBB_HH_sauve(d_J_r_epsBB_HH);
	  Tenseur3HH d_I_B_epsBB_HH_sauve=d_I_B_epsBB_HH;
	  Tenseur3HH d_II_B_epsBB_HH_sauve = d_II_B_epsBB_HH;
	  Tenseur3HH d_III_B_epsBB_HH_sauve = d_III_B_epsBB_HH;
	  // dimensionnement
	  Tableau <Tenseur3HH> d_J_r_epsBB_HH_num(d_J_r_epsBB_HH);
	  // on va boucler sur les composantes de gijBB
	  for (int i=1;i<=3;i++) 
	    for (int j=1;j<=3;j++)
	     { gijBBtdt_N = gijBB_tdt;
	       gijBBtdt_N.Coor(i,j) += delta;
	       // en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i)
	       // car le tenseur est symétrique
	       // on a donc en variation numérique la somme des deux dérivées
	       // on définit un coeff multiplicatif qui vaut 1 ou 0.5
	       double coef=1.; if (i != j) coef = 0.5;
	       gijHHtdt_N = gijBBtdt_N.Inverse();
	       double jacobien_N=sqrt(gijBBtdt_N.Det());
	       // calcul des invariants
	       InvariantsEtVar1(gijBB_0_,gijHH_0_,(const TenseurBB &)gijBBtdt_N,(const TenseurHH &)gijHHtdt_N
	                         ,jacobien_0,jacobien_N);
	       for (int r=1;r<=3;r++)
	         {double& der_num = d_J_r_epsBB_HH_num(r).Coor(i,j);
	          // comme on a la somme des dérivées (i,j) et (j,i) qui doivent être égale
	          // on divise par 2 lorsque i != j	          
	          der_num=coef*2.*(J_r(r)-J_r_sauve(r))*unSurDelta; // dérivée numérique
	          double& der = d_J_r_epsBB_HH(r).Coor(i,j); // dérivée analytique
           // comparaison avec les valeurs de dérivées analytiques
           bool erreur = false;
 //          if (r!=2)
           if (diffpourcent(der_num,der,MaX(Dabs(der_num),Dabs(der)),0.1))
            if (MaX(Dabs(der_num),Dabs(der)) > 1.e-5) //1000.*delta)
              {if (MiN(Dabs(der_num),Dabs(der)) == 0.)
                {if ( MaX(Dabs(der_num),Dabs(der)) > 5.e-6 ) erreur = true;} //50.*delta) erreur = true;}
               else erreur = true;
              };
           if (erreur)
             { 
               // calcul des dérivées des I pour voir
               double derI_B_num = coef * 2. * (I_B-I_B_sauve)*unSurDelta;
               double derII_B_num = coef * 2. * (II_B-II_B_sauve)*unSurDelta;
               double derIII_B_num = coef * 2. * (III_B-III_B_sauve)*unSurDelta;
               double& derI_B = d_I_B_epsBB_HH_sauve.Coor(i,j);
               double& derII_B = d_II_B_epsBB_HH_sauve.Coor(i,j);
               double& derIII_B = d_III_B_epsBB_HH_sauve.Coor(i,j);
             
               // calcul de dérivées partielles pour voir
//	           double dernumpartielle = pow(III_B,-1./3.)*2.*(I_B-I_B_sauve)*unSurDelta;
//	           double deranapartielle = 2.*J_r(3)*B_HH(i,j);
//	           double dernumpartielle = I_B * 2. *(pow(III_B,-1./3.)-pow(III_B_sauve,-1./3.))*unSurDelta;
//	           double deranapartielle = -2./3.*pow(III_B_sauve,-1./3.)*gijHH_tdt(i,j);
               double dernumpartielle = 2. *(III_B-III_B_sauve)*unSurDelta;
               double deranapartielle = 2.*III_B_sauve*gijHH_tdt(i,j);
               cout << "\n dernumpartielle= " << dernumpartielle << ", deranapartielle= " << deranapartielle; 
               //
               cout << "\n erreur dans le calcul analytique des derivees de l'invariants: " << r
                    << "\n der_num= " << der_num << " der= " << der <<" der_N=" << d_J_r_epsBB_HH(r)(i,j)
                    << " ietj= " << i << " " << j
                    << "\n derI_B_num= " << derI_B_num << " derI_B= " << derI_B 
                    << " derII_B_num= " << derII_B_num << " derII_B= " << derII_B 
                    << " derIII_B_num= " << derIII_B_num << " derIII_B= " << derIII_B ;
               cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
               cout << "\n un caractere ";
               string toto;
               toto=lect_chaine();
//       Sortie(1);
             };
	         };
	     };
	  // retour aux valeurs sauvegardées
	  J_r = J_r_sauve; 
	  I_B=I_B_sauve; II_B = II_B_sauve; III_B = III_B_sauve;
	  for (int e=1;e<=3;e++)
	    d_J_r_epsBB_HH(e)=d_J_r_epsBB_HH_sauve(e);	
	};

// calcul des dérivées numériques premières et seconde
void   Hyper_W_gene_3D::Calcul_derivee_numerique2(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
		                       ,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
		                       ,const  double& jacobien_0,const double& jacobien)
	{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
   const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); //  "
   const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); //  "
   const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); //  "
	
   Tenseur3BB  gijBBtdt_N; // tenseur  modifié
   Tenseur3HH  gijHHtdt_N; // idem_0 
   const Tenseur3HH & B_HH = gijHH_0;	  
	  double delta = ConstMath::unpeupetit*10.;
	  double unSurDelta = 1./delta;
	  // sauvegarde des invariants actuelles et de leur variations
	  Vecteur J_r_sauve(J_r) ; 
	  double I_B_sauve=I_B; double II_B_sauve=II_B; double III_B_sauve=III_B;
	  Tableau <Tenseur3HH> d_J_r_epsBB_HH_sauve(d_J_r_epsBB_HH);
      // idem pour la partie variation seconde
	  Tenseur3HHHH d_J_1_eps2BB_HHHH_sauve = d_J_1_eps2BB_HHHH;
	  Tenseur3HHHH d_J_2_eps2BB_HHHH_sauve = d_J_2_eps2BB_HHHH; 
	  Tenseur3HHHH d_J_3_eps2BB_HHHH_sauve = d_J_3_eps2BB_HHHH;
      
	  // dimensionnement
	  Tableau <Tenseur3HH> d_J_r_epsBB_HH_num(d_J_r_epsBB_HH);
	  // on va boucler sur les composantes de gijBB
	  for (int i=1;i<=3;i++) 
	    for (int j=1;j<=3;j++)
	     { gijBBtdt_N = gijBB_tdt;
	       gijBBtdt_N.Coor(i,j) += delta;
	       // en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i)
	       // car le tenseur est symétrique
	       // on a donc en variation numérique la somme des deux dérivées
	       // on définit un coeff multiplicatif qui vaut 1 ou 0.5
	       double coef=1.; if (i != j) coef = 0.5;
	       gijHHtdt_N = gijBBtdt_N.Inverse();
	       double jacobien_N=sqrt(gijBBtdt_N.Det());
	       // calcul des invariants
	       InvariantsEtVar1(gijBB_0_,gijHH_0_,(const TenseurBB &)gijBBtdt_N,(const TenseurHH &)gijHHtdt_N
	                         ,jacobien_0,jacobien_N);
	       for (int r=1;r<=3;r++)
	         {double& der_num = d_J_r_epsBB_HH_num(r).Coor(i,j);
	          der_num=coef*2.*(J_r(r)-J_r_sauve(r))*unSurDelta; // dérivée numérique
	          double& der = d_J_r_epsBB_HH(r).Coor(i,j); // dérivée analytique
           // comparaison avec les valeurs de dérivées analytiques
           bool erreur = false;
 //          if (r!=2)
           if (diffpourcent(der_num,der,MaX(Dabs(der_num),Dabs(der)),0.01))
            if (MaX(Dabs(der_num),Dabs(der)) > 1000.*delta)
             {if (MiN(Dabs(der_num),Dabs(der)) == 0.)
                {if ( MaX(Dabs(der_num),Dabs(der)) > 50.*delta) erreur = true;}
              else erreur = true;
             };
           if (erreur)
            { 
              // calcul des dérivées des I pour voir
              cout << "\n erreur dans le calcul analytique des derivees de l'invariants: " << r
                   << "\n der_num= " << der_num << " der= " << der << " ietj= " << i << " " << j
                   << " J_r= " << J_r(r);
              cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
              cout << "\n un caractere ";
              string toto;
              toto=lect_chaine();
            };
	         };
	       // cas des variations secondes
	         // calcul des dérivées secondes analytiques
	       InvariantsEtVar2(gijBB_0_,gijHH_0_,(const TenseurBB &)gijBBtdt_N,(const TenseurHH &)gijHHtdt_N
	                                ,jacobien_0,jacobien_N);
	       // calcul des dérivées numériques et comparaisons
	       for (int k=1;k<=3;k++)
	         for (int l=1;l<=3;l++)
	           { // cas de d_J_1_eps2BB_HHHH
	             double der2num = coef * 2.*(d_J_r_epsBB_HH(1)(k,l) - d_J_r_epsBB_HH_sauve(1)(k,l) )*unSurDelta;
	             double der2ana = 0.5*(d_J_1_eps2BB_HHHH(k,l,i,j) + d_J_1_eps2BB_HHHH(k,l,j,i));
              bool erreur = false;
              if ((Dabs(d_J_r_epsBB_HH(1)(k,l))> 0.8) && (Dabs(der2ana)<30))
                if (diffpourcent(der2num,der2ana,MaX(Dabs(der2num),Dabs(der2ana)),0.01))
                  if (MaX(Dabs(der2num),Dabs(der2ana)) > 1000.*delta) 
                    {if (MiN(Dabs(der2num),Dabs(der2ana)) == 0.)
                       {if ( MaX(Dabs(der2num),Dabs(der2ana)) > 50.*delta) erreur = true;}
                     else erreur = true;
                    };
              if (erreur)
                { 
                  // calcul des dérivées d'éléments intermédiaires pour voir
                  double der2anapartielle=-4./3.*pow(III_B,-1./3.)*(B_HH(k,l)-1./3.*I_B*gijHH_tdt(k,l))*gijHH_tdt(i,j)
                      +2.*pow(III_B,-1./3.)*(-2./3.*gijHH_tdt(k,l)*B_HH(i,j)+2./3.*I_B*gijHH_tdt(k,i)*gijHH_tdt(l,j));
                  double der2numpartielle=coef*2.*(2.*pow(III_B,-1./3.)*(B_HH(k,l)-1./3.*I_B*gijHHtdt_N(k,l))
                                       -(2.*pow(III_B_sauve,-1./3.)*(B_HH(k,l)-1./3.*I_B_sauve*gijHH_tdt(k,l)) ))*unSurDelta;
                  cout << "\n der2numpartielle= " << der2numpartielle << " der2anapartielle= " << der2anapartielle ;
                  
                  
                  
                  //
                  cout << "\n erreur dans le calcul analytique des derivees secondes de l'invariants: " << 1
                       << "\n der2num= " << der2num << " der2ana= " << der2ana 
                       << " klij= " << k << " " << l << " " << i << " " << j
                       << " d_J_r_epsBB_HH(1)(k,l)= " << d_J_r_epsBB_HH(1)(k,l);
                  cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
                  cout << "\n un caractere ";
                  string toto;
                  toto=lect_chaine();
                };
 //   d_J_1_eps2BB_HHHH.Change(k,l,i,j,der2num);     // a virer        
                 // cas de d_J_2_eps2BB_HHHH
	             der2num = coef * 2.*(d_J_r_epsBB_HH(2)(k,l) - d_J_r_epsBB_HH_sauve(2)(k,l) )*unSurDelta;
	             der2ana = 0.5*(d_J_2_eps2BB_HHHH(k,l,i,j) + d_J_2_eps2BB_HHHH(k,l,j,i));
              erreur = false;
              if (diffpourcent(der2num,der2ana,MaX(Dabs(der2num),Dabs(der2ana)),0.01))
                if (MaX(Dabs(der2num),Dabs(der2ana)) > 1000.*delta)
                  {if (MiN(Dabs(der2num),Dabs(der2ana)) == 0.)
                    {if ( MaX(Dabs(der2num),Dabs(der2ana)) > 50.*delta) erreur = true;}
                   else erreur = true;
                  };
              if (erreur)
                { 
                  // calcul des dérivées des I pour voir
                  cout << "\n erreur dans le calcul analytique des derivees  secondes de l'invariants: " << 2
                       << "\n der2num= " << der2num << " der2ana= " << der2ana 
                       << " klij= " << k << " " << l << " " << i << " " << j
                       << " d_J_r_epsBB_HH(2)(k,l)= " << d_J_r_epsBB_HH(2)(k,l);
                  cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
                  cout << "\n un caractere ";
                  string toto;
    //              cin >> toto;
                };
 //   d_J_2_eps2BB_HHHH.Change(k,l,i,j,der2num);     // a virer        
                 // cas de d_J_3_eps2BB_HHHH
	             der2num = coef * 2.*(d_J_r_epsBB_HH(3)(k,l) - d_J_r_epsBB_HH_sauve(3)(k,l) )*unSurDelta;
	             der2ana = 0.5*(d_J_3_eps2BB_HHHH(k,l,i,j) + d_J_3_eps2BB_HHHH(k,l,j,i));
              erreur = false;
              if (diffpourcent(der2num,der2ana,MaX(Dabs(der2num),Dabs(der2ana)),0.1))
                if (MaX(Dabs(der2num),Dabs(der2ana)) > 1000.*delta)
                  {if (MiN(Dabs(der2num),Dabs(der2ana)) == 0.)
                    {if ( MaX(Dabs(der2num),Dabs(der2ana)) > 50.*delta) erreur = true;}
                   else erreur = true;
                  };
              if (erreur)
                 { 
                   // calcul des dérivées des I pour voir
                   cout << "\n erreur dans le calcul analytique des derivees  secondes de l'invariants: " << 3
                        << "\n der2num= " << der2num << " der2ana= " << der2ana 
                        << " klij= " << k << " " << l << " " << i << " " << j
                        << " d_J_r_epsBB_HH(3)(k,l)= " << d_J_r_epsBB_HH(3)(k,l);
                   cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
                   cout << "\n un caractere ";
                   string toto;
       //            cin >> toto;
                 };
  //  d_J_3_eps2BB_HHHH.Change(k,l,i,j,der2num);     // a virer        
  
	           };
	     };
	  // retour aux valeurs sauvegardées
	  J_r = J_r_sauve;
	  I_B=I_B_sauve; II_B = II_B_sauve; III_B = III_B_sauve;
	  for (int e=1;e<=3;e++)
	    d_J_r_epsBB_HH(e)=d_J_r_epsBB_HH_sauve(e);
	  // pour les dérivées secondes  	
	  d_J_1_eps2BB_HHHH = d_J_1_eps2BB_HHHH_sauve;
	  d_J_2_eps2BB_HHHH = d_J_2_eps2BB_HHHH_sauve; 
	  d_J_3_eps2BB_HHHH = d_J_3_eps2BB_HHHH_sauve;
	};