// 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 "Deformation.h"
#include "ConstMath.h"
#include "ParaGlob.h"
#include "MathUtil.h"
#include "Tenseur3.h"
#include "Tenseur2.h"
#include "Tenseur1.h"
# include "TypeConsTens.h"
#include <math.h>

// -----------------------------------------------------------------
//      cas de la déformation Logarithmique
// -----------------------------------------------------------------

  	  // cas implicite : tous les parametres sont de resultats	
const Met_abstraite::Impli& Deformation::Cal_implicit_Logarithmique (bool 
           ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
           ,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul
           ,const Met_abstraite::Impli& ex)
{
	 const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
	 double deltat=vartemps.IncreTempsCourant();
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//	 // dans le cas où le pas de temps est trop petit on met un message d'erreur
//	 if (Dabs(deltat) < ConstMath::trespetit)
//	 	{ cout << "\n erreur le pas de temps " << deltat 
//	 	       << " est trop faible pour un calcul correcte de la vitesse moyenne";
//	 	  if (ParaGlob::NiveauImpression() > 3)      
//	 	       cout << "\n Deformation::Cal_implicit_Logarithmique( .... ";
//	 	  Sortie(1);     
//	 	};
  // dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
  if (!premier_calcul)
   { Deformation::Stmet& a_0 = saveDefResul->meti_00;  // pour simplifier
     Deformation::Stmet& a_t = saveDefResul->meti_t;  // pour simplifier     
     metrique->Recup_grandeur_0_t(*a_0.giB_,*a_0.giH_,*a_t.giB_,*a_t.giH_
                                  ,*a_0.gijBB_,*a_0.gijHH_,*a_t.gijBB_,*a_t.gijHH_
                                  ,*a_t.gradVmoyBB_,*a_t.jacobien_,*a_0.jacobien_);
   };
  bool variation = true; // pour la variation de B et Palpha_BH
  // appel pour le calcul de la déformation logarithmique et de sa variation
  Cal_Logarithmique (*ex.gijBB_0,*ex.gijHH_0,d_epsBB_tdt,*ex.gijBB_tdt,*ex.gijHH_tdt,epsBB_tdt
                     ,*ex.d_gijBB_tdt,variation);
  // -- calcul de l'acroissement de déformation et de la vitesse de déformation
  delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
  // -- calcul de la vitesse de déformation: a priori on utilise 
  // la variation du tenseur déformation sur delta t ce qui correspond à l'utilisation du gradient
  // de vitesse de déformation moyen calculé à t+(delta t)/2
	 if (Dabs(deltat) >= ConstMath::trespetit)
	    { DepsBB = delta_epsBB/deltat; }
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//	 else // dans le cas où l'incrément de temps est nul pour l'instant on met la vitesse de déformation à 0
//	    { DepsBB.Inita(0.);}
     else
        // on utilise la vitesse précédente pour la vitesse de déformation
        { DepsBB = (*saveDefResul->D_BB_t);};
 
//	    //on utilise le gradient de vitesse moyen -> pas bon car le gradient de vitesse moyen sera également null donc calcul inutile
//	      // pour éviter une division par 0
//	   { DepsBB = 0.5 * ((*ex.gradVmoyBB_tdt) + ex.gradVmoyBB_tdt->Transpose()); }

     // calcul éventuelle de la dérivée de D* c'est-à-dire la dérivée seconde de la déformation d'almansi !!
     // et qui sert pour la dérivée complète de la puissance interne (c'est différent de la dérivée seconde de la def log)
     // qui elle ne sert ensuite à rien !!
//     if (cal_derSeconde || premier_calcul)  
//	  { // recup des variations secondes de la déformation
//	   int d2_epsBB_tdtTaille1 = d2_epsBB_tdt.Taille1();
//       for (int i=1; i<= d2_epsBB_tdtTaille1; i++)
//        for (int j=1; j<= i; j++)  // symétrie du  tableau et du tenseur
//	      { *(d2_epsBB_tdt(i,j)) = 0.5 * (*((*(ex.d2_gijBB_tdt))(i,j))) ;
//	        *(d2_epsBB_tdt(j,i)) = *(d2_epsBB_tdt(i,j)) ; 
//	       }
//	   }    


     // vérification éventuelle
  //   VerifCal_deflog(gradV_instantane,ex,epsBB_tdt,d_epsBB_tdt);
     return ex;
};
  
  	  // cas explicite : tous les parametres sont de resultats	
const Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (bool 
           ,TenseurBB & epsBB_t,Tableau <TenseurBB *> & d_epsBB,TenseurBB& DepsBB
           ,TenseurBB& delta_epsBB,bool premier_calcul,const Met_abstraite::Expli& ex)
  {  
     // récup des infos sur le temps 
	 const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
	 double deltat=vartemps.IncreTempsCourant();
	 double temps =vartemps.TempsCourant();
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//	 // dans le cas où le pas de temps est trop petit on met un message d'erreur
//	 if (Dabs(deltat) < ConstMath::trespetit)
//	 	{ cout << "\n erreur le pas de temps " << deltat 
//	 	       << " est trop faible pour un calcul correcte de la vitesse moyenne";
//	 	  if (ParaGlob::NiveauImpression() > 3)      
//	 	      cout  << "\n Deformation::Cal_explicit_Almansi( .... ";
//	 	  Sortie(1);     
//	 	};

     // dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
     if (!premier_calcul)
      { Deformation::Stmet& a_0 = saveDefResul->meti_00;  // pour simplifier
        metrique->Recup_grandeur_0(*a_0.giB_,*a_0.giH_,*a_0.gijBB_,*a_0.gijHH_,*a_0.jacobien_);
      };
     bool variation = false; // pour la variation de B et Palpha_BH
     // appel pour le calcul de la déformation logarithmique et de sa variation
     Cal_Logarithmique (*ex.gijBB_0,*ex.gijHH_0,d_epsBB,*ex.gijBB_t,*ex.gijHH_t,epsBB_t
                        ,*ex.d_gijBB_t,variation);
	 
     // -- calcul de l'acroissement de déformation et de la vitesse de déformation
     // ici on considère que delta_epsBB est petit donc de l'ordre de delta_t * DijBB
     // ******* a reconsidérer par la suite ***********
     if (Dabs(temps) > ConstMath::trespetit)
     // on peut utiliser le temps en dénominateur
       {
            // normalement on considère que le delta_epsBB est proportionnel à epsBB_t ! c-a-d le delta de 0 à t
         if (Dabs(deltat) >= ConstMath::trespetit)
              { delta_epsBB = (deltat/temps) * epsBB_t;
                DepsBB = delta_epsBB/deltat;}
         else// dans le cas où l'incrément de temps est nul on garde l'incrément total pour deltaeps
           { delta_epsBB =  epsBB_t;
             // pour la vitesse de déformation, on utilise une approximation linéaire sur le temps complet
             if (Dabs(temps) >= ConstMath::trespetit)
              { DepsBB = delta_epsBB/temps;}
             else // dans le cas où l'incrément de temps est nul on utilise le gradient de vitesse moyen
           // pour éviter une division par 0
       // non ça ne marche pas car le calcul du gradient moyen utilise lui aussi le delta t, mais la question qu'il faut
       // se poser c'est pourquoi delta t est si petit plutôt que d'essayer d'y palier
       // normalement passe jamais ici vu le test plus haut !!	      
 // modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//                { DepsBB = 0.5 * ((*ex.gradVmoyBB_t) + ex.gradVmoyBB_t->Transpose()); }
              {DepsBB = (*saveDefResul->D_BB_t);};
            };
       }
     else // sinon le temps est nul
       { delta_epsBB.Inita(0.); // on met à 0 l'acroissement de déformation de t à tdt
         // on utilise la vitesse précédente pour la vitesse de déformation
         DepsBB = (*saveDefResul->D_BB_t);
       };
  
	 // sauvegarde des infos à 0 éventuellement
     if (premier_calcul) saveDefResul->MiseAJourGrandeurs_a_0(metrique);
     // sauvegarde des infos à t à chaque passage
     saveDefResul->MiseAJourGrandeurs_a_tdt(metrique,DepsBB);

     return ex;
  };
  
  	  // cas explicite à t et tdt: tous les parametres sont de resultats	
const Met_abstraite::Expli_t_tdt& Deformation::Cal_explicit_logarithmique_tdt (bool 
           ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB,TenseurBB& DepsBB
           ,TenseurBB& delta_epsBB,bool premier_calcul,const Met_abstraite::Expli_t_tdt& ex)
  {  
	 const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
	 double deltat=vartemps.IncreTempsCourant();
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//	 // dans le cas où le pas de temps est trop petit on met un message d'erreur
//	 if (Dabs(deltat) < ConstMath::trespetit)
//	 	{ cout << "\n erreur le pas de temps " << deltat 
//	 	       << " est trop faible pour un calcul correcte de la vitesse moyenne";
//	 	  if (ParaGlob::NiveauImpression() > 3)      
//	 	       cout << "\n Deformation::Cal_explicit_logarithmique_tdt( .... ";
//	 	  Sortie(1);     
//	 	};
     // dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
     if (!premier_calcul)
      { Deformation::Stmet& a_0 = saveDefResul->meti_00;  // pour simplifier
        Deformation::Stmet& a_t = saveDefResul->meti_t;  // pour simplifier     
        metrique->Recup_grandeur_0_t(*a_0.giB_,*a_0.giH_,*a_t.giB_,*a_t.giH_
                                     ,*a_0.gijBB_,*a_0.gijHH_,*a_t.gijBB_,*a_t.gijHH_
                                     ,*a_t.gradVmoyBB_,*a_t.jacobien_,*a_0.jacobien_);
      };
     bool variation = false; // pour la variation de B et Palpha_BH
     // appel pour le calcul de la déformation logarithmique et de sa variation
     Cal_Logarithmique (*ex.gijBB_0,*ex.gijHH_0,d_epsBB,*ex.gijBB_tdt,*ex.gijHH_tdt,epsBB_tdt
                        ,*ex.d_gijBB_tdt,variation);
     // -- calcul de l'acroissement de déformation et de la vitesse de déformation
     delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
     // -- calcul de la vitesse de déformation: a priori on utilise 
     // la variation du tenseur déformation sur delta t ce qui correspond à l'utilisation du gradient
     // de vitesse de déformation moyen calculé à t+(delta t)/2
	 if (Dabs(deltat) >= ConstMath::trespetit)
	    { DepsBB = delta_epsBB/deltat; }
// modif du 3 nov 2014 -> introduction de la sauvegarde de D : D_BB_t et D_BB_tdt
//	 else // dans le cas où l'incrément de temps est nul on utilise le gradient de vitesse moyen
//	      // pour éviter une division par 0
//	   { DepsBB = 0.5 * ((*ex.gradVmoyBB_tdt) + ex.gradVmoyBB_tdt->Transpose()); }
     else
        // on utilise la vitesse précédente pour la vitesse de déformation
        { DepsBB = (*saveDefResul->D_BB_t);};

     return ex;
  };

    // fonction qui calcule le tenseur logarithmique et sa variation éventuelle
