

// 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/>.


/***********************************************************************
*     DATE:        23/01/97                                            *
*                                                                $     *
*     AUTEUR:      G RIO   (mailto:gerardrio56@free.fr)                *
*                                                                $     *
*     PROJET:      Herezh++                                            *
*                                                                $     *
************************************************************************
*     BUT:   Definir les tenseurs de differentes composantes.          *
*     Les classes sont virtuelles pures.
*     Elles se declinent en fonction de la dimension du probleme.
*     L'objectif principal est de surcharger les differentes operations
*     classiques.
*                                                                $     *
*     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *
************************************************************************/

#ifndef Tenseur_H
#define Tenseur_H

//#include "Debug.h"
#include <iostream>
#include <stdlib.h>
#include "Sortie.h"
#include "ParaGlob.h"
#include "Mat_pleine.h"
#include "Base.h"
#include <iomanip>
#include "Algo_zero.h"
#include "Tableau2_T.h"

/** @defgroup Les_classes_tenseurs_virtuelles_ordre2 Les_classes_tenseurs_virtuelles_ordre2
*
*     BUT:   Definir les tenseurs d'ordre 2 de differentes composantes.
*     Les classes sont virtuelles pures.
*     Elles se declinent en fonction de la dimension du probleme.
*     L'objectif principal est de surcharger les differentes operations
*     classiques.
*
*   concernant le produit contracte un fois, en particulier pour les tenseurs mixtes
*   il y a contraction du 2ieme indice du premier tenseur avec le premier indice du second
*   tenseur : Aij * Bjk = C ik <-> A * B = C
*   le tenseur inverse par rapport au produit contracte est defini de la maniere suivante
*   Inverse(A) * A = Id, ainsi l'inverse d'un tenseur BH est un BH idem pour les HB
*   mais l'inverse d'un BB est un HH, et l'inverse d'un HH est un BB
*
*   NB: lorsque les tenseurs mixtes sont issues de tenseurs HH ou BB symetrique, l'ordre
*   de contraction des indices n'a pas d'importance sur le resultat !!
*
*   le produit contracte de deux tenseurs symetriques quelconques ne donne pas un tenseur
*   symetrique !!, donc par exemple la contraction d'un tenseur HB avec HH n'est pas forcement
*   symetrique. Le resultat est symetrique SEULEMENT lorsque ces operations sont effectues
*   avec le tenseur metrique.
*
* \author    Gérard Rio
* \version   1.0
* \date       23/01/97
* \brief       Définition des classes virtuelles pures de type Tenseur d'ordre 2, en coordonnées avec différentes variances.
*
*/




         class TenseurBB;   // def anticipee
         class TenseurHB;   // pour l'utilisation dans les produits contractes
         class TenseurBH;   //
         
/// @addtogroup Les_classes_tenseurs_virtuelles_ordre2
///  @{
///

//------------------------------------------------------------------
//!      TenseurHH:  cas des composantes deux fois contravariantes
//------------------------------------------------------------------
/// \author    Gérard Rio
/// \version   1.0
/// \date       23/01/97

class TenseurHH
{ friend TenseurHH & operator * (double  r, const  TenseurHH & t) 
                { return ( t*r); } ;
  /// produit contracte avec un vecteur
  friend CoordonneeH operator * ( const CoordonneeB & v ,  const TenseurHH & t) 
         { return (t*v); };
 
  public :
   ///  DESTRUCTEUR :
    virtual ~TenseurHH() {};
    
   // METHODES PUBLIQUES :
//1)   non virtuelles
    ///  retourne la dimension du tenseur
    int Dimension() const{ return dimension;};
//2)    virtuelles 
    ///  initialise toutes les composantes à val
    virtual void Inita(double val) = 0;   
//    //fonctions définissant le produit tensoriel normal de deux vecteurs
//    // *this=aH(i).bH(j) gBi \otimes gBj
//    // le résultat est a priori non symétrique
//    static TenseurHH &  Prod_tensoriel(const CoordonneeH & aH, const CoordonneeH & bH) ;
//     //fonctions définissant le produit tensoriel d'un vecteur avec lui-même
//    // *this=aH(i).aH(j) gBi \otimes gBj
//    // le résultat est  symétrique
//    static TenseurHH &  Prod_tensoriel(const CoordonneeH & aH) ;
    ///  operations +
    virtual TenseurHH & operator + ( const TenseurHH &) const  = 0;
    ///  operations +=
    virtual void operator += ( const TenseurHH &) = 0;
    ///  operations -
    virtual TenseurHH & operator - () const  = 0;  //  oppose du tenseur
    ///  operations - tens
    virtual TenseurHH & operator - ( const TenseurHH &) const  = 0;
    ///  operations -=
    virtual void  operator -= ( const TenseurHH &) = 0;
    ///  operations =
    virtual TenseurHH & operator = ( const TenseurHH &) = 0;
    ///  operations * double
    virtual TenseurHH & operator * (const double &) const  = 0  ;
    ///  operations *= double
    virtual void  operator *= ( const double &) = 0;
    ///  operations /
    virtual TenseurHH & operator / ( const double &) const  = 0;
    ///  operations +/=
    virtual void  operator /= ( const double &) = 0;
    
    ///  produit contracte avec un vecteur
    virtual CoordonneeH operator * ( const CoordonneeB & )  const =0;
    
    ///  produit contracte contracté une fois A(i,j)*B(j,k)=A.B
    ///  -> donc c'est l'indice du milieu qui est contracté
    /// avec BH
    virtual TenseurHH & operator * ( const TenseurBH &) const  = 0;
    /// idem avec BB
    virtual TenseurHB & operator * ( const TenseurBB &)  const = 0;
    ///  produit  contracté deux fois A(i,j)*B(j,i)=A..B
    ///  -> on contracte d'abord l'indice du milieu puis l'indice externe
    virtual double  operator && ( const TenseurBB &) const  = 0;
    
