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

// idem que Dsig_d_, mais pour le cas où l'on calcul l'avancement à l'aide d'un runge_Kutta
// -> calcul de la matrice M
void Hysteresis3D::Dsig_d_Runge(MatLapack& MPC)
 {  // 1)   on expédie le cas du trajet neutre
    //  ************ 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 !!!
    double const untiers=1./3.;
    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));
            MPC(ig,jg) = -2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j));
           };
        };
       return; // retour  pour faire court
     };
    // 2)   ----  maintenant la suite c'est le cas d'un trajet non-neutre  -----
    // si c'est un trajet non neutre, on utilise l'ensemble de la loi constitutive
    // récupération du vecteur delta_sigma
//    rdelta_sigma_barre_BH_Ratdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH - sigma_barre_BH_R;
    // -- calcul de QdeltaSigmaRatdt
    QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
    // -- calculs relatifs aux angles de phases, en particulier les variations  w_prime
//**** il faudrait regarder si ce calcul est utile vue que wprime et les grandeurs associées ont déjà été calculé
//**** lors de l'équation d'avancement temporel     
    if (avecVarWprimeS)
     { Cal_wprime_et_der(QdeltaSigmaRatdt,delta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR
                      ,wBase,delta_sigma_barre_tdt_BH,wprime,d_wprime,wprime_point,d_wprime_point,trajet_neutre);
     }
    else 
     {Cal_wprime(QdeltaSigmaRatdt,delta_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.);          
     };
    // 3) on reregarde si on est arrivé sur un trajet neutre (ce qui vient d'être déterminé par cal_wprime...)
    //  ************ 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));
            MPC(ig,jg) = -2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j));
           };
        };
       return; // retour  pour faire court
     };
  
    // --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
// en fait on utilise ici les derniers paramètres calculs lors de l'algorithme d'avancement, ce n'est pas
// la peine de les recalculer (perte de temps), idem pour la dépendance à la def équivalente   
//    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 Beta
    // normalement arrivée ici wprime n'est pas nulle car on n'est pas sur un trajet neutre
    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 des termes de la matrice H_BHBH (i,j,f,e)
    double betaQdeltaSig2Sur2muSurWprime = beta * QdeltaSigmaRatdt * QdeltaSigmaRatdt / (2. * xmu*wprime);
    if (avecVarWprimeS) // cas le plus complexe (trajet non radial)
     { 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 e=Tenseur3BH::idx_i(jg); int f=Tenseur3BH::idx_j(jg);
            matriceH(ig,jg) = IdBH3(i,e)*IdBH3(f,j) 
                            + betaQdeltaSig2Sur2muSurWprime * d_wprime(e,f);
           };
        };
       // calcul des termes de la matrice M_BHBH (i,j,l,k)
       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);
            // le moins c'est pour le calcul de - 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 !!!
            matriceM(ig,jg) = -(2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j))
                                + beta * delta_sigma_barre_BH_Ratdt(i,j) * delta_sigma_barre_BH_Ratdt(l,k));
           };
        };
       // -- on calcul l'inverse de H
//    matriceH.Affiche();
       matriceH.Change_Choix_resolution(GAUSS,RIEN_PRECONDITIONNEMENT);
       matriceH.Inverse(matHmoinsun);
       // -- puis on calcul le produit des deux matrices 
       MPC =  matHmoinsun * matriceM;
     }
    else 
     // cas simplifié où le trajet est radial ou bien cas ou on ne prend pas en compte les variations
     // du au cos(alpha) et cos(alpha)Point   
     { // calcul directe des termes de la matrice MPC (i,j,l,k) = M_BHBH (i,j,l,k)
       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);
            // le moins c'est pour le calcul de - 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 !!!
            MPC(ig,jg) = -(2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j))
                                + beta * delta_sigma_barre_BH_Ratdt(i,j) * delta_sigma_barre_BH_Ratdt(l,k));
           };
        };
     };
