// FICHIER : Hysteresis3D.cc
// CLASSE : Hysteresis3D


// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.


//#include "Debug.h"

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

#include "Hysteresis3D.h"

 // ========== codage des METHODES VIRTUELLES  protegees:================
        // calcul des contraintes 
void Hysteresis3D::Calcul_SigmaHH (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl,
             TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB&  epsBB_,
             TenseurBB& delta_epsBB_,TenseurBB & gijBB_,
		  	 TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,double& ,double& ,
		  	  TenseurHH & sigHH_
		  	,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Expli_t_tdt& )
 {   
    #ifdef MISE_AU_POINT
    if (epsBB_.Dimension() != 3)
       { cout << "\nErreur : la dimension devrait etre 3 !\n";
         cout << " Hysteresis3D::Calcul_SigmaHH\n";
         Sortie(1);
       };
    if (tab_ddl.NbDdl() != d_gijBB_.Taille())
       { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_  !\n";
         cout << " Hysteresis3D::Calcul_SigmaHH\n";
         Sortie(1);
       };
    #endif

    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_); // passage en dim 3
    const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3
    const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_); //   "      "  "  "
    const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_); //   "      "  "  "
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_);  //   "      "  "  "
//    Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t);  //   "      "  "  "
    SaveResulHysteresis3D & save_resul = *((SaveResulHysteresis3D*) saveResul);
    // initialisation des variables de travail
    save_resul.Init_debut_calcul();
    // initialisation éventuelle des variables thermo-dépendantes
    Init_thermo_dependance();
    
    // cas d'une dépendance à la phase de la déformation totale
    if (avec_phase_paradefEqui || avec_phaseEps)
      { xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
	       ParaMateriauxAvecPhaseEpsilon(gijHH,epsBB,xnp,xmu,Qzero);
      };
    
    // === 1 == calcul de l'avancement temporel sur 1 pas,
    int cas = 1; // cas normal avec dérivée de Jauman
//   try { // débug
    Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
// débug---------
/*      }
    catch (ErrNonConvergence_loiDeComportement excep) 
       { // initialisation des variables de travail
         save_resul.Init_debut_calcul();
         // initialisation éventuelle des variables thermo-dépendantes
         Init_thermo_dependance();
         Avancement_temporel(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ);
         }
    catch ( ... )  
       { // initialisation des variables de travail
         save_resul.Init_debut_calcul();
         // initialisation éventuelle des variables thermo-dépendantes
         Init_thermo_dependance();
         Avancement_temporel(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ);
         }
*/			 
//    Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// fin débug -------------
    // on libère les tenseurs intermédiaires        
    LibereTenseur();LibereTenseurQ();        
 };
        
        // calcul des contraintes a t+dt et de ses variations  
void Hysteresis3D::Calcul_DsigmaHH_tdt (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl
            ,BaseB& ,TenseurBB & ,TenseurHH &
            ,BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> &
            ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,TenseurBB & delta_epsBB_
		          ,TenseurBB & gijBB_tdt ,TenseurHH & gijHH_tdt
		          ,Tableau <TenseurBB *>& d_gijBB_tdt
		  	       ,Tableau <TenseurHH *>& d_gijHH_tdt,double& ,double&
		  	       ,Vecteur& ,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
		  	       ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite
            ,double&  module_cisaillement,const Met_abstraite::Impli& )
 {
    #ifdef MISE_AU_POINT
     if (epsBB_tdt.Dimension() != 3)
        { cout << "\nErreur : la dimension devrait etre 3 !\n";
       cout << " Hysteresis3D::Calcul_DsigmaHH_tdt\n";
       Sortie(1);
     };
     if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille())
        { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt  !\n";
       cout << " Hysteresis3D::Calcul_DsigmaHH_tdt\n";
       Sortie(1);
     };
    #endif

    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
    const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3
    const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); //   "      "  "  "
    const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_tdt); //   "      "  "  "
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt);  //   "      "  "  "
//    Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t);  //   "      "  "  "
    SaveResulHysteresis3D & save_resul = *((SaveResulHysteresis3D*) saveResul);
    // initialisation des variables de travail
    save_resul.Init_debut_calcul(); 
    // initialisation éventuelle des variables thermo-dépendantes
    Init_thermo_dependance();
    
    // cas d'une dépendance à la phase de la déformation totale
    if (avec_phase_paradefEqui || avec_phaseEps)
      { xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
	       ParaMateriauxAvecPhaseEpsilon(gijHH,epsBB,xnp,xmu,Qzero);
      };
    
    // === 1 == calcul de l'avancement temporel sur 1 pas,
    int cas = 1; // cas normal avec dérivée de Jauman
    Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
          
    //  ==== 2 ===  calcul de l'opérateur tangent =============
    Hysteresis3D::Dsig_d_ddl(d_sigHH,epsBB_tdt,gijHH_tdt,d_gijBB_tdt,d_gijHH_tdt,d_epsBB );
     
// débug---------
//    Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// fin débug -------------

    // on libère les tenseurs intermédiaires        
    LibereTenseur();LibereTenseurQ();        
 }; 
	 
		  	
// calcul des contraintes et ses variations  par rapport aux déformations a t+dt
// en_base_orthonormee:  le tenseur de contrainte en entrée est  en orthonormee
//                  le tenseur de déformation et son incrémentsont également en orthonormees
//                 si = false: les bases transmises sont utilisées
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void Hysteresis3D::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB&
            ,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB_,double& ,double& 
		  	      ,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps
		  	,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double&  module_cisaillement
		  	,const Met_abstraite::Umat_cont& ex)

 {
   #ifdef MISE_AU_POINT	 	 
	  if (epsBB_tdt.Dimension() != 3)
	    { cout << "\nErreur : la dimension devrait etre 3 !\n";
       cout << " Hysteresis3D::Calcul_dsigma_deps\n";
		     Sortie(1);
     };
    #endif    

    const Tenseur3HH & gijHH = *((Tenseur3HH*) ex.gijHH_tdt); //   "      "  "  "
    const Tenseur3BB & gijBB = *((Tenseur3BB*) ex.gijBB_tdt); //   "      "  "  "
    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
    const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3
    Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt);  //   "      "  "  "
//    Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t);  //   "      "  "  "
    SaveResulHysteresis3D & save_resul = *((SaveResulHysteresis3D*) saveResul);
    // initialisation des variables de travail
    save_resul.Init_debut_calcul(); 
    // initialisation éventuelle des variables thermo-dépendantes
    Init_thermo_dependance();
    
    // cas d'une dépendance à la phase de la déformation totale
    if (avec_phase_paradefEqui || avec_phaseEps)
      { xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
        if (en_base_orthonormee)
	          ParaMateriauxAvecPhaseEpsilon(IdHH3,epsBB,xnp,xmu,Qzero);
        else
           ParaMateriauxAvecPhaseEpsilon(gijHH,epsBB,xnp,xmu,Qzero);
      };
    
    // === 1 == calcul de l'avancement temporel sur 1 pas,
    // ici il n'y a pas de différence entre orthonormeee ou pas au niveau du type de calcul
    int cas = 1; // donc cas avec dérivée de Jauman en mixte
    if (en_base_orthonormee)
       Avancement_temporel(delta_epsBB,IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ);
    else
       Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
      
    //  ==== 2 ===  calcul de l'opérateur tangent =============
//    Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps);

//    try { // débug
    Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps,sigHH);
// débug---------
//      }
//    catch (ErrNonConvergence_loiDeComportement excep)
//       { // initialisation des variables de travail
//         save_resul.Init_debut_calcul();
//         // initialisation éventuelle des variables thermo-dépendantes
//         Init_thermo_dependance();
//         if (en_base_orthonormee)
//            Avancement_temporel(delta_epsBB,IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ);
//         else
//            Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
//         Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps);
//         }
//    catch ( ... )
//       { // initialisation des variables de travail
//         racine.Change_taille(9);
//         save_resul.Init_debut_calcul();
//         // initialisation éventuelle des variables thermo-dépendantes
//         Init_thermo_dependance();
//         if (en_base_orthonormee)
//            Avancement_temporel(delta_epsBB,IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ);
//         else
//            Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
//         Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps);
//       };

//    Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// fin débug -------------



    // on libère les tenseurs intermédiaires        
    LibereTenseur();LibereTenseurQ();        
 }; 
	 
		  	

// =============== fonction protegee ============

		  	
// calcul de la fonction résidu de la résolution de l'équation constitutive
// l'argument test ramène
//     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
//       fatal, qui invalide le calcul du résidu
Vecteur& Hysteresis3D::Residu_constitutif (const double & alpha,const Vecteur & x, int& test)
  { // récupération du vecteur delta_sigma 
    for (int i=1;i<=9;i++) // accroissement de t à tdt
      rdelta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = x(i);
    test = 1; // init par défaut
    // variation de sigma de R à tdt
    rdelta_sigma_barre_BH_Ratdt = sigma_i_barre_BH + rdelta_sigma_barre_tdt_BH - sigma_barre_BH_R;
    delta_barre_alpha_epsBH = alpha * delta_barre_epsBH; // limitation de la charge
    // -- calcul de QdeltaSigma
    QdeltaSigmaRatdt = sqrt(rdelta_sigma_barre_BH_Ratdt.II());
    // -- calculs relatifs aux angles de phases
    Cal_wprime(QdeltaSigmaRatdt,rdelta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR 
               ,wBase,rdelta_sigma_barre_tdt_BH,wprime,wprime_point,trajet_neutre);
    // --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
    if (avec_phaseDeltaSig)
     {xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
      ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,rdelta_sigma_barre_BH_Ratdt,xmu,Qzero);
     };

	   // -- modif 	7 mai 2012 ---
	   // dans le cas où le phidt est négatif, l'équation constitutive n'est plus valide, dans ce cas on considère
	   // que phi est nulle, c-a-d que l'on évolue suivant un trajet neutre ce qui n'est sans doute pas l'idéal, mais permet d'éviter l'explosion
	       // -- calcul de phidt
	   double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH) ; // init partie radial   // -- modif 16 juin 2012 --
	   if (!trajet_neutre)																									   // -- modif 16 juin 2012 --
	     phidt -=  QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime);			// -- modif 16 juin 2012 --

// 	  if ((phidt < 0) && force_phi_zero_si_negatif)  // -- modif 16 juin 2012 --
 	  if (phidt < 0)                                      // -- modif 9 fev 2015 --
     {if (Permet_affichage() > 8)
        cout << "\n  phidt= " << phidt << "\n Hysteresis3D::Residu_constitutif(..." << endl;
      if (  Dabs(phidt) < force_phi_zero_si_negatif)   // -- modif 9 fev 2015 --
       { if (Permet_affichage() > 5)
          {	cout << "\n attention on utilise une valeur de phi negative dans pour le calcul de l'equation constitutive "
                 << " phidt= " << phidt;
            if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
            cout  << "\n Hysteresis3D::Residu_constitutif(..." << endl;
          };
         phidt = 0.;trajet_neutre = true;
         test = -1; // c'est un pb mais pas fatal
       };
     };
	   // fin modif modif 	7 mai 2012 ---

    if (trajet_neutre)
     { // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
       // -- calcul du résidu
       residuBH = rdelta_sigma_barre_tdt_BH - 2. * xmu * delta_barre_alpha_epsBH;
       // --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
       if (avec_defEqui)
        { if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
          double phitau = 0.;
		        if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
          ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
         };
     }
    else
     { // cas d'un trajet non neutre
		     // modif 	7 mai 2012 --- on ne recalcule plus le phi qui a été calculé plus haut pour savoir s'il était négatif
//       // -- calcul de phidt
//       double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH)
//                      - QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime);

       // --- si dépendance à la déformation équivalente 
       if (avec_defEqui)
        { if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
		        if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
          ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phidt,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
        };
       // -- calcul de Beta
       double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
       double beta=0.; // init dans le cas où ne le calcul pas après
       if (QdeltaSigmaRatdt >= ConstMath::petit)
        // on ne calcul beta que si la norme de deltaSigmaRatdt est significative
        {if (xnp == 2.) { beta = - 2. * xmu *  unsurwprimeQ0_puiss_np;}
         else {beta = - 2. * xmu *  unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRatdt,xnp-2.);};
        };
       // -- calcul du résidu
       residuBH = rdelta_sigma_barre_tdt_BH 
                 - 2. * xmu * delta_barre_alpha_epsBH - beta * phidt * rdelta_sigma_barre_BH_Ratdt;
     };
    // -- retour du résidu
    for (int i=1;i<=9;i++) 
       residu(i) = residuBH(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i));
    return residu;
    };  	

// calcul de la matrice tangente de la résolution de l'équation constitutive
// l'argument test ramène
//     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
//       fatal, qui invalide le calcul du résidu et de la dérivée
Mat_abstraite& Hysteresis3D::Mat_tangente_constitutif
                              (const double & alphap,const Vecteur & x, Vecteur& residu, int& test)
  { // récupération du vecteur delta_sigma
    for (int i=1;i<=9;i++) // accroissement de t à tdt
      rdelta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = x(i);
    test = 1; // init par défaut
    rdelta_sigma_barre_BH_Ratdt = sigma_i_barre_BH + rdelta_sigma_barre_tdt_BH - sigma_barre_BH_R;
    delta_barre_alpha_epsBH = alphap * delta_barre_epsBH; // limitation de la charge
    // -- calcul de QdeltaSigmaRatdt
    QdeltaSigmaRatdt = sqrt(rdelta_sigma_barre_BH_Ratdt.II());
    // -- calculs relatifs aux angles de phases
    if (avecVarWprimeS)
     {Cal_wprime_et_der(QdeltaSigmaRatdt,rdelta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR
                      ,wBase,rdelta_sigma_barre_tdt_BH,wprime,d_wprime,wprime_point,d_wprime_point,trajet_neutre);
     }
    else 
     {Cal_wprime(QdeltaSigmaRatdt,rdelta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR 
               ,wBase,rdelta_sigma_barre_tdt_BH,wprime,wprime_point,trajet_neutre);
      d_wprime.Inita(0.); 
      d_wprime_point.Inita(0.);        
     };
    // --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
    if (avec_phaseDeltaSig)
     {xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
      ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,rdelta_sigma_barre_BH_Ratdt,xmu,Qzero);
     };

//---- débug
//    wprime = 1.;wprime_point=0.;trajet_neutre = false;
// --- fin débug 

    // -- modif 	7 mai 2012 ---
    // dans le cas où le phidt est négatif, l'équation constitutive n'est plus valide, dans ce cas on considère
    // que phi est nulle, c-a-d que l'on évolue suivant un trajet neutre ce qui n'est sans doute pas l'idéal, mais permet d'éviter l'explosion
        // calcul de phidt
    double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH) ; // init partie radial   // -- modif 16 juin 2012 -- 
    if (!trajet_neutre)																									   // -- modif 16 juin 2012 --
	        phidt -=  QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime);			// -- modif 16 juin 2012 -- 