    ///  test
    virtual int operator == ( const TenseurHH &) const  = 0;
    /// test
    virtual int operator != ( const TenseurHH &)  const = 0;
    ///  determinant de la matrice des coordonnees
    virtual double  Det() const  = 0;     
    
    ///  calcul du tenseur inverse par rapport au produit contracte
    virtual TenseurBB & Inverse() const  = 0;
    
    ///  changement de base (cf. théorie)
    ///
    
    ///  changement de base (cf. théorie) : la matrice beta est telle que:
    ///  \n par défaut inverse = false : (c'est l'utilisation historique)
    ///  \n    beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
    ///  \n    gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
    ///  \n gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    ///  \n et on a la matrice gamma telle que:
    ///  \n gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    ///  \n gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    ///  \n   c-a-d= gp^i = gamma^i_j * g^j
    ///  \n rappel des différentes relations entre beta et gamma
    ///  \n [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    ///  \n [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    ///  \n changement de base pour de deux fois contravariants:
    ///  \n [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
    ///  \n     donc le changement de base s'effectue directement à l'aide de gamma

    ///  \n si inverse = true:
    ///  \n     beta(i,j) représente les coordonnées de l'ancienne base gB dans la nouvelle gpB
    ///  \n     gB(i) = beta(i,j) * gpB(j), i indice de ligne, j indice de colonne
    ///  \n     la formule de changement de base est alors:
    ///  \n [Ap^kl] = [beta]^T * [A^ij] * [beta]
    ///
    ///  \n   nécessite d'inverser la matrice gamma pour calculer beta
    ///
    void  ChBase( const Mat_pleine& gamma, bool inverse = false) ;
 
    ///  il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///
    ///  \n il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///  \n connaissant sa variation dans la base actuelle
    ///  \n  this       : le tenseur
    ///  \n var_tensBB  : en entrée: la variation du tenseur dans la base initiale qu'on appelle g_i
    ///  \n  var_tensHH : en sortie: la variation du tenseur dans la base finale qu'on appelle gp_i
    ///  \n gamma       : en entrée gpH(i) = gamma(i,j) * gH(j)
    ///  \n var_gamma   : en entrée : la variation de gamma
    ///  \n [Ap^kl] = [gamma] * [A_ij] * [gamma]^T
    void Var_tenseur_dans_nouvelle_base
              (const Mat_pleine& gamma,TenseurHH& var_tensHH, const Mat_pleine& var_gamma);
 
    ///  Affectation_trans_dimension
    ///
    ///  \n affectation de B dans this,  plusZero = false: les données manquantes sont inchangées,
    ///  \n plusZero = true: les données manquantes sont mises à 0
    ///  \n si au contraire la dimension de B est plus grande que *this, il y a uniquement affectation
    ///  \n des données possibles
    virtual void Affectation_trans_dimension(const TenseurHH & B,bool plusZero) = 0;
    
    /// calcul des composantes du tenseur dans la base absolue: cas HH
    ///
    ///  \n calcul des composantes du tenseur dans la base absolue
    ///  \n la variance du résultat peut-être  quelconque d'où quatre possibilités
    ///  \n en fonction de l'argument A. Dans tous les cas les composantes sont identiques
    ///  \n car la sortie est en absolue
    ///  \n en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n différente du tenseur courant suivant que la dimension absolue et la dimension locale
    ///  \n sont égales ou différentes,  retour d'une reference sur A
    
    TenseurHH & BaseAbsolue(TenseurHH & A,const BaseB & gi) const;
    /// idem: cas BB
    TenseurBB & BaseAbsolue(TenseurBB & A,const BaseB & gi) const;
    /// idem: cas BH
    TenseurBH & BaseAbsolue(TenseurBH & A,const BaseB & gi) const;
    /// idem: cas HB
    TenseurHB & BaseAbsolue(TenseurHB& A,const BaseB & gi) const;
    
    /// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///
    ///  \n calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///  \n en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n différente du tenseur courant suivant que la dimension absolue et la dimension locale
    ///  \n sont égales ou différentes ,  retour d'une reference sur A
    TenseurHH & Baselocale(TenseurHH & A,const BaseH & gi) const;
    
    ///  ATTENTION creation d'un tenseur transpose qui est supprime par Libere
    virtual TenseurHH & Transpose() const  = 0; 
 
    ///  calcul du maximum en valeur absolu des composantes du tenseur
    virtual double MaxiComposante() const = 0;        
    
    ///  ---- manipulation d'indice ---- -> création de nouveaux tenseurs
    virtual TenseurBB& Baisse2Indices() const = 0;
    ///  ---- manipulation d'indice ---- -> création de nouveaux tenseurs
    virtual TenseurBH& BaissePremierIndice() const = 0;
    ///  ---- manipulation d'indice ---- -> création de nouveaux tenseurs
    virtual TenseurHB& BaisseDernierIndice() const = 0;

//    // conversion de type
//     operator TenseurHH & (void)
//      {cout << " appel de la conversion de type\n";
//       return *this;}; 
 
    ///  Retourne la composante i,j du tenseur
    ///  acces en lecture et en ecriture
    virtual double& Coor(int i, int j) = 0;
 
    ///  retourne la matrice contenant les composantes du tenseur
    Mat_pleine& Matrice_composante(Mat_pleine& a) const;
    ///  opération inverse: affectation en fonction d'une matrice
    ///  donc sans vérif de variance !
    void Affectation(const Mat_pleine& a);