// ---- vérif débug
/*if (wprime_point==0.)
	{
    double QdeltaSigmaRat = sqrt(delta_sigma_barre_BH_Rat.II());
    double beta_dep=0.; // init dans le cas où ne le calcul pas après
    double unsurQdeltaSigma_dep= ConstMath::grand; // on limite la valeur par exemple pour le démarrage
    double puiss_dep = 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_dep=1./QdeltaSigmaRat;
          // pour l'instant on utilise le w' final
         if (xnp == 2.) { beta_dep = - 2. * xmu *  unsurwprimeQ0_puiss_np;}
         else {puiss = pow(QdeltaSigmaRat,xnp-2.);
               beta_dep = - 2. * xmu *  unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRat,xnp-2.);
               };
        };
	
    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);
         // le moins c'est pour le calcul de - 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 !!!
         MPC(ig,jg) = -(2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j))
                          + beta * delta_sigma_barre_BH_Ratdt(i,j) * delta_sigma_barre_BH_Ratdt(l,k) 
                         );
//                          + 0.5*(beta * rdelta_sigma_barre_BH_Ratdt(i,j) * rdelta_sigma_barre_BH_Ratdt(l,k)
//                       +beta_dep * delta_sigma_barre_BH_Rat(i,j) * delta_sigma_barre_BH_Rat(l,k)));
        };
     };
		
	}
else
 	{
 		cout << "\n ********************************cas avec wprime_point non nul " << endl;
 	}*/	

// ----- fin vérif    
 };

// calcul des paramètres matériaux dépendants de la phase de delta sigma
// en entrée: les paramètres sans phase = paramètres lues ou fonction de la température
// en sortie: les paramètres en tenant compte de la phase de delta_sigma_Oi_tdt
void Hysteresis3D::ParaMateriauxAvecPhaseDeltaSig(double & xnp,const Tenseur3BH&  delta_sigma_barre_BH_OiaR
		                           ,const Tenseur3BH& delta_sig_barre_BH_Ratdt
		                           , double & xmu, double& Qzero )
	{// on calcul delta_sigma_Oi_tdt
	 Tenseur3BH delta_sigmaBarre_BH_Oi_tdt(delta_sigma_barre_BH_OiaR + delta_sig_barre_BH_Ratdt);
	 // on regarde si le tenseur n'est pas trop petit
	 double Qsig =  sqrt(Dabs(delta_sigmaBarre_BH_Oi_tdt.II())); // Dabs: pour les pb de transports éventuels
	 double Qsig3 = Qsig * Qsig * Qsig;
  
  if (mini_Qsig_pour_phase_sigma_Oi_tdt > 0.) // cas historique
	  {if (Qsig > mini_Qsig_pour_phase_sigma_Oi_tdt )
     { // on peut calculer un cos3phi pas débile
       double bIIIb = delta_sigmaBarre_BH_Oi_tdt.III() / 3.;
       double cos3phi = 3. * sqrt(6.)  * bIIIb/ Qsig3;
  //--debug
  //if (abs(cos3phi) > 1.) cout << " cos3phi= " << cos3phi;
  //-- fin débug
          // on limite les valeurs au normal
          if (cos3phi > 1.) { cos3phi = 1.;}
          else if (cos3phi < -1.) { cos3phi = -1.;};
       
       // on calcul les différentes grandeurs dépendantes de la phase
       if (Qzero_phase  != 0)
        { double coefQ  = Qzero_phase->Valeur(cos3phi);
          Qzero *= coefQ;
  //--debug
  //cout << " Qzero= " << Qzero;
  //-- fin débug
        };
       if (xmu_phase  != 0)
        { double coefxmu  = xmu_phase->Valeur(cos3phi);
          xmu *= coefxmu;
        };
       if (xnp_phase  != 0)
        { double coefxnp = xnp_phase->Valeur(cos3phi);
          xnp *= coefxnp;
        };
  //--debug
  //cout << " Qzero= " << Qzero << " xmu= " << xmu << " xnp= " << xnp ;
  //-- fin débug
     };
   }
  else // ici il s'agit du cas avec régularisation
     { // on peut calculer un cos3phi pas débile grace à la régularisation
       double bIIIb = delta_sigmaBarre_BH_Oi_tdt.III() / 3.;
       double cos3phi = 3. * sqrt(6.)  * bIIIb/ (Qsig3-mini_Qsig_pour_phase_sigma_Oi_tdt);
  //--debug
  //if (abs(cos3phi) > 1.) cout << " cos3phi= " << cos3phi;
  //-- fin débug
       // on limite les valeurs au normal
       if (cos3phi > 1.) { cos3phi = 1.;}
       else if (cos3phi < -1.) { cos3phi = -1.;};
       
       // on calcul les différentes grandeurs dépendantes de la phase
       if (Qzero_phase  != 0)
        { double coefQ  = Qzero_phase->Valeur(cos3phi);
          Qzero *= coefQ;
  //--debug
  //cout << " Qzero= " << Qzero;
  //-- fin débug
        };
       if (xmu_phase  != 0)
        { double coefxmu  = xmu_phase->Valeur(cos3phi);
          xmu *= coefxmu;
        };
       if (xnp_phase  != 0)
        { double coefxnp = xnp_phase->Valeur(cos3phi);
          xnp *= coefxnp;
        };
  //--debug
  //cout << " Qzero= " << Qzero << " xmu= " << xmu << " xnp= " << xnp ;
  //-- fin débug
     };
	 // sinon on ne fait rien
	};