// 	  if ((phidt < 0) && force_phi_zero_si_negatif)  // -- modif 16 juin 2012 --
 	  if (phidt < 0)                                      // -- modif 9 fev 2015 --
     {if (Permet_affichage() > 8)
        cout << "\n  phidt= " << phidt << "\n Hysteresis3D::Mat_tangente_constitutif(..." << endl;
      if (  Dabs(phidt) < force_phi_zero_si_negatif)   // -- modif 9 fev 2015 --
       { if (Permet_affichage() > 5)
           {	cout << "\n attention on utilise une valeur de phi negative dans pour le calcul de l'equation constitutive "
                       << " phidt= " << phidt;
             if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
             cout << "\n Hysteresis3D::Mat_tangente_constitutif(..." << endl;
           };
         phidt = 0.;	trajet_neutre=true;
         test = -1; // c'est un pb mais pas fatal
       };
     };
    // fin modif modif 	7 mai 2012 ---

    if (trajet_neutre)
     { // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
       // -- calcul du résidu
       residuBH = rdelta_sigma_barre_tdt_BH - 2. * xmu * delta_barre_alpha_epsBH;
       //  variation du résidu et retour du résidu
       for (int ig=1;ig<=9;ig++)
       {int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
        residu(ig) = residuBH(i,j);
        for (int jg=1;jg<=9;jg++)
         {int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
          derResidu(ig,jg) = IdBH3(i,k)*IdBH3(l,j);
         };
       };
       // --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
       if (avec_defEqui)
        { if (!avec_phaseDeltaSig  && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
          double phitau = 0.;
		    if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
          ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
         };
     }
    else
     { // cas d'un trajet non neutre
		 // modif 	7 mai 2012 --- on ne recalcule plus le phi qui a été calculé plus haut pour savoir s'il était négatif
//       // calcul de phidt
//       double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH)
//                      - QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime);
							 
       // --- si dépendance à la déformation équivalente 
       if (avec_defEqui)
        { if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
		        if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
          ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phidt,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
         };
       // -- calcul de Beta
       double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
       double beta=0.; // init dans le cas où ne le calcul pas après
       double unsurQdeltaSigma= ConstMath::grand; // on limite la valeur par exemple pour le démarrage
       double puiss = 1.; // variable intermédiaire qui resert par la suite
       if (QdeltaSigmaRatdt >= ConstMath::petit)
        // on ne calcul beta que si la norme de deltaSigmaRatdt est significative
        // idem pour unsurQdeltaSigma
        {unsurQdeltaSigma=1./QdeltaSigmaRatdt;
         if (xnp == 2.) { beta = - 2. * xmu *  unsurwprimeQ0_puiss_np;}
         else {puiss = pow(QdeltaSigmaRatdt,xnp-2.);
               beta = - 2. * xmu *  unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRatdt,xnp-2.);
              };
        };
       // calcul du résidu
       residuBH = rdelta_sigma_barre_tdt_BH 
                 - 2. * xmu * delta_barre_alpha_epsBH - beta * phidt * rdelta_sigma_barre_BH_Ratdt;
       
       // -- calcul de la variation du résidu
       // 1 # variation de QdeltaSigmaRatdt ;
       // -- calcul de la variation de QdeltaSigmaRatdt par rapport à sigma_barre        
       Tenseur3BH d_QdeltaSigmaRatdt(unsurQdeltaSigma * rdelta_sigma_barre_BH_Ratdt);
       d_QdeltaSigmaRatdt.PermuteHautBas();
       // 2 # variation de beta
       Tenseur3BH d_beta; // initialisé par défaut à 0
       if (QdeltaSigmaRatdt >= ConstMath::petit)
        // dans le cas ou la norme de deltaSigmaRatdt n'est pas significative, beta est mis à zero
        // il en est de même pour le cas où xnp = 2, pour lequel beta est constant ainsi que sa variation. 
        //  A part ces deux cas  on peut calculer d_beta
        {if (xnp != 2.)
           {double a = 2. * xmu * unsurwprimeQ0_puiss_np * puiss;
            d_beta = (a * (2.-xnp) * unsurQdeltaSigma)  * d_QdeltaSigmaRatdt
                             + ( a *xnp / wprime) * d_wprime;
           }; // si xnp==2: cas où d_beta est nul, déjà fait dans l'initialisation de d_beta
        };
       // 3 # variation de phidt
       Tenseur3BH d_phidt(delta_barre_alpha_epsBH);
       d_phidt.PermuteHautBas();
       d_phidt -= (wprime_point*QdeltaSigmaRatdt/(xmu*wprime))*d_QdeltaSigmaRatdt;
       d_phidt -= (QdeltaSigmaRatdt*QdeltaSigmaRatdt/2. * xmu)*
                   (d_wprime_point/wprime - (wprime_point/(wprime*wprime))*d_wprime);
       // 4 # variation du résidu et retour du résidu
       for (int ig=1;ig<=9;ig++)
        {int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
         residu(ig) = residuBH(i,j);
         for (int jg=1;jg<=9;jg++)
           {int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
            derResidu(ig,jg) = (1.-beta * phidt)*IdBH3(i,k)*IdBH3(l,j)
                             - (d_beta(k,l) * phidt + beta * d_phidt(k,l))
                                *  rdelta_sigma_barre_BH_Ratdt(i,j);
           };
        };
     };// fin du cas non neutre
    calcul_dResi_dsig = true; // on indique que la matrice est calculée 
    return derResidu;
   };	

// calcul de l'opérateur tangent : dsigma/depsilon
TenseurHHHH& Hysteresis3D::Dsig_depsilon
      (bool en_base_orthonormee,const Tenseur3HH & gijHH_tdt, TenseurHHHH& dsig_deps_
       ,TenseurHH & sigHH)
 {  // -- calcul de la matrice MPC
    Dsig_d_(MPC);
    // -- calcul de dsigma / d_eps
    // comme le tenseur des contraintes est symétrique, idem pour le tenseur de déformation
    // on doit avoir un tenseur d'ordre 4 symétrique par rapport aux deux premiers indices
    // et symétrique par rapport aux deux derniers indices
    // donc modif à faire -> on devrait itérer de 1 à 6 et non de 1 à 9 <<<<< !!!!!
// ===> non ! car en fait il y a de la non symétrie qui traîne... voir modif plus bas
    Tenseur3HHHH & dsig_deps = *((Tenseur3HHHH*) &dsig_deps_);
//    TenseurQ3geneHHHH & dsig_deps = *((TenseurQ3geneHHHH*) &dsig_deps_);
// dans le cas d'une base orthonormeee, la position des indices importe peu, donc l'opérateur
    // déjà calculé est ok
    if (en_base_orthonormee)
     {for (int ig=1;ig<=9;ig++)
       {int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
        for (int jg=1;jg<=9;jg++)
          {int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
           dsig_deps.Change(i,j,k,l,-MPC(ig,jg));                       
          };
       };
     }
    else // dans le cas d'une base quelconque il faut intervenir sur le premier et dernier indice
     {
     
//     for (int i=1;i<4;i++)
//       for (int j=1;j<4;j++)
//        for (int k=1;k<4;k++)
//         for (int l=1;l<4;l++)
//          {double inter=0.;
//           for (int e=1;e<4;e++)
//             for (int f=1;f<4;f++)
//               {int ej = Tenseur3BH::OdVect(e,j);
//                int kf = Tenseur3HB::OdVect(k,f); // en fait idem mixte dans un sens ou l'autre
////                inter += gijHH_tdt(i,e) * MPC(ej,kf) * gijHH_tdt(f,l);
//// essai (a virer si pas concluant ) 1 oct 2020
//// a priori c'est concluant, cela fct nettement mieux, en particulier quand on l'utilise
//// en implicit : LoiAdditiveEnSigma::Calcul_DsigmaHH_via_eps_tdt
//// Je ne comprends pas pourquoi la matrice ne conduit pas à une symétrie, mais ...
//// forcé de constater qu'il faut intégrer tous les termes... et donc on force une symétrie !!
//                int je = Tenseur3BH::OdVect(j,e);
//                int fk = Tenseur3HB::OdVect(f,k); // en fait idem mixte dans un sens ou l'autre
//                inter += gijHH_tdt(i,e) *
//                          0.25* (MPC(ej,kf)+MPC(je,kf)+MPC(ej,fk)+MPC(je,fk))
//                          * gijHH_tdt(f,l);
////$$$$$ il manque les variation de gijHH_tdt par rapport à epsilon !!
//
//// fin essai: on garde car a priori concluant
//               }
//           dsig_deps.Change(i,j,k,l,-inter);
//          };
////$$$$ à voir la prise en compte de la partie der de Jauman
//// essai en suivant exactement la même procédure que dsigma/d_ddl


          // nouveau calcul car il manquait des termes
          for (int oj=1;oj<=9;oj++)
           {int o=Tenseur3BH::idx_i(oj); int j=Tenseur3BH::idx_j(oj);
            for (int ef=1;ef<=9;ef++)
             {int e=Tenseur3BH::idx_i(ef); int f=Tenseur3BH::idx_j(ef);
              double inter=0.;
              for (int i=1;i<4;i++)
               for (int k=1;k<4;k++)
                for (int l=1;l<4;l++)
                 { int ij = Tenseur3BH::OdVect(i,j);
                   int kl = Tenseur3BH::OdVect(k,l);
                   // ici un - car MPC(ij,kl) est = - d sig_i^j / d eps_k^l
                   inter -= gijHH_tdt(o,i) * MPC(ij,kl) * IdBH3(k,e) * gijHH_tdt(f,l);
                 };
              // ici - car c'est un - dans la formule
              inter -= 2.* gijHH_tdt(o,e) * sigHH(f,j);
              dsig_ojef_HHHH.Change(o,j,e,f,inter);
             };
           };
          // on symétrise
          dsig_deps.TransfertDunTenseurGeneral(dsig_ojef_HHHH.Symetrise1et2_3et4());
      };

    // retour
    return dsig_deps;
  };	
     
// calcul de l'opérateur tangent : dsigma/d_ddl
void Hysteresis3D::Dsig_d_ddl(Tableau <TenseurHH *>& d_sigHH,TenseurBB & epsBB_tdt
               ,TenseurHH & gijHH_tdt,Tableau <TenseurBB *>& d_gijBB_tdt
               ,Tableau <TenseurHH *>& d_gijHH_tdt,Tableau <TenseurBB *>& d_epsBB )
 {  // -- calcul de la matrice MPC
    Dsig_d_(MPC);
////----- debug
//cout << "\n debug Hysteresis3D::Dsig_d_ddl( ";
//MPC.Affiche();
//cout << endl;
////--- fin debug
    // -- calcul de dsigma / d_ddl
    // nombre de ddl pour la variation du tenseur des contraintes
    int nbddl = d_gijBB_tdt.Taille();
    Tenseur3BH depsBH_dll;Tenseur3BH dsig3BH; Tenseur_ns3HH dsig1_HH;
    const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); //   "      "  "  "
    const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
    for (int i = 1; i<= nbddl; i++)
      { // on fait  uniquement une égalité d'adresse de manière à ne pas utiliser
        // le constructeur d'ou la profusion d'* et de ()
        Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i)));  // passage en dim 1
        const Tenseur3BB & d_gijBB =  *((Tenseur3BB*)(d_gijBB_tdt(i)));  // passage en dim 1
        const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
        const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); //    "
        depsBH_dll = (depsBB * gijHH +   epsBB  * dgijHH );
        
        for (int ig=1;ig<=9;ig++)
         {int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
          dsig3BH.Coor(i,j) = 0.;
          for (int jg=1;jg<=9;jg++)
           {int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
//            dsig3BH(i,j) += -(dRdeps(ig,jg) * depsBH_dll(k,l)); 
            dsig3BH.Coor(i,j) += -(MPC(ig,jg) * depsBH_dll(k,l));
           };
         };
        // et en deux fois contra
        dsig1_HH = gijHH * dsig3BH + dgijHH * sigma_t_barre_tdt ;
        // et on symétrise
        dsigHH = 0.5*(dsig1_HH + dsig1_HH.Transpose());         
      }
 };	

// calcul de la première phase de l'opérateur tangent, est utilisé par Dsig_d_ddl et
// Dsig_depsilon
// calcul de la matrice MPC
void Hysteresis3D::Dsig_d_(MatLapack& MPC)
 {  // dans le cas où la résolution est faite avec du Runge on appelle une routine spécifique
    if (type_calcul_comportement_tangent==2)
      { Dsig_d_Runge(MPC);
    	   return;	
      };
    // ----  sinon il s'agit de l'algo avec linéarisation 
    double const untiers=1./3.;
    double beta=0.; // init dans le cas où ne le calcul pas après
    if (!trajet_neutre)
     { // si c'est un trajet non neutre, on utilise l'ensemble de la loi constitutive
       // -- calcul de Beta
       //     remarque: le xnp que l'on utilise ici est le dernier calculé dans l'algo
       //     d'avancement temporel (on ne le met pas à jour)
       double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
       if (QdeltaSigmaRatdt >= ConstMath::petit)
        // on ne calcul beta que si la norme de deltaSigmaRatdt est significative
        {if (xnp == 2.)
           { beta = - 2. * xmu *  unsurwprimeQ0_puiss_np ;}
         else
          	{ beta = - 2. * xmu *  unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRatdt,xnp-2.);}; 
        };
     }; // fin du calcul de beta pour le cas d'un trajet non neutre
    // normalement la matrice der_at_racine (qui vient d'être calculée
    // pour résoudre l'équation constitutive, contient  d_R_i^j / d_sigbarre_k^l)
// mais ce serait peut-être mieux de systématiquement la recalculer pour alpha = 1
// on aurait ainsi une matrice sécante ???    
    if (!calcul_dResi_dsig)
     {// cas où la matrice tangente n'a pas été calculée
     // calcul de la matrice tangente de la résolution de l'équation constitutive
// essai d'une nouvelle méthode

/* a prioiri ne fonctionne pas
     // on utilise une petite portion de tàtdt
     MatLapack der_at_racine_inter(9,9); 
     int nb_decoup = 10;
     double alphap = 1./nb_decoup; //00; // pour voir 1/100
     
     // sigma_i_barre_BH représente le point de ref sur la même courbe d'évolution 
     // delta_sigma_barre_tdt_BH représente le delta sigma depuis le point de ref
     Tenseur3BH sigma_courant_BH(sigma_i_barre_BH);
     Tenseur3BH sigma_i_barre_sauve_BH(sigma_i_barre_BH); // sauvegarde
     Tenseur3BH incr_sig_BH(delta_sigma_barre_tdt_BH);incr_sig_BH /= nb_decoup;
     for (int ig=1;ig<=9;ig++)
         racine(ig) = incr_sig_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
     for (int i=1;i<=nb_decoup;i++)
     	{ // on prend une nouvelle  référence ad Hoc
          sigma_i_barre_BH = sigma_courant_BH; 
          sigma_courant_BH += incr_sig_BH;
          // on calcul l'opérateur tangent pour le sigma_courant_BH
          der_at_racine_inter = Mat_tangente_constitutif(alphap,racine,residu);
          // on cumule l'opérateur tangent
     	  der_at_racine += der_at_racine_inter; 
     	};
     // calcul de l'opérateur tangent finale
     der_at_racine *= alphap; 
     
     */    
     
     
  // ici une solution qui fonctionne bien pour alpha = 10 !! incompréhensible !! 
/*     double alphap = 10;       
     // récup de la solution
     for (int ig=1;ig<=9;ig++)
//        racine(ig) = alphap * delta_sigma_barre_tdt_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
        racine(ig) =  delta_sigma_barre_tdt_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));

     // on prend une nouvelle  référence, proche de l'arrivée
//     sigma_i_barre_BH = sigma_t_barre_tdt-alphap * delta_sigma_barre_tdt_BH;  
     sigma_i_barre_BH = sigma_t_barre_tdt-delta_sigma_barre_tdt_BH;  
     der_at_racine = Mat_tangente_constitutif(alphap,racine,residu); */
     
// fin essai de la nouvelle méthode 
     // la méthode normale !    
     double alphap = 1.;
     // récup de la solution
     for (int ig=1;ig<=9;ig++)
        racine(ig) = delta_sigma_barre_tdt_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
     int test; // ne sert pas pour l'instant
     der_at_racine = Mat_tangente_constitutif(alphap,racine,residu,test);

// !!!!!!!! pour le débug     
//cout << " norme residu= " << residu.Max_val_abs()/delta_sigma_barre_tdt_BH.MaxiComposante() ;
//if (residu.Max_val_abs() <= ConstMath::petit)
//	   cout << " la norme est nulle !";
// !!!!!!!! pour le débug     

      };
    
    // 1-- on calcul son inverse
    der_at_racine.Inverse(dRdsigMoinsun);
    
// ------ pour le débug
//     MatLapack produit(9,9);
//     produit = dRdsigMoinsun * der_at_racine;
//     for (int i=1;i<=9;i++) 
//      { if (Dabs(produit(i,i)-1. ) > ConstMath::unpeupetit)
//          cout << "\n diag("<<i<<","<<i<<") n'est pas unitaire =" << produit(i,i);
//        for (int j=1;j<= 9; j++)
//       	{ if (i!=j)
//       	   if (Dabs(produit(i,j)) > ConstMath::unpeupetit)
//              cout << "\n  terme ("<<i<<","<<j<<") non nul =" << produit(i,j);
//       	  if ((produit(i,j) != produit(j,i)) && (Dabs(produit(i,j)) > ConstMath::unpeupetit))
//             cout << "\n  non symetrique ("<<i<<","<<j<<") =" << produit(i,j) << " et " << produit(j,i);
//      		
//       	}
//      }
//

//------ fin pour le débug    
    
    
    
    // 2-- maintenant il faut calculer la variation du résidu constitutif / aux déformations
    dRdeps.Initialise(0.);
    