    ///  Retourne la composante i,j du tenseur
    ///  acces en lecture seule
	   virtual double  operator () (int i, int j) const = 0;
	 
	   ///  lecture et écriture de données
    virtual istream & Lecture(istream & entree) = 0;
    ///  lecture et écriture de données
    virtual ostream & Ecriture(ostream & sort) const = 0;

//  protected :
	
   // variables
     int dimension; // dimension du tenseur
     double * t;  // pointeur des données
       
  protected :
	
	///  sortie d'un message standard
	///  dim = dimension du tenseur argument
	void Message(int dim, string mes) const ;

   
 };
/// @}  // end of group

 
 /** @defgroup Les_classes_Maillons_tenseurs
 *
 *     BUT:    On gère les tenseurs intermédiaires au travers d'une classe  "LesMaillonsXX"
 *  qui stocke les pointeurs sur les tenseurs intermédiaire
 *
 * \author    Gérard Rio
 * \version   1.0
 * \date       23/01/97
 * \brief       Définition des classes qui stockent les pointeurs sur les tenseurs intermédiaire.
 *
 */

 // on gère les tenseurs intermédiaires au travers d'une classe  "LesMaillonsHH"
 // qui stocke les pointeurs sur les tenseurs intermédiaire

 class PtTenseurHH;    // def de la liste chainee des tenseurs intermediaires
    
 /// @addtogroup Les_classes_Maillons_tenseurs
 ///  @{

 /// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class LesMaillonsHH 
   { public :
	     ///liberation de la place occupee par des tenseurs crees pour les
	     /// operations intermediaires
      static void Libere();
      /// enregistrement de l'ajout d'un tenseur
      static void NouveauMaillon(const TenseurHH *);
      
      /// dernier maillon
      static PtTenseurHH * maille ;
      #ifdef MISE_AU_POINT
        /// nombre de maillon courant sauvegarde
        static int nbmailHH;
      #endif
    };
/// @}  // end of group
         
/// @addtogroup Les_classes_tenseurs_virtuelles_ordre2
///  @{///

//
//------------------------------------------------------------------
//!     TenseurBB:   cas des composantes deux fois covariantes
//------------------------------------------------------------------
      /// \author    Gérard Rio
      /// \version   1.0
      /// \date       23/01/97

class TenseurBB
{ friend TenseurBB & operator * (double  r, const  TenseurBB & t)
                { return ( t*r); } ;
  /// produit contracte avec un vecteur
  friend CoordonneeB operator * ( const CoordonneeH & v , const TenseurBB & t)
         { return (t*v); };
  public :
    /// DESTRUCTEUR :
    virtual ~TenseurBB() {};
      
    // METHODES PUBLIQUES :
//1)   non virtuelles
    ///  retourne la dimension du tenseur
    int Dimension() const{ return dimension;}; // retourne la dimension du tenseur
//2)    virtuelles        
    /// initialise toutes les composantes à val
    virtual void Inita(double val) = 0;   
    ///  operations +
    virtual TenseurBB & operator + ( const TenseurBB &)  const = 0;
    ///  operations +=
    virtual void operator += (const TenseurBB &) = 0;
    ///  operations -
    virtual TenseurBB & operator - () const  = 0;  // oppose du tenseur
    ///  operations - tens
    virtual TenseurBB & operator - ( const TenseurBB &) const  = 0;
    ///  operations -=
    virtual void operator -= ( const TenseurBB &) = 0;
    ///  operations =
    virtual TenseurBB & operator = ( const TenseurBB &) = 0;
    ///  operations * double
    virtual TenseurBB & operator * ( const double &) const = 0;
    ///  operations *= double
    virtual void operator *= ( const double &) = 0;
    ///  operations /
    virtual TenseurBB & operator / (const double &) const = 0;
    ///  operations /=
    virtual void operator /= ( const double &) = 0;
    
    /// produit contracte avec un vecteur
    virtual CoordonneeB operator * ( const CoordonneeH & ) const  =0;
    
    /// produit contracte contracté une fois A(i,j)*B(j,k)=A.B
    /// -> donc c'est l'indice du milieu qui est contracté
    virtual TenseurBB & operator * ( const TenseurHB &) const  = 0;
    /// idem en BH
    virtual TenseurBH & operator * ( const TenseurHH &) const  = 0;
    /// produit contracte contracté deux fois A(i,j)*B(j,i)=A..B
    /// -> on contracte d'abord l'indice du milieu puis l'indice externe
    virtual double  operator && ( const TenseurHH &)  const = 0;
    
    /// test
    virtual int operator == ( const TenseurBB &) const  = 0;
    /// test
    virtual int operator != ( const TenseurBB &) const  = 0;
    ///  determinant de la matrice des coordonnees
    virtual double  Det() const  = 0;     // determinant de la matrice des coordonnees
    
    ///  changement de base (cf. théorie)
    ///
    
    /// changement de base (cf. théorie) : la matrice beta est telle que:
    ///  \n  gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    ///  \n  et la matrice gamma telle que:
    ///  \n  gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    ///  \n  gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    ///  \n    c-a-d= gp^i = gamma^i_j * g^j
    ///  \n  rappel des différentes relations entre beta et gamma
    ///  \n  [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    ///  \n  [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    ///  \n  changement de base pour de deux fois covariants:
    ///  \n  [Ap_kl] = [beta] * [A_ij] * [beta]^T
    void ChBase( const Mat_pleine& beta);
 