// calcul des paramètres matériaux dépendants de la déformation équivalente
// en entrée: les paramètres sans dépendance = paramètres lues ou fonction de la température ou phase
// en sortie: les paramètres en tenant compte de la dépendance 
void Hysteresis3D::ParaMateriauxAvecDefEqui(double & xnp,const double& QdeltaSigmaRatdt,const double& phi_2_dt
		                           ,const double& wprime,const double& def_equi_R
		                           ,const bool& trajet_neutre,const double& def_i_equi,const double& defEqui_maxi)
{ // calcul de la def équivalente
  double defEqui = 0.; // init

  // la déformation équivalente calculée ici est approchée. Ce n'est pas celle sauvegardée, cette dernière est calculée en fait au niveau
  // du calcul de l'énergie (mais suivant une formule équivalente)
  if (trajet_neutre)
    {// calcul de la déformation équivalente 
	    defEqui = def_i_equi ; // le long d'un trajet neutre la def équivalente est constante
    }
  else
    {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 ) && (phi_2_dt > 0.)) 
	    // il peut arriver dans des itérations non correctes que phi_2_dt soit négatif 
	 	   { double unsurQdeltaSigma=1./QdeltaSigmaRatdt;
	 	     defEqui = def_i_equi + unsurwprime * phi_2_dt * unsurQdeltaSigma;
      }
	    else
	 	   {defEqui = def_i_equi ;};
    };
  // calcul de l'évolution de la def équivalente du point de référence à la situation actuelle
  double defEqui_Ratdt = defEqui-def_equi_R;
  
  // --- on s'occupe éventuellement de Qzero
  // la valeur de Qero est tout d'abord fonction du fait que l'on est sur la courbe de première charge ou non
  if (Qzero_defEqui != NULL)
   {if (wBase == 1)
     // cas où l'on est sur la courbe de première charge
	    {// on utilise la déformation équivalente courante
        double coefQzero = Qzero_defEqui->Valeur(defEqui_Ratdt);
        Qzero *= coefQzero;		  
	    }
    else
     // cas de secondes charges	: on garde pour Qzero le maxi que l'on a obtenue sur la première charge
	    { double   coefQzero = Qzero_defEqui->Valeur(defEqui_maxi);
	      Qzero *= coefQzero;
	    };
//cout <<"d Q " <<  defEqui_maxi << " " << Qzero << " " ; 
   }; 
  
  
  // calcul éventuellement du nouveau paramètre xnp
  if (xnp_defEqui != NULL)
   {// on définit un paramètre para qui correspond grosso-modo au paramètre: W * Qzero^2 / mu , d'une part 
    // de la thèse de pégon annexe II.2 formule A.2.1 ou thèse de laurent, chatpitre II.55 formule II-38
    double para =  defEqui_Ratdt * (2. * Qzero * Qzero / (wprime * xmu )) * QdeltaSigmaRatdt;
// --- débug ---  
//   double para =  defEqui_Ratdt ;  // ---------------- pour le débug
//  cout << " " << para << " ";
// --- fin débug ----
//para = sqrt(para);  // !!!!!!! pour voir

    // --- on s'occupe tout d'abord de xnp  
    double coefxnp = xnp_defEqui->Valeur(para);
//--debug
//cout << " val_coef_paradefEqui= " << val_coef_paradefEqui << " coefxnp= "<<coefxnp;
//-- fin débug
    xnp = (2.- xnp * val_coef_paradefEqui ) * coefxnp + xnp * val_coef_paradefEqui;
  
// --- débug ---  
//  cout << "  xnp= " << xnp << " ";
// --- fin débug ----
   };
  

};