//  ************ important ********************
    // en fait dans la suite on calcul - dSigma/ dEps et non + 
    // ensuite dans l'utilisation on met un moins (ce serait peut-être plus judicieux de mettre moins ici
    // et pas après !!!
    
    
    if (trajet_neutre)
     { // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
       for (int ig=1;ig<=9;ig++)
        {int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
         for (int jg=1;jg<=9;jg++)
           {int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
            dRdeps(ig,jg) = -2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j));
           };
        };
     }
    else
    { // cas non neutre 
       for (int ig=1;ig<=9;ig++)
        {int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
         for (int jg=1;jg<=9;jg++)
           {int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
             dRdeps(ig,jg) = -2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j)) // originale
//            dRdeps(ig,jg) = 2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j))    // essai
// //                         + beta * rdelta_sigma_barre_BH_Ratdt(i,j)
// //                                * rdelta_sigma_barre_BH_Ratdt(l,k);
//                            + beta * delta_sigma_barre_BH_Ratdt(i,j) // originale
                            // normalement suivant les formules, on devrait avoir un - ???? incompréhensible
                            // car avec un + ça marche mieux ??
                            - beta * delta_sigma_barre_BH_Ratdt(i,j)  // essai
                                   * delta_sigma_barre_BH_Ratdt(l,k);
           };
        };
    };
    // 3-- puis on calcul le produit des deux matrices 
    MPC =  dRdsigMoinsun * dRdeps;  // originale 
//    MPC =   -dRdeps * dRdsigMoinsun;   // essai
    
 };
 