    ///  il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///
    ///  \n  il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///  \n  connaissant sa variation dans la base actuelle
    ///  \n   this       : le tenseur
    ///  \n  var_tensBB  : en entrée: la variation du tenseur dans la base initiale qu'on appelle g^i
    ///  \n  var_tensBB : en sortie: la variation du tenseur dans la base finale qu'on appelle gp^i
    ///  \n  beta       : en entrée gpB(i) = beta(i,j) * gB(j)
    ///  \n  var_beta   : en entrée : la variation de beta
    void Var_tenseur_dans_nouvelle_base
              (const Mat_pleine& beta,TenseurBB& var_tensBB, const Mat_pleine& var_beta);

    ///  Affectation_trans_dimension
    ///
    ///  \n  affectation de B dans this,  plusZero = false: les données manquantes sont inchangées,
    ///  \n  plusZero = true: les données manquantes sont mises à 0
    ///  \n  si au contraire la dimension de B est plus grande que *this, il y a uniquement affectation
    ///  \n  des données possibles
    virtual void Affectation_trans_dimension(const TenseurBB & B,bool plusZero) = 0;
    
    /// calcul des composantes du tenseur dans la base absolue: cas BB
    ///
    ///  \n  calcul des composantes du tenseur dans la base absolue
    ///  \n  la variance du résultat peut-être  quelconque d'où quatre possibilités
    ///  \n  en fonction de l'argument A. Dans tous les cas les composantes sont identiques
    ///  \n  car la sortie est en absolue
    ///  \n  en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n  différente du tenseur courant suivant que la dimension absolue et la dimension locale
    ///  \n  sont égales ou différentes,  retour d'une reference sur A
    TenseurBB & BaseAbsolue(TenseurBB & A,const BaseH & gi) const ;
    /// idem pour HH
    TenseurHH & BaseAbsolue(TenseurHH & A,const BaseH & gi) const ;
    /// idem pour BH
    TenseurBH & BaseAbsolue(TenseurBH & A,const BaseH & gi) const ;
    /// idem pour HB
    TenseurHB & BaseAbsolue(TenseurHB & A,const BaseH & gi) const ;
    
    /// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///
    ///  \n  calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///  \n  en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n  différente du tenseur courant suivant que la dimension absolue et la dimension locale
    ///  \n  sont égales ou différentes ,  retour d'une reference sur A
    TenseurBB & Baselocale(TenseurBB & A,const BaseB & gi) const;
    
    /// ATTENTION creation d'un tenseur transpose qui est supprime par Libere
    virtual TenseurBB & Transpose()  const = 0; 
    
    /// ---- manipulation d'indice ---- -> création de nouveaux tenseurs
    virtual TenseurHH& Monte2Indices() const = 0;
    /// ---- manipulation d'indice ---- -> création de nouveaux tenseurs
    virtual TenseurHB& MontePremierIndice() const = 0;
    /// ---- manipulation d'indice ---- -> création de nouveaux tenseurs
    virtual TenseurBH& MonteDernierIndice() const = 0;

    /// calcul du maximum en valeur absolu des composantes du tenseur
    virtual double MaxiComposante() const = 0;

    /// calcul du tenseur inverse par rapport au produit contracte
    virtual TenseurHH & Inverse() const  = 0;
        
//    // conversion de type
//     operator TenseurBB & (void)
//      {cout << " appel de la conversion de type\n";
//       return *this;}; 
 
    /// Retourne la composante i,j du tenseur
    /// acces en lecture et en ecriture
    virtual double& Coor(int i, int j) = 0;
 
    /// retourne la matrice contenant les composantes du tenseur
    Mat_pleine& Matrice_composante(Mat_pleine& a) const;
    /// opération inverse: affectation en fonction d'une matrice
    /// donc sans vérif de variance !
    void Affectation(const Mat_pleine& a);

    /// Retourne la composante i,j du tenseur
    /// acces en lecture seulement
	   virtual double operator () (int i, int j) const  = 0;
	 
	   /// lecture et écriture de données
    virtual istream & Lecture(istream & entree) = 0;
    /// lecture et écriture de données
    virtual ostream & Ecriture(ostream & sort) const = 0;
 
//  protected :
	
    // variables
    int dimension; // dimension du tenseur
    double * t;  // pointeur des données
    
  protected :
	
	/// sortie d'un message standard
	/// dim = dimension du tenseur argument
	void Message(int dim, string mes)  const ;   
};
/// @}  // end of group

// on gère les tenseurs intermédiaires au travers d'une classe  "LesMaillonsBB"
// qui stocke les pointeurs sur les tenseurs intermédiaire

class PtTenseurBB;    // pour la liste chainee des tenseurs intermediaires
    
/// @addtogroup Les_classes_Maillons_tenseurs
///  @{

 /// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class LesMaillonsBB
   { public :
      ///liberation de la place occupee par des tenseurs crees pour les
      /// operations intermediaires
      static void Libere();
      /// enregistrement de l'ajout d'un tenseur
      static void NouveauMaillon(const TenseurBB *); // enregistrement de l'ajout d'un tenseur
      
      /// dernier maillon
      static PtTenseurBB * maille ;// dernier maillon
      #ifdef MISE_AU_POINT
        /// nombre de maillon courant sauvegarde
        static int nbmailBB;  // nombre de maillon courant sauvegarde
      #endif
    }; 
/// @}  // end of group
         
/// @addtogroup Les_classes_tenseurs_virtuelles_ordre2
///  @{///

//------------------------------------------------------------------
//!      TenseurBH: cas des composantes mixtes BH
//------------------------------------------------------------------
/// \author    Gérard Rio
/// \version   1.0
/// \date       23/01/97
class TenseurBH
{ friend TenseurBH & operator * (double r, const TenseurBH & t)
                { return ( t*r); } ;
 /// produit contracte avec un vecteur
  friend CoordonneeH operator * ( const CoordonneeH & v , const TenseurBH & t);
  