// calcul des paramètres matériaux dépendants de la phase de la déformation totale
// en entrée: les paramètres sans phase = paramètres lues ou fonction de la température
// en sortie: les paramètres en tenant compte de la phase de epsBH
void Hysteresis3D::ParaMateriauxAvecPhaseEpsilon(const Tenseur3HH&  gijHH,const Tenseur3BB&  epsBB
                                                 ,double & xnp, double & xmu, double& Qzero )
	{// on calcul la déformation en mixte
	 Tenseur3BH epsBH(epsBB * gijHH);
	 // on regarde si le tenseur n'est pas trop petit
	 double Qepsilon =  sqrt(epsBH.II()); 
	 double Qepsilon3 = Qepsilon * Qepsilon * Qepsilon;
	 val_coef_paradefEqui = 1. ; // init par défaut
  if (mini_QepsilonBH > 0.) // cas historique
     {if (Qepsilon > mini_QepsilonBH )
       { // on peut calculer un cos3phi pas débile
         double bIIIb = epsBH.III() / 3.;
         double cos3phi = 3. * sqrt(6.)  * bIIIb/ Qepsilon3;
    //--debug
    //if (abs(cos3phi) > 1.) cout << " cos3phi= " << cos3phi;
    //-- fin débug
            // on limite les valeurs au normal
            if (cos3phi > 1.) { cos3phi = 1.;}
            else if (cos3phi < -1.) { cos3phi = -1.;};
         
         // on calcul les différentes grandeurs dépendantes de la phase
         if (avec_phase_paradefEqui && (coef_paradefEqui  != NULL))
           val_coef_paradefEqui = coef_paradefEqui->Valeur(cos3phi);
    //--debug
    //cout << " val_coef_paradefEqui= " << val_coef_paradefEqui;
    //-- fin débug
            if (avec_phaseEps)
          {if (Qzero_phase  != 0)
          { double coefQ  = Qzero_phase->Valeur(cos3phi);
            Qzero *= coefQ;
    //--debug
    //cout << " Qzero= " << Qzero;
    //-- fin débug
          };
           if (xmu_phase  != 0)
          { double coefxmu  = xmu_phase->Valeur(cos3phi);
            xmu *= coefxmu;
          };
           if (xnp_phase  != 0)
          { double coefxnp = xnp_phase->Valeur(cos3phi);
            xnp *= coefxnp;
          };
          };
    //--debug
    //cout << " Qzero= " << Qzero << " xmu= " << xmu << " xnp= " << xnp ;
    //-- fin débug
       };
     }
  else // cas on mini_QepsilonBH est un facteur de régularisation
     { // on peut calculer un cos3phi pas débile
       double bIIIb = epsBH.III() / 3.;
       double cos3phi = 3. * sqrt(6.)  * bIIIb/ (Qepsilon3-mini_QepsilonBH);
  //--debug
  //if (abs(cos3phi) > 1.) cout << " cos3phi= " << cos3phi;
  //-- fin débug
          // on limite les valeurs au normal
          if (cos3phi > 1.) { cos3phi = 1.;}
          else if (cos3phi < -1.) { cos3phi = -1.;};
       
       // on calcul les différentes grandeurs dépendantes de la phase
       if (avec_phase_paradefEqui && (coef_paradefEqui  != NULL))
           val_coef_paradefEqui = coef_paradefEqui->Valeur(cos3phi);
  //--debug
  //cout << " val_coef_paradefEqui= " << val_coef_paradefEqui;
  //-- fin débug
       if (avec_phaseEps)
        {if (Qzero_phase  != 0)
           { double coefQ  = Qzero_phase->Valeur(cos3phi);
             Qzero *= coefQ;
     //--debug
     //cout << " Qzero= " << Qzero;
     //-- fin débug
           };
         if (xmu_phase  != 0)
           { double coefxmu  = xmu_phase->Valeur(cos3phi);
             xmu *= coefxmu;
           };
         if (xnp_phase  != 0)
           { double coefxnp = xnp_phase->Valeur(cos3phi);
             xnp *= coefxnp;
           };
        };
  //--debug
  //cout << " Qzero= " << Qzero << " xmu= " << xmu << " xnp= " << xnp ;
  //-- fin débug
     };
  
	 // sinon on ne fait rien
	};

// calcul de la déformation cumulée, et de la déformation maxi à prendre en compte 
void  Hysteresis3D::DefEqui_et_maxi(double& QdeltaSigmaRatdt_,const bool& trajet_neutre_,double& defEquiMaxi_
                    ,const double& def_equi_R_,const double& def_i_equi_, double& def_equi_,const double & phi_2_dt_
						  ,const double& wprime_)
{ // calcul de la def équivalente
  def_equi_  = 0.; // init

  // la déformation équivalente calculée ici est approchée. 
  if (trajet_neutre_)
    // calcul de la déformation équivalente 
	 {def_equi_ = def_i_equi_ ;} // le long d'un trajet neutre la def équivalente est constante
  else
    {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 ) && (phi_2_dt_ > 0.)) 
	 // il peut arriver dans des itérations non correctes que phi_2_dt soit négatif 
	 	{ double unsurQdeltaSigma=1./QdeltaSigmaRatdt_;
	 	  def_equi_ = def_i_equi_ + unsurwprime * phi_2_dt_ * unsurQdeltaSigma;
	 	}
	 else
	 	{def_equi_ = def_i_equi_ ;};
    };
};

