// 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-2021 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 "Tenseur.h"
#include "MathUtil.h"
#include "MathUtil2.h"
#include "ParaGlob.h"

#ifndef  Tenseur_H_deja_inclus

//-----------------------------------------------------------------------
// def de templates qui servent à diminuer les temps de calcul pour
// les méthodes de calcul des tenseurs dans la base absolue
//-----------------------------------------------------------------------
// 1) pour le passage repère local en repère absolu
//  T1 et T2 sont des tenseurs, Base représente une base
// pour des tenseurs T1 et T2 symétriques HH ou BB      
template <class T1,class T2,class Base>  
 T1&  produit1(T1& A,const T2& Nous, const Base & gi)
    {  int dim = abs(A.Dimension());// dimension du tenseur résultat
       int dimint = abs(Nous.Dimension());// dimension du tenseur en local
       for (int al=1;al<= dim; al++)
          for (int be=1; be<= al; be++)
             { A.Coor(al,be) = 0.;
               for (int i=1; i<= dimint; i++)
                  for (int j=1; j<= dimint; j++)
                      A.Coor(al,be) += Nous(i,j)*gi(i)(al)*gi(j)(be);
             }
       return A;
    };
// pour un tenseurs T1 non symétrique (HB ou BH) et T2 symétriques HH ou BB      
template <class T1,class T2,class Base>  
 T1&  produit2(T1& A,const T2& Nous, const Base & gi)
    {  int dim = abs(A.Dimension());
        int dimint = abs(Nous.Dimension());
        for (int al=1;al<= dim; al++)
             for (int be=1; be<= dim; be++)
               { A.Coor(al,be) = 0.;
                  for (int i=1; i<= dimint; i++)
                     for (int j=1; j<= dimint; j++)
                        A.Coor(al,be) += Nous(i,j)*gi(i)(al)*gi(j)(be);
                 }       
        return A;      
      };
// pour des tenseurs non symétriques en variance T1 et T2
// il y a deux bases, une de chaque variance     
template <class T1,class T2,class Base2,class Base1>  
T1 & produit3(T1 & A,const T2& Nous, const Base2 & gi2, const Base1 & gi1)
    {  int dim = abs(A.Dimension());
        int dimint = abs(Nous.Dimension());
        for (int i=1;i<= dim; i++)
             for (int j=1; j<= dim; j++)
               { A.Coor(i,j) = 0.;
                  for (int al=1; al<= dimint; al++)
                     for (int be=1; be<= dimint; be++)
                        A.Coor(i,j) += Nous(al,be)*gi1(al)(i)*gi2(be)(j);
                 }       
        return A;      
      };
// 2) pour le passage repère local en repère absolu
//  T1 et T2 sont des tenseurs, Base représente une base
// pour des tenseurs T1 et T2 symétriques HH ou BB      
template <class T1,class T2,class Base>  
 T1&  produit4(T1& A,const T2& Nous, const Base & gi)
    {  int dim = abs(A.Dimension());
        int dimint = abs(Nous.Dimension());
        for (int i=1;i<= dim; i++)
             for (int j=1; j<= i; j++)
               { A.Coor(i,j) = 0.;
                  for (int al=1; al<= dimint; al++)
                     for (int be=1; be<= dimint; be++)
                        A.Coor(i,j) += Nous(al,be)*gi(i)(al)*gi(j)(be);
                 }       
        return A;      
      };
//  la même chose que produit4, mais lorsque le tenseur A est de stockage
// non symétrique, il faut donc à ce moment là calculer A(i,j) et A(j,i)      
template <class T1,class T2,class Base>  
 T1&  produit6(T1& A,const T2& Nous, const Base & gi)
    {  int dim = abs(A.Dimension());
        int dimint = abs(Nous.Dimension());
        for (int i=1;i<= dim; i++)
             for (int j=1; j<= dim; j++)
               { A.Coor(i,j) = 0.;
                  for (int al=1; al<= dimint; al++)
                     for (int be=1; be<= dimint; be++)
                        A.Coor(i,j) += Nous(al,be)*gi(i)(al)*gi(j)(be);
                 }       
        return A;      
      };
// pour des tenseurs non symétriques en variance T1 et T2
// il y a deux bases, une de chaque variance     
template <class T1,class T2,class Base2,class Base1>  
T1 & produit5(T1 & A,const T2& Nous, const Base2 & gi2, const Base1 & gi1)
    {  int dim = abs(A.Dimension());
        int dimint = abs(Nous.Dimension());
        for (int i=1;i<= dim; i++)
             for (int j=1; j<= dim; j++)
               { A.Coor(i,j) = 0.;
                  for (int al=1; al<= dimint; al++)
                     for (int be=1; be<= dimint; be++)
 //                       A(i,j) += Nous(al,be)*gi2(al)(i)*gi1(be)(j);
                        A.Coor(i,j) += Nous(al,be)*gi2(i)(al)*gi1(j)(be);
                 }       
        return A;      
      };

//----------- fin de la définition de template --------------------------

// produit contracte avec un vecteur
#ifndef MISE_AU_POINT
  inline 
#endif
 CoordonneeH operator * (const CoordonneeH & v ,const TenseurBH & t)
      { return (t.Transpose() * v) ; };

// ------------------tenseur HH ----------------------------------
//PtTenseurHH * LesMaillonsHH::maille = NULL; // initialisation de la première maille
 #ifdef MISE_AU_POINT
   int LesMaillonsHH::nbmailHH =0; // initialisation du nombre de tenseur intermediaire
 #endif

// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurHH 
   { public :
       PtTenseurHH * t1;     // adresse du maillon precedent
       TenseurHH * sauvePtr; // sauvegarde du pointeur de tenseur
       PtTenseurHH (const PtTenseurHH * x1,const TenseurHH * x2) 
               // ici on ruse pour endormir le compilateur avec un cast
               // par la suite on utilisera t1 et x2 pour effacer le tenseur !!
               // ce qui n'est pas possible avec const !! normalement
             { t1 = (PtTenseurHH *) x1;sauvePtr = (TenseurHH *) x2;};
       ~PtTenseurHH () 
         { delete sauvePtr;};
    }; 