  public :
    /// DESTRUCTEUR :
    virtual ~TenseurBH() {};
      
    // METHODES PUBLIQUES :
//1)   non virtuelles
    ///  retourne la dimension du tenseur
    int Dimension()  const { return dimension;}; // retourne la dimension du tenseur
//2)    virtuelles        
    /// initialise toutes les composantes à val
    virtual void Inita(double val) = 0;   
    ///  operations +
    virtual TenseurBH & operator + ( const TenseurBH &) const = 0;
    ///  operations +=
    virtual void operator += ( const TenseurBH &) = 0;
    ///  operations -
    virtual TenseurBH & operator - ( const TenseurBH &) const = 0;
    ///  operations opposé
    virtual TenseurBH & operator - () const  = 0;  // oppose du tenseur
    ///  operations -=
    virtual void operator -= ( const TenseurBH &) = 0;
    ///  operations =
    virtual TenseurBH & operator = ( const TenseurBH &) = 0;
    ///  operations *
    virtual TenseurBH & operator * ( const double &) const  = 0;
    ///  operations *=
    virtual void operator *= ( const double &) = 0;
    ///  operations /
    virtual TenseurBH & operator / ( const double &) const  = 0;
    ///  operations /=
    virtual void operator /= ( const double &) = 0;
    
    /// produit contracte avec un vecteur
    virtual CoordonneeB operator * ( const CoordonneeB & )  const =0;
    
    /// produit contracte contracté une fois A(i,j)*B(j,k)=A.B
    /// -> donc c'est l'indice du milieu qui est contracté
    virtual TenseurBB & operator * ( const TenseurBB &) const  = 0;
    virtual TenseurBH & operator * ( const TenseurBH &)  const = 0; // produit une fois contracte
    /// produit contracte contracté deux fois A(i,j)*B(j,i)=A..B
    /// -> on contracte d'abord l'indice du milieu puis l'indice externe
    virtual double  operator && ( const TenseurBH &)  const = 0; // produit deux fois contracte
        
    /// test
    virtual int operator == ( const TenseurBH &) const  = 0;
    virtual int operator != ( const TenseurBH &) const  = 0;
        
    /// calcul du tenseur inverse par rapport au produit contracte
    virtual TenseurBH & Inverse() const  = 0;
    
    /// trace du tenseur ou premier invariant
    virtual double  Trace()  const = 0;
    /// second invariant = trace (A*A)
    virtual double  II() const  = 0;
    /// troisieme invariant = trace ((A*A)*A)
    virtual double  III()  const = 0;
    /// determinant de la matrice des coordonnees
    virtual double  Det()  const = 0;
    
    /// calcul des valeurs propres
    ///
    ///  \n valeurs propre  dans le vecteur  de retour, classée par ordres "décroissants"
    ///  \n cas indique le cas de valeur propre:
    ///  \n quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
    ///  \n                            dans ce cas les valeurs propres de retour sont nulles par défaut
    ///  \n dim = 1, cas=1 pour une valeur propre;
    ///  \n dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
    ///  \n dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
    ///  \n         , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
    ///  \n         , cas = 2 si val(1)=val(2) ( val(1), et val(3)  dans les 2 premières composantes du retour)
    ///  \n         , cas = 3 si val(2)=val(3) ( val(1), et val(2)  dans les 2 premières composantes du retour)
    virtual Coordonnee ValPropre(int& cas) const =0 ;

    /// calcul des valeurs propres
    ///
    ///  \n idem met en retour la matrice mat contiend par colonne les vecteurs propre
    ///  \n elle doit avoir la dimension du tenseur
    ///  \n les vecteurs propre sont exprime dans le repere dual (contrairement au HB) pour les tenseurs dim 3
    ///  \n pour dim=2:le premier vecteur propre est exprime dans le repere dual
    ///  \n            le second vecteur propre est exprimé dans le repère naturel
    ///  \n pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
    virtual Coordonnee ValPropre(int& cas, Mat_pleine& mat) const = 0 ;

    ///  \n calcul des vecteurs propres, les valeurs propres
    ///  \n étant déjà connues
    ///
    ///  \n ici il s'agit uniquement de calculer les vecteurs propres, les valeurs propres
    ///  \n étant déjà connues
    ///  \n en retour VP les vecteurs propre : doivent avoir la dimension du tenseur
    ///  \n les vecteurs propre sont exprime dans le repere naturel (pour les tenseurs dim 3
    ///  \n pour dim=2:le premier vecteur propre est exprime dans le repere naturel
    ///  \n            le second vecteur propre est exprimé dans le repère dual
    ///  \n pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
    ///  \n   en sortie cas = -1 s'il y a eu un problème, dans ce cas, V_P est quelconque
    ///  \n             sinon si tout est ok, cas est identique en sortie avec l'entrée
    virtual void VecteursPropres(const Coordonnee& Val_P,int& cas, Tableau <Coordonnee>& V_P) const = 0;
 
    ///  changement de base (cf. théorie)
    ///
    
    ///  \n changement de base (cf. théorie) : la matrice beta est telle que:
    ///  \n gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    ///  \n et la matrice gamma telle que:
    ///  \n gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    ///  \n gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    ///  \n   c-a-d= gp^i = gamma^i_j * g^j
    ///  \n rappel des différentes relations entre beta et gamma
    ///  \n [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    ///  \n [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    ///  \n formule de changement de base
    ///  \n [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
    ///  \n     beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
    ///  \n     gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
    void ChBase( const Mat_pleine& beta,const Mat_pleine& gamma);
 