// calcul de l'avancement temporel sur 1 pas,
// utilisé par les 3 programmes principaux: 
// Calcul_SigmaHH, Calcul_DsigmaHH_tdt, Calcul_dsigma_deps, 
void Hysteresis3D::Avancement_temporel(const Tenseur3BB & delta_epsBB,const Tenseur3HH & gijHH
                                 ,const Tenseur3BB & gijBB,int cas,SaveResulHysteresis3D & save_resul
                                 ,Tenseur3HH & sigHH,const EnergieMeca & energ_t,EnergieMeca & energ)
{   
    // . le tenseur des contraintes initiale en mixte, qui sera considéré comme la contrainte de début
    // de calcul, pour toute la suite, 
    // . s'il y a plusieurs passage dans la boucle while qui suit,  sigma_i_barre_BH variera, à chaque passage
    // . on utilise le sigmaBH sauvegardé plutôt que le sigmaHH ce qui signifie que l'on effectue un transport
    // en mixte plutôt qu'en 2 fois contravariants !!
 
    // mise en place d'un transport éventuelle sur toutes les grandeurs
    Transport_grandeur(gijBB,gijHH,save_resul);
 
 
    sigma_i_barre_BH = save_resul.sigma_barre_BH_t; // récup des info
    // même procédé pour la def équivalente de début de calcul, initié à t qui sert pour calculer la def équi de R à tdt
    def_i_equi = save_resul.def_equi_at;
    // init énergie
    energ = energ_t;
    double W_a = save_resul.fonction_aide_t; // valeur courante de la fonction d'aide
    double delta_W_a = 0.; // valeur courante du delta W
 
    
    // -- récup et initialisation des variables de travail
    // iatens_princ: pointeur sur la liste déjà enregistrée, iatens: pointeur courant
    List_io <Tenseur3BH>::iterator iatens_princ = save_resul.sigma_barre_BH_R.begin();
    List_io <Tenseur3BH>::iterator iatens = iatens_princ;
    // idem pour la fonction d'aide
    List_io <double>::iterator iafct_princ = save_resul.fct_aide.begin();
    List_io <double>::iterator iafct = iafct_princ;
    aucun_pt_inversion = true; // pour gérer le fait qu'au début il n'y a pas de pt d'inversion !
    if (iatens_princ != save_resul.sigma_barre_BH_R.end())
       aucun_pt_inversion = false; // cas où on a déjà des inversions enregistrées
    // idem pour les déformations équivalentes, fonctionnent avec l'indicateur "aucun_pt_inversion"
    List_io <double>::iterator iadef_equi_princ = save_resul.def_equi.begin();
    List_io <double>::iterator iadef_equi = iadef_equi_princ;
    // idem pour les centres de surfaces de coïncidence
    List_io <Tenseur3BH>::iterator iaOi_princ = save_resul.oc_BH_R.begin();
    List_io <Tenseur3BH>::iterator iaOi = iaOi_princ;
    // on récupère la def maxi de première charge, c'est celle qui correspond au premier point d'inversion
    // si on est sur la première charge on prend la dernière valeur enregistrée sinon on récupère
    // la valeur du premier point d'inversion
    if (aucun_pt_inversion) { defEqui_maxi = def_i_equi;}
    else { List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();
        		 iil++;defEqui_maxi = *iil; // la dernière valeur est nulle, c'est l'avant dernière qui est la bonne
         };
    centre_initial = true; // pour gérer le fait qu'au début le centre initial = tenseur nul
    if (iaOi_princ != save_resul.oc_BH_R.end())
       centre_initial = false; // cas où on a déjà de nouveaux centres enregistrés
    // un pointeur qui indique si les pointeurs courants sont sur les listes principales ou
    // sur les nouvelles listes
    bool pt_sur_principal = true; bool oi_sur_principal = true;

    wBase = save_resul.wBase_t; // initialisation du paramètre de masing
    double unSur_wBaseCarre= 1./(wBase*wBase);
    // -- la déformation
    delta_epsBH = delta_epsBB * gijHH; // incrément de def en mixte
    double trace = delta_epsBH.Trace();
    delta_barre_epsBH = delta_epsBH - (trace/3.)*IdBH3;
    //-- partie spécifique gestion via fct aide
    Tenseur3BH delta_epsBH_totale(delta_epsBH); // sauvegarde delta epsilon totale
	   Tenseur3BH delta_barre_epsBH_totale(delta_barre_epsBH);  // sauvegarde delta_barre epsilon totale
    // a priori, on fait un pas de calcul correspondant à la totalité de la déformation
    // sauf cas inversion ou coïncidence
    double ratio_deformation_en_cours = 0.; // utilisée pour l'instant uniquement avec l'algo fct aide
    double ratio_multiplicatif_coincidence = 1.; // ""
    int phase_coincidence = 0; // pour la gestion de la coincidence
    //-- fin partie spécifique gestion via fct aide
 
    // on définit un indicateur pour le type de gestion de la mémorisation discrète
    bool bascule_fct_aide=false; // init
    bool premier_passage = true; // init
    // indicateur pour gérer la fin du traitement
    bool fin_traitement= false; bool force_fin_traitement = false;
//--------débug- traçage
    #ifdef MISE_AU_POINT
    if (Permet_affichage() == -77)
	      { cout << "\n$$$ === Hysteresis3D::Avancement_temporel: debut de l'algo " << endl ;
       };
    #endif
//----- fin débug
    //  ==== 1 ===  calcul de la contrainte à tdt =============
    // faire tant que fin_traitement n'est pas bon
    int nb_coin_inver=0; // compteur pour éviter une boucle infinie
    while (!fin_traitement)
    {// initialisation des variables de travail qui peuvent varier à chaque passage du while
     // -- 1) sigma_barre_BH_R : initialisé sur la liste principal, mais ensuite peut évoluer sur liste secondaire
     if (aucun_pt_inversion) {sigma_barre_BH_R.Inita(0.);def_equi_R=0.;}
     else                    {sigma_barre_BH_R = *iatens; // la contrainte de référence
                              def_equi_R = *iadef_equi; }; // def équi au pt de référence
     delta_sigma_barre_BH_Rat = sigma_i_barre_BH - sigma_barre_BH_R;
     // -- 2) centre et rayon de la première surface de coïncidence
     if (centre_initial) {oc_BH_R.Inita(0.);}
     else                {oc_BH_R = *iaOi;};

     delta_sigma_barre_BH_OiaR = sigma_barre_BH_R-oc_BH_R; // le tenseur rayon
     // -- 3) calcul du rayon: la première fois le rayon est infini
     if (aucun_pt_inversion) {Q_OiaR = ConstMath::tresgrand;}
     else                    {Q_OiaR = sqrt(delta_sigma_barre_BH_OiaR.II());};
     
     // -- détermination du type de gestion de la mémorisation discrète au premier passage
     // restera constant pour tout le reste du traitement
     // le choix est donc fait par rapport aux valeurs de départ
     if (premier_passage)
      {Tenseur3BH delta_sigma_oi_a_t_BH = sigma_i_barre_BH - oc_BH_R; // OA_t
       double rayon_t = sqrt(delta_sigma_oi_a_t_BH.II());
       if (rayon_t > niveau_bascule_fct_aide * Qzero)
         bascule_fct_aide = true;
       premier_passage = false;
      };

     // le calcul de la contrainte s'effectue par la résolution de l'équation différentielle
     // du schéma constitutif en BH, et toutes les grandeurs sont relatives à BH
     // bien voir que que ce n'est qu'une partie de l'équation, et que pour l'intégration de 
     // la dérivée de Jauman, il faut aussi la partie HB qui est différente, mais on a une symétrie
     // formelle (cf théorie) 
     //  sigma_point = 2*mu*D_barre + beta*phi*(delta_barre de t à R de Sigma) 
	  
// ---7 mai 2012 ---	  
// point nouveau: cette équation n'est valide "que" si phi est positif !!	
// en effet, si phi est négatif, le sigma_point au lieu de saturé, a tendance a exploser (c'est normal car la tendance est inverse)  
   
     // méthode permettant le calcul de sigma à tdt par différente méthodes: linéarisation
     // ou kutta
     // en sortie calcul de : 
     //    sigma_t_barre_tdt, delta_sigma_barre_tdt_BH, delta_sigma_barre_BH_Ratdt
     //    également les grandeurs: QdeltaSigmaRatdt, et wprime,wprime_point, calculées à tdt
//--------débug- traçage
     #ifdef MISE_AU_POINT
     if (Permet_affichage() == -77)
	       { cout << " \nsigR_sur_principal= " << pt_sur_principal << ",      oi_sur_principal= " << oi_sur_principal
               << " save_resul.modif= " << save_resul.modif
               << " \n nb_Oi= "<<save_resul.oc_BH_R.size() << " nb_sigR= "<<save_resul.sigma_barre_BH_R.size()
               << " || nb_Oi_sec= " << save_resul.oc_BH_t_a_tdt.size()
               << " nb_sigR_sec= " << save_resul.sigma_barre_BH_R_t_a_tdt.size()
               << " || wBase= " << wBase;
        };
     #endif
//----- fin débug
     double phi_2_dt = 0.; // on initialise avec un trajet neutre par défaut
     CalculContrainte_tdt(phi_2_dt,save_resul.indicateurs_resolution);

//--------débug- traçage
     #ifdef MISE_AU_POINT
     if (Permet_affichage() == -77)
	       { cout << "\nCalculContrainte_tdt(phi_2_dt...";
          cout << " phi_2_dt = " << setprecision(ParaGlob::NbdigdoEC()) <<  phi_2_dt
               << " || wBase= " << wBase
          << endl ;
        };
     #endif
//----- fin débug
     
     // 2==== gestion de la mémoire discrète                        
     Tenseur3BH delta_sigma_oi_a_tdt_BH = sigma_t_barre_tdt - oc_BH_R; // OA
	    double rayon_tdt = sqrt(delta_sigma_oi_a_tdt_BH.II());
     double d_W_a_essai = 0.; // init
     if (wprime != 0.) // cas d'un trajet non neutre
       {d_W_a_essai = 2. * phi_2_dt / (wprime * wprime);};
//--------débug- traçage
     #ifdef MISE_AU_POINT
     if (Permet_affichage() == -77)
        { cout << "\n rayon_tdt = " << setprecision(ParaGlob::NbdigdoEC()) <<  rayon_tdt
               << " Q_OiaR= " << setprecision(ParaGlob::NbdigdoEC()) <<  Q_OiaR << endl ;
        };
     #endif
//----- fin débug

     if (!bascule_fct_aide) // cas d'une gestion à l'aide des rayons
      { if ((rayon_tdt < Q_OiaR)&&(!force_fin_traitement))
         // cas où le nouveau rayon est inférieur au rayon de référence (= infini si courbe de première charge)
         {int pb_calcul_centre_ok = 0; // init pour une gestion locale du cas où le calcul du centre cercle pose pb,
          
          if (phi_2_dt < 0) //===========================================================================
          // ce qui correspond également au cas d_W_a_essai < O
          // a priori on a détecté un point d'inversion
          // le calcul effectué n'est pas valide et on doit refaire une itération en initialisant les grandeurs ad hoc
          // ------------------ cas d'une inversion ----------
           {// -- on vérifie la cohérence de la valeur de la fonction d'aide (qui n'a rien à voir ici
            // et donc doit être inférieur au maxi)
            if (!aucun_pt_inversion) // si on a déjà des pt d'inversion, on vérifie les niveaux de la fonction d'aide
              if (W_a > *iafct)
                {if (Permet_affichage() > 5)
                  {cout << "\n warning dans l'algorithme de gestion de la memoire discrete,  "
                      << " bien que l'on ait detecte un point d'inversion, il se trouve que la fonction d'aide "
                      << " est superieur au dernier maxi, ce qui signifie que le calcul de l'integrale de la puissance "
                      << " non reversible (phi(t)) n'a pas ete correcte (c'est une erreur cumulee), du sans doute a des pas de temps trop grand, "
                      << " le calcul aura peut-etre des pb si on gere avec la fct d'aide, essayer de le reprendre avec un pas de temps plus petit ! "
                      << " \n W_a= " << W_a  << " iafct= " << *iafct
                      << "\n Hysteresis3D::Avancement_temporel(...";
                   if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
                  };
                };
       //--------débug- traçage
            #ifdef MISE_AU_POINT
            if (Permet_affichage() == -77)
               { cout << "\ninversion " << endl ;
               };
            #endif
       //----- fin débug
            // -- calcul des centre de cercle
            if (!aucun_pt_inversion) // s'il y a au moins déjà un point d'inversion valide on peut traiter
             {Tenseur3BH sig0i_BH;
              pb_calcul_centre_ok=CentreCercle(sig0i_BH,oc_BH_R,delta_sigma_barre_BH_OiaR
                          ,delta_sigma_barre_BH_Rat);
              if (!pb_calcul_centre_ok) // si le calcul a été ok
                {save_resul.oc_BH_t_a_tdt.push_front(sig0i_BH);
                 iaOi = save_resul.oc_BH_t_a_tdt.begin(); // init du centre courant pour la suite
                 oi_sur_principal = false;                // va être mis à jour au début de la prochaine boucle
                 centre_initial = false;                 // du while (!fin_traitement)
                };
             }
            else // si auparavant il n'y avait pas de pt d'inversion, maintenant il y en a un
             {aucun_pt_inversion = false;  // on le signal pour la suite
              // point d'inversion
              wBase = 2; unSur_wBaseCarre = 0.25; //=1./(wBase*wBase);
              defEqui_maxi = def_i_equi;
             };
            // -- s'il n'y a pas eu de pb de calcul de centre on met tout à jour ==> fonctionnement normal
            if (!pb_calcul_centre_ok)
              { // gestion du paramètre "modif" de saveresult
                bool inversion = true;
                Gestion_para_Saveresult_Modif(pt_sur_principal,save_resul,inversion);
           //--------débug- traçage
                #ifdef MISE_AU_POINT
                if (Permet_affichage() == -77)
                   { cout << "\npas de probleme de calcul de centre " ;
                   };
                #endif
           //----- fin débug
                // calcul du nouveau centre d'inversion et sauvegarde, sauf dans le cas du premier point
                // d'inversion, pour lequel le centre reste  encore le tenseur nul pour la suite
                
                // ajout de la contrainte de référence = sigma at
                save_resul.sigma_barre_BH_R_t_a_tdt.push_front(sigma_i_barre_BH);
                save_resul.fct_aide_t_a_tdt.push_front(W_a);
                iafct = save_resul.fct_aide_t_a_tdt.begin();
                W_a=0.; delta_W_a = 0.; // on ré-initialise pour les futures coincidences éventuelles
                iatens = save_resul.sigma_barre_BH_R_t_a_tdt.begin(); // pour l'init du début du while générale
                pt_sur_principal = false; // on indique que maintenant on pointe sur la nouvelle liste
                // on indique le point d'inversion pour les futurs traitements
                sigma_barre_BH_R =save_resul.sigma_barre_BH_R_t_a_tdt.front();
                // -- idem pour la déformation équivalente 
                save_resul.def_equi_t_a_tdt.push_front(def_i_equi);
                iadef_equi = save_resul.def_equi_t_a_tdt.begin(); // pour l'init du début du while générale
                 // on indique la déformation équivalente d'inversion pour les futurs traitements
                def_equi_R  = save_resul.def_equi_t_a_tdt.front();
                // on met à jour la def de début  de calcul (valant t au début, puis tau= pt d'inver ensuite ) pour def_equi_atdt
                def_i_equi =  def_equi_R;                    
                // le calcul d'énergie se fera ensuite dans le cas normale ou avec coïncidence
                // car la valeur actuelle de sigma est erroné, et phi est négatif

                fin_traitement = false; // -- il faut recalculer avec le nouveau point d'inversion
           //--------débug- traçage
                #ifdef MISE_AU_POINT
                if (Permet_affichage() == -77)
                   { cout << " fin_traitement = " << fin_traitement
                          << " save_resul.modif= " << save_resul.modif
                          << endl ;
                   };
                #endif
           //----- fin débug
              };
           // sinon on est dans le cas où il y a eu un pb dans le calcul du centre --> ce qui invalide le point d'inverion
           // on indique qu'il faut traiter la suite comme si il n'y avait pas de point d'inversion
           // ceci se fait au travers de l'indicateur: pb_calcul_centre_ok
            
           };//-- fin du cas de l'inversion
          
          // --- maintenant on traite le cas courant ou il n'y a pas de point d'inversion ---
          if ((phi_2_dt >= 0) || pb_calcul_centre_ok)
           //======================================================================================
           // c-a-d       si phi_2_dt = 0 : on se promène sur un cercle
           //             si phi_2_dt > 0 : progression normale
           //             si pb_calcul_centre_ok != 0 : on est dans un cas limite pour lequel il est préférable de ne pas
           //                                           introduire de point d'inversion
           { // --- cas d'une évolution normale sans inversion ni coincidence ---
             // mais cette évolution peut-être la fin après une coincidence par exemple, donc
             // on ne modifie pas save_resul.modif
             fin_traitement = true;
             // dans l'expression suivante, on utilise une intégration incrémentale ...
             // dans l'expression suivante, on utilise une intégration incrémentale ...
             // s'il y a eu un pb de calcul de centre, cela veut dire que phi_2_dt était légèrement
             // négatif à la prec près, il ne faut pas le compter
             if (pb_calcul_centre_ok) {delta_W_a = 0;}
             else if (wprime != 0.) // cas d'un trajet non neutre
               { delta_W_a =  2. * phi_2_dt / (wprime * wprime);}
             else {delta_W_a = 0.;};
            
             W_a += delta_W_a;
             save_resul.fonction_aide_tdt = W_a;
             // -- uptdate des énergies
             UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
                         ,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
                         ,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R,def_i_equi,def_equi_atdt);
             if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
             else  // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
             // non c'est pas bon
   /*				   defEqui_maxi = 0; // init
        // pour la liste principale, la dernière valeur est nulle, c'est l'avant dernière qui est la bonne
        if (save_resul.def_equi.size())
          { List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
                    // pour la liste secondaire c'est directement la dernière valeur
        if (save_resul.def_equi_t_a_tdt.size())
          { List_io <double>::reverse_iterator iil = save_resul.def_equi_t_a_tdt.rbegin();
            defEqui_maxi = MaX(*iil,defEqui_maxi);
         };*/
                 {   List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
       //--------débug- traçage
            #ifdef MISE_AU_POINT
            if (Permet_affichage() == -77)
               { cout << "\ncas courant sans point d'inversion " ;
                 cout << " fin_traitement = " << fin_traitement
                      << " save_resul.modif= " << save_resul.modif
                      << endl ;
               };
            #endif
       //----- fin débug
            
           };
         } // -- fin du cas où le nouveau rayon est inférieur au rayon de référence
        else if (!force_fin_traitement) // sinon cas d'une coincidence (rayon_tdt >= Q_OiaR )
         { // -- on calcul exactement le point de coïncidence
             // auparavant on sauvegarde le delta eps pour connaître la partie concernée pour le 
           Tenseur3BH delta_eps_sur_deltaTau_BH = delta_barre_epsBH;  // calcul d'énergie
           double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas 
           // directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie sur une valeur =  def_equi_R)
           // mais c'est pour bien marquer le fait que à la sortie de la coïncidence, le pointeur iadef_equi de def_equi à changé
           fin_traitement = Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB,save_resul
                                        ,W_a,iatens,iadef_equi
                                        ,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ,iadef_equi_princ
                                        ,iaOi_princ,iafct_princ,rayon_tdt,false);
       //--------débug- traçage
           #ifdef MISE_AU_POINT
           if (Permet_affichage() == -77)
              { cout << "\ncas d'une coïncidence ";
                cout << " fin_traitement = " << fin_traitement
                     << " save_resul.modif= " << save_resul.modif
                     << " save_resul.nb_coincidence= " << save_resul.nb_coincidence
                     << endl ;
              };
           #endif
       //----- fin débug

           // dans le cas d'un trajet strictement radial, on a deux coïncidences à suivre
           // mais on laisse l'algorithme passer deux fois, car sinon il y a des pb dans le cas 
           // de grand incrément
          
           // -- uptdate des énergies
           // récup du delta eps concerné: delta_barre_epsBH est modifié dans Coïncidence, maintenant
           // il contient le reste de la déformation
           delta_eps_sur_deltaTau_BH -= delta_barre_epsBH; 
    
           // on recalcule la valeur de phi car ici on est après le calcul d'une coïncidence, la contrainte qu'il faut considérer est 
           // celle de la coïncidence, la deformation à considérer est delta_eps_sur_deltaTau_BH, il faut donc calculer le phi qui
           // correspond à ces deux grandeurs.
           delta_sigma_barre_BH_Ra_i = sigma_i_barre_BH - sigma_barre_BH_R;
           double QdeltaSigmaRai = sqrt(delta_sigma_barre_BH_Ra_i.II());
           phi_2_dt = (delta_sigma_barre_BH_Ra_i && delta_eps_sur_deltaTau_BH)
                          - QdeltaSigmaRai * QdeltaSigmaRai * wprime_point / (2.*xmu * wprime);
           // si le phi est négatif, on le met à 0 pour limiter les ennuis dans les calculs d'énergies (ici phi ne sert
           // ne sert que là par la suite), mais c'est quand même pas normal !! mais ça peut être due à une légère erreur
           // au niveau de la coïncidence 
           if (phi_2_dt < 0.) 
              phi_2_dt = 0.;              
                          
                          
           UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
                         ,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_eps_sur_deltaTau_BH
                         ,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
           if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
           else  // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
             // non c'est pas bon
                 { List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
           // on met à jour la def de début  de calcul (valant t au début )
           // il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
           // au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
           def_i_equi =  def_equi_R_ancien;  // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
           // NB: 1)on pourrait également prendre  def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
           // 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
         }
        else // cas où on force la fin du traitement: force_fin_traitement = true
         { fin_traitement = true;
           // comme (rayon_tdt >= Q_OiaR ) on force la coincidence
           double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas 
           // directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie sur une valeur =  def_equi_R)
           // mais c'est pour bien marquer le fait que à la sortie de la coïncidence, le pointeur iadef_equi de def_equi à changé
           Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB,save_resul,W_a,iatens,iadef_equi
                       ,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ,iadef_equi_princ
                       ,iaOi_princ,iafct_princ,rayon_tdt,true);
       //--------débug- traçage
           #ifdef MISE_AU_POINT
           if (Permet_affichage() == -77)
              { cout << "\non force la fin du traitement et la coincidence ";
                cout << " fin_traitement = " << fin_traitement
                     << " save_resul.modif= " << save_resul.modif
                     << endl ;
              };
           #endif
       //----- fin débug
          
           // si le phi est négatif, on le met à 0 pour limiter les ennuis dans les calculs d'énergies (ici phi ne sert
           // ne sert que là par la suite)
           if (phi_2_dt < 0.) 
              phi_2_dt = 0.;              
                          
           // -- uptdate des énergies
           UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
                         ,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
                         ,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
           if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
           else  // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
             // non c'est pas bon
   /*				   defEqui_maxi = 0; // init
        // pour la liste principale, la dernière valeur est nulle, c'est l'avant dernière qui est la bonne
        if (save_resul.def_equi.size())
          { List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
                    // pour la liste secondaire c'est directement la dernière valeur
        if (save_resul.def_equi_t_a_tdt.size())
          { List_io <double>::reverse_iterator iil = save_resul.def_equi_t_a_tdt.rbegin();
            defEqui_maxi = MaX(*iil,defEqui_maxi);
         };*/
                 {   List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;
                 };
           // on met à jour la def de début  de calcul (valant t au début )
           // il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
           // au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
           def_i_equi =  def_equi_R_ancien;  // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
           // NB: 1)on pourrait également prendre  def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
           // 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
         };
      }
	// --- fin de la gestion par rapport au rayon
     else // sinon cas d'une gestion avec la fonction d'aide
      { int pb_calcul_centre_ok = 0; // init pour une gestion locale du cas où le calcul du centre cercle pose pb,
        if (((phase_coincidence==0) && (d_W_a_essai < 0.)) && (!force_fin_traitement))
          { // on a un point d'inversion
            // -- on vérifie la cohérence de la valeur de la fonction d'aide
            if (!aucun_pt_inversion) // si on a déjà des pt d'inversion, on vérifie les niveaux de la fonction d'aide
              if (W_a > *iafct)
                {cout << "\n erreur dans l'algorithme de gestion de la memoire discrete,  "
                      << " bien que l'on ait detecte un point d'inversion, il se trouve que la fonction d'aide "
                      << " est superieur au dernier maxi, ce qui signifie que le calcul de l'integrale de la puissance "
                      << " non reversible (phi(t)) n'a pas ete correcte (c'est une erreur cumulee), du sans doute a des pas de temps trop grand, "
                      << " il n'est pas possible de continuer le calcul, essayer de le reprendre avec un pas de temps plus petit ! "
                      << " \n W_a= " << W_a  << " iafct= " << *iafct
                      << "\n Hysteresis3D::Avancement_temporel(...";
                 if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
                 Sortie(1);
                };
            // -- calcul des centre de cercle
            if (!aucun_pt_inversion) // s'il y a au moins déjà un point d'inversion valide on peut traiter
             {Tenseur3BH sig0i_BH;
              pb_calcul_centre_ok=CentreCercle(sig0i_BH,oc_BH_R,delta_sigma_barre_BH_OiaR
                          ,delta_sigma_barre_BH_Rat);
              if (!pb_calcul_centre_ok) // si le calcul a été ok
                {save_resul.oc_BH_t_a_tdt.push_front(sig0i_BH);
                 iaOi = save_resul.oc_BH_t_a_tdt.begin(); // init du centre courant pour la suite
                 oi_sur_principal = false;                // va être mis à jour au début de la prochaine boucle
                 centre_initial = false;                 // du while (!fin_traitement)
                };
             }
            else // si auparavant il n'y avait pas de pt d'inversion, maintenant il y en a un
             {aucun_pt_inversion = false;  // on le signal pour la suite
              // point d'inversion
              wBase = 2; unSur_wBaseCarre = 0.25; //=1./(wBase*wBase);
              defEqui_maxi = def_i_equi;
             };
            // -- s'il n'y a pas eu de pb de calcul de centre on met tout à jour ==> fonctionnement normal
            if (!pb_calcul_centre_ok)
              { // calcul du nouveau centre d'inversion et sauvegarde, sauf dans le cas du premier point
                // d'inversion, pour lequel le centre reste  encore le tenseur nul pour la suite

                // gestion du paramètre "modif" de saveresult
                bool inversion = true;
                Gestion_para_Saveresult_Modif(pt_sur_principal,save_resul,inversion);
               
                // ajout de la contrainte de référence = sigma at
                save_resul.sigma_barre_BH_R_t_a_tdt.push_front(sigma_i_barre_BH);
                // on ajoute le W_a et non le W_a_essai, car ce dernier n'est pas correcte
                // par contre le W_a est de dernier maxi correcte
                save_resul.fct_aide_t_a_tdt.push_front(W_a);
                W_a=0.; delta_W_a = 0.; // on ré-initialise pour les futures coincidences éventuelles
                iatens = save_resul.sigma_barre_BH_R_t_a_tdt.begin(); // pour l'init du début du while générale
                iafct = save_resul.fct_aide_t_a_tdt.begin();
                pt_sur_principal = false; // on indique que maintenant on pointe sur la nouvelle liste
                // on indique le point d'inversion pour les futurs traitements
                sigma_barre_BH_R =save_resul.sigma_barre_BH_R_t_a_tdt.front();
                // -- idem pour la déformation équivalente 
                save_resul.def_equi_t_a_tdt.push_front(def_i_equi);
                iadef_equi = save_resul.def_equi_t_a_tdt.begin(); // pour l'init du début du while générale
                 // on indique la déformation équivalente d'inversion pour les futurs traitements
                def_equi_R  = save_resul.def_equi_t_a_tdt.front();
                // on met à jour la def de début  de calcul (valant t au début, puis tau= pt d'inver ensuite ) pour def_equi_atdt
                def_i_equi =  def_equi_R;                    
                // le calcul d'énergie se fera ensuite dans le cas normale ou avec coïncidence
                // car la valeur actuelle de sigma est erroné, et phi est négatif

                fin_traitement = false; // -- il faut recalculer avec le nouveau point d'inversion
              };
           // sinon on est dans le cas où il y a eu un pb dans le calcul du centre --> ce qui invalide le point d'inverion
           // on indique qu'il faut traiter la suite comme si il n'y avait pas de point d'inversion
           // ceci se fait au travers de l'indicateur: pb_calcul_centre_ok
             
          };//-- fin du cas de l'inversion
        //-------------------------------------------------
        // arrivée ici on peu éventuellement calculer une valeur de la fonction d'aide correcte
        // dans l'expression suivante, on utilise une intégration incrémentale ...
        // s'il y a eu un pb de calcul de centre, cela veut dire que phi_2_dt était légèrement
        // négatif à la prec près, il ne faut pas le compter
        if (pb_calcul_centre_ok) {delta_W_a = 0;}
        else { delta_W_a =  d_W_a_essai;};
        if (phase_coincidence != 1)
           W_a += delta_W_a;
       
        if (((phase_coincidence==1) || (delta_W_a > 0.)) || force_fin_traitement)
          {if (W_a > *iafct) // on regarde si on est en coïncidence
              phase_coincidence=1;
          
           
           if (((phase_coincidence==0)&&(W_a <= *iafct)) && (!force_fin_traitement))// cas normal après plusieurs inversion
            //======================================================================================
            // c-a-d       si phi_2_dt = 0 : on se promène sur un cercle
            //             si phi_2_dt > 0 : progression normale
            //             si pb_calcul_centre_ok != 0 : on est dans un cas limite pour lequel il est préférable de ne pas
            //                                           introduire de point d'inversion
             { // --- cas d'une évolution normale sans inversion ni coincidence ---
                // mais cette évolution peut-être la fin après une coincidence par exemple, donc
                // on ne modifie pas save_resul.modif
                fin_traitement = true;
                save_resul.fonction_aide_tdt = W_a;
                // -- uptdate des énergies
                UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
                            ,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
                            ,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R,def_i_equi,def_equi_atdt);
                if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
                else  // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
                // non c'est pas bon
                    {   List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
              }
           else if ((phase_coincidence==1) && (!force_fin_traitement))
            //======================================================================================
            // cas où on est en coïncidence: on gère toutes les coïncidence en même temps c-a-d courbe secondaire
            // ou retour sur la courbe principale
            { // -- on calcul exactement le point de coïncidence
              // auparavant on sauvegarde le delta eps pour connaître la partie concernée pour le
              Tenseur3BH delta_eps_sur_deltaTau_BH = delta_barre_epsBH;  // calcul d'énergie
              double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas 
              // directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie
              // sur une valeur =  def_equi_R)
              // mais c'est pour bien marquer le fait que à la sortie de la coïncidence,
              // le pointeur iadef_equi de def_equi à changé
              double ratio_deformation_pour_coincidence=0;
              bool force_coincidence = false; // init par défaut
              if (phase_coincidence ==1)  // si on suit une phase de coincidence qui vient d'être repéré, à la boucle précédente
                force_coincidence = true; // on force au deuxième tour la coincidence
//--- debug


//-- fin debug
                          
              bool coincidence_exacte = Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB
                               ,save_resul,W_a
                               ,iatens,iadef_equi,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ
                               ,iadef_equi_princ,ratio_deformation_pour_coincidence,iaOi_princ,iafct_princ,delta_W_a,force_coincidence);
             
              // si la coïncidence est parfaite, dans ce cas seulement elle a été validée dans l'algo  Coïncidence
              // et on peut updater les énergies, sinon on refait un passage avec une déformation diminuée pour
              // obtenir une coïncidence exacte
              if (coincidence_exacte)
               { // -- cas d'une coïncidence validée
                 // dans le cas d'un trajet strictement radial, on a deux coïncidences à suivre
                 // mais on laisse l'algorithme passer deux fois, car sinon il y a des pb dans le cas
                 // de grand incrément
                
                 // -- uptdate des énergies
                
                 UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
                               ,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
                               ,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
                 if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
                 else  // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
                   // non c'est pas bon
                       { List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
                 // on met à jour la def de début  de calcul (valant t au début )
                 // il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
                 // au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
                 def_i_equi =  def_equi_R_ancien;  // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
                 // NB: 1)on pourrait également prendre  def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
                 // 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
                
                 // Mise en place de la déformation restante pour le calcul qui suit
                 delta_epsBH_totale -= delta_epsBH; delta_epsBH = delta_epsBH_totale;
                 delta_barre_epsBH_totale -= delta_barre_epsBH; delta_barre_epsBH = delta_barre_epsBH_totale;
                 // idem pour la contrainte de départ
                 sigma_i_barre_BH = sigma_t_barre_tdt;
//                 ratio_deformation_en_cours += ratio_multiplicatif_coincidence;
//                 ratio_multiplicatif_coincidence = 1.;
                 phase_coincidence = 0; // on réinitialise pour les autres passages
               }
              else
               { // -- cas non valide, on prépare pour une passe suivante
                 phase_coincidence = 1;
//                 ratio_multiplicatif_coincidence *= ratio_deformation_pour_coincidence;
//                 if (ratio_deformation_en_cours + ratio_multiplicatif_coincidence > 1.)
//                     // on ne peut pas dépasser la déformation imposée on limite donc
//                     ratio_multiplicatif_coincidence = 1.-ratio_deformation_en_cours;
                 delta_barre_epsBH *= ratio_deformation_pour_coincidence;
                 delta_epsBH *= ratio_deformation_pour_coincidence;
              }
            }
           else // cas où on force la fin du traitement: force_fin_traitement = true
            { // ici comme on force une coïncidence exacte, il n'y a pas de calcul de point intermédiaire
              // donc on reprend l'appel du cas des rayons, mais c'est équivalent à l'autre appel de Coincidence
              // avec la fonction d'aide
              fin_traitement = true;
              // comme (rayon_tdt >= Q_OiaR ) on force la coincidence
              double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas 
              // directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie sur une valeur =  def_equi_R)
              // mais c'est pour bien marquer le fait que à la sortie de la coïncidence, le pointeur iadef_equi de def_equi à changé
              Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB,save_resul,W_a,iatens,iadef_equi
                          ,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ,iadef_equi_princ
                          ,iaOi_princ,iafct_princ,rayon_tdt,true);
             
              // si le phi est négatif, on le met à 0 pour limiter les ennuis dans les calculs d'énergies (ici phi ne sert
              // ne sert que là par la suite)
              if (phi_2_dt < 0.) 
                 phi_2_dt = 0.;              
                             
              // -- uptdate des énergies
              UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
                            ,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
                            ,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
              if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
              else  // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
                // non c'est pas bon
                    {   List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;
                    };
              // on met à jour la def de début  de calcul (valant t au début )
              // il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
              // au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
              def_i_equi =  def_equi_R_ancien;  // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
              // NB: 1)on pourrait également prendre  def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
              // 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
            };
          };
      };
        
     // calcul de la contrainte de retour si le traitement est fini
     // et mise à jour a doc
     if (fin_traitement)
        {save_resul.sigma_barre_BH_tdt = sigma_t_barre_tdt;
         save_resul.wBase_tdt = wBase;
         save_resul.def_equi_atdt = def_equi_atdt;
         switch (cas)
          {case 1: 
            	{ // cas de la dérivée de jauman
                  Tenseur_ns3HH sig1HH(gijHH * sigma_t_barre_tdt);
                  sigHH = 0.5*(sig1HH+sig1HH.Transpose()); // jauman
            	  break;	
            	}
           case 2:
            	{ // pas de symétrisation du tenseur
                  sigHH = gijHH * sigma_t_barre_tdt; 
            	  break;	
            	}
           default:
             cout << "\n erreur: cas= " << cas << " non prevu pour le type de derivee de la contrainte ";
             if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
             cout << "\n Hysteresis3D::Avancement_temporel(...)" << endl;
             Sortie(1);      	
          };
       //--------débug- traçage
         #ifdef MISE_AU_POINT
         if (Permet_affichage() == -77)
            { cout << "\n$$$ === fin de l'algo d'avancement : calcul de la contrainte " << endl ;
            };
         #endif
       //----- fin débug
        };
     // on regarde si l'on n'a pas dépassé le nombre de boucle maxi permis: ce qui permet
     // d'éviter les boucles infinies
     nb_coin_inver++;
     if (nb_coin_inver > Abs(nb_maxInvCoinSurUnPas))
     	  { // dans ce cas cela signifie qu'il y a pb et on arrête en générant éventuellement
          // une interruption d'erreur de convergence dans une loi de comportement
          if (Permet_affichage() > 3)
             {cout << "\n warning dans le calcul de la contrainte a t+dt, le nombre de coincidence et/ou"
                   << " d'inversion (" << nb_coin_inver << ") est superieur a la limite fixee: = " 
                   << nb_maxInvCoinSurUnPas ;
              if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
             };
                   
// débug---------
//    Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
//    Sortie(1);
// fin débug -------------
                   
          if (save_resul.nb_coincidence == 0) 
           { // on est dans un cas particulier où l'algo fait du sur-place: un coup il trouve une coïncidence
             // puis le coup après il trouve une inversion --> il revient au même endroit --> boucle infini
             // dans ce cas on décide de manière arbitraire a imposer la coïncidence ce qui garantie la cohérence de la
             // succession des rayons on contraire d'imposer l'inversion, qui conduira à un rayon incohérent (ce qui conduit
             // normalement tout de suite après à une coïncidence)
             force_fin_traitement = true;
             if (Permet_affichage() > 3)
               cout << " on est dans le cas d'aller-retour identiques !! on impose la coincidence " << endl;
           }
          else // on regarde s'il faut imposer la valeur finale obtenue ou générer une erreur
           { if (nb_maxInvCoinSurUnPas > 0)
                { throw ErrNonConvergence_loiDeComportement();break;}  // envoie d'un signal a traiter en dehors
             else  
                { force_fin_traitement = true;
                  if (Permet_affichage() > 3)
                    cout << " on impose la coincidence " << endl;
                }; // cas où on impose coûte que coûte (résultat pas sûr !!)
           };
     	  };
     
    };// -- fin du While (!fin_traitement)
		
};
 
