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