void Deformation::Cal_Logarithmique 
                     (const TenseurBB & gijBB_0,const TenseurHH& _gijHH_0,Tableau <TenseurBB *> & d_epsBB
                      ,const TenseurBB& _gijBB,const TenseurHH& _gijHH,TenseurBB & _epsBB
                      ,const Tableau <TenseurBB *>& _d_gijBB,bool variation)
{// dimensions                                                     
 int dima = gijBB_0.Dimension(); // dimension des tenseurs
 int  dimddl = d_epsBB.Taille(); // nombre de ddl                                         
 // def des grandeurs de travail
 DimensionementVarLog(dima,variation,dimddl);                  
 // choix suivant la dimension des tenseurs de manière à optimiser la vitesse et les calculs
 switch (dima)
 {case 1:    // $$$$$ dimension 1 avec les tenseurs $$$$$$
		{// simplification et transformation en tenseur de dim 3 (optimisation des calculs)
   const Tenseur1HH & gijHH_0 = *((const Tenseur1HH*) & _gijHH_0);
   const Tenseur1BB & gijBB = *((const Tenseur1BB *) & _gijBB);
   Tenseur1BB & epsBB = *((Tenseur1BB *) & _epsBB);
   const Tenseur1HH & gijHH = *((const Tenseur1HH *) & _gijHH);
   const Tableau <Tenseur1BB *>& d_gijBB =*((Tableau <Tenseur1BB *>*) & _d_gijBB);
   // def du tenseur B_BH (cauchy_green gauche)
   (*B_BH_tr) = gijBB * gijHH_0; Tenseur1BH&  B_BH= *((Tenseur1BH*) B_BH_tr);
   // calcul des valeurs propres et projection propres du tenseur B_BH 
   Tableau <Tenseur1BH * >& d_B_BH=*((Tableau <Tenseur1BH * >*) d_B_BH_tr);
   if (variation) // dans le cas de variation, calcul des variation de B_BH
   // en fait ici il s'agit des variations des composantes et non de tout le tenseur
         { for (int ic=1;ic<=dimddl;ic++) {*(d_B_BH(ic)) = (*(d_gijBB(ic))) * gijHH_0;};};
   // def des grandeurs de travail
   Tableau <Tenseur1BH* >& Palpha_BH=*((Tableau <Tenseur1BH* >*) Palpha_BH_tr); // les projections propres
   Tableau <Tableau <Tenseur1BH* > >& d_Palpha_BH
                 =*((Tableau <Tableau <Tenseur1BH* > >*) d_Palpha_BH_tr); // variation des projections propres
   // appel de la méthode calculant les grandeurs de travail
   int cas_ki; // indique le cas des valeurs propres traitées (cf TenseurBH)
   Val_et_projection_prop_tenseur(*B_BH_tr,*d_B_BH_tr,ki,*Palpha_BH_tr,variation,*d_ki,*d_Palpha_BH_tr,cas_ki);
   // calcul des composantes du tenseur logarithmique dans la base naturelle
   switch (cas_ki)
    { case 1: // cas d'un calcul correcte de la valeur propre
        { // vérif du signe de la valeur propre
          if (ki(1) < 0.) 
           { cout << "\n erreur (8) on trouve une valeur propre negative pour le tenseur B (Cauchy Green droit)"
                  << "\n on ne peut pas calculer le tenseur log correspondant : val= " << ki;
             if (ParaGlob::NiveauImpression() >= 4) 
                 cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
             Sortie(1);
            };
          // calcul du tenseur log
          double lambda1=sqrt((double) (ki(1)));
          double log_lambda1 = log(lambda1); 
          Tenseur1BH& epsBH = *( (Tenseur1BH*)  NevezTenseurBH(*(Palpha_BH(1)))); // un tenseur de travail
          epsBH = log_lambda1 * (*(Palpha_BH(1)));
          epsBB = epsBH * gijBB;
          // calcul de la variation des composantes du tenseur log
          if (variation)
           {for (int i=1; i<= dimddl; i++)
            { Tenseur1BB & depsBB = *((Tenseur1BB *) (d_epsBB(i))); //    "
              // variation des composantes du log dans sa base
              double d_log_lambda1 = 0.5 * (*d_ki)(i)(1)/ki(1);
              // calcul de la variation du tenseur log, en fait des composantes et non du tenseur totale
              depsBB = (d_log_lambda1 * (*(Palpha_BH(1))) 
                     + log_lambda1 * (*(d_Palpha_BH(i)(1)))) * gijBB + epsBH * (*d_gijBB(i));
             };
            };                
          break;
        }                    
      default:
        cout << "\n erreur (9) il y a une erreur dans le calcul des valeurs propres pour le tenseur "
             << " B (Cauchy Green droit), B= ";   
        B_BH_tr->Ecriture(cout);                        
        if (ParaGlob::NiveauImpression() >= 4) 
            cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
        Sortie(1);
    };// -- fin du switch sur les différents cas de calcul de valeurs propres      
   // -- fin du cas de la dimension 1 pour les tenseurs 
   break;       
  }
  case 2:    // $$$$$ dimension 2 avec les tenseurs $$$$$$
		{// simplification et transformation en tenseur de dim 2 (optimisation des calculs)
         const Tenseur2HH & gijHH_0 = *((const Tenseur2HH *) & _gijHH_0);
         const Tenseur2BB & gijBB = *((const Tenseur2BB *) & _gijBB);
         Tenseur2BB & epsBB = *(( Tenseur2BB *) & _epsBB);
         const Tenseur2HH & gijHH = *((const Tenseur2HH *) & _gijHH);
         const Tableau <Tenseur2BB *>& d_gijBB =*((Tableau <Tenseur2BB *>*)& _d_gijBB);
         // def du tenseur B_BH (cauchy_green gauche)
         (*B_BH_tr) = gijBB * gijHH_0; Tenseur2BH&  B_BH= *((Tenseur2BH*) B_BH_tr);
         // calcul des valeurs propres et projection propres du tenseur B_BH 
         Tableau <Tenseur2BH * >& d_B_BH=*((Tableau <Tenseur2BH * >*) d_B_BH_tr);
         if (variation) // dans le cas de variation, calcul des variation de B_BH
               { for (int ic=1;ic<=dimddl;ic++) {*(d_B_BH(ic)) = (*(d_gijBB(ic))) * gijHH_0;};};
         // def des grandeurs de travail
         Tableau <Tenseur2BH* >& Palpha_BH=*((Tableau <Tenseur2BH* >*) Palpha_BH_tr); // les projections propres
         Tableau <Tableau <Tenseur2BH* > >& d_Palpha_BH
                       =*((Tableau <Tableau <Tenseur2BH* > >*) d_Palpha_BH_tr); // variation des projections propres
         // appel de la méthode calculant les grandeurs de travail
         int cas_ki; // indique le cas des valeurs propres traitées (cf TenseurBH)
         Val_et_projection_prop_tenseur(*B_BH_tr,*d_B_BH_tr,ki,*Palpha_BH_tr,variation,*d_ki,*d_Palpha_BH_tr,cas_ki);
         // calcul des composantes du tenseur logarithmique dans la base naturelle
         switch (cas_ki)
          { case 1: // cas de deux valeurs propres distinctes
            { // vérif du signe des valeurs propres
              if ((ki(1) < 0.) || (ki(2) < 0.))
               { cout << "\n erreur (5) on trouve une valeur propre negative pour le tenseur B (Cauchy Green droit)"
                      << "\n on ne peut pas calculer le tenseur log correspondant : val= " << ki;
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
                }
              // calcul du tenseur log
              double lambda1=sqrt((double)(ki(1)));double lambda2=sqrt(((double)ki(2)));
              double log_lambda1 = log(lambda1);double log_lambda2 = log(lambda2); 
              Tenseur2BH& epsBH = *( (Tenseur2BH*)  NevezTenseurBH(*(Palpha_BH(1)))); // un tenseur de travail
              epsBH = (log_lambda1 * (*(Palpha_BH(1))) + log_lambda2 * (*(Palpha_BH(2))));
              epsBB = epsBH * gijBB;
              // calcul de la variation des composantes du tenseur log
//              Tenseur2BB totoBB;Tenseur2BH totoBH;Tenseur_ns2BB titiBB;
              if (variation)
               {for (int i=1; i<= dimddl; i++)
                { Tenseur2BB & depsBB = *((Tenseur2BB *) (d_epsBB(i))); //    "
                  // variation des composantes du log dans sa base
                  // *** bizarre la formule fausse suivante donne de meilleures résultats que la formule juste ****
                  double d_log_lambda1 = 0.5 * (*d_ki)(i)(1)/lambda1;//ki(1);
                  double d_log_lambda2 = 0.5 * (*d_ki)(i)(2)/lambda2; //ki(2);
                  // calcul de la variation du tenseur log, en fait des composantes et non du tenseur totale
  /*                totoBH = ((*d_ki)(i)(1) * (*(Palpha_BH(1))) + ki(1) * (*(d_Palpha_BH(i)(1)))
                           + (*d_ki)(i)(2) * (*(Palpha_BH(2))) + ki(2) * (*(d_Palpha_BH(i)(2))));
                  titiBB = totoBH * gijBB + (*B_BH_tr) * (*d_gijBB(i));
                  
                  totoBB = ((*d_ki)(i)(1) * (*(Palpha_BH(1))) + ki(1) * (*(d_Palpha_BH(i)(1)))
                           + (*d_ki)(i)(2) * (*(Palpha_BH(2))) + ki(2) * (*(d_Palpha_BH(i)(2)))) * gijBB 
                           + (*B_BH_tr) * (*d_gijBB(i));*/
                  
                  depsBB = (d_log_lambda1 * (*(Palpha_BH(1))) + log_lambda1 * (*(d_Palpha_BH(i)(1)))
                          +d_log_lambda2 * (*(Palpha_BH(2))) + log_lambda2 * (*(d_Palpha_BH(i)(2)))) * gijBB 
                          + epsBH * (*d_gijBB(i));
            // **** on ajoute la matrice tangente d'Almansi, cela marche mieux mais on ne sait pas pourquoi ****
                  depsBB +=  0.5* (*d_gijBB(i));       
                  //depsBB(1,2)+=0.5* (*d_gijBB(i))(1,2);depsBB(2,1)+=0.5* (*d_gijBB(i))(2,1);  
                 }
                };                
              break;
              }                    
            case 0: // cas de deux valeurs propres identiques
            { // vérif du signe de la valeur propre
              if (ki(1) < 0.)
               { cout << "\n erreur (6) on trouve une valeur propre negative pour le tenseur B (Cauchy Green droit)"
                      << "\n on ne peut pas calculer le tenseur log correspondant : val= " << ki(1);
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
                }
              // calcul du tenseur log
              double lambda1=sqrt((double) (ki(1)));double log_lambda1 = log(lambda1);
              Tenseur2BH& epsBH = *( (Tenseur2BH*)  NevezTenseurBH(*(Palpha_BH(1)))); // un tenseur de travail
              epsBH = log_lambda1 * (*(Palpha_BH(1)));
              epsBB = epsBH * gijBB;
              // calcul de la variation des composantes du tenseur log
//              Tenseur2BB totoBB;Tenseur2BH B_BH_nBH;Tenseur_ns2BB titiBB;
              Tenseur2BH depsBH;
//              Vecteur ki_n(2);Vecteur dki_n(2);ki(2)=ki(1);
//              Tableau <Tenseur2BH > Palpha_s_BH(2);
//              Palpha_s_BH(1) = *(Palpha_BH(1));Palpha_s_BH(2) = *(Palpha_BH(2));
              if (variation)
               {for (int i=1; i<= dimddl; i++)
                { Tenseur2BB & depsBB = *((Tenseur2BB *) (d_epsBB(i))); //    "
                  // dans le cas oû le tenseur est sphérique on utilise d'autres relations
                  // B_1^1=ki(1) et B_2^2=ki(2), 
                  double d_log_lambda1 = 0.5 * (*d_B_BH(i))(1,1)/ki(1);
                  double d_log_lambda2 = 0.5 * (*d_B_BH(i))(2,2)/ki(1);
                  depsBH.Inita(0.);
                  depsBH.Coor(1,1)=d_log_lambda1; depsBH.Coor(2,2)=d_log_lambda2;
                  // pour les deux autres composantes on utilise les éléments du tenseurs d'almansi
                  // ce qui évite d'avoir une raideur nulle en rotation
                  depsBB = depsBH *  gijBB + epsBH * (*d_gijBB(i));
                  depsBB.Coor(1,2)+=0.5* (*d_gijBB(i))(1,2);//depsBB(2,1)+=0.5* (*d_gijBB(i))(2,1);
             // **** on ajoute la matrice tangente d'Almansi, cela marche mieux mais on ne sait pas pourquoi ****
              //    depsBB +=  0.5* (*d_gijBB(i));       
                  depsBB.Coor(1,2)+=0.5* (*d_gijBB(i))(1,2);//depsBB(2,1)+=0.5* (*d_gijBB(i))(2,1);
                 
  
 /* fonctionne mais ne me plait pas car pb de perturbation 
                  // on cherche à lever l'indétermination par différence finie
                  double delta = ConstMath::unpeupetit;
                  B_BH_nBH = B_BH + (*(d_B_BH(i))) * delta;
                  // on ne calcul que les valeurs propres et projection propre                  
                  Val_et_projection_prop_tenseur
                     (B_BH_nBH,*d_B_BH_tr,ki,*Palpha_BH_tr,false,*d_ki,*d_Palpha_BH_tr,cas_ki);
                  if (cas_ki == 0)
                   { dki_n= ( ki_n - ki)/ delta;
                     double d_log_lambda1 = 0.5 * (dki_n(1))/ki(1);
                     depsBB = (d_log_lambda1 * (*(Palpha_BH(1)))
                            + (log_lambda1 / delta) * ((*(Palpha_BH(1))) - Palpha_s_BH(1))
                             ) * gijBB 
                            + epsBH * (*d_gijBB(i));
                    }
                  else
                   { dki_n= ( ki_n - ki)/ delta;
                     double d_log_lambda1 = 0.5 * (dki_n(1))/ki(1);
                     double d_log_lambda2 = 0.5 * (dki_n(2))/ki(1);
                     depsBB = (d_log_lambda1 * (*(Palpha_BH(1)))
                            + (log_lambda1 / delta) * ((*(Palpha_BH(1))) - Palpha_s_BH(1))
                          +d_log_lambda2 * (*(Palpha_BH(2))) 
                            + (log_lambda1 / delta) * ((*(Palpha_BH(2))) - Palpha_s_BH(2))
                             ) * gijBB 
                          + epsBH * (*d_gijBB(i));
                    };      
     */
                 }
                };                
              break;
              }                    
            default:
              cout << "\n erreur (7) il y a une erreur dans le calcul des valeurs propres pour le tenseur "
                   << " B (Cauchy Green droit), B= ";   
              B_BH_tr->Ecriture(cout);                        
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
           };// -- fin du switch sur les différents cas de calcul de valeurs propres      
         // -- fin du cas de la dimension 2 pour les tenseurs 
         break;       
         }
  case 3:   // $$$$$$ dimension 3 pour les tenseurs $$$$$$$
		{// simplification et transformation en tenseur de dim 3 (optimisation des calculs)
         const Tenseur3HH & gijHH_0 = *((const Tenseur3HH *) & _gijHH_0);
         const Tenseur3BB & gijBB = *((const Tenseur3BB *) & _gijBB);
         Tenseur3BB & epsBB = *(( Tenseur3BB *) & _epsBB);
         const Tenseur3HH & gijHH = *((const Tenseur3HH *) & _gijHH);
         const Tableau <Tenseur3BB *>& d_gijBB =*((Tableau <Tenseur3BB *>*) & _d_gijBB);
         // def du tenseur B_BH (cauchy_green gauche)
         (*B_BH_tr) = gijBB * gijHH_0; Tenseur3BH&  B_BH= *((Tenseur3BH*) B_BH_tr);
         // calcul des valeurs propres et projection propres du tenseur B_BH 
         Tableau <Tenseur3BH * >& d_B_BH=*((Tableau <Tenseur3BH * >*) d_B_BH_tr);
         if (variation) // dans le cas de variation, calcul des variation de B_BH
               { for (int ic=1;ic<=dimddl;ic++) {*(d_B_BH(ic)) = (*(d_gijBB(ic))) * gijHH_0;};};
         // def des grandeurs de travail
         Tableau <Tenseur3BH* >& Palpha_BH=*((Tableau <Tenseur3BH* >*) Palpha_BH_tr); // les projections propres
         Tableau <Tableau <Tenseur3BH* > >& d_Palpha_BH
                       =*((Tableau <Tableau <Tenseur3BH* > >*) d_Palpha_BH_tr); // variation des projections propres
         // appel de la méthode calculant les grandeurs de travail
         int cas_ki; // indique le cas des valeurs propres traitées (cf TenseurBH)
         Val_et_projection_prop_tenseur(*B_BH_tr,*d_B_BH_tr,ki,*Palpha_BH_tr,variation,*d_ki,*d_Palpha_BH_tr,cas_ki);
         // calcul des composantes du tenseur logarithmique dans la base naturelle
         switch (cas_ki)
          { case 1: // cas de trois valeurs propres distinctes
            { // vérif du signe des valeurs propres
              if ((ki(1) < 0.) || (ki(2) < 0.)  || (ki(3) <= 0.))
               { cout << "\n erreur (1) on trouve une valeur propre negative pour le tenseur B (Cauchy Green droit)"
                      << "\n on ne peut pas calculer le tenseur log correspondant : val= " << ki;
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
                }
              // calcul du tenseur log
              double lambda1=sqrt(ki(1));double log_lambda1 = log(lambda1);
              double lambda2=sqrt(ki(2));double log_lambda2 = log(lambda2);
              double lambda3=sqrt(ki(3));double log_lambda3 = log(lambda3);
              Tenseur3BH& epsBH = *( (Tenseur3BH*)  NevezTenseurBH(*(Palpha_BH(1)))); // un tenseur de travail
              epsBH = (log_lambda1 * (*(Palpha_BH(1))) + log_lambda2 * (*(Palpha_BH(2))) 
                       + log_lambda3 * (*(Palpha_BH(3))));
              epsBB = epsBH * gijBB;
              // calcul de la variation des composantes du tenseur log
              if (variation)
               {for (int i=1; i<= dimddl; i++)
                { Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); //    "
                  // variation des composantes du log dans sa base
                  double d_log_lambda1 = 0.5 * (*d_ki)(i)(1)/ki(1);
                  double d_log_lambda2 = 0.5 * (*d_ki)(i)(2)/ki(2);
                  double d_log_lambda3 = 0.5 * (*d_ki)(i)(3)/ki(3);
                  // calcul de la variation du tenseur log, en fait des composantes et non du tenseur totale
                  depsBB = (d_log_lambda1 * (*(Palpha_BH(1))) + log_lambda1 * (*(d_Palpha_BH(i)(1)))
                                  +d_log_lambda2 * (*(Palpha_BH(2))) + log_lambda2 * (*(d_Palpha_BH(i)(2)))
                                  +d_log_lambda3 * (*(Palpha_BH(3))) + log_lambda3 * (*(d_Palpha_BH(i)(3)))) * gijBB
                            + epsBH * (*d_gijBB(i));
            // **** on ajoute la matrice tangente d'Almansi, cela marche mieux mais on ne sait pas pourquoi ****
         //         depsBB +=  0.5* (*d_gijBB(i));       
                 }
                };                
              break;
              }                    
            case 0: // cas de trois valeurs propres identiques
            { // vérif du signe de la valeur propre
              if (ki(1) < 0.)
               { cout << "\n erreur (2) on trouve une valeur propre negative pour le tenseur B (Cauchy Green droit)"
                      << "\n on ne peut pas calculer le tenseur log correspondant : val= " << ki(1);
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
                }
              // calcul du tenseur log
              double lambda1=sqrt((double) (ki(1))); 
              Tenseur3BH& epsBH = *( (Tenseur3BH*)  NevezTenseurBH(*(Palpha_BH(1)))); // un tenseur de travail
              epsBH = (log(lambda1) * (*(Palpha_BH(1))));
              epsBB = epsBH * gijBB;
              // calcul de la variation des composantes du tenseur log
              if (variation)
               {Tenseur3BH depsBH; // tenseur de travail
                for (int i=1; i<= dimddl; i++)
                { Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); //    "
                  // dans le cas oû le tenseur est sphérique on utilise d'autres relations
                  // B_1^1=ki(1) et B_2^2=ki(2), B_3^3=ki(3) 
                  double d_log_lambda1 = 0.5 * (*d_B_BH(i))(1,1)/ki(1);
                  double d_log_lambda2 = 0.5 * (*d_B_BH(i))(2,2)/ki(1);
                  double d_log_lambda3 = 0.5 * (*d_B_BH(i))(3,3)/ki(1);
                  depsBH.Inita(0.);
                  depsBH.Coor(1,1)=d_log_lambda1; depsBH.Coor(2,2)=d_log_lambda2;depsBH.Coor(3,3)=d_log_lambda3;
                  // pour les 3 autres composantes on utilise les éléments du tenseurs d'almansi
                  // ce qui évite d'avoir une raideur nulle en rotation
                  depsBB = depsBH *  gijBB + epsBH * (*d_gijBB(i));
                  depsBB.Coor(1,2)+=0.5* (*d_gijBB(i))(1,2);depsBB.Coor(1,3)+=0.5* (*d_gijBB(i))(1,3);
                  depsBB.Coor(2,3)+=0.5* (*d_gijBB(i))(2,3);
                 }
                };                
              break;
              }                    
            case 2: case 3: // cas de deux valeurs propres identiques
            { // vérif du signe des  valeurs propres distinctes
              if ((ki(1) < 0.) || (ki(2) < 0.))
               { cout << "\n erreur (3) on trouve une valeur propre negative pour le tenseur B (Cauchy Green droit)"
                      << "\n on ne peut pas calculer le tenseur log correspondant : val= " << ki(1) << " " << ki(2);
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
                }
              // calcul du tenseur log
              double lambda1=sqrt(ki(1));double log_lambda1 = log(lambda1);
              double lambda2=sqrt(ki(2));double log_lambda2 = log(lambda2);
              Tenseur3BH& epsBH = *( (Tenseur3BH*)  NevezTenseurBH(*(Palpha_BH(1)))); // un tenseur de travail
              epsBH = (log_lambda1 * (*(Palpha_BH(1))) + log_lambda2 * (*(Palpha_BH(2)))); 
              epsBB = epsBH * gijBB;
              // calcul de la variation des composantes du tenseur log
              Tenseur3BB totoBB;Tenseur_ns3BB interBB;
              if (variation)
               {for (int i=1; i<= dimddl; i++)
                { Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); //    "
                  // variation des composantes du log dans sa base
                  double d_log_lambda1 = 0.5 * (*d_ki)(i)(1)/ki(1);
                  double d_log_lambda2 = 0.5 * (*d_ki)(i)(2)/ki(2);
                  // calcul de la variation du tenseur log, en fait des composantes et non du tenseur totale
        /*          totoBB = ((*d_ki)(i)(1) * (*(Palpha_BH(1))) + ki(1) * (*(d_Palpha_BH(i)(1)))
                           + (*d_ki)(i)(2) * (*(Palpha_BH(2))) + ki(2) * (*(d_Palpha_BH(i)(2)))) * gijBB 
                           + (*B_BH_tr) * (*d_gijBB(i));  */
                  depsBB = (d_log_lambda1 * (*(Palpha_BH(1))) + log_lambda1 * (*(d_Palpha_BH(i)(1)))
                                  +d_log_lambda2 * (*(Palpha_BH(2))) + log_lambda2 * (*(d_Palpha_BH(i)(2)))) * gijBB
                            + epsBH * (*d_gijBB(i));
            //     depsBB =  (interBB + interBB.Transpose());          
            // **** on ajoute la matrice tangente d'Almansi, cela marche mieux mais on ne sait pas pourquoi ****
          //        depsBB +=  0.5* (*d_gijBB(i));       
                 }
                };                
              break;
              }
            default:
              cout << "\n erreur (4) il y a une erreur dans le calcul des valeurs propres pour le tenseur "
                   << " B (Cauchy Green droit), B= ";   
              B_BH_tr->Ecriture(cout);                        
                 if (ParaGlob::NiveauImpression() >= 4) 
                     cout << "\n Met_abstraite::Expli& Deformation::Cal_explicit_Logarithmique (...";
                 Sortie(1);
           };// -- fin du switch sur les différents cas de calcul de valeurs propres      
         // -- fin du cas de la dimension 3 pour les tenseurs 
         break;       
         }
 }; // -- fin du switch sur la dimension des tenseurs    
};

    // calcul de la deformation au temps donné  dans le cas d'une déformation d'Almansi
    // on suppose que les métriques en HH ou BB sont définis à 0 et au temps