    ///  il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///
    ///  \n il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///  \n connaissant sa variation dans la base actuelle
    ///  \n  this       : le tenseur
    ///  \n var_tensBB  : en entrée: la variation du tenseur dans les bases initiales qu'on appelle g_i et g^i
    ///  \n var_tensBH : en sortie: la variation du tenseur dans les bases finales qu'on appelle gp_i et gp^j
    ///  \n beta       : en entrée gpB(i) = beta(i,j) * gB(j)
    ///  \n var_beta   : en entrée : la variation de beta
    ///  \n gamma       : en entrée gpH(i) = gamma(i,j) * gH(j)
    ///  \n var_gamma   : en entrée : la variation de gamma
    ///  \n [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
    void 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);

    
    ///  Affectation_trans_dimension
    ///
    ///  \n affectation de B dans this,  plusZero = false: les données manquantes sont inchangées,
    ///  \n plusZero = true: les données manquantes sont mises à 0
    ///  \n si au contraire la dimension de B est plus grande que *this, il y a uniquement affectation
    ///  \n des données possibles
    virtual void Affectation_trans_dimension(const TenseurBH & B,bool plusZero) = 0;
    
    /// calcul des composantes du tenseur dans la base absolue: cas HH
    ///
    ///  \n calcul des composantes du tenseur dans la base absolue
    ///  \n la variance du résultat peut-être  quelconque d'où quatre possibilités
    ///  \n mais on n'en conserve que les non symétriques qui sont le cas général
    ///  \n se qui évite de symétriser accidentellement un tenseur non symétrique
    ///  \n en fonction de l'argument A.Dans tous les cas les composantes sont identiques
    ///  \n car la sortie est en absolue.
    ///  \n en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n différente du tenseur courant,  retour d'une reference sur A
    TenseurBH & BaseAbsolue(TenseurBH & A,const BaseB & giB,const BaseH & giH) const;
    /// idem cas HB
    TenseurHB & BaseAbsolue(TenseurHB & A,const BaseB & giB,const BaseH & giH) const;
    
    /// calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///
    ///  \n calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///  \n en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n différente du tenseur courant suivant que la dimension absolue et la dimension locale
    ///  \n sont égales ou différentes ,  retour d'une reference sur A
    TenseurBH & Baselocale(TenseurBH & A,const BaseH & gih,const BaseB & gib) const;

    /// ATTENTION creation d'un tenseur transpose qui est supprime par Libere
    virtual TenseurHB & Transpose() const  = 0; 
    /// permute Bas Haut, mais reste dans le même tenseur
    virtual void PermuteHautBas() =0;
    
    /// calcul du maximum en valeur absolu des composantes du tenseur
    virtual double MaxiComposante() const = 0;
    
//    // conversion de type
//     operator TenseurBH & (void)
//      { return *this;};
 
    /// Retourne la composante i,j du tenseur
    /// acces en lecture et en ecriture
    virtual double& Coor(int i, int j) = 0;
 
    /// retourne la matrice contenant les composantes du tenseur
    Mat_pleine& Matrice_composante(Mat_pleine& a) const;
    /// opération inverse: affectation en fonction d'une matrice
    /// donc sans vérif de variance !
    void Affectation(const Mat_pleine& a);

    /// Retourne la composante i,j du tenseur
    /// acces en lecture seulement
	   virtual double operator () (int i, int j) const = 0;
	 
	   /// lecture et écriture de données
    virtual istream & Lecture(istream & entree) = 0;
    /// lecture et écriture de données
    virtual ostream & Ecriture(ostream & sort) const = 0;

//  protected :
	
   // variables
    int dimension; // dimension du tenseur
    double * t;  // pointeur des données
        
  protected :
	
    static Algo_zero  alg_zero; // algo pour la recherche de zero	
    /// sortie d'un message standard
    /// dim = dimension du tenseur argument
    void Message(int dim, string mes) const ;
 
    /// méthode interne pour le calcul de vecteurs propres
    Coordonnee ValPropre_int(int& cas, Mat_pleine& mat) const;
};
/// @}  // end of group

// on gère les tenseurs intermédiaires au travers d'une classe  "LesMaillonsBH"
// qui stocke les pointeurs sur les tenseurs intermédiaire

class PtTenseurBH;    // pour la liste chainee des tenseurs intermediaires
    
/// @addtogroup Les_classes_Maillons_tenseurs
///  @{

/// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class LesMaillonsBH
   { public :
      ///liberation de la place occupee par des tenseurs crees pour les
      /// operations intermediaires
      static void Libere();
      /// enregistrement de l'ajout d'un tenseur
      static void NouveauMaillon(const TenseurBH *); // enregistrement de l'ajout d'un tenseur
      
      /// dernier maillon
      static PtTenseurBH * maille ;// dernier maillon
      #ifdef MISE_AU_POINT
        /// nombre de maillon courant sauvegarde
        static int nbmailBH;  // nombre de maillon courant sauvegarde
      #endif
    }; 
/// @}  // end of group
         
/// @addtogroup Les_classes_tenseurs_virtuelles_ordre2
///  @{///

//------------------------------------------------------------------
//!   TenseurHB:  cas des composantes mixtes HB
//------------------------------------------------------------------
/// \author    Gérard Rio
/// \version   1.0
/// \date       23/01/97
class TenseurHB
{ friend TenseurHB & operator * (double  r,const TenseurHB & t)
                { return ( t*r); } ;
 /// produit contracte avec un vecteur
  friend CoordonneeB operator * (const CoordonneeB & v ,const TenseurHB & t)
        { return (t.Transpose() * v) ; };