// Transport des différentes grandeurs en fonction de la variable : type_de_transport_memorisation
// coordonnées initiales du tenseur en mixte: aBH,
// coordonnées finales: rBH : donc correspondant au tenseur transporté en l'état actuel
void Hysteresis3D::Transport_grandeur(const Tenseur3BB&  gijBB, const Tenseur3HH&  gijHH
                        ,const Tenseur3BH& aBH, Tenseur3BH & rBH)
  { switch (type_de_transport_memorisation)
     { case 0 : // cas historique on ne fait rien
        rBH = aBH;
        break;
       case -1 : // cas d'un transport cohérent avec Jauman
        {// on crée un tenseur transposé
         Tenseur3HB  cHB = aBH.Transpose();
         // le transporté
         rBH = gijBB * cHB * gijHH;
         rBH = 0.5*(aBH+rBH);
         break;
        }
       default: break;
     };
  };

// transport éventuel de toutes les grandeurs tensorielles mémorisées dans save_resul
void Hysteresis3D::Transport_grandeur(const Tenseur3BB&  gijBB, const Tenseur3HH&  gijHH,SaveResulHysteresis3D & save_resul)
  { switch (type_de_transport_memorisation)
     { case 0 : // cas historique on ne fait rien
        break;
       case -1 : // cas d'un transport cohérent avec Jauman
        {// on crée un tenseur transposé
         { // sigma_barre_BH_R
           List_io <Tenseur3BH>::iterator itens,itens_fin = save_resul.sigma_barre_BH_R.end();
           Tenseur3BH interBH; // tenseur intermédiaire de travail
           for (itens=save_resul.sigma_barre_BH_R.begin();itens != itens_fin;itens++)
            { interBH = *itens;
              Transport_grandeur(gijBB,gijHH,interBH,*itens);
            };
           // oc_BH_R
           List_io <Tenseur3BH>::iterator ioc_i,ioc_i_fin = save_resul.oc_BH_R.end();
           for (ioc_i=save_resul.oc_BH_R.begin();ioc_i != ioc_i_fin;ioc_i++)
            { interBH = *ioc_i;
              Transport_grandeur(gijBB,gijHH,interBH,*ioc_i);
            };
         };
         break;
        }
       default: break;
     };
  };

  // gestion du paramètre "modif" de saveresult
  // inversion = true: on met à jour après une inversion
  //           = false : on met à jour après une coïncidence
  void Hysteresis3D::Gestion_para_Saveresult_Modif
                 (const bool& pt_sur_principal,SaveResulHysteresis3D & save_resul
                  ,const bool& inversion
                 )
  {
    // --- gestion de l'indicateur modif
    switch (save_resul.modif)
     { case 0: // le cas le plus simple: on avait rien avant
        {if (inversion) save_resul.modif = 2; // cas d'une inversion
         else           save_resul.modif = 1; // cas d'une coïncidence
         break;
        }
       case 1: // on a déjà une ou plusieurs coïncidence(s)
        {if (inversion) save_resul.modif = 3; // cas d'une inversion
         else           save_resul.modif = 1; // cas d'une coïncidence
         break;
        }
       case 2: // on avait une ou plusieurs inversions
        {if (inversion) save_resul.modif = 2; // cas d'une inversion sup
         else // cas d'une coïncidence
          {if (pt_sur_principal) save_resul.modif = 0; // on est revenu sur la courbe principal
           else save_resul.modif = 2; // il reste encore des pt secondaires
          };
         break;
        }
       case 3: // on avait un mixte d'inversion et de coïncidence
        {if (inversion) save_resul.modif = 3; // cas d'une inversion sup
         else // cas d'une coïncidence
          {if (pt_sur_principal) save_resul.modif = 1; // on est revenu en coïncidence résiduelle
           else save_resul.modif = 3; // il reste encore des pt secondaires
          };
         break;
        }
       default:
         cout << "\n Erreur *** Hysteresis3D::Gestion_para_Saveresult_Modif: cas non prevu !!"
              << "  erreur save_resul.modif= " << save_resul.modif << endl;
         Sortie(1);
       break;
     };
  };