// nouveau maillon
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsHH::NouveauMaillon(const TenseurHH * ptr)
   { PtTenseurHH * poi = maille;
     maille = new PtTenseurHH ( poi,ptr);
     #ifdef MISE_AU_POINT
        nbmailHH++;
     #endif
   };
// liberation de l'espace
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsHH::Libere()
   { while (maille != NULL)
      { PtTenseurHH * poi = maille->t1; // recup de l'adresse precedente
        delete maille; // liberation de l'espace de tenseur
        maille = poi;
      }; 
     #ifdef MISE_AU_POINT
           nbmailHH = 0;
     #endif
   };
	
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurHH::Message(int dim, string mes) const 
  { // on sort le message en fonction du niveau de commentaires
    if (ParaGlob::NiveauImpression() != 0)
      { cout << " \n erreur de calcul sur les tenseurs HH ***** ";
        cout << " dimension de l'argument = " << dim ;
        cout << '\n' << mes << endl;
       };
    // dans le cas où le niveau est supérieur à une valeur on s'arrête   
    if (ParaGlob::NiveauImpression() > 5)
       Sortie(1);
  };

    // retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
  inline
#endif
Mat_pleine& TenseurHH::Matrice_composante(Mat_pleine& res)const
   { int dim = Abs (this->dimension);
     if  ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
       res.Initialise(dim,dim,0.);
     if (this->dimension < 0)
      // cas non symétrique
      {for (int i=1;i<=dim;i++)
        for (int j=1;j<=dim;j++)
         res(i,j) = (*this)(i,j);
      }
     else // cas symétrique
      {for (int i=1;i<=dim;i++)    // le résultat peut ne pas être symétrique
        for (int j=1;j<=i;j++)
         res(j,i) = res(i,j) = (*this)(i,j);
      };
     return res;
   };

    // opération inverse: affectation en fonction d'une matrice
    // donc sans vérif de variance !
#ifndef MISE_AU_POINT
  inline
#endif
void TenseurHH::Affectation(const Mat_pleine& a)
   { int dim = Abs (this->dimension);
     #ifdef MISE_AU_POINT
      if  ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
      { cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
             << a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
             << " a une dimension differente du tenseur = "
             << dim << " !\n";
        cout << "TenseurHH::Affectation(const Mat_pleine& a)  \n";
        Sortie(1);
      };
     #endif

     if (this->dimension < 0)
      // cas non symétrique
      {for (int i=1;i<=dim;i++)
        for (int j=1;j<=dim;j++)
          (*this).Coor(i,j) = a(i,j);
      }
     else // cas symétrique
      {
       #ifdef MISE_AU_POINT
       if (!a.Symetrie())
         {cout << "\n attention la matrice est non symetrique et on veut l'affecter a"
               << " un tenseur symetrique !";
          a.Affiche();
          cout << "TenseurHH::Affectation(const Mat_pleine& a)  \n";
         };
       #endif

       for (int i=1;i<=dim;i++)    // le résultat peut ne pas être symétrique
        for (int j=1;j<=i;j++)     // à un chouat près
          (*this).Coor(i,j) = 0.5* (a(j,i) + a(i,j));
      };
   };


    // changement de base (cf. théorie) : la matrice beta est telle que:
    // par défaut inverse = false : (c'est l'utilisation historique)
    //     beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
    //     gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
    // gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    // et on a la matrice gamma telle que:
    // gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    // gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    //   c-a-d= gp^i = gamma^i_j * g^j
    // rappel des différentes relations entre beta et gamma
    // [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    // [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    // changement de base pour de deux fois contravariants:
    // [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
    //     donc le changement de base s'effectue directement à l'aide de gamma


    // si inverse = true:
    //     beta(i,j) représente les coordonnées de l'ancienne base gB dans la nouvelle gpB
    //     gB(i) = beta(i,j) * gpB(j), i indice de ligne, j indice de colonne
    //     la formule de changement de base est alors:
    // [Ap^kl] = [beta]^T * [A^ij] * [beta]
    //
    //   nécessite d'inverser la matrice gamma pour calculer beta
    //
    //
    //
#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurHH::ChBase(const Mat_pleine& gamma,bool inverse)
  { // marche également si le tenseur n'est pas symétrique
    int dim = Abs (this->dimension); 
    #ifdef MISE_AU_POINT	 	 
     if  ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
     { cout << "\nErreur : la nouvelle base doit etre de dimension " 
            << dim << " !\n";
       cout << "TenseurHH::ChBase(Mat_pleine& gamma)  \n";
       Sortie(1);
     };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)    // le résultat peut ne pas être symétrique
      for (int j=1;j<=dim;j++)
        res(i,j) = (*this)(i,j);
    if (inverse)
     {// calcul de la matrice inverse transposée
      Mat_pleine beta = (gamma.Inverse()).Transpose();
////--- debug
//cout << "\n debug TenseurHH::ChBase( ";
//cout << "\n beta = "; beta.Affiche();
//cout << "\n res avant calcul = "; res.Affiche();
    // [Ap^kl] = [beta]^T * [A^ij] * [beta]

     res = ((beta.Transpose()) * res) * beta;
//cout << "\n res après calcul = "; res.Affiche();
//// fin debug ---
     }
    else
     {
////--- debug
//cout << "\n debug TenseurHH::ChBase( ";
//cout << "\n beta = "; beta.Affiche();
//cout << "\n gamma = "; gamma.Affiche();
//cout << "\n res avant calcul = "; res.Affiche();
      
    // [Ap^kl] = [gamma] * [A^ij] * [gamma]^T

      // calcul de la transformation
//      res = ((gamma.Transpose()) * res) * gamma; ==> non  c'est faux voir théorie
      res = (gamma * res) * gamma.Transpose();
      
//cout << "\n res après calcul = "; res.Affiche();
//// fin debug ---

     };
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        this->Coor(i,j) = res(i,j) ;
  };                

// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
//  this       : le tenseur
//  var_tensHH : en entrée: la variation du tenseur this dans la base initiale qu'on appelle g_i
//  var_tensHH : en sortie: la variation du tenseur dans la base finale qu'on appelle gp_i
// gamma       : en entrée gpH(i) = gamma(i,j) * gH(j)
// var_gamma   : en entrée : la variation de gamma
// [Ap^kl] = [gamma] * [A^ij] * [gamma]^T