// initialisation éventuelle des variables thermo-dépendantes
void Hysteresis3D::Init_thermo_dependance()
 { // cas d'une thermo dépendance, on calcul les grandeurs en fonction de la température
   if (xnp_temperature != NULL) xnp = xnp_lue = xnp_temperature->Valeur(*temperature);
   if (Qzero_temperature != NULL) Qzero = Qzero_lue = Qzero_temperature->Valeur(*temperature);
   if (xmu_temperature != NULL) xmu = xmu_lue = xmu_temperature->Valeur(*temperature);
//--- pour le debug
//cout << "\n np=" << xnp << " Q0=" << Qzero <<	" xmu=" << xmu  << " T=" << *temperature;
//--- fin debug
   		
 };		  	
 

// méthode permettant le calcul de sigma à tdt par différente méthodes: linéarisation
// ou kutta
// en sortie calcul de : 
//    sigma_t_barre_tdt, delta_sigma_barre_tdt_BH, delta_sigma_barre_BH_Ratdt
//    également les grandeurs: QdeltaSigmaRatdt, et wprime,wprime_point, calculées à tdt
// ramène en retour la valeur de phi_2_dt

void Hysteresis3D::CalculContrainte_tdt(double & phi_2_dt,Tableau<double>& indicateurs_resolution)
	{ // le calcul de la contrainte s'effectue par la résolution de l'équation différentielle
      // du schéma constitutif en BH, et toutes les grandeurs sont relatives à BH
      // bien voir que que ce n'est qu'une partie de l'équation, et que pour l'intégration de 
      // la dérivée de Jauman, il faut aussi la partie HB qui est différente, mais on a une symétrie
      // formelle (cf théorie) 
      //  sigma_point = 2*mu*D_barre + beta*phi*(delta_barre de t à R de Sigma)
  
      if ((sortie_post)&&(indicateurs_resolution.Taille()!= 5)) // dimensionnement éventuelle de la sortie d'indicateurs
        indicateurs_resolution.Change_taille(5); 
      
      calcul_dResi_dsig=false; // on indique que pour l'instant la matrice tangente n'est pas calculée
      // comme la matrice n'est pas forcément définie positive on utilise GAUSS avec MatLapack     
      der_at_racine.Change_Choix_resolution(GAUSS,RIEN_PRECONDITIONNEMENT); 
      bool passage_par_newton = false; // permettra de savoir si on est passé par du newton dans le cas d'un basculement de méthode
      //-- choix de la méthode 
      switch (type_resolution_equa_constitutive)
       {case 1: // cas de la linéarisation de l'équation
         {// ----- pour ce faire on appelle une methode de recherche de zero
          val_initiale.Zero(); // on démarre la recherche à la valeur à t 
          racine.Zero();            // init du résultat à 0.
          der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
          // on met à jour éventuellement le tableau d'indicateurs
          // 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
          int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
          bool conver = false; // init par défaut
          try // on met le bloc sous surveillance
           {
             conver=alg_zero.Newton_raphson
                   (*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
                    ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                    ,maxi_delta_var_sig_sur_iter_pour_Newton);
           }
          catch (ErrSortieFinale)
               // cas d'une direction voulue vers la sortie
               // on relance l'interuption pour le niveau supérieur
             { ErrSortieFinale toto;
               throw (toto);
             }
          catch ( ... ) //(ErrNonConvergence_Newton erreur)
           { double absracinemax=racine.Max_val_abs();
             if (Permet_affichage() > 0)
                {cout << "\n non convergence 2 sur l'algo de la resolution du schema constitutif"
                      << " abs_racine_max " << absracinemax;
                 if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
                 cout << "\n Hysteresis3D::CalculContrainte_tdt(.." << endl;
                };
             LibereTenseur();
             LibereTenseurQ();
             // on génère une exception
             throw ErrNonConvergence_loiDeComportement();break;
           };
          
          if(sortie_post) // sauvegarde éventuelle des indicateurs
            {indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;
////---- debug
//               cout << "\n Hysteresis3D::CalculContrainte_tdt: indicateurs_resolution(1)= "<< indicateurs_resolution(1)
//                    << "  indicateurs_resolution(2)=" << indicateurs_resolution(2)
//                    << " nb_incr_total "<< nb_incr_total << " nb_iter_total "<< nb_iter_total<< endl;
////---- fin debug
            };
//          if (!conver)
          double absracinemax=racine.Max_val_abs();
          if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
            { if (Permet_affichage() > 0)
                {cout << "\n non convergence sur l'algo de la resolution du schema constitutif"
                      << " abs_racine_max " << absracinemax;
                 if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
                 cout << "\n Hysteresis3D::CalculContrainte_tdt(.." << endl;
                };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
            }; 
          // récup de la solution
          for (int ig=1;ig<=9;ig++)
            delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
           
          sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH;
          if (Permet_affichage() > 5)
            // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
            { double sigma_tau_II = sigma_tauBH.II();
              if ( sigma_tau_II > Qzero*Qzero)
                cout << "\n sigma_tau_II= "<< sigma_tau_II << "  Hysteresis3D::CalculContrainte_tdt" << endl;
            };

          // ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
          // car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et 
          // pour t+dt donc pas de pb 
          delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
          passage_par_newton = true; // c'est vrai compte tenu de la méthode
          
          //modif : 11 juin 2014
          // on vérifie que l'amplitude de la contrainte calculée
          {double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
           if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
             {if (Permet_affichage() > 0)
               {cout << "\n erreur10 fatale dans  l'algo de la resolution du schema constitutif"
                     << " concernant le niveau de saturation de la contrainte  "
                     << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase);
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
                cout << "\n Hysteresis3D::CalculContrainte_tdt() (..." << endl;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
             };
          };
          break;
         }
        case 2: // cas d'une résolution par intégration explicite par du kutta
                // si cela ne converge pas on bascule sur de la linéarisation
         {//initialisation  des variables de calcul 
          for (int ig=1;ig<=9;ig++)
            sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
          double tdeb=0.,tfi=1.;
          // calcul de la dérivée initiale
          int erreur=0; //init d'une erreur de calcul de Sigma_point
          Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
//// --- debug ---
//cout << "\n Hysteresis3D::CalculContrainte_tdt avant kutta"
//     << "\n sig_initiale= ";sig_initiale.Affiche();
//cout << "\n dersig_initiale= "; dersig_initiale.Affiche();
//cout <<  endl;      double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());

//// --- fin debug ---
          if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
            { if (Permet_affichage() > 0)
               {cout << "\n erreur0 fatale dans  l'algo de la resolution du schema constitutif"
                     << " au niveau du calcul de la derivee initiale "
                     << "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II())
                     << " au lieu de "<< (Qzero*depassement_Q0*wBase);
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
                cout << "\n Hysteresis3D::CalculContrainte_tdt() (..." << endl;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
            }; 
          double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
          int nombreAppelF=0,nb_step=0;            //  "          "
          double erreur_maxi_global=0.;            //   "
          
          // appel de la fonction kutta
          int conver=alg_edp.Pilotage_kutta
                      (cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                       ,sig_initiale,dersig_initiale
                       ,tdeb,tfi,erreurAbsolue,erreurRelative
                       ,sig_finale,dersig_finale,dernierTemps,dernierdeltat
                       ,nombreAppelF,nb_step,erreur_maxi_global);
          if(sortie_post) // sauvegarde éventuelle des indicateurs
            {indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
             indicateurs_resolution(5)+=erreur_maxi_global;
            };
           // gestion de l'erreur de retour
          if (Permet_affichage() > 5)
            cout << "\n Hysteresis3D::CalculContrainte_tdt "
                 << "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
                 << erreur_maxi_global << endl;

           double abssigmax=sig_finale.Max_val_abs();
           if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
           	{// dans le cas où le Runge ne marche pas on essaie avec la première méthode
             // linéarisation de l'équation
             // ----- pour ce faire on appelle une methode de recherche de zero             
             val_initiale.Zero(); // on démarre la recherche à la valeur à t 
             racine.Zero();            // init du résultat à 0.
             der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
             int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
             // 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
             bool converge=alg_zero.Newton_raphson
                   (*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
                    ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                    ,maxi_delta_var_sig_sur_iter_pour_Newton);
             if(sortie_post) // sauvegarde éventuelle des indicateurs
               {indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;};
             double absracinemax=racine.Max_val_abs();
             if ((converge) && (isfinite(absracinemax)) && (!isnan(absracinemax)))
             	{ // on a convergé, donc on termine le calcul et on sort du switch
                  // récup de la solution
                  for (int ig=1;ig<=9;ig++)
                    delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
           
                  sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH; 
                  // ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
                  // car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et 
                  // pour t+dt donc pas de pb 
                  delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
                  passage_par_newton = true; // on indique donc que l'on est passé par du newton
                  if (Permet_affichage() > 3)
                  	 { cout << "\n probleme dans la resolution de l'equation constitutive avec Runge Kutta";
                      if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
                      if (conver != 2)
                         alg_edp.Affiche_Explication_erreur_RG(conver);
                      if (!isfinite(abssigmax))
                         cout  << " abssigmax n'est pas fini = " << abssigmax ;
                      if (isnan(abssigmax))
                         cout  << " abssigmax = nan = " << abssigmax ;
                  	   if (Permet_affichage() > 3)
                       { cout << "\n  abs_racine_max " << absracinemax;}
                      if (Permet_affichage() > 4)
                       {cout <<  " (2) utilisation avec succes de la linearisation ";};
                      if (Permet_affichage() > 5)
                          cout << "\n Hysteresis3D::CalculContrainte_tdt (..."<<endl;
                    };
                  // puis on met conver à 2 pour signaler que pour la suite que ça a convergé
                  conver=2;        
                  break;
             	}
             else // sinon on n'a pas convergé également avec la linéarisation
              // on indique à toute fin utile que la première n'a pas convergé au cas ou cela n'aurait pas été noté
              { if (conver == 2)
                   conver = 0 ;
              };
           	};
           // arrivé ici cela signifie que: soit kutta a convergé -> conver=2, 
           // soit kutta n'a pas convergé et la méthode de linéarisation n'a pas également convergé -> conver != 2
           // --- si on n'a pas convergé:
           if (conver !=2) // cas d'une non convergence
            { // on appel du kutta adoc sans gestion d'erreur !!
              double deltat=tfi-tdeb;
              switch (cas_kutta)
              	{ case 3:
              	  	{ alg_edp.Runge_Kutta_step23(*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                                               ,sig_initiale,dersig_initiale,tdeb,deltat,sig_finale,estime_erreur);
                    break;                           
              	  	}
                 case 4:
              	  	{ alg_edp.Runge_Kutta_step34(*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                                               ,sig_initiale,dersig_initiale,tdeb,deltat,sig_finale,estime_erreur);
                    break;                           
              	  	}
                 case 5:
              	  	{ alg_edp.Runge_Kutta_step45(*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                                               ,sig_initiale,dersig_initiale,tdeb,deltat,sig_finale,estime_erreur);
                    break;                           
              	  	}
              	  default:
              	   cout << "\n erreur , cas Runge kutta non prevu ! "
              	        << "\n Hysteresis3D::CalculContrainte_tdt(..";
                  if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
              	   Sortie(1);
              	};
              if (Permet_affichage() > 3)
              	{cout << "\n probleme dans la resolution de l'equation constitutive avec Runge Kutta";
                alg_edp.Affiche_Explication_erreur_RG(conver);
                cout <<  ", appel direct de kutta-> erreur estimee= " << estime_erreur(1)
                      << "\n Hysteresis3D::CalculContrainte_tdt (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               //  Sortie(1);
               // pour le debug
     //          if ((Dabs(estime_erreur(1)) > 10) || (Dabs(sig_finale(1)) > 10))
     //             alg_edp.Runge_Kutta_step45(*this,& Hysteresis3D::Sigma_point,sig_initiale,dersig_initiale
     //                    ,tdeb,deltat,sig_finale,estime_erreur);
                 
              	};
              double abssigmax=sig_finale.Max_val_abs();
              if ((!isfinite(abssigmax)) || (isnan(abssigmax)) )
               // dans le cas où la valeur de sig_finale est infinie, ce n'est pas la peine d'aller plus loin
               {
              	   cout << "\n erreur , l'appel directe de Kutta donne une valeur infinie de contrainte ! "
              	        << "\n Hysteresis3D::CalculContrainte_tdt(..";
                  if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                  // on génère une exception
                  throw ErrNonConvergence_loiDeComportement();
              	   Sortie(1);     	
               };
            }
          
           // -- arrivée ici, cela veut dire que: 1) soit le RG de départ a marché correctement,
           // 2) soit la linéarisation a marché correctement, 3) soit si les deux autres n'ont pas marché,
           // on a eu un appel directe à un Rg d'un ordre donné, avec un résultat non contrôlé
           else if (!passage_par_newton)// sinon cela veut dire que c'est le résultat avec Runge que l'on exploite
            {// récup des résultats
             for (int ig=1;ig<=9;ig++)
               sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
             delta_sigma_barre_tdt_BH =  sigma_t_barre_tdt - sigma_i_barre_BH;   
             // ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
             // car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
             // d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
             delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;          
             QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
            };
          // le dernier cas possible avec convergence est celui du passage par newton, et dans ce cas
          // la sauvegarde a déjà été faite
          
          //modif : 11 juin 2014
          // on vérifie que l'amplitude de la contrainte calculée
          {double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
           if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
             {if (Permet_affichage() > 0)
               {cout << "\n erreur11 fatale dans  l'algo de la resolution du schema constitutif"
                     << " concernant le niveau de saturation de la contrainte  "
                     << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
                     << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
             };
          };

          break;
         }
        case 3: // cas d'une résolution par intégration explicite par du kutta
                // mais ici si ça ne converge pas on génère une exception c'est tout , on ne cherche pas
                // a utiliser une autre méthode de secours
         {//initialisation  des variables de calcul 
          for (int ig=1;ig<=9;ig++)
            sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
          double tdeb=0.,tfi=1.;
          int erreur=0; //init d'une erreur de calcul de Sigma_point
          // calcul de la dérivée initiale
          Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
          if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
            { if (Permet_affichage() > 0)
               {cout << "\n erreur1 fatale dans  l'algo de la resolution du schema constitutif"
                     << " au niveau du calcul de la derivee initiale "
                     << "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "<< Qzero
                     << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
            }; 
          
          double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
          int nombreAppelF=0,nb_step=0;            //  "          "
          double erreur_maxi_global=0.;            //   "
          
          // appel de la fonction kutta
          int conver=alg_edp.Pilotage_kutta
                      (cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                       ,sig_initiale,dersig_initiale
                       ,tdeb,tfi,erreurAbsolue,erreurRelative
                       ,sig_finale,dersig_finale,dernierTemps,dernierdeltat
                       ,nombreAppelF,nb_step,erreur_maxi_global);
          if(sortie_post) // sauvegarde éventuelle des indicateurs
            {indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
             indicateurs_resolution(5)+=erreur_maxi_global;
            };
          if (Permet_affichage() > 5)
            cout << "\n Hysteresis3D::CalculContrainte_tdt "
                 << "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
                 << erreur_maxi_global << endl;
           // gestion de l'erreur de retour
           double abssigmax=sig_finale.Max_val_abs();
           if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
             { if (Permet_affichage() > 0)
                 {cout << "\n non convergence sur l'algo de la resolution du schema constitutif"
                       << "\n Hysteresis3D::CalculContrainte_tdt(.." ;
                  if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                 };
                // on génère une exception
                throw ErrNonConvergence_loiDeComportement();break;
             }; 
          // récup des résultats
          for (int ig=1;ig<=9;ig++)
            sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
          delta_sigma_barre_tdt_BH =  sigma_t_barre_tdt - sigma_i_barre_BH;   
          // ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
          // car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
          // d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
          delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;          
          QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
          
          //modif : 11 juin 2014
          // on vérifie que l'amplitude de la contrainte calculée
          {if ( QdeltaSigmaRatdt > (Qzero*depassement_Q0*wBase))
             {if (Permet_affichage() > 0)
               {cout << "\n erreur12 fatale dans  l'algo de la resolution du schema constitutif"
                     << " concernant le niveau de saturation de la contrainte  "
                     << "\n QdeltaSigma =" << QdeltaSigmaRatdt << " au lieu de "<< (Qzero*depassement_Q0*wBase)
                     << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
             };
          };
          
          break;
         }
        case 4: // cas d'une résolution par intégration explicite par du kutta
         {//initialisation  des variables de calcul 
          for (int ig=1;ig<=9;ig++)
            sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
          double tdeb=0.,tfi=1.;
          // calcul de la dérivée initiale
          int erreur=0; //init d'une erreur de calcul de Sigma_point
          Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
          if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
            { if (Permet_affichage() > 0)
               {cout << "\n erreur2 fatale dans  l'algo de la resolution du schema constitutif"
                     << " au niveau du calcul de la derivee initiale "
                     << "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "<< (Qzero*depassement_Q0*wBase)
                     << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
            }; 
          double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
          int nombreAppelF=0,nb_step=0;            //  "          "
          double erreur_maxi_global=0.;            //   "
          
          // appel de la fonction kutta
          int conver=alg_edp.Pilotage_kutta
                      (cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                       ,sig_initiale,dersig_initiale
                       ,tdeb,tfi,erreurAbsolue,erreurRelative
                       ,sig_finale,dersig_finale,dernierTemps,dernierdeltat
                       ,nombreAppelF,nb_step,erreur_maxi_global);
          if(sortie_post) // sauvegarde éventuelle des indicateurs
            {indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
             indicateurs_resolution(5)+=erreur_maxi_global;
            };
          if (Permet_affichage() > 5)
           {cout << "\n Hysteresis3D::CalculContrainte_tdt "
                 << "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
                 << erreur_maxi_global;
            if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
           };
           // gestion de l'erreur de retour
           double abssigmax=sig_finale.Max_val_abs();
           if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
           	{// dans le cas où le Runge ne marche pas on essaie avec la première méthode
             // linéarisation de l'équation
             // ----- pour ce faire on appelle une methode de recherche de zero
             val_initiale.Zero(); // on démarre la recherche à la valeur à t 
             racine.Zero();            // init du résultat à 0.
             der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
             int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
             // 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
             bool converge=alg_zero.Newton_raphson
                   (*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
                    ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                    ,maxi_delta_var_sig_sur_iter_pour_Newton);
             if(sortie_post) // sauvegarde éventuelle des indicateurs
               {indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;};
             double absracinemax=racine.Max_val_abs();
             if ((converge) && (isfinite(absracinemax)) && (!isnan(absracinemax)))
             	{ // on a convergé, donc on termine le calcul et on sort du switch
                  // récup de la solution
                  for (int ig=1;ig<=9;ig++)
                    delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
                  sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH;
                  // ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
                  // car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et 
                  // pour t+dt donc pas de pb 
                  delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
                  passage_par_newton = true; // on indique donc que l'on est passé par du newton
                  if (Permet_affichage() > 3)
                  	 { if (Permet_affichage() > 4)
                       {cout << " -> probleme dans la resolution de l'equation constitutive avec Runge Kutta";
                        alg_edp.Affiche_Explication_erreur_RG(conver);
                        cout <<  ", (1) utilisation avec succes de la linearisation ";
                        cout << "\n  abs_racine_max " << absracinemax;
                        if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;                      
                       };
                      if (Permet_affichage() > 5)
                          cout << "\n Hysteresis3D::CalculContrainte_tdt (..."<<endl;
                    };
                  // puis on met conver à 2 pour signaler que pour la suite que ça a convergé
                  conver=2;        
                  break;
             	}
             else
           	  { if (Permet_affichage() > 3)
                  {cout << "\n non convergence sur l'algo de la resolution du schema constitutif"
                       << " abs_racine_max " << absracinemax << " conver= " << conver << " abssigmax= " << abssigmax;
                   if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                  };
                
                if (Permet_affichage() > 7)
                  cout   << "\n Hysteresis3D::CalculContrainte_tdt(.." << endl;
                // on génère une exception
                throw ErrNonConvergence_loiDeComportement();break;
             }; 
           	}
           else // sinon cela veut dire que le calcul avec Runge est correcte on sauvegarde les résutats
            {// récup des résultats
             for (int ig=1;ig<=9;ig++)
               sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
             delta_sigma_barre_tdt_BH =  sigma_t_barre_tdt - sigma_i_barre_BH;   
             // ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
             // car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
             // d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
             delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;          
             QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
            };
          
          //modif : 11 juin 2014
          // on vérifie que l'amplitude de la contrainte calculée
          {double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
           if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
             {if (Permet_affichage() > 0)
               {cout << "\n erreur13 fatale dans  l'algo de la resolution du schema constitutif"
                     << " concernant le niveau de saturation de la contrainte  "
                     << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
                     << "\n Hysteresis3D::CalculContrainte_tdt() (..." ;
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
             };
          };
          
          break;
         }
        case 5: case 6:// cas d'une résolution par intégration explicite par du kutta
         // l'idée est ici de diminuer la précision jusqu'à un résultat si on ne converge pas normalement
         {//initialisation  des variables de calcul 
          for (int ig=1;ig<=9;ig++)
            sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
          double tdeb=0.,tfi=1.;
          // calcul de la dérivée initiale
          int erreur=0; //init d'une erreur de calcul de Sigma_point
          Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
          if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
            { if (Permet_affichage() > 0)
               {cout << "\n erreur3 fatale dans  l'algo de la resolution du schema constitutif"
                     << " au niveau du calcul de la derivee initiale "
                     << "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "
                     << (Qzero*depassement_Q0*wBase)
                     << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
            }; 
          double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
          int nombreAppelF=0,nb_step=0;            //  "          "
          double erreur_maxi_global=0.;            //   "
          
          // appel de la fonction kutta
          int conver=alg_edp.Pilotage_kutta
                      (cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                       ,sig_initiale,dersig_initiale
                       ,tdeb,tfi,erreurAbsolue,erreurRelative
                       ,sig_finale,dersig_finale,dernierTemps,dernierdeltat
                       ,nombreAppelF,nb_step,erreur_maxi_global);
          if(sortie_post) // sauvegarde éventuelle des indicateurs
            {indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
             indicateurs_resolution(5)+=erreur_maxi_global;
            };
          if (Permet_affichage() > 5)
            cout << "\n Hysteresis3D::CalculContrainte_tdt "
                 << "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
                 << erreur_maxi_global << endl;
          // gestion de l'erreur de retour
          double abssigmax=sig_finale.Max_val_abs();
          if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
           	{// dans le cas où le Runge ne marche pas on essaie avec la première méthode
             // linéarisation de l'équation
             // ----- pour ce faire on appelle une methode de recherche de zero
             val_initiale.Zero(); // on démarre la recherche à la valeur à t 
             racine.Zero();            // init du résultat à 0.
             der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
             int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
             // 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
             bool converge=alg_zero.Newton_raphson
                   (*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
                    ,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
                    ,maxi_delta_var_sig_sur_iter_pour_Newton);
             if(sortie_post) // sauvegarde éventuelle des indicateurs
               {indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;};
//             if (converge)
             double absracinemax=racine.Max_val_abs();
             if ((converge) && (isfinite(absracinemax)) && (!isnan(absracinemax)))
             	{ // on a convergé, donc on termine le calcul et on sort du switch
                  // récup de la solution
                  for (int ig=1;ig<=9;ig++)
                    delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
           
                  sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH; 
                  // ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
                  // car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et 
                  // pour t+dt donc pas de pb 
                  delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
                  passage_par_newton = true; // on indique donc que l'on est passé par du newton
                  if (Permet_affichage() > 3)
                  	 { cout << "\n  estime_erreur= " << estime_erreur(1);
                      if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                      if (Permet_affichage() > 4)
                       {cout << " -> probleme dans la resolution de l'equation constitutive avec Runge Kutta";
                        alg_edp.Affiche_Explication_erreur_RG(conver);
                        cout << ", (3) utilisation avec succes de la linearisation ";
                       };
                      if (Permet_affichage() > 5)
                          cout << "\n Hysteresis3D::CalculContrainte_tdt (..."<<endl;
                    };
                  // puis on met conver à 2 pour signaler que pour la suite que ça a convergé
                  conver=2;        
                  break;
             	}
             else // sinon on n'a pas convergé également avec la linéarisation
              // on indique à toute fin utile que la première n'a pas convergé au cas ou cela n'aurait pas été noté
              { if (conver == 2)
                   conver = 0 ;
              };
           	};
          // arrivé ici cela signifie que: soit kutta a convergé -> conver=2, 
          // soit kutta n'a pas convergé et la méthode de linéarisation n'a pas également convergé -> conver != 2
          if (conver !=2) // cas d'une non convergence
            { // l'idée est de diminuer la précision jusqu'à une solution
              int conver_interne = 0;
              double erreurAsolue_inter=erreurAbsolue;
              double erreurRelative_inter = erreurRelative;
              int nb_boucle=1;
              double abssigmax0; // init
              do
              	{  erreurAsolue_inter *=10;
              	   erreurRelative_inter *=10;
              	   nb_boucle ++;
                  for (int ig=1;ig<=9;ig++)
                     sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
                  double tdeb=0.,tfi=1.;
                  // calcul de la dérivée initiale
                  int erreur=0; //init d'une erreur de calcul de Sigma_point
                  Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
                  if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
                   { if (Permet_affichage() > 0)
                      {cout << "\n erreur4 fatale dans  l'algo de la resolution du schema constitutif"
                            << " au niveau du calcul de la derivee initiale "
                            << "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "<< Qzero
                            << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                       if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                      };
                     // on génère une exception
                     throw ErrNonConvergence_loiDeComportement();break;
                   }; 
                  double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
                  int nombreAppelF=0,nb_step=0;            //  "          "
                  erreur_maxi_global=0.;           
                  // appel de la fonction kutta
                  conver_interne=alg_edp.Pilotage_kutta
                      (cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
                       ,sig_initiale,dersig_initiale
                       ,tdeb,tfi,erreurAsolue_inter,erreurRelative_inter
                       ,sig_finale,dersig_finale,dernierTemps,dernierdeltat
                       ,nombreAppelF,nb_step,erreur_maxi_global);
                  if(sortie_post) // sauvegarde éventuelle des indicateurs
                    {indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
                     indicateurs_resolution(5)+=erreur_maxi_global;
                    };
                  if (Permet_affichage() > 5)
                    cout << "\n Hysteresis3D::CalculContrainte_tdt "
                         << "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
                         << erreur_maxi_global << endl;
              	   if (nb_boucle > 10) // si en * par 10^10 on n'a pas réussi on arrête
              	    { if (Permet_affichage() > 0)
                       {cout << "\n non convergence extreme !! sur l'algo de la resolution du schema constitutif"
                             << "\n Hysteresis3D::CalculContrainte_tdt(.." ;
                        if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                       };
                      // on génère une exception
                      throw ErrNonConvergence_loiDeComportement();break;
                   }; 
                  abssigmax=sig_finale.Max_val_abs();
              	}
              while ((conver_interne != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) );
             
             
              if (Permet_affichage() > 3)
              	{cout << "\n probleme dans la resolution de l'equation constitutive avec Runge Kutta";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                alg_edp.Affiche_Explication_erreur_RG(conver);
                cout  << ", decroissance de la precision -> erreur estimee= "
                      << erreur_maxi_global << " nb cycle= " << nb_boucle
                      << "\n Hysteresis3D::CalculContrainte_tdt (..."<< endl;
              	};
             }
           else if (!passage_par_newton)// sinon cela veut dire que le calcul avec Runge est correcte on sauvegarde les résutats
            {// récup des résultats
             for (int ig=1;ig<=9;ig++)
               sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
             delta_sigma_barre_tdt_BH =  sigma_t_barre_tdt - sigma_i_barre_BH;   
             // ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
             // car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
             // d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
             delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;          
             QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
            };
          // le dernier cas possible avec convergence est celui du passage par newton, et dans ce cas
          // la sauvegarde a déjà été faite
                    //modif : 11 juin 2014
          // on vérifie que l'amplitude de la contrainte calculée
          {double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
           if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
             {if (Permet_affichage() > 0)
               {cout << "\n erreur14 fatale dans  l'algo de la resolution du schema constitutif"
                     << " concernant le niveau de saturation de la contrainte  "
                     << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
                     << "\n Hysteresis3D::CalculContrainte_tdt() (...";
                if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
              // on génère une exception
              throw ErrNonConvergence_loiDeComportement();break;
             };
          };

          break;
         }
       };
      //--- // fin du switch sur les différentes méthodes pour calculer la contraintes à tdt 
      // --- on calcul certaines grandeurs globales et locales
//------ debug------
//cout << "\n dd " << delta_sigma_barre_BH_Ratdt.MaxiComposante() << ", "
//     << " fini ? " << isfinite(delta_sigma_barre_BH_Ratdt.MaxiComposante())
//     << " isnan . " << isnan(delta_sigma_barre_BH_Ratdt.MaxiComposante());
//cout << "\n delta_sigma_barre_BH_Ratdt=== ";
//delta_sigma_barre_BH_Ratdt.Ecriture(cout);
//if  ((!isfinite(delta_sigma_barre_BH_Ratdt.MaxiComposante()))||(isnan(delta_sigma_barre_BH_Ratdt.MaxiComposante())))
// { cout << " delta_sigma_barre_BH_Ratdt= ";
//   delta_sigma_barre_BH_Ratdt.Ecriture(cout);
// };
//----- fin debug -----
      if ((type_resolution_equa_constitutive != 1)&& (!passage_par_newton))
      	{ // sauf le cas 1, tous les autres c'est des kutta  on doit calculer wprime,wprime_point
          // car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
          // d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
          Cal_wprime(QdeltaSigmaRatdt,delta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR 
                     ,wBase,delta_sigma_barre_tdt_BH,wprime,wprime_point,trajet_neutre);      
          // --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
          if (avec_phaseDeltaSig)
             {xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
              ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,delta_sigma_barre_BH_Ratdt,xmu,Qzero);
             };
      		
      	};
      // --<< calcul de phi >>---	
      if (trajet_neutre) // cas où l'on n'a pas un trajet neutre
          { phi_2_dt = 0.;} // pas de variation sur un trajet neutre
      else 
          { phi_2_dt = (delta_sigma_barre_BH_Ratdt && delta_barre_epsBH)
                        - QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2.*xmu * wprime);
//------ debug------
//if ((!isfinite(phi_2_dt))||(isnan(phi_2_dt)))
// { cout << "\n";
//   delta_sigma_barre_BH_Ratdt.Ecriture(cout);
//   delta_barre_epsBH.Ecriture(cout);
//   cout << " q= " << QdeltaSigmaRatdt;
////   Sortie(1);
// };
//// if (phi_2_dt < 0.) cout << "\n phi_2_dt = " << phi_2_dt;
// cout << " " << phi_2_dt;
//----- fin debug -----
           
           
          };
      // --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
      if (avec_defEqui)
          { if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
		          if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
            ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phi_2_dt,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
          };
	};
 
// calcul de l'expression permettant d'obtenir la dérivée temporelle de la contrainte
// en fait il s'agit de l'équation constitutive
// utilisée dans la résolution explicite (runge par exemple) de l'équation constitutive
// erreur : indique si le calcul est valide =0 dans ce cas
//          sinon (diff de 0) indique que les valeurs calculées ne sont pas licite
//          =1:  la norme de sigma est supérieure à la valeur limite de saturation
Vecteur& Hysteresis3D::Sigma_point(const double & tau, const Vecteur & sigma_tau, Vecteur& sig_point,int& erreur)
  {
    double titi= tau; // sert à rien, c'est pour taire le compilo car tau ne sert pas directement
    // récup de la contrainte à tau
    for (int i=1;i<=9;i++)
      sigma_tauBH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = sigma_tau(i);
    //modif : 1 juin 2014
// non on fait la vérif sur le delta cf. plus bas
//    // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
//    double sigma_tau_II = sigma_tauBH.II();
//    if ( sigma_tau_II > Qzero*Qzero)
//     {erreur = 1;
//  //--- pour le debug
//      if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
//        cout << "\n sigma_tau_II =" << sigma_tau_II << " au lieu de "<< Qzero*Qzero
//             << " Hysteresis3D::Sigma_point " << endl ;
//  //--- fin debug
//     }
//    else // sinon on continue 
//     { erreur = 0;
   
  //--- pour le debug
  //cout << "\n tt np=" << xnp << " Q0=" << Qzero <<	" xmu=" << xmu  << " T=" << *temperature;
  //--- fin debug
      // variation de sigma de R à tau
      delta_sigma_barre_R_a_tauBH = sigma_tauBH - sigma_barre_BH_R; // ne sert que dans Sigma_point
      delta_sigma_barre_t_a_tauBH = sigma_tauBH - sigma_i_barre_BH; // ne sert que dans Sigma_point
      // calcul de QdeltaSigma: on prend la valeur absolue car dans les petits nombres on peut avoir des presques 0 !!
      // négatifs qui conduisent à une racine infinie !!
      double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
      
      //modif : 1 juin 2014
      // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
      // on fait la vérification sur le delta, car si on est sur une branche non initiale, on peut très bien
      // dépasser le Qzero, localement, mais ce sera invalidé ensuite avec l'algo de coincidence
      if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
       {erreur = 1;
        if (Permet_affichage() > 5)
          cout << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
               << " Hysteresis3D::Sigma_point " << endl ;
       }
      else // sinon on continue 
       { erreur = 0;
      
      // -- calculs relatifs aux angles de phases : calcul de -> wprime,wprime_point,trajet_neutre
      Cal_wprime(QdeltaSigma,delta_sigma_barre_R_a_tauBH,Q_OiaR,delta_sigma_barre_BH_OiaR 
                 ,wBase,delta_sigma_barre_t_a_tauBH,wprime,wprime_point,trajet_neutre);
      // --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
      if (avec_phaseDeltaSig)
       {xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
        // calcul de xnp, xmu et Qzero
        ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,delta_sigma_barre_R_a_tauBH,xmu,Qzero);
       };

      double phitau = (delta_sigma_barre_R_a_tauBH && delta_barre_epsBH) ; // init partie radial // -- modif 16 juin 2012 --
      if (!trajet_neutre)																							 // -- modif 16 juin 2012 --
        phitau -=  QdeltaSigma * QdeltaSigma * wprime_point / (2. * xmu * wprime);				// -- modif 16 juin 2012 -- 

      // -- modif 	7 mai 2012 ---
      // dans le cas où le phidt est négatif, l'équation constitutive n'est plus valide, dans ce cas on considère
      // que phi est nulle, c-a-d que l'on évolue suivant un trajet neutre ce qui n'est sans doute pas l'idéal, mais permet d'éviter l'explosion
          // -- calcul de phidt
//      if ((phitau < 0) && force_phi_zero_si_negatif)
      if (phitau < 0)                                      // -- modif 9 fev 2015 --
       {if (Permet_affichage() > 8)
          cout << "\n  phidt= " << phitau << "\n Hysteresis3D::Sigma_point(..." << endl;
         if (  Dabs(phitau) < force_phi_zero_si_negatif)   // -- modif 9 fev 2015 --
           { if (Permet_affichage() > 5)
               {	cout << "\n attention on utilise une valeur de phi negative dans pour le calcul de l'equation constitutive "
                      << " phitau= " << phitau
                      << "\n Hysteresis3D::Sigma_point(..." ;
                 if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
               };
             phitau = 0.;trajet_neutre = true;
           };
        };
      // fin modif modif 	7 mai 2012 ---


      if (trajet_neutre)						// -- modif 16 juin 2012 --
       { // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
         // variables intermédiaires pour les tests
         betaphideltasigHB.Inita(0.);
         deuxmudeltaepsHB = (2. * xmu ) * delta_barre_epsBH;
         // formule pour un trajet neutre :   sigma_pointBH = 2. * xmu * delta_barre_epsBH ;
         sigma_pointBH = deuxmudeltaepsHB;
         // --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
         if (avec_defEqui)
          { if (!avec_phaseDeltaSig  && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
            double phitau = 0.;
            if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
            ParaMateriauxAvecDefEqui(xnp,QdeltaSigma,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
          };
       }
      else
       { // cas d'un trajet non neutre
     // modif 	7 mai 2012 --- on ne recalcule plus le phi qui a été calculé plus haut pour savoir s'il était négatif
  //       // --- calcul de phi
  //       double phitau = (delta_sigma_barre_R_a_tauBH && delta_barre_epsBH)
  //                      - QdeltaSigma * QdeltaSigma * wprime_point / (2. * xmu * wprime);
         // --- si dépendance à la déformation équivalente 
         if (avec_defEqui)
          { if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
            if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
              // calcul de xnp (qui est le seul paramètre de retour de la méthode)
              ParaMateriauxAvecDefEqui(xnp,QdeltaSigma,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
          };
         // --- calcul de Beta
         double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
         double beta=0.; // init dans le cas où ne le calcul pas après
         if (QdeltaSigma >= ConstMath::petit)
          // on ne calcul beta que si la norme de QdeltaSigma est significative
          {if (xnp == 2.) { beta = - 2. * xmu *  unsurwprimeQ0_puiss_np;}
           else {beta = - 2. * xmu *  unsurwprimeQ0_puiss_np * pow(QdeltaSigma,xnp-2.);};
          };
         // -- calcul de la dérivée temporelle de la contrainte
         // on considère que D est constant = delta_barre_epsBH/deltat et comme deltat = 1 
         //  ==> = delta_barre_epsBH
         // variables intermédiaires pour les tests
         betaphideltasigHB = beta * phitau * delta_sigma_barre_R_a_tauBH;
         deuxmudeltaepsHB = (2. * xmu ) * delta_barre_epsBH;
         // formule normale :   sigma_pointBH = 2. * xmu * delta_barre_epsBH + beta * phitau * delta_sigma_barre_R_a_tauBH;
         sigma_pointBH = deuxmudeltaepsHB + betaphideltasigHB;
       }; // fin du test sur  trajet_neutre

      
  // si pb : a faire sur toutes les composantes: limitations des dérivées: spécifiques  à la loi d'évolution
      // on limite la variation de la dérivée de la contrainte 
  /*    for (int i=1;i<=3;i++) for (int j=1;j<=3;j++)
        if (deuxmudeltaepsHB(i,j) > 0.)
         {// cas d'une contrainte positive et normalement la dérivée doit évoluer entre 2mu et 0 (saturation)
          if (sigma_pointBH(i,j) < 0. ) sigma_pointBH(i,j) = 0.;
          if (sigma_pointBH(i,j) > deuxmudeltaepsHB(i,j)) sigma_pointBH(i,j) = deuxmudeltaepsHB(i,j);
         }
        else
         {// cas d'une contrainte négative et normalement la dérivée doit évoluer entre -2mu et 0 (saturation)
          if (sigma_pointBH(i,j) < deuxmudeltaepsHB(i,j) ) sigma_pointBH(i,j) = deuxmudeltaepsHB(i,j);
          if (sigma_pointBH(i,j) > 0.) sigma_pointBH(i,j) = 0.;
         }; */
       
  // si pb : a faire sur toutes les composantes: limitations des dérivées: spécifiques  à la loi d'évolution
      // on limite la variation de la dérivée de la contrainte au niveau de la dérivée à l'origine
  // en fait cela introduit une erreur systématique "dramatique !!" sur le résultat qui est alors erroné !!!
  // car il y a accumulation, et il n'y a pas de détection de cette erreur
  // donc au lieu de regarder si l'on est juste supérieur, on met une petite marge ???
  // non pas bon non plus, donc on supprime !!! à voir par la suite
  // sans doute que s'il y a une erreur il faut récupérer une exception à un autre endroit
      if (type_resolution_equa_constitutive==6)  // cas exploratoire
        { double marge = 100.;    
  ////    for (int i=1;i<=3;i++) for (int j=1;j<=3;j++)
  ////        if (Dabs(sigma_pointBH(i,j)) > marge*Dabs(deuxmudeltaepsHB(i,j))) sigma_pointBH(i,j) = deuxmudeltaepsHB(i,j);
  // on tente une autre stratégie: on regarde l'intensité et on limite, parcontre on garde la même direction que celle calculée
          // calcul de l'intensité de cisaillement de la vitesse de déformation
          double deuxmuQdeltaeps = sqrt(deuxmudeltaepsHB.II());
          // on n'intervient que si deuxmuQdeltaeps a une valeur non nulle
          if (deuxmuQdeltaeps > ConstMath::petit)
            if (QdeltaSigma > marge * deuxmuQdeltaeps) // normalement pb
              {sigma_pointBH *= marge*deuxmuQdeltaeps/QdeltaSigma;
               if (Permet_affichage() > 3)
                    {cout << "\n limitation de sigma_point ";
                     if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
                    };
              };
        };
          
      // retour de la dérivée temporelle de la contrainte
      for (int i=1;i<=9;i++) 
         sig_point(i) = sigma_pointBH(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i));
     };
    // retour de la dérivée, qui n'a été calculée que si l'erreur = 0
	   return sig_point;
  };

// vérification de l'intégrité du sigma calculé
// erreur : =0: le calcul est licite, si diff de 0, indique qu'il y a eu une erreur
//          =1:  la norme de sigma est supérieure à la valeur limite de saturation
void Hysteresis3D::Verif_integrite_Sigma(const double & , const Vecteur & sigma_tau,int & erreur)
  { // récup de la contrainte à tau
    for (int i=1;i<=9;i++)
      sigma_tauBH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = sigma_tau(i);
    //modif : 1 juin 2014
// non on fait la vérif sur le delta cf. plus bas
//    // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
//    double sigma_tau_II = sigma_tauBH.II();
//    if ( sigma_tau_II > Qzero*Qzero)
//     {erreur = 1;
//  //--- pour le debug
//      if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
//        cout << "\n sigma_tau_II =" << sigma_tau_II << " au lieu de "<< Qzero*Qzero
//             << " Hysteresis3D::Sigma_point " << endl ;
//  //--- fin debug
//     }
//    else // sinon on continue 
//     { erreur = 0;
   
  //--- pour le debug
  //cout << "\n tt np=" << xnp << " Q0=" << Qzero <<	" xmu=" << xmu  << " T=" << *temperature;
  //--- fin debug
      // variation de sigma de R à tau
      delta_sigma_barre_R_a_tauBH = sigma_tauBH - sigma_barre_BH_R; // ne sert que dans Sigma_point
      delta_sigma_barre_t_a_tauBH = sigma_tauBH - sigma_i_barre_BH; // ne sert que dans Sigma_point
      // calcul de QdeltaSigma: on prend la valeur absolue car dans les petits nombres on peut avoir des presques 0 !!
      // négatifs qui conduisent à une racine infinie !!
      double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
      
      //modif : 1 juin 2014
      // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
      // on fait la vérification sur le delta, car si on est sur une branche non initiale, on peut très bien
      // dépasser le Qzero, localement, mais ce sera invalidé ensuite avec l'algo de coincidence
      if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
       {erreur = 1;
        if (Permet_affichage() > 5)
          cout << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
               << " Hysteresis3D::Verif_integrite_Sigma " << endl ;
       }
      else // sinon on continue 
       { erreur = 0;
     };
  };



// update des différentes énergies sur le pas de temps ou la partie de pas de temps effectué
// et calcul de defEqui finalisée
void Hysteresis3D::UpdateEnergieHyste(EnergieMeca & energ
		                              ,const Tenseur3BH & sigma_deb_BH,const Tenseur3BH & sigma_fin_BH
		                              ,const Tenseur3BH & delta_sigma_sur_tau_BH
		                              ,const Tenseur3BH & delta_sigma_R_tau_BH
		                              ,const double& phi_2_dt
		                              ,const Tenseur3BH & delta_eps_sur_deltaTau_BH,double& QdeltaSigmaRatdt
		                              ,const double& wprime,const double& w_prime_point,const bool& trajet_neutre
		                              ,const double& def_equi_R,const double& def_i_equi, double& def_equi)
{ // delta tau = 1 ici, en fait c'est un paramètre formel
  // double delta_ener_interne = - 0.5*((sigma_deb_BH + sigma_fin_BH) && delta_eps_sur_deltaTau_BH);
  double delta_ener_interne = - (sigma_fin_BH && delta_eps_sur_deltaTau_BH);
  // -------- débug au cas où -----------
  //if (phi_2_dt < 0)
  //  cout << "\n attention, phi_2_dt = " << phi_2_dt << " est negatif ??, normalement, a l'issue d'une convergence"
  //       << " phi est positif ou nul , a regarder en debug "
  //       << "\n Hysteresis3D::UpdateEnergieHyste(.... ";
  // ------- fin débug --------       	

	
  if (trajet_neutre)
    {// -- calcul du delta énergie interne
     // a priori ici tout est élastique: phi = 0, tau_C=0, 
     // -- calcul du taux fondamental C (formule thèse Mayssa Chache p. 71: formule 2.51, avec les formules
     //     2.52, )
     // calcul de la déformation équivalente 
     def_equi = def_i_equi ;
     // -- delta énergies : on considère un tau C constant sur l'incrément de temps ainsi que phi
     double delta_energie_elastique = - delta_ener_interne ;
     double delta_energie_frotteurs = 0.;
     EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
     energ += delta_ener;
    }
  else
    {// -- calcul du delta énergie interne
     // -- calcul du taux fondamental C (formule thèse Mayssa Chache p. 71: formule 2.51, avec les formules
     //     2.52, )
     double unsurwprime = 1./wprime;
     // si QdeltaSigmaRatdt est très petit on considère qu'il n'y a qu'un calcul élastique, pas de dissipation
     // on plus simple on considère que la def équivalente ne bouge pas 
     if (QdeltaSigmaRatdt >= ConstMath::petit )
      { double unsurQdeltaSigma=1./QdeltaSigmaRatdt;
	 	     def_equi = def_i_equi + unsurwprime * phi_2_dt * unsurQdeltaSigma;
	 	     double def_equi_Ratdt = def_equi-def_equi_R;
        // calcul de QpointDeltaSigma
        double QpointDeltaSigma = unsurQdeltaSigma * (delta_sigma_R_tau_BH && delta_sigma_sur_tau_BH);
        // calcul du tau_C qui représente en fait : delta_R^tdt eps : sigma_point
        double tau_C = (wprime * QpointDeltaSigma - w_prime_point * QdeltaSigmaRatdt) * def_equi_Ratdt;
        // -- delta énergies : on considère un tau C constant sur l'incrément de temps ainsi que phi
        double delta_energie_elastique = - delta_ener_interne - unsurwprime * phi_2_dt + unsurwprime * tau_C;
        double delta_energie_frotteurs = unsurwprime * ( phi_2_dt - tau_C);
        // si l'énergie des frotteurs est négative, on considère que c'est due aux différentes imprécisions, 
        // comme je ne vois pas de solution, donc
        // on la met à 0, et on modifie le tau_C en conséquence, et donc
        // on recalcule l'énergie élastique en conséquence, ce qui conduit, tout calcul fait à un trajet 
        // purement élastique ( finalement idem que neutre ??)
        if (delta_energie_frotteurs < 0.)
         { delta_energie_frotteurs = 0.;
            delta_energie_elastique = - delta_ener_interne;
         }
//-----pour le debug
//   if (delta_energie_frotteurs < 0.)
//      cout << "\n (1) delta_energie_frotteurs= " << delta_energie_frotteurs << endl; 
// ----- fin débug

	       // si l'énergie frotteur est négative, on considère que c'est due à une erreur d'arrondie, on met 0 à la place
//	      if (delta_energie_frotteurs < 0.)
//	        delta_energie_frotteurs = 0.;
	       EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
	       energ += delta_ener;
	 	   }
	    else
	 	   { // cas où l'accroissement de sigma est trop faible donc considéré nul
        // on a sigma_point = delta_sigma/delta_t, avec ici delta_t=1 donc delta_sigma = 0 => sigma_point=0
        // comme en fait : tau_C= delta_R^(t+deltat) eps  : sigma_point, on a tau_C = 0 
        // on considère donc tau_C = 0
        double delta_energie_elastique = - delta_ener_interne - unsurwprime * phi_2_dt ;
        double delta_energie_frotteurs = unsurwprime * ( phi_2_dt );
        if (Permet_affichage() > 5)
          if (delta_energie_frotteurs < 0)
            cout << "\n (2) delta_energie_frotteurs= " << delta_energie_frotteurs << endl; 
        EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
        energ += delta_ener;
      }; 
//-------- ancien calcul en commentaire -------	 
//    double unsurQdeltaSigma= ConstMath::grand; // on limite la valeur par exemple pour le démarrage
//     if (QdeltaSigmaRatdt >= ConstMath::petit) unsurQdeltaSigma=1./QdeltaSigmaRatdt;
//	   // calcul de la déformation équivalente 
//	   def_equi = def_equi_R + unsurwprime * phi_2_dt * unsurQdeltaSigma;

//	 // calcul de QpointDeltaSigma
//	 double QpointDeltaSigma = unsurQdeltaSigma * (delta_sigma_R_tau_BH && delta_sigma_sur_tau_BH);
	 // calcul du tau_C qui représente en fait : delta_R^tdt eps : sigma_point
//	 double tau_C = (wprime * QpointDeltaSigma - w_prime_point * QdeltaSigmaRatdt) * def_equi;
	 // -- delta énergies : on considère un tau C constant sur l'incrément de temps ainsi que phi
//	 double delta_energie_elastique = - delta_ener_interne - unsurwprime * phi_2_dt + unsurwprime * tau_C;
//	 double delta_energie_frotteurs = unsurwprime * ( phi_2_dt - tau_C);
//	 EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
//	 energ += delta_ener;
//-------- fin ancien calcul en commentaire -------	 
    };
};

// vérif de la coïncidence: prog de mise au point
// premiere_fois: comme il y a une récursion, pour arrêter la boucle infinie        
void Hysteresis3D::Verif_coincidence(const Tenseur3BB & delta_epsBB,const Tenseur3HH & gijHH
                                 ,const Tenseur3BB & gijBB,int cas,SaveResulHysteresis3D & save_resul
                                 ,Tenseur3HH & sigHH,const EnergieMeca & energ_t,EnergieMeca & energ
                                 ,bool premiere_fois)
{   
    // l'objectif est ici de vérifier que l'on a bien la position de la nouvelle contrainte cohérente avec
    // le rayon de référence
    // --- 1ere étape, récupérer les infos de position de référence qui vont-être enregistrées
    Tenseur3BH sigBH_de_ref; Tenseur3BH centreBH_de_ref;  // init à 0 par défaut
    int aucun_pt_inv=false;
    switch (save_resul.modif)
      { case 0: // rien ne change sur le pas en terme de coincidence et d'inversion 
         { if(save_resul.sigma_barre_BH_R.begin() != save_resul.sigma_barre_BH_R.end())
              {sigBH_de_ref = *(save_resul.sigma_barre_BH_R.begin());}
           else aucun_pt_inv=true;   
           if(save_resul.oc_BH_R.begin() != save_resul.oc_BH_R.end())
              centreBH_de_ref = *(save_resul.oc_BH_R.begin());
           break;
          }
        case 1: // cas où l'on a une ou plusieurs coincidences 
         { // a chaque coincidence il faudra supprimer 1 pt de ref et un centre de cercle
           // sauf si l'on arrive au centre nul ou au point de référence null
	       List_io <Tenseur3BH>::iterator itens;List_io <Tenseur3BH>::iterator ioc_i;
	       itens = save_resul.sigma_barre_BH_R.begin(); ioc_i = save_resul.oc_BH_R.begin();
	       for (int i=1;i<=save_resul.nb_coincidence;i++)
	         { if (itens != save_resul.sigma_barre_BH_R.end()) // si itens == sigma_barre_BH_R.end(), cela signifie
	                itens++;  // que la référence est le tenseur nulle, on arrête de dépiler
	           if (ioc_i != save_resul.oc_BH_R.end())  // si ioc_i == oc_BH_R.end(), cela signifie que le centre de référence
	                ioc_i++;   // est le tenseur nulle 
	          };
	       if(itens != save_resul.sigma_barre_BH_R.end())   
              {sigBH_de_ref = *(itens);}
           else aucun_pt_inv=true;   
	       if (ioc_i != save_resul.oc_BH_R.end())
              centreBH_de_ref = *(ioc_i);
          break;
          }
        case 2: // cas où l'on a une ou plusieurs inversions 
         { // on traite les centres différement des points d'inversions
           // --- cas des centres ---
	       if (save_resul.oc_BH_t_a_tdt.size() != 0)
	          {centreBH_de_ref = *(save_resul.oc_BH_t_a_tdt.begin());} // on récupère le premier
	       else if (save_resul.oc_BH_R.size() != 0)   
	          {centreBH_de_ref = *(save_resul.oc_BH_R.begin());}; // on récupère le premier
	       // sinon tous les centres sont nuls on utilise le centre de départ initialisé à 0   
	       // --- cas des ref ---   
	       if (save_resul.sigma_barre_BH_R_t_a_tdt.size() != 0)
	          {sigBH_de_ref = *(save_resul.sigma_barre_BH_R_t_a_tdt.begin());}
	       else if (save_resul.sigma_barre_BH_R.size() != 0)
	          {sigBH_de_ref = *(save_resul.sigma_barre_BH_R.begin());}
	       else  aucun_pt_inv=true;   
	       
          break;
          }
        case 3: // cas où l'on a une ou plusieurs coincidence(s) et inversion (s)
         { // on traite les centres différement des points d'inversions
           // --- cas des centres ---
	       if (save_resul.oc_BH_t_a_tdt.size() != 0) // cas on ajoutera des inversions 
	          {centreBH_de_ref = *(save_resul.oc_BH_t_a_tdt.begin());} // on récupère le premier
	       else // cas où on a déjà une coincidence et pas d'inversion
	       	{ List_io <Tenseur3BH>::iterator ioc_i = save_resul.oc_BH_R.begin();
	          for (int i=1;i<=save_resul.nb_coincidence;i++)
	             if (ioc_i != save_resul.oc_BH_R.end())  // si ioc_i == oc_BH_R.end(), cela signifie que le centre de référence
	                  ioc_i++;   // est le tenseur nulle 
	          if (ioc_i != save_resul.oc_BH_R.end())
                  centreBH_de_ref = *(ioc_i);
	       	};
	       // --- cas des références ----
	       if (save_resul.sigma_barre_BH_R_t_a_tdt.size() != 0)
	          {sigBH_de_ref = *(save_resul.sigma_barre_BH_R_t_a_tdt.begin());}
	       else // cas où on a déjà une coincidence et pas d'inversion
	       	{ List_io <Tenseur3BH>::iterator itens = save_resul.sigma_barre_BH_R.begin();
	          for (int i=1;i<=save_resul.nb_coincidence;i++)
	            if (itens != save_resul.sigma_barre_BH_R.end()) // si itens == sigma_barre_BH_R.end(), cela signifie
	                itens++;  // que la référence est le tenseur nulle, on arrête de dépiler
	          if(itens != save_resul.sigma_barre_BH_R.end())   
                  {sigBH_de_ref = *(itens);}
              else aucun_pt_inv=true;   
	       	};
           break;
          }         
        default:
         cout << "\n erreur  modif= " << save_resul.modif
              << "\n Hysteresis3D::Verif_coincidence(... ";
         // on génère une exception
         throw ErrNonConvergence_loiDeComportement();
         Sortie(1);
        }       
      
      // --- calcul du rayon de Oi à R  sigBH_de_ref
      Tenseur3BH ddelta_sigma_barre_BH_OiaR = sigBH_de_ref - centreBH_de_ref; // le tenseur rayon
      // calcul du rayon: la première fois le rayon est infini
      double QQ_OiaR=0.;
      if (aucun_pt_inv) {QQ_OiaR = ConstMath::tresgrand;}
      else              {QQ_OiaR = sqrt(ddelta_sigma_barre_BH_OiaR.II());};  

      // --- calcul du nouveau rayon ----
      Tenseur3BH delta_sigma_oi_a_tdt_BH = save_resul.sigma_barre_BH_tdt - centreBH_de_ref; // OA
      double rayon_tdt = sqrt(delta_sigma_oi_a_tdt_BH.II());
      // ---- vérif ---
      if (rayon_tdt > QQ_OiaR) // cas d'un inchoérence
      	{ cout << "\n *********** inchoerence de coincidence ";
         if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
         // re-appel de l'algo pour passer en débug dedans
         // initialisation des variables de travail
         save_resul.Init_debut_calcul();
         // initialisation éventuelle des variables thermo-dépendantes
         Init_thermo_dependance();
         // vérif de l'emboitement
         save_resul.Verif_centre_reference(Permet_affichage());
         // re-calcul pour vérif pas à pas
         Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
         if (premiere_fois)
          {premiere_fois = false;
           Verif_coincidence(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ,premiere_fois);
          };
         // on génère une exception
         throw ErrNonConvergence_loiDeComportement();
         Sortie(1);
      };
};

// vérif de la cohérence des centres et références: les rayons associés doivent être de taille
// décroissante
void Hysteresis3D::SaveResulHysteresis3D::Verif_centre_reference(const int &permet_affichage)
	{ // on va balayer les différents centres et points de références, pour construire les différents rayons
	  // qui devront être de tailles croissantes à partir du premier
	  double rayon_precedent = 0.; List_io <double> les_rayons;
	  List_io <Tenseur3BH>::iterator itens,itens_fin = sigma_barre_BH_R.end();
	  List_io <Tenseur3BH>::iterator ioc_i,ioc_i_fin = oc_BH_R.end();
	  for (itens=sigma_barre_BH_R.begin(), ioc_i = oc_BH_R.begin();
	       itens != itens_fin; itens++,ioc_i++)
	  	{double rayon = 0;
	  	 if (ioc_i != ioc_i_fin)
	  	 	{ Tenseur3BH dsigmaBH = *itens - *ioc_i;
	  	 	  rayon = sqrt(dsigmaBH.II());
	  	 	}
	  	 else 	
	  	 	{ Tenseur3BH dsigmaBH = *itens;
	  	 	  rayon = sqrt(dsigmaBH.II());
	  	 	}
	  	 les_rayons.push_front(rayon);	
	  	 if (rayon >= rayon_precedent)
	  	 	{ rayon_precedent = rayon;}
	  	 else
	  	 	{ cout << "\n ---- ******* attention la croissance des rayons n'est pas verifiee !! " ;
        if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5)) || (permet_affichage > 5))
          {// les listes
           cout << "\n save_resul.sigma_barre_BH_R_t_a_tdt= " << sigma_barre_BH_R_t_a_tdt;
           cout << "\n save_resul.sigma_barre_BH_R= " << sigma_barre_BH_R;
           cout << "\n save_resul.oc_BH_t_a_tdt= " << oc_BH_t_a_tdt;
           cout << "\n save_resul.oc_BH_R= " << oc_BH_R;
           cout << "\n les rayons= " << les_rayons;
           // on vide le buffer
           cout << endl;
          };
       
//	  	      Sortie(1);
	  	 	}; 	
	  	};
	};