  public :
    /// DESTRUCTEUR :
    virtual ~TenseurHB() { };
    
    // METHODES PUBLIQUES :
//1)   non virtuelles
    ///  retourne la dimension du tenseur
    int Dimension() const { return dimension;}; // retourne la dimension du tenseur
//2)    virtuelles        
    /// initialise toutes les composantes à val
    virtual void Inita(double val) = 0;   
    ///  operations +
    virtual TenseurHB & operator + ( const TenseurHB &) const  = 0;
    ///  operations +=
    virtual void operator += ( const TenseurHB &) = 0;
    ///  operations opposé
    virtual TenseurHB & operator - () const  = 0;  // oppose du tenseur
    ///  operations -
    virtual TenseurHB & operator - ( const TenseurHB &) const  = 0;
    ///  operations -=
    virtual void operator -= ( const TenseurHB &) = 0;
    ///  operations =
    virtual TenseurHB & operator = ( const TenseurHB &) = 0;
    ///  operations *
    virtual TenseurHB & operator * ( const double &) const  = 0;
    ///  operations *=
    virtual void operator *= ( const double &) = 0;
    ///  operations /
    virtual TenseurHB & operator / ( const double &)  const = 0;
    ///  operations /=
    virtual void operator /= ( const double &) = 0;
    
    /// produit contracte avec un vecteur
    virtual CoordonneeH operator * ( const CoordonneeH & ) const  =0;
    
    /// produit contracte contracté une fois A(i,j)*B(j,k)=A.B
    /// -> donc c'est l'indice du milieu qui est contracté
    virtual TenseurHH & operator * ( const TenseurHH &) const  = 0;
    virtual TenseurHB & operator * ( const TenseurHB &) const  = 0; // produit une fois contracte
    /// produit contracte contracté deux fois A(i,j)*B(j,i)=A..B
    /// -> on contracte d'abord l'indice du milieu puis l'indice externe
    virtual double  operator && ( const TenseurHB &) const  = 0; // produit deux fois contracte
        
    /// test
    virtual int operator == ( const TenseurHB &) const = 0;
    virtual int operator != ( const TenseurHB &)  const = 0;    
    
    /// calcul du tenseur inverse par rapport au produit contracte
    virtual TenseurHB & Inverse()  const = 0;
    
    /// trace du tenseur ou premier invariant
    virtual double  Trace() const  = 0;
    /// second invariant = trace (A*A)
    virtual double  II()  const = 0;
    /// troisieme invariant = trace ((A*A)*A)
    virtual double  III() const  = 0;
    /// determinant de la matrice des coordonnees
    virtual double  Det()  const = 0;
    
    /// calcul des valeurs propres
    ///

    ///  \n valeurs propre  dans le vecteur  de retour, classée par ordres "décroissants"
    ///  \n cas indique le cas de valeur propre:
    ///  \n quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
    ///  \n                            dans ce cas les valeurs propres de retour sont nulles par défaut
    ///  \n dim = 1, cas=1 pour une valeur propre;
    ///  \n dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
    ///  \n dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
    ///  \n         , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
    ///  \n         , cas = 2 si val(1)=val(2) ( val(1), et val(3)  dans les 2 premières composantes du retour)
    ///  \n         , cas = 3 si val(2)=val(3) ( val(1), et val(2)  dans les 2 premières composantes du retour)
    virtual Coordonnee ValPropre(int& cas) const =0 ;
    
    /// calcul des valeurs propres
    ///
    ///  \n idem met en retour la matrice mat contiend par colonne les vecteurs propre
    ///  \n elle doit avoir la dimension du tenseur
    ///  \n les vecteurs propre sont exprime dans le repere naturel (pour les tenseurs dim 3
    ///  \n pour dim=2:le premier vecteur propre est exprime dans le repere naturel
    ///  \n            le second vecteur propre est exprimé dans le repère dual
    ///  \n pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
    virtual Coordonnee ValPropre(int& cas, Mat_pleine& mat) const = 0;

    
    ///  \n calcul des vecteurs propres, les valeurs propres
    ///  \n étant déjà connues
    ///
    ///  \n ici il s'agit uniquement de calculer les vecteurs propres, les valeurs propres
    ///  \n étant déjà connues
    ///  \n en retour VP les vecteurs propre : doivent avoir la dimension du tenseur
    ///  \n les vecteurs propre sont exprime dans le repere naturel (pour les tenseurs dim 3
    ///  \n pour dim=2:le premier vecteur propre est exprime dans le repere naturel
    ///  \n            le second vecteur propre est exprimé dans le repère dual
    ///  \n pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
    ///  \n   en sortie cas = -1 s'il y a eu un problème, dans ce cas, V_P est quelconque
    ///  \n             sinon si tout est ok, cas est identique en sortie avec l'entrée
    virtual void VecteursPropres(const Coordonnee& Val_P,int& cas, Tableau <Coordonnee>& V_P) const = 0;
 
//    // conversion de type
//    operator TenseurHB & (void)
//       {return *this;};        
 
 
    
    ///  changement de base (cf. théorie)
    ///
       