void Deformation::Cal_logarithmique_auTemps (Enum_dure temps, TenseurBB & epsBB)
  {  switch (temps)
	  { case TEMPS_t: 
	     { // récup des infos
	       const Met_abstraite::InfoExp_t& ex = metrique->Recup_InfoExp_t(); 
        // appel pour le calcul de la déformation logarithmique et de sa variation
        bool variation = false; // pour la variation de B et Palpha_BH
        Tableau <TenseurBB *>  bidon_d_epsBB; // un tableau qui ne sert pas ici
         Cal_Logarithmique (*ex.gijBB_0,*ex.gijHH_0,bidon_d_epsBB,*ex.gijBB_t,*ex.gijHH_t,epsBB
                            ,bidon_d_epsBB,variation);
        break;
      }
	    case TEMPS_tdt: 
	     { // récup des infos
	       const Met_abstraite::InfoExp_tdt& ex = metrique->Recup_InfoExp_tdt(); 
        // appel pour le calcul de la déformation logarithmique et de sa variation
        bool variation = false; // pour la variation de B et Palpha_BH
        Tableau <TenseurBB *>  bidon_d_epsBB; // un tableau qui ne sert pas ici
         Cal_Logarithmique (*ex.gijBB_0,*ex.gijHH_0,bidon_d_epsBB,*ex.gijBB_tdt,*ex.gijHH_tdt,epsBB
                            ,bidon_d_epsBB,variation);
        break;
      }
     case TEMPS_0:
      { // pb on ne peut pas calculer
        cout << "\n erreur: on demande la def a t=0 , elle n'est pas calculee ! ";
        if (ParaGlob::NiveauImpression() > 0)
          cout << "\n Deformation::Cal_logarithmique_auTemps (.. ";
        cout << endl;
        Sortie(1);
        break;
      }
   };
  };