#ifndef MISE_AU_POINT
  inline
#endif
void TenseurHH::Var_tenseur_dans_nouvelle_base
          (const Mat_pleine& gamma,TenseurHH& var_tensHH, const Mat_pleine& var_gamma)
 {  int dim = Abs (this->dimension);
    #ifdef MISE_AU_POINT
    if  ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
     { cout << "\nErreur : la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurHH::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    if  ((var_gamma.Nb_ligne() != dim) || (var_gamma.Nb_colonne() != dim))
     { cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurHH::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    Mat_pleine res1(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j)= (*this)(i,j);
    // [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
    // donc d[Ap^kl] = d[gamma] * [A_^ij] * [gamma]^T
    //                 + [gamma] * d[A^ij] * [gamma]^T
    //                 + [gamma] * [Aij] * d[gamma]^T
    res = var_gamma * res1 * gamma.Transpose();
    res += gamma * res1 * var_gamma.Transpose();
    // récup des variations
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = var_tensHH(i,j);
    // ajout de la ligne restante: res maintenant correspond à d[A_ij]
    res += gamma * res1 * gamma.Transpose();
  
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        var_tensHH.Coor(i,j) = res(i,j) ;
 };


// calcul des composantes du tenseur dans la base absolue
// la variance du résultat peut-être  quelconque d'où quatre possibilités
// en fonction de l'argument A. Dans tous les cas les composantes sont identiques
// car la sortie est en absolue 
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHH & TenseurHH::BaseAbsolue(TenseurHH & A,const BaseB & gi) const
    {    // on traite en fonction également de la symétrie de stockage ou pas
         if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
         else                                        { return produit2(A,(*this),gi);};
     };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBB & TenseurHH::BaseAbsolue(TenseurBB & A,const BaseB & gi) const
    {    // on traite en fonction également de la symétrie de stockage ou pas
         if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
         else                                        { return produit2(A,(*this),gi);};      
      };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBH & TenseurHH::BaseAbsolue(TenseurBH & A,const BaseB & gi) const
    {    return produit2(A,(*this),gi);      
      };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHB & TenseurHH::BaseAbsolue(TenseurHB& A,const BaseB & gi) const
    {    return produit2(A,(*this),gi);      
      };
// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes ,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHH & TenseurHH::Baselocale(TenseurHH & A,const BaseH & gi) const
 {  // on traite en fonction également de la symétrie de stockage ou pas
    if ((A.Dimension() > 0) && (dimension > 0)) { return produit4(A,(*this),gi);}
    else                                        { return produit6(A,(*this),gi);};
  };
     
// ------------------tenseur BB ----------------------------------
//PtTenseurBB * LesMaillonsBB::maille = NULL; // initialisation de la première maille
 #ifdef MISE_AU_POINT
   int LesMaillonsBB::nbmailBB =0; // initialisation du nombre de tenseur intermediaire
 #endif

// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurBB 
   { public :
       PtTenseurBB * t1;     // adresse du maillon precedent
       TenseurBB * sauvePtr; // sauvegarde du pointeur de tenseur
       PtTenseurBB (const PtTenseurBB * x1,const TenseurBB * x2) 
               // ici on ruse pour endormir le compilateur avec un cast
               // par la suite on utilisera t1 et x2 pour effacer le tenseur !!
               // ce qui n'est pas possible avec const !! normalement
             { t1 = (PtTenseurBB *) x1;sauvePtr = (TenseurBB *) x2;};
       ~PtTenseurBB () 
         { delete sauvePtr;};
    }; 
// nouveau maillon
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsBB::NouveauMaillon(const TenseurBB * ptr)
   { PtTenseurBB * poi = maille;
     maille = new PtTenseurBB ( poi,ptr);
     #ifdef MISE_AU_POINT
        nbmailBB++;
     #endif
   };
// liberation de l'espace
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsBB::Libere()
   { while (maille != NULL)
      { PtTenseurBB * poi = maille->t1; // recup de l'adresse precedente
        delete maille; // liberation de l'espace de tenseur
        maille = poi;
      }; 
     #ifdef MISE_AU_POINT
           nbmailBB = 0;
     #endif
   };
	
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurBB::Message(int dim, string mes) const 
  { cout << " \n erreur de calcul sur les tenseurs BB ***** ";
    cout << " dimension de l'argument = " << dim ;
    cout << '\n' << mes << endl;
    Sortie(1);
  }  

    // retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
  inline
#endif
Mat_pleine& TenseurBB::Matrice_composante(Mat_pleine& res)const
   { int dim = Abs (this->dimension);
     if  ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
       res.Initialise(dim,dim,0.);
     if (this->dimension < 0)
      // cas non symétrique
      {for (int i=1;i<=dim;i++)
        for (int j=1;j<=dim;j++)
         res(i,j) = (*this)(i,j);
      }
     else // cas symétrique
      {for (int i=1;i<=dim;i++)    // le résultat peut ne pas être symétrique
        for (int j=1;j<=i;j++)
         res(j,i) = res(i,j) = (*this)(i,j);
      };
     return res;
   };

    // opération inverse: affectation en fonction d'une matrice
    // donc sans vérif de variance !
#ifndef MISE_AU_POINT
  inline
#endif
void TenseurBB::Affectation(const Mat_pleine& a)
   { int dim = Abs (this->dimension);
     #ifdef MISE_AU_POINT
      if  ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
      { cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
             << a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
             << " a une dimension differente du tenseur = "
             << dim << " !\n";
        cout << "TenseurBB::Affectation(const Mat_pleine& a)  \n";
        Sortie(1);
      };
     #endif

     if (this->dimension < 0)
      // cas non symétrique
      {for (int i=1;i<=dim;i++)
        for (int j=1;j<=dim;j++)
          (*this).Coor(i,j) = a(i,j);
      }
     else // cas symétrique
      {
       #ifdef MISE_AU_POINT
       if (!a.Symetrie())
         {cout << "\n attention la matrice est non symetrique et on veut l'affecter a"
               << " un tenseur symetrique !";
          a.Affiche();
          cout << "TenseurBB::Affectation(const Mat_pleine& a)  \n";
         };
       #endif

       for (int i=1;i<=dim;i++)    // le résultat peut ne pas être symétrique
        for (int j=1;j<=i;j++)     // à un chouat près
          (*this).Coor(i,j) = 0.5* (a(j,i) + a(i,j));
      };
   };

    // changement de base (cf. théorie) : la matrice beta est telle que:
    // gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    // et la matrice gamma telle que:
    // gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    // gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    //   c-a-d= gp^i = gamma^i_j * g^j
    // rappel des différentes relations entre beta et gamma
    // [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    // [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    // changement de base pour de deux fois covariants:
    // [Ap_kl] = [beta] * [A_ij] * [beta]^T

#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurBB::ChBase(const Mat_pleine& beta)
  { // marche également si le tenseur n'est pas symétrique
    int dim = Abs (this->dimension); 
    #ifdef MISE_AU_POINT	 	 
    if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
     { cout << "\nErreur : la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurBB::ChBase(Mat_pleine& beta)  \n";
       Sortie(1);
     };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res(i,j) = (*this)(i,j); 
    // calcul de la matrice inverse
         // on n'en a pas besoin mais comme cela le determinant est verifie
// non on vire car c'est débile de faire autant de calcul pour un déterminant         
//    Mat_pleine beta = beta.Inverse();
    // calcul de la transformation
//    res = ((beta.Transpose()) * res) * beta; ==> vieux calcul erreur
    res = (beta * res) * beta.Transpose();
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        this->Coor(i,j) = res(i,j) ;
  };

// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
//  this       : le tenseur
// var_tensBB : en entrée: la variation du tenseur dans la base initiale qu'on appelle g^i
// var_tensBB : en sortie: la variation du tenseur dans la base finale qu'on appelle gp^i
// beta       : en entrée gpB(i) = beta(i,j) * gB(j)
// var_beta   : en entrée : la variation de beta

#ifndef MISE_AU_POINT
  inline
#endif
void TenseurBB::Var_tenseur_dans_nouvelle_base
          (const Mat_pleine& beta,TenseurBB& var_tensBB, const Mat_pleine& var_beta)
 {  int dim = Abs (this->dimension);
    #ifdef MISE_AU_POINT
    if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
     { cout << "\nErreur : la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurBB::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    if  ((var_beta.Nb_ligne() != dim) || (var_beta.Nb_colonne() != dim))
     { cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurBB::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    Mat_pleine res1(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = (*this)(i,j);
    // [Ap_kl] = [beta] * [A_ij] * [beta]^T
    // donc d[Ap_kl] = d[beta] * [A_ij] * [beta]^T
    //                 + [beta] * d[A_ij] * [beta]^T
    //                 + [beta] * [A_ij] * d[beta]^T
    res = var_beta * res1 * beta.Transpose();
    res += beta * res1 * var_beta.Transpose();
    // récup des variations
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = var_tensBB(i,j);
    // ajout de la ligne restante: res maintenant correspond à d[A_ij]
    res += beta * res1 * beta.Transpose();
  
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        var_tensBB.Coor(i,j) = res(i,j) ;
 };

// calcul des composantes du tenseur dans la base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBB & TenseurBB::BaseAbsolue(TenseurBB & A,const BaseH & gi)  const
    {    // on traite en fonction également de la symétrie de stockage ou pas
         if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
         else                                        { return produit2(A,(*this),gi);};
     };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHH & TenseurBB::BaseAbsolue(TenseurHH & A,const BaseH & gi)  const
    {    // on traite en fonction également de la symétrie de stockage ou pas
         if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
         else                                        { return produit2(A,(*this),gi);};
     };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBH & TenseurBB::BaseAbsolue(TenseurBH & A,const BaseH & gi) const
    {    return produit2(A,(*this),gi);      
      };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHB & TenseurBB::BaseAbsolue(TenseurHB & A,const BaseH & gi)  const
    {    return produit2(A,(*this),gi);      
      };
// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes ,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBB & TenseurBB::Baselocale(TenseurBB & A,const BaseB & gi)  const
 {  // on traite en fonction également de la symétrie de stockage ou pas
    if ((A.Dimension() > 0) && (dimension > 0)) { return produit4(A,(*this),gi);}
    else                                        { return produit6(A,(*this),gi);};
  };
    
    
// ------------------tenseur BH ----------------------------------
//PtTenseurBH * LesMaillonsBH::maille = NULL; // initialisation de la première maille
 #ifdef MISE_AU_POINT
   int LesMaillonsBH::nbmailBH =0; // initialisation du nombre de tenseur intermediaire
 #endif

// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurBH 
   { public :
       PtTenseurBH * t1;     // adresse du maillon precedent
       TenseurBH * sauvePtr; // sauvegarde du pointeur de tenseur
       PtTenseurBH (const PtTenseurBH * x1,const TenseurBH * x2) 
               // ici on ruse pour endormir le compilateur avec un cast
               // par la suite on utilisera t1 et x2 pour effacer le tenseur !!
               // ce qui n'est pas possible avec const !! normalement
             { t1 = (PtTenseurBH *) x1;sauvePtr = (TenseurBH *) x2;};
       ~PtTenseurBH () 
         { delete sauvePtr;};
    }; 
// nouveau maillon
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsBH::NouveauMaillon(const TenseurBH * ptr)
   { PtTenseurBH * poi = maille;
     maille = new PtTenseurBH ( poi,ptr);
     #ifdef MISE_AU_POINT
        nbmailBH++;
     #endif
   };
// liberation de l'espace
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsBH::Libere()
   { while (maille != NULL)
      { PtTenseurBH * poi = maille->t1; // recup de l'adresse precedente
        delete maille; // liberation de l'espace de tenseur
        maille = poi;
      }; 
     #ifdef MISE_AU_POINT
           nbmailBH = 0;
     #endif
   };
	
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurBH::Message(int dim, string mes) const 
  { cout << " \n erreur de calcul sur les tenseurs BH ***** ";
    cout << " dimension de l'argument = " << dim ;
    cout << '\n' << mes << endl;
    Sortie(1);
  }  

    // retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
  inline
#endif
Mat_pleine& TenseurBH::Matrice_composante(Mat_pleine& res)const
   { int dim = Abs (this->dimension);
     if  ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
       res.Initialise(dim,dim,0.);
     // toujours non symétrique
     for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
         res(i,j) = (*this)(i,j);
     return res;
   };

    // opération inverse: affectation en fonction d'une matrice
    // donc sans vérif de variance !
#ifndef MISE_AU_POINT
  inline
#endif
void TenseurBH::Affectation(const Mat_pleine& a)
   { int dim = Abs (this->dimension);
     #ifdef MISE_AU_POINT
      if  ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
      { cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
             << a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
             << " a une dimension differente du tenseur = "
             << dim << " !\n";
        cout << "TenseurBH::Affectation(const Mat_pleine& a)  \n";
        Sortie(1);
      };
     #endif

     // on est toujours en non symétrique
     for (int i=1;i<=dim;i++)
        for (int j=1;j<=dim;j++)
          (*this).Coor(i,j) = a(i,j);

   };

/* on vire la méthode générale qui utilise NRC, car il y a une implantation spécifique bien plus
  performante !!

    // valeurs propre  dans le vecteur  de retour, classée par ordres "décroissants"
    // cas indique le cas de valeur propre: 
    // quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
    //                            dans ce cas les valeurs propres de retour sont nulles par défaut
    // dim = 1, cas=1 pour une valeur propre; 
    // dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
    // dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
    //         , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
    //         , cas = 2 si val(1)=val(2) ( val(1), et val(3)  dans les 2 premières composantes du retour)
    //         , cas = 3 si val(2)=val(3) ( val(1), et val(2)  dans les 2 premières composantes du retour)
#ifndef MISE_AU_POINT
  inline 
#endif
Coordonnee TenseurBH::ValPropre(int& cas) const 
  { int dim = this->dimension;
    Mat_pleine mat(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        mat(i,j) = (*this)(i,j);        
    int cass; Coordonnee vec =  VectValPropre(mat,cass);
    // classement des valeurs propres
    // cas ou dim == 1 on ne fait rien
    Tableau <int> ii(dim); // indicateur de l'ordre final
    switch (dim)
     { case 1: ii(1)=1; break;
       case 2: 
       { if (vec(2) > vec(1))  { ii(1)=2;ii(2)=1;}
         else                  { ii(1)=1;ii(2)=2;};
         break;
        }
       case 3:
       { if (vec(1) >= vec(2)) 
           { if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         else  
           { if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         break;
         }   
      };    
    // maintenant on met dans l'ordre
    Vecteur vec_bis(vec);
    for (int i=1;i<=dim;i++)  vec(i)=vec_bis(ii(i));     
    // définition du cas de valeurs propres
    if ( difftrespetit(vec(1),vec(2)) && difftrespetit(vec(2),vec(3))) {cas = 1;}
    else if ( !(difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 0;}
    else if ( (difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 2;}
    else if ( !(difftrespetit(vec(1),vec(2))) && difftrespetit(vec(2),vec(3))) {cas = 3;}
    // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
    if (cass <= 0) cas=-1;
    // retour de vec et de cas
    return vec;
   };      
*/
// idem met en retour la matrice mat contiend par colonne les vecteurs propre
// elle doit avoir la dimension du tenseur
// les vecteurs propre sont exprime dans le repere dual (contrairement au HB)
    // valeurs propre  dans le vecteur  de retour, classée par ordres "décroissants"
    // cas indique le cas de valeur propre: 
    // quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
    //                            dans ce cas les valeurs propres de retour sont nulles par défaut
    // dim = 1, cas=1 pour une valeur propre; 
    // dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
    // dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
    //         , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
    //         , cas = 2 si val(1)=val(2) ( val(1), et val(3)  dans les 2 premières composantes du retour)
    //         , cas = 3 si val(2)=val(3) ( val(1), et val(2)  dans les 2 premières composantes du retour)
    // les vecteurs propre sont exprime dans le repere dual (contrairement au HB) pour les tenseurs dim 3
    // pour dim=2:le premier vecteur propre est exprime dans le repere dual 
    //            le second vecteur propre est exprimé dans le repère naturel
    // pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas) 

#ifndef MISE_AU_POINT
  inline 
#endif
Coordonnee TenseurBH::ValPropre_int(int& cas, Mat_pleine& mat) const   
  { int dim = this->dimension;
  
//---- a priori c'est méthode n'est plus utilisée par les tenseurs donc à virer par la suite ------
  
    #ifdef MISE_AU_POINT	 	 
    if  ((mat.Nb_ligne() != dim) || (mat.Nb_colonne() != dim))
      { cout << "\nErreur : la matrice en parametre doit etre de dimension "
             << dim << " !\n";
        cout << "TenseurBH::ValPropre(Mat_pleine& mat)  \n";
        Sortie(1);
      };
    #endif
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        mat(i,j) = (*this)(i,j);
    int cass; 
    Vecteur vec(MathUtil2::VectValPropre(mat,cass)); // calcul des vecteurs et valeurs propres
    // classement des valeurs propres
    // cas ou dim == 1 on ne fait rien
    Tableau <int> ii(dim); // indicateur de l'ordre final
    switch (dim)
     { case 1: ii(1)=1; break;
       case 2: 
       { if (vec(2) > vec(1))  { ii(1)=2;ii(2)=1;}
         else                  { ii(1)=1;ii(2)=2;};
         break;
        }
       case 3:
       { if (vec(1) >= vec(2)) // d'abord le plus grand
           { if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         else  
           { if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         // puis les deux derniers
         if (vec(ii(2)) < vec(ii(3))) {int ix=ii(2); ii(2)=ii(3);ii(3)=ix;};   
         break;
         }   
      };    
    // maintenant on met dans l'ordre
    Mat_pleine mat_bis(mat);
    Coordonnee coor(vec.Taille());  
    for (int i=1;i<=dim;i++)  
      {coor(i)=vec(ii(i)); mat.RemplaceColonneSpe(i,mat_bis.Colonne (ii(i)));};
   
    // définition du cas de valeurs propres
    switch (dim)
     { case 1:  break; // il n'y a rien à faire
       case 2: 
        { if (coor(2) > coor(1))
           { double x = coor(1); coor(1) = coor(2); coor(2)=x;};
          // def du cas
          if ( difftrespetit(coor(1),coor(2))) { cas = 1;} else {cas = 0;};
          break;
        }
       case 3:
        { if ( difftrespetit(coor(1),coor(2)) && difftrespetit(coor(2),coor(3))) {cas = 1;}
          else if ( !(difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 0;}
          else if ( (difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 2;}
          else if ( !(difftrespetit(coor(1),coor(2))) && difftrespetit(coor(2),coor(3))) {cas = 3;}
            {cout << "\n classement des valeurs propres : non realise !! "
                  << "\n valeurs propres: coor(1)= "<<coor(1) << ", coor(2)= "<<coor(2) << ", coor(3)= "<<coor(3)
                  << "\n TenseurBH::ValPropre_int(... "<<endl;
            };
          break;
        };
       default: ; // sinon rien (ne devrait pas arriver !
     };
   
    // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
    if (cass <= 0) cas=-1;
    // retour de vec et de cas
    return coor;
   };      


    // changement de base (cf. théorie) : la matrice beta est telle que:
    // gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    // et la matrice gamma telle que:
    // gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    // gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    //   c-a-d= gp^i = gamma^i_j * g^j
    // rappel des différentes relations entre beta et gamma
    // [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    // [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    // formule de changement de base
    // [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
    //     beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
    //     gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne

#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurBH::ChBase( const Mat_pleine& beta,const Mat_pleine& gamma)
  { // marche également si le tenseur n'est pas symétrique
    int dim = Abs (this->dimension); 
    #ifdef MISE_AU_POINT	 	 
	    if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
       { cout << "\nErreur : la nouvelle base doit etre de dimension "
              << dim << " !\n";
         cout << "TenseurBH::ChBase(Mat_pleine& beta...  \n";
         Sortie(1);
       };
     if  ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
       { cout << "\nErreur : la nouvelle base doit etre de dimension "
              << dim << " !\n";
         cout << "TenseurBH::ChBase(Mat_pleine& beta, Mat_pleine& gamma)  \n";
         Sortie(1);
       };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res(i,j) = (*this)(i,j); 
    // calcul de la matrice inverse transposée
//    Mat_pleine gamma = (beta.Inverse()).Transpose();
    // calcul de la transformation
//    res = (gamma.Transpose() * res) * beta; ===> vieux calcul erreur voir théorie
    // [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
    res = (beta * res) * gamma.Transpose();
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        this->Coor(i,j) = res(i,j) ;
  };                

// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
//  this       : le tenseur
// var_tensBB : en entrée: la variation du tenseur dans les bases initiales qu'on appelle g_i et g^i
// var_tensBH : en sortie: la variation du tenseur dans les bases finales qu'on appelle gp_i et gp^j
// beta       : en entrée gpB(i) = beta(i,j) * gB(j)
// var_beta   : en entrée : la variation de beta
// gamma       : en entrée gpH(i) = gamma(i,j) * gH(j)
// var_gamma   : en entrée : la variation de gamma
    // [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
#ifndef MISE_AU_POINT
  inline
#endif
void TenseurBH::Var_tenseur_dans_nouvelle_base
          (const Mat_pleine& beta,TenseurBH& var_tensBH, const Mat_pleine& var_beta
           ,const Mat_pleine& gamma,const Mat_pleine& var_gamma)
 {  int dim = Abs (this->dimension);
    #ifdef MISE_AU_POINT
    if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
     { cout << "\nErreur : la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurBH::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    if  ((var_beta.Nb_ligne() != dim) || (var_beta.Nb_colonne() != dim))
     { cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurBH::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    Mat_pleine res1(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = (*this)(i,j);
    // [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
    // donc d[Ap_k^l] = d[beta] * [A_i^j] * [gamma]^T
    //                 + [beta] * [A_i^j] * d[gamma]^T
    //                 + [beta] * d[A_i^j] * [gamma]^T
    res = var_beta * res1 * gamma.Transpose();
    res += beta * res1 * var_gamma.Transpose();
    // récup des variations
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = var_tensBH(i,j);
    // ajout de la ligne restante: res maintenant correspond à d[A_i^j]
    res += beta * res1 * gamma.Transpose();
  
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        var_tensBH.Coor(i,j) = res(i,j) ;
 };


// calcul des composantes du tenseur dans la base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBH & TenseurBH::BaseAbsolue(TenseurBH & A,const BaseB & giB,const BaseH & giH)  const
    {    return produit3(A,(*this),giB,giH);      
      };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHB & TenseurBH::BaseAbsolue(TenseurHB & A,const BaseB & giB,const BaseH & giH)  const
    {    return produit3(A,(*this),giB,giH);      
      };
// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes ,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBH & TenseurBH::Baselocale(TenseurBH & A,const BaseH & giH,const BaseB & giB)  const
 {  return produit5(A,(*this),giB,giH);
  };
       
    
// ------------------tenseur HB ----------------------------------
//PtTenseurHB * LesMaillonsHB::maille = NULL; // initialisation de la première maille
 #ifdef MISE_AU_POINT
   int LesMaillonsHB::nbmailHB =0; // initialisation du nombre de tenseur intermediaire
 #endif

// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurHB 
   { public :
       PtTenseurHB * t1;     // adresse du maillon precedent
       TenseurHB * sauvePtr; // sauvegarde du pointeur de tenseur
       PtTenseurHB (const PtTenseurHB * x1,const TenseurHB * x2) 
               // ici on ruse pour endormir le compilateur avec un cast
               // par la suite on utilisera t1 et x2 pour effacer le tenseur !!
               // ce qui n'est pas possible avec const !! normalement
             { t1 = (PtTenseurHB *) x1;sauvePtr = (TenseurHB *) x2;};
       ~PtTenseurHB () 
         { delete sauvePtr;};
    }; 
// nouveau maillon
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsHB::NouveauMaillon(const TenseurHB * ptr)
   { PtTenseurHB * poi = maille;
     maille = new PtTenseurHB ( poi,ptr);
     #ifdef MISE_AU_POINT
        nbmailHB++;
     #endif
   };
// liberation de l'espace
#ifndef MISE_AU_POINT
  inline 
#endif
void LesMaillonsHB::Libere()
   { while (maille != NULL)
      { PtTenseurHB * poi = maille->t1; // recup de l'adresse precedente
        delete maille; // liberation de l'espace de tenseur
        maille = poi;
      }; 
     #ifdef MISE_AU_POINT
           nbmailHB = 0;
     #endif
   };
	
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurHB::Message(int dim, string mes) const 
  { cout << " \n erreur de calcul sur les tenseurs HB ***** ";
    cout << " dimension de l'argument = " << dim ;
    cout << '\n' << mes << endl;
    Sortie(1);
  }  

    // retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
  inline
#endif
Mat_pleine& TenseurHB::Matrice_composante(Mat_pleine& res)const
   { int dim = Abs (this->dimension);
     if  ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
       res.Initialise(dim,dim,0.);
     // toujours non symétrique
     for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
         res(i,j) = (*this)(i,j);
     return res;
   };

    // opération inverse: affectation en fonction d'une matrice
    // donc sans vérif de variance !
#ifndef MISE_AU_POINT
  inline
#endif
void TenseurHB::Affectation(const Mat_pleine& a)
   { int dim = Abs (this->dimension);
     #ifdef MISE_AU_POINT
      if  ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
      { cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
             << a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
             << " a une dimension differente du tenseur = "
             << dim << " !\n";
        cout << "TenseurHB::Affectation(const Mat_pleine& a)  \n";
        Sortie(1);
      };
     #endif

     // on est toujours en non symétrique
     for (int i=1;i<=dim;i++)
        for (int j=1;j<=dim;j++)
          (*this).Coor(i,j) = a(i,j);

   };

/* on vire la méthodes générale qui utilise NRC, car il y a une implantation spécifique bien plus
  performante !!

    // valeurs propre  dans le vecteur  de retour, classée par ordres "décroissants"
    // cas indique le cas de valeur propre: 
    // quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
    //                            dans ce cas les valeurs propres de retour sont nulles par défaut
    // dim = 1, cas=1 pour une valeur propre; 
    // dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
    // dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
    //         , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
    //         , cas = 2 si val(1)=val(2) ( val(1), et val(3)  dans les 2 premières composantes du retour)
    //         , cas = 3 si val(2)=val(3) ( val(1), et val(2)  dans les 2 premières composantes du retour)
#ifndef MISE_AU_POINT
  inline 
#endif
Coordonnee TenseurHB::ValPropre(int& cas) const 
  { int dim = this->dimension;
    Mat_pleine mat(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        mat(i,j) = (*this)(i,j);
    int caas;            
    Coordonnee vec =  VectValPropre(mat,caas);
    // classement des valeurs propres
    // cas ou dim == 1 on ne fait rien
    Tableau <int> ii(dim); // indicateur de l'ordre final
    switch (dim)
     { case 1: ii(1)=1; break;
       case 2: 
       { if (vec(2) > vec(1))  { ii(1)=2;ii(2)=1;}
         else                  { ii(1)=1;ii(2)=2;};
         break;
        }
       case 3:
       { if (vec(1) >= vec(2)) 
           { if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         else  
           { if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         break;
         }   
      };    
    // maintenant on met dans l'ordre
    Vecteur vec_bis(vec);
    for (int i=1;i<=dim;i++)  vec(i)=vec_bis(ii(i));     
    // définition du cas de valeurs propres
    if ( difftrespetit(vec(1),vec(2)) && difftrespetit(vec(2),vec(3))) {cas = 1;}
    else if ( !(difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 0;}
    else if ( (difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 2;}
    else if ( !(difftrespetit(vec(1),vec(2))) && difftrespetit(vec(2),vec(3))) {cas = 3;}
    // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
    if (caas <= 0) cas=-1;
    // retour de vec et de cas
    return vec;
   };      
*/
// idem met en retour la matrice mat contiend par colonne les vecteurs propre
// elle doit avoir la dimension du tenseur
// les vecteurs propre sont exprime dans le repere naturel  
#ifndef MISE_AU_POINT
  inline 
#endif
Coordonnee TenseurHB::ValPropre_int(int& cas, Mat_pleine& mat) const
  { int dim = this->dimension;
  
//---- a priori c'est méthode n'est plus utilisée par les tenseurs donc à virer par la suite ------
  
    #ifdef MISE_AU_POINT	 	 
	  if  ((mat.Nb_ligne() != dim) || (mat.Nb_colonne() != dim))
			{ cout << "\nErreur : la matrice en parametre doit etre de dimension " 
			       << dim << " !\n";
			  cout << "TenseurBH::ValPropre(Mat_pleine& mat)  \n";
			  Sortie(1);
			};
    #endif
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        mat(i,j) = (*this)(i,j);
    int caas; 
    Vecteur vec(MathUtil2::VectValPropre(mat,caas)); // calcul des vecteurs et valeurs propres
    // classement des valeurs propres
    // cas ou dim == 1 on ne fait rien
    Tableau <int> ii(dim); // indicateur de l'ordre final
    switch (dim)
     { case 1: ii(1)=1; break;
       case 2: 
       { if (vec(2) > vec(1))  { ii(1)=2;ii(2)=1;}
         else                  { ii(1)=1;ii(2)=2;};
         break;
        }
       case 3:
       { if (vec(1) >= vec(2)) // d'abord le plus grand
           { if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         else  
           { if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
             else                  { ii(1)=3; ii(2)=1; ii(3)=2;};
            }
         // puis les deux derniers
         if (vec(ii(2)) < vec(ii(3))) {int ix=ii(2); ii(2)=ii(3);ii(3)=ix;};   
         break;
         }   
      };    
    // maintenant on met dans l'ordre
    Mat_pleine mat_bis(mat);
    Coordonnee coor(vec.Taille());  
    for (int i=1;i<=dim;i++)  
      {coor(i)=vec(ii(i)); mat.RemplaceColonneSpe(i,mat_bis.Colonne (ii(i)));};     
    // définition du cas de valeurs propres
   
    // définition du cas de valeurs propres
    switch (dim)
     { case 1:  break; // il n'y a rien à faire
       case 2: 
        { if (coor(2) > coor(1))
           { double x = coor(1); coor(1) = coor(2); coor(2)=x;};
          // def du cas
          if ( difftrespetit(coor(1),coor(2))) { cas = 1;} else {cas = 0;};
          break;
        }
       case 3:
        { if ( difftrespetit(coor(1),coor(2)) && difftrespetit(coor(2),coor(3))) {cas = 1;}
          else if ( !(difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 0;}
          else if ( (difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 2;}
          else if ( !(difftrespetit(coor(1),coor(2))) && difftrespetit(coor(2),coor(3))) {cas = 3;}
          else
            {cout << "\n classement des valeurs propres : non realise !! "
                  << "\n valeurs propres: coor(1)= "<<coor(1) << ", coor(2)= "<<coor(2) << ", coor(3)= "<<coor(3)
                  << "\n TenseurHB::ValPropre_int(... "<<endl;
            };
          break;
        };
       default: ; // sinon rien (ne devrait pas arriver !
     };
   
    // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
    if (caas <= 0) cas=-1;
    // retour de coor et de cas
    return coor;
   };
   

    // changement de base (cf. théorie) : la matrice beta est telle que:
    // gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    // et la matrice gamma telle que:
    // gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    // gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    //   c-a-d= gp^i = gamma^i_j * g^j
    // rappel des différentes relations entre beta et gamma
    // [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    // [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    // formule de changement de base
    // [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
    //     beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
    //     gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne

#ifndef MISE_AU_POINT
  inline 
#endif
void TenseurHB::ChBase( const Mat_pleine& beta,const Mat_pleine& gamma)
  { // marche également si le tenseur n'est pas symétrique
    int dim = Abs (this->dimension); 
    #ifdef MISE_AU_POINT	 	 
	  if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
			{ cout << "\nErreur : la nouvelle base doit etre de dimension " 
			       << dim << " !\n";
			  cout << "TenseurHB::ChBase(Mat_pleine& beta)  \n";
			  Sortie(1);
			};
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res(i,j) = (*this)(i,j); 
    // calcul de la matrice inverse transposée
//    Mat_pleine gamma = (beta.Inverse()).Transpose();    // calcul de l'inverse de beta
    // calcul de la transformation
//    res = (beta.Transpose() * res) * gamma;  ===> vieux calcul faux
    // [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
    res = (gamma * res) * beta.Transpose();
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        this->Coor(i,j) = res(i,j) ;
  };                

// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
//  this       : le tenseur
// var_tensBB : en entrée: la variation du tenseur dans les bases initiales qu'on appelle g_i et g^i
// var_tensHB : en sortie: la variation du tenseur dans les bases finales qu'on appelle gp_i et gp^j
// beta       : en entrée gpB(i) = beta(i,j) * gB(j)
// var_beta   : en entrée : la variation de beta
// gamma       : en entrée gpH(i) = gamma(i,j) * gH(j)
// var_gamma   : en entrée : la variation de gamma
// [Ap^k_l] = [gamma] * [A^i^_j] * [beta]^T
#ifndef MISE_AU_POINT
  inline
#endif
void TenseurHB::Var_tenseur_dans_nouvelle_base
          (const Mat_pleine& beta,TenseurBH& var_tensHB, const Mat_pleine& var_beta
           ,const Mat_pleine& gamma,const Mat_pleine& var_gamma)
 {  int dim = Abs (this->dimension);
    #ifdef MISE_AU_POINT
    if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
     { cout << "\nErreur : la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurHB::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    if  ((var_beta.Nb_ligne() != dim) || (var_beta.Nb_colonne() != dim))
     { cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
            << dim << " !\n";
       cout << "TenseurHB::Var_tenseur_dans_nouvelle_base(..  \n";
       Sortie(1);
     };
    #endif
    Mat_pleine res(dim,dim);  // matrice intermediaire
    Mat_pleine res1(dim,dim);  // matrice intermediaire
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = (*this)(i,j);
    // [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
    // donc d[Ap^k_l] = d[gamma] * [A^i_j] * [beta]^T
    //                 + [gamma] * [A^i_j] * d[beta]^T
    //                 + [gamma] * d[A^i_j] * [beta]^T
    res = var_gamma * res1 * beta.Transpose();
    res += gamma * res1 * var_beta.Transpose();
    // récup des variations
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        res1(i,j) = var_tensHB(i,j);
    // ajout de la ligne restante: res maintenant correspond à d[A^i_j]
    res += gamma * res1 * beta.Transpose();
  
    // retour au tenseur
    for (int i=1;i<=dim;i++)
      for (int j=1;j<=dim;j++)
        var_tensHB.Coor(i,j) = res(i,j) ;
 };


// calcul des composantes du tenseur dans la base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHB & TenseurHB::BaseAbsolue(TenseurHB & A,const BaseH & giH,const BaseB & giB) const
    {    return produit3(A,(*this),giH,giB);      
      };
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurBH & TenseurHB::BaseAbsolue(TenseurBH & A,const BaseH & giH,const BaseB & giB) const
    {    return produit3(A,(*this),giH,giB);      
      };    
// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes ,  retour d'une reference sur A
#ifndef MISE_AU_POINT
  inline 
#endif
TenseurHB & TenseurHB::Baselocale(TenseurHB& A,const BaseB & giB,const BaseH & giH) const
 {  return produit5(A,(*this),giH,giB);
  };
  
// ----------------- fonction externe globale -----------------------
   
// liberation des tenseurs intermediaires
#ifndef MISE_AU_POINT
  inline 
#endif
void LibereTenseur() 
   {  // appel du membre static  Libere des differents tenseurs
     LesMaillonsHH::Libere();
     LesMaillonsBB::Libere();
     LesMaillonsHB::Libere();
     LesMaillonsBH::Libere();
   };  
                
                
#endif  // fin pour l'insertion dans le .h