    ///  \n changement de base (cf. théorie) : la matrice beta est telle que:
    ///  \n gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    ///  \n et la matrice gamma telle que:
    ///  \n gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    ///  \n gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    ///  \n   c-a-d= gp^i = gamma^i_j * g^j
    ///  \n rappel des différentes relations entre beta et gamma
    ///  \n [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    ///  \n [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    ///  \n formule de changement de base
    ///  \n [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
    ///  \n     beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
    ///  \n     gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
    void ChBase( const Mat_pleine& beta,const Mat_pleine& gamma);
 
    
    ///  il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///
    ///  \n il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
    ///  \n connaissant sa variation dans la base actuelle
    ///  \n  this       : le tenseur
    ///  \n var_tensBB  : en entrée: la variation du tenseur dans les bases initiales qu'on appelle g_i et g^i
    ///  \n var_tensHB : en sortie: la variation du tenseur dans les bases finales qu'on appelle gp_i et gp^j
    ///  \n beta       : en entrée gpB(i) = beta(i,j) * gB(j)
    ///  \n var_beta   : en entrée : la variation de beta
    ///  \n gamma       : en entrée gpH(i) = gamma(i,j) * gH(j)
    ///  \n var_gamma   : en entrée : la variation de gamma
    ///  \n [Ap^k_l] = [gamma] * [A^i^_j] * [beta]^T
    void 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);

 
    ///  Affectation_trans_dimension
    ///
    ///  \n affectation de B dans this,  plusZero = false: les données manquantes sont inchangées,
    ///  \n plusZero = true: les données manquantes sont mises à 0
    ///  \n si au contraire la dimension de B est plus grande que *this, il y a uniquement affectation
    ///  \n des données possibles
    virtual void Affectation_trans_dimension(const TenseurHB & B,bool plusZero) = 0;
    
    /// calcul des composantes du tenseur dans la base absolue: cas HH
    ///
    ///  \n calcul des composantes du tenseur dans la base absolue
    ///  \n la variance du résultat peut-être  quelconque d'où quatre possibilités
    ///  \n mais on n'en conserve que les non symétriques qui sont le cas général
    ///  \n se qui évite de symétriser accidentellement un tenseur non symétrique
    ///  \n en fonction de l'argument A. Dans tous les cas les composantes sont identiques
    ///  \n car la sortie est en absolue.
    ///  \n en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n différente du tenseur courant,  retour d'une reference sur A
    TenseurHB & BaseAbsolue(TenseurHB & A,const BaseH & giH,const BaseB & giB) const;
    TenseurBH & BaseAbsolue(TenseurBH & A,const BaseH & giH,const BaseB & giB) const;
    ///  \n calcul des composantes locales du tenseur considéré s'exprimé dans une  base absolue
    ///  \n en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
    ///  \n différente du tenseur courant suivant que la dimension absolue et la dimension locale
    ///  \n sont égales ou différentes ,  retour d'une reference sur A
    TenseurHB & Baselocale(TenseurHB& A,const BaseB & gib,const BaseH & gih) const;
    
    ///  ATTENTION creation d'un tenseur transpose qui est supprime par Libere
    virtual TenseurBH & Transpose() const  = 0;
    /// permute Bas Haut, mais reste dans le même tenseur
    virtual void PermuteHautBas() =0;
    
    /// calcul du maximum en valeur absolu des composantes du tenseur
    virtual double MaxiComposante() const = 0;
   
    /// Retourne la composante i,j du tenseur
    /// acces en lecture et en ecriture
    virtual double& Coor(int i, int j) = 0;
 
    /// retourne la matrice contenant les composantes du tenseur
    Mat_pleine& Matrice_composante(Mat_pleine& a) const;
    /// opération inverse: affectation en fonction d'une matrice
    /// donc sans vérif de variance !
    void Affectation(const Mat_pleine& a);

    /// Retourne la composante i,j du tenseur
    /// acces en lecture seulement
	   virtual double operator () (int i, int j) const = 0;
	 
	   /// lecture et écriture de données
    virtual istream & Lecture(istream & entree) = 0;
    /// lecture et écriture de données
    virtual ostream & Ecriture(ostream & sort) const = 0;

//  protected :
	
     // variables
    int dimension; // dimension du tenseur
    double * t;  // pointeur des données
    
  protected :
	
    static Algo_zero  alg_zero; // algo pour la recherche de zero	
    /// sortie d'un message standard
    /// dim = dimension du tenseur argument
    void Message(int dim, string mes) const ;
 
    /// méthode interne pour le calcul de vecteurs propres
    Coordonnee ValPropre_int(int& cas, Mat_pleine& mat) const;

};
/// @}  // end of group

 // on gère les tenseurs intermédiaires au travers d'une classe  "LesMaillonsHB"
 // qui stocke les pointeurs sur les tenseurs intermédiaire

 class PtTenseurHB;    // pour  la liste chainee des tenseurs intermediaires
    
/// @addtogroup Les_classes_Maillons_tenseurs
///  @{

/// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class LesMaillonsHB
   { public :
      ///liberation de la place occupee par des tenseurs crees pour les
      /// operations intermediaires
      static void Libere();
      /// enregistrement de l'ajout d'un tenseur
      static void NouveauMaillon(const TenseurHB *); // enregistrement de l'ajout d'un tenseur
      
      /// dernier maillon
      static PtTenseurHB * maille ;// dernier maillon
      #ifdef MISE_AU_POINT
        /// nombre de maillon courant sauvegarde
        static int nbmailHB;  // nombre de maillon courant sauvegarde
      #endif
    }; 
/// @}  // end of group



//==============================================================================
// probleme de gestion de la definition de tenseurs supplementaires lors 
// d'ecriture de grandes expressions. Il est necessaire de liberer l'espace
// apres les calculs
// la fonction libereTenseur libere tout l'espace de tous les types de tenseurs
//==============================================================================

 void LibereTenseur() ;


#ifndef MISE_AU_POINT
  #include "Tenseur.cc"
  #define  Tenseur_H_deja_inclus
#endif
 
#endif  