// -----------------   méthodes de vérifications------- ----
    
void Deformation::VerifCal_deflog(bool gradV_instantane,const Met_abstraite::Impli & ,TenseurBB& epsBB_tdt
                                  ,Tableau <TenseurBB *> & d_epsBB_tdt)
 { // l'idée est de faire une vérification des dérivées à l'aide d'une méthode de différence finie
   int dim = ParaGlob::Dimension();
   // dans le cas du premier passage on indique qu'il y a vérification
   if (indic_VerifCal_implicit == 0)
     { cout << "\n ****vérification du calcul de la déformation et des éléments de la métrique associé****";
       cout << "\n Deformation::VerifCal_implicit \n";
      } 
   indic_VerifCal_implicit++;
  
   // on cré une seconde métrique pour éviter de détruire la métrique originale
   Met_abstraite metrique_bis(*metrique);
   // ici on considère que l'on a même nombre de ddl par noeud = dim
   // on va modifier chaque ddl de chaque noeud systématiquement
   int nbnoeud = tabnoeud->Taille();
   // le deltat pour les différences finis
   double delta = ConstMath::unpeupetit;
   double mini_val = ConstMath::pasmalpetit;
   int numddl = 1; // le compteur de ddl
   bool erreur = false; // indicateur d'erreur
   bool premier_calcul=true;
   for (int inoeud=1;inoeud<=nbnoeud;inoeud++)
    {// on récupère les coordonnées du noeud
     Coordonnee coordtdt = (*tabnoeud)(inoeud)->Coord2();
     for (int ix= 1;ix<=dim;ix++,numddl++)
     { Coordonnee X(dim); X(ix) += delta;
       (*tabnoeud)(inoeud)->Ajout_coord2(X);
       
       // appel de la metrique
       const Met_abstraite::Expli_t_tdt& ex_n = 
          metrique_bis.Cal_explicit_tdt(*tabnoeud,premier_calcul,(*tabDphi)(numInteg),
                                         nbNoeud,(*tabPhi)(numInteg),gradV_instantane);
     TenseurBB * epsBB_tdt_n = NevezTenseurBB(epsBB_tdt);
     //  calcul de la base duale et de la métrique à t=0
     const Met_abstraite::gijHH_0_et_giH_0& ex1_n = metrique_bis.Cal_gijHH_0_et_giH_0_apres_im_expli();         
     bool variation = false; // pour la variation de B et Palpha_BH
     Tableau <TenseurBB *>  d_epsBB_tdt_n(d_epsBB_tdt);
     for (int i=1;i<= d_epsBB_tdt_n.Taille();i++) d_epsBB_tdt_n(i)=NevezTenseurBB(*d_epsBB_tdt(1));
     // appel pour le calcul de la déformation logarithmique et de sa variation
     Cal_Logarithmique (*ex_n.gijBB_0,*ex1_n.gijHH_0,d_epsBB_tdt_n,*ex_n.gijBB_tdt,*ex_n.gijHH_tdt,*epsBB_tdt_n
                        ,*ex_n.d_gijBB_tdt,variation);
     // calcul des dérivées numériques et  vérification
/*     for (int j=1;j<=dim;j++)
        { // variation des vecteurs giB_tdt
          CoordonneeB dgiB = ((*ex_n.giB_tdt)(j) -(*ex.giB_tdt)(j))/delta;
          for (int i=1;i<=dim;i++)
             if (diffpourcent(dgiB(i),(*ex.d_giB_tdt)(numddl)(j)(i),MaX(Dabs(dgiB(i)),Dabs((*ex.d_giB_tdt)(numddl)(j)(i))),0.05))
                if (MiN(Dabs(dgiB(i)),Dabs((*ex.d_giB_tdt)(numddl)(j)(i))) <= mini_val)
                   {if ( MaX(Dabs(dgiB(i)),Dabs((*ex.d_giB_tdt)(numddl)(j)(i))) > 50.*delta) erreur = true;}
                else erreur = true;
          // variation des vecteurs giH_tdt
          CoordonneeH dgiH = ((*ex_n.giH_tdt)(j) - (*ex.giH_tdt)(j))/delta;
          for (int i=1;i<=dim;i++)
             if (diffpourcent(dgiH(i),(*ex.d_giH_tdt)(numddl)(j)(i),MaX(Dabs(dgiH(i)),Dabs((*ex.d_giH_tdt)(numddl)(j)(i))),0.05))
                if (MiN(Dabs(dgiH(i)),Dabs((*ex.d_giH_tdt)(numddl)(j)(i))) <= mini_val)
                   {if ( MaX(Dabs(dgiH(i)),Dabs((*ex.d_giH_tdt)(numddl)(j)(i))) > 50.*delta) erreur = true;}
                else erreur = true;
         }    */   
       // variation du tenseur epsBB_tdt
       TenseurBB * ddepsBB = NevezTenseurBB(epsBB_tdt);
       *ddepsBB = (*epsBB_tdt_n - epsBB_tdt) / delta; // la dérivée numérique
       for (int i1=1;i1<=dim;i1++)
          for (int j1=1;j1<=dim;j1++)
             if (diffpourcent((*ddepsBB)(i1,j1),(*(d_epsBB_tdt(numddl)))(i1,j1),
                MaX(Dabs((*ddepsBB)(i1,j1)),Dabs((*(d_epsBB_tdt(numddl)))(i1,j1))),0.05))
                {if (MiN(Dabs((*ddepsBB)(i1,j1)),Dabs((*(d_epsBB_tdt(numddl)))(i1,j1))) <= mini_val)
                   {if( MaX(Dabs((*ddepsBB)(i1,j1)),Dabs((*(d_epsBB_tdt(numddl)))(i1,j1))) > 50.*delta) erreur = true;}
                 else erreur = true;
                };
  /*     // variation du tenseur gijHH_tdt        
       TenseurHH * gijHH = NevezTenseurHH(*ex_n.gijHH_tdt);
       *gijHH = (*ex_n.gijHH_tdt - *ex.gijHH_tdt) / delta;
       for (int i1=1;i1<=dim;i1++)
           for (int j1=1;j1<=dim;j1++)
             if (diffpourcent((*gijHH)(i1,j1),(*(*ex. d_gijHH_tdt)(numddl))(i1,j1),
                MaX(Dabs((*gijHH)(i1,j1)),Dabs((*(*ex. d_gijHH_tdt)(numddl))(i1,j1))),0.05))
                if (MiN(Dabs((*gijHH)(i1,j1)),Dabs((*(*ex. d_gijHH_tdt)(numddl))(i1,j1))) <= mini_val)
                   {if( MaX(Dabs((*gijHH)(i1,j1)),Dabs((*(*ex. d_gijHH_tdt)(numddl))(i1,j1))) > 50.*delta) erreur = true;}
                else erreur = true;
       // variation du jacobien 
       double djaco = ((*ex_n.jacobien_tdt) - (*ex.jacobien_tdt))/delta;
       if (diffpourcent(djaco,(*ex.d_jacobien_tdt)(numddl),MaX(Dabs(djaco),Dabs((*ex.d_jacobien_tdt)(numddl))),0.05))
            if (MiN(Dabs(djaco),Dabs((*ex.d_jacobien_tdt)(numddl))) <= mini_val)
              {if( MaX(Dabs(djaco),Dabs((*ex.d_jacobien_tdt)(numddl))) > 50.*delta) erreur = true;}
                else erreur = true;
       // effacement des tenseurs intermédiaires
       delete gijBB; delete  gijHH;        
       // maintenant on remet les coordonnées du noeud à l'état initial         
       (*tabnoeud)(inoeud)->Change_coord2(coordtdt); */       
      } // fin de boucle sur la dimension de coordonnée
    }  // fin de boucle sur les noeuds
  // message d'erreur si besoin
  if (erreur)
     { cout << "\n erreur dans le calcul analytique des dérivees de la metrique";
       cout << "\n Deformation::VerifCal_implicit(.."
            << " , numero d'increment = " << indic_VerifCal_implicit;
       Sortie(1);
      } 
 };