rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/tenseurs_mai99/Tenseur/Tenseur.cc

1491 lines
58 KiB
C++
Raw Normal View History

ajout des fichiers .h et .cc du répertoire tenseurs_mai99
2021-09-07 09:39:21 +02:00
// 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.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
ajout des fichiers .h et .cc du répertoire tenseurs_mai99
2021-09-07 09:39:21 +02:00
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
//#include "Debug.h"
# include "Tenseur.h"
#include "MathUtil.h"
#include "MathUtil2.h"
#include "ParaGlob.h"
#ifndef Tenseur_H_deja_inclus
//-----------------------------------------------------------------------
// def de templates qui servent à diminuer les temps de calcul pour
// les méthodes de calcul des tenseurs dans la base absolue
//-----------------------------------------------------------------------
// 1) pour le passage repère local en repère absolu
// T1 et T2 sont des tenseurs, Base représente une base
// pour des tenseurs T1 et T2 symétriques HH ou BB
template <class T1,class T2,class Base>
T1& produit1(T1& A,const T2& Nous, const Base & gi)
{ int dim = abs(A.Dimension());// dimension du tenseur résultat
int dimint = abs(Nous.Dimension());// dimension du tenseur en local
for (int al=1;al<= dim; al++)
for (int be=1; be<= al; be++)
{ A.Coor(al,be) = 0.;
for (int i=1; i<= dimint; i++)
for (int j=1; j<= dimint; j++)
A.Coor(al,be) += Nous(i,j)*gi(i)(al)*gi(j)(be);
}
return A;
};
// pour un tenseurs T1 non symétrique (HB ou BH) et T2 symétriques HH ou BB
template <class T1,class T2,class Base>
T1& produit2(T1& A,const T2& Nous, const Base & gi)
{ int dim = abs(A.Dimension());
int dimint = abs(Nous.Dimension());
for (int al=1;al<= dim; al++)
for (int be=1; be<= dim; be++)
{ A.Coor(al,be) = 0.;
for (int i=1; i<= dimint; i++)
for (int j=1; j<= dimint; j++)
A.Coor(al,be) += Nous(i,j)*gi(i)(al)*gi(j)(be);
}
return A;
};
// pour des tenseurs non symétriques en variance T1 et T2
// il y a deux bases, une de chaque variance
template <class T1,class T2,class Base2,class Base1>
T1 & produit3(T1 & A,const T2& Nous, const Base2 & gi2, const Base1 & gi1)
{ int dim = abs(A.Dimension());
int dimint = abs(Nous.Dimension());
for (int i=1;i<= dim; i++)
for (int j=1; j<= dim; j++)
{ A.Coor(i,j) = 0.;
for (int al=1; al<= dimint; al++)
for (int be=1; be<= dimint; be++)
A.Coor(i,j) += Nous(al,be)*gi1(al)(i)*gi2(be)(j);
}
return A;
};
// 2) pour le passage repère local en repère absolu
// T1 et T2 sont des tenseurs, Base représente une base
// pour des tenseurs T1 et T2 symétriques HH ou BB
template <class T1,class T2,class Base>
T1& produit4(T1& A,const T2& Nous, const Base & gi)
{ int dim = abs(A.Dimension());
int dimint = abs(Nous.Dimension());
for (int i=1;i<= dim; i++)
for (int j=1; j<= i; j++)
{ A.Coor(i,j) = 0.;
for (int al=1; al<= dimint; al++)
for (int be=1; be<= dimint; be++)
A.Coor(i,j) += Nous(al,be)*gi(i)(al)*gi(j)(be);
}
return A;
};
// la même chose que produit4, mais lorsque le tenseur A est de stockage
// non symétrique, il faut donc à ce moment là calculer A(i,j) et A(j,i)
template <class T1,class T2,class Base>
T1& produit6(T1& A,const T2& Nous, const Base & gi)
{ int dim = abs(A.Dimension());
int dimint = abs(Nous.Dimension());
for (int i=1;i<= dim; i++)
for (int j=1; j<= dim; j++)
{ A.Coor(i,j) = 0.;
for (int al=1; al<= dimint; al++)
for (int be=1; be<= dimint; be++)
A.Coor(i,j) += Nous(al,be)*gi(i)(al)*gi(j)(be);
}
return A;
};
// pour des tenseurs non symétriques en variance T1 et T2
// il y a deux bases, une de chaque variance
template <class T1,class T2,class Base2,class Base1>
T1 & produit5(T1 & A,const T2& Nous, const Base2 & gi2, const Base1 & gi1)
{ int dim = abs(A.Dimension());
int dimint = abs(Nous.Dimension());
for (int i=1;i<= dim; i++)
for (int j=1; j<= dim; j++)
{ A.Coor(i,j) = 0.;
for (int al=1; al<= dimint; al++)
for (int be=1; be<= dimint; be++)
// A(i,j) += Nous(al,be)*gi2(al)(i)*gi1(be)(j);
A.Coor(i,j) += Nous(al,be)*gi2(i)(al)*gi1(j)(be);
}
return A;
};
//----------- fin de la définition de template --------------------------
// produit contracte avec un vecteur
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeH operator * (const CoordonneeH & v ,const TenseurBH & t)
{ return (t.Transpose() * v) ; };
// ------------------tenseur HH ----------------------------------
//PtTenseurHH * LesMaillonsHH::maille = NULL; // initialisation de la première maille
#ifdef MISE_AU_POINT
int LesMaillonsHH::nbmailHH =0; // initialisation du nombre de tenseur intermediaire
#endif
// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurHH
{ public :
PtTenseurHH * t1; // adresse du maillon precedent
TenseurHH * sauvePtr; // sauvegarde du pointeur de tenseur
PtTenseurHH (const PtTenseurHH * x1,const TenseurHH * x2)
// ici on ruse pour endormir le compilateur avec un cast
// par la suite on utilisera t1 et x2 pour effacer le tenseur !!
// ce qui n'est pas possible avec const !! normalement
{ t1 = (PtTenseurHH *) x1;sauvePtr = (TenseurHH *) x2;};
~PtTenseurHH ()
{ delete sauvePtr;};
};
// nouveau maillon
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsHH::NouveauMaillon(const TenseurHH * ptr)
{ PtTenseurHH * poi = maille;
maille = new PtTenseurHH ( poi,ptr);
#ifdef MISE_AU_POINT
nbmailHH++;
#endif
};
// liberation de l'espace
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsHH::Libere()
{ while (maille != NULL)
{ PtTenseurHH * poi = maille->t1; // recup de l'adresse precedente
delete maille; // liberation de l'espace de tenseur
maille = poi;
};
#ifdef MISE_AU_POINT
nbmailHH = 0;
#endif
};
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHH::Message(int dim, string mes) const
{ // on sort le message en fonction du niveau de commentaires
if (ParaGlob::NiveauImpression() != 0)
{ cout << " \n erreur de calcul sur les tenseurs HH ***** ";
cout << " dimension de l'argument = " << dim ;
cout << '\n' << mes << endl;
};
// dans le cas où le niveau est supérieur à une valeur on s'arrête
if (ParaGlob::NiveauImpression() > 5)
Sortie(1);
};
// retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
Mat_pleine& TenseurHH::Matrice_composante(Mat_pleine& res)const
{ int dim = Abs (this->dimension);
if ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
res.Initialise(dim,dim,0.);
if (this->dimension < 0)
// cas non symétrique
{for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
}
else // cas symétrique
{for (int i=1;i<=dim;i++) // le résultat peut ne pas être symétrique
for (int j=1;j<=i;j++)
res(j,i) = res(i,j) = (*this)(i,j);
};
return res;
};
// opération inverse: affectation en fonction d'une matrice
// donc sans vérif de variance !
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHH::Affectation(const Mat_pleine& a)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
{ cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
<< a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
<< " a une dimension differente du tenseur = "
<< dim << " !\n";
cout << "TenseurHH::Affectation(const Mat_pleine& a) \n";
Sortie(1);
};
#endif
if (this->dimension < 0)
// cas non symétrique
{for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
(*this).Coor(i,j) = a(i,j);
}
else // cas symétrique
{
#ifdef MISE_AU_POINT
if (!a.Symetrie())
{cout << "\n attention la matrice est non symetrique et on veut l'affecter a"
<< " un tenseur symetrique !";
a.Affiche();
cout << "TenseurHH::Affectation(const Mat_pleine& a) \n";
};
#endif
for (int i=1;i<=dim;i++) // le résultat peut ne pas être symétrique
for (int j=1;j<=i;j++) // à un chouat près
(*this).Coor(i,j) = 0.5* (a(j,i) + a(i,j));
};
};
// changement de base (cf. théorie) : la matrice beta est telle que:
// par défaut inverse = false : (c'est l'utilisation historique)
// beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
// gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// et on a la matrice gamma telle que:
// gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
// c-a-d= gp^i = gamma^i_j * g^j
// rappel des différentes relations entre beta et gamma
// [beta]^{-1} = [gamma]^T ; [beta]^{-1T} = [gamma]
// [beta] = [gamma]^{-1T} ; [beta]^{T} = [gamma]^{-1}
// changement de base pour de deux fois contravariants:
// [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
// donc le changement de base s'effectue directement à l'aide de gamma
// si inverse = true:
// beta(i,j) représente les coordonnées de l'ancienne base gB dans la nouvelle gpB
// gB(i) = beta(i,j) * gpB(j), i indice de ligne, j indice de colonne
// la formule de changement de base est alors:
// [Ap^kl] = [beta]^T * [A^ij] * [beta]
//
// nécessite d'inverser la matrice gamma pour calculer beta
//
//
//
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHH::ChBase(const Mat_pleine& gamma,bool inverse)
{ // marche également si le tenseur n'est pas symétrique
int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurHH::ChBase(Mat_pleine& gamma) \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++) // le résultat peut ne pas être symétrique
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
if (inverse)
{// calcul de la matrice inverse transposée
Mat_pleine beta = (gamma.Inverse()).Transpose();
////--- debug
//cout << "\n debug TenseurHH::ChBase( ";
//cout << "\n beta = "; beta.Affiche();
//cout << "\n res avant calcul = "; res.Affiche();
// [Ap^kl] = [beta]^T * [A^ij] * [beta]
res = ((beta.Transpose()) * res) * beta;
//cout << "\n res après calcul = "; res.Affiche();
//// fin debug ---
}
else
{
////--- debug
//cout << "\n debug TenseurHH::ChBase( ";
//cout << "\n beta = "; beta.Affiche();
//cout << "\n gamma = "; gamma.Affiche();
//cout << "\n res avant calcul = "; res.Affiche();
// [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
// calcul de la transformation
// res = ((gamma.Transpose()) * res) * gamma; ==> non c'est faux voir théorie
res = (gamma * res) * gamma.Transpose();
//cout << "\n res après calcul = "; res.Affiche();
//// fin debug ---
};
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
this->Coor(i,j) = res(i,j) ;
};
// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
// this : le tenseur
// var_tensHH : en entrée: la variation du tenseur this dans la base initiale qu'on appelle g_i
// var_tensHH : en sortie: la variation du tenseur dans la base finale qu'on appelle gp_i
// gamma : en entrée gpH(i) = gamma(i,j) * gH(j)
// var_gamma : en entrée : la variation de gamma
// [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHH::Var_tenseur_dans_nouvelle_base
(const Mat_pleine& gamma,TenseurHH& var_tensHH, const Mat_pleine& var_gamma)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurHH::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
if ((var_gamma.Nb_ligne() != dim) || (var_gamma.Nb_colonne() != dim))
{ cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurHH::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
Mat_pleine res1(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j)= (*this)(i,j);
// [Ap^kl] = [gamma] * [A^ij] * [gamma]^T
// donc d[Ap^kl] = d[gamma] * [A_^ij] * [gamma]^T
// + [gamma] * d[A^ij] * [gamma]^T
// + [gamma] * [Aij] * d[gamma]^T
res = var_gamma * res1 * gamma.Transpose();
res += gamma * res1 * var_gamma.Transpose();
// récup des variations
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = var_tensHH(i,j);
// ajout de la ligne restante: res maintenant correspond à d[A_ij]
res += gamma * res1 * gamma.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
var_tensHH.Coor(i,j) = res(i,j) ;
};
// calcul des composantes du tenseur dans la base absolue
// la variance du résultat peut-être quelconque d'où quatre possibilités
// en fonction de l'argument A. Dans tous les cas les composantes sont identiques
// car la sortie est en absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes, retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & TenseurHH::BaseAbsolue(TenseurHH & A,const BaseB & gi) const
{ // on traite en fonction également de la symétrie de stockage ou pas
if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
else { return produit2(A,(*this),gi);};
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & TenseurHH::BaseAbsolue(TenseurBB & A,const BaseB & gi) const
{ // on traite en fonction également de la symétrie de stockage ou pas
if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
else { return produit2(A,(*this),gi);};
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & TenseurHH::BaseAbsolue(TenseurBH & A,const BaseB & gi) const
{ return produit2(A,(*this),gi);
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & TenseurHH::BaseAbsolue(TenseurHB& A,const BaseB & gi) const
{ return produit2(A,(*this),gi);
};
// calcul des composantes locales du tenseur considéré s'exprimé dans une base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes , retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & TenseurHH::Baselocale(TenseurHH & A,const BaseH & gi) const
{ // on traite en fonction également de la symétrie de stockage ou pas
if ((A.Dimension() > 0) && (dimension > 0)) { return produit4(A,(*this),gi);}
else { return produit6(A,(*this),gi);};
};
// ------------------tenseur BB ----------------------------------
//PtTenseurBB * LesMaillonsBB::maille = NULL; // initialisation de la première maille
#ifdef MISE_AU_POINT
int LesMaillonsBB::nbmailBB =0; // initialisation du nombre de tenseur intermediaire
#endif
// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurBB
{ public :
PtTenseurBB * t1; // adresse du maillon precedent
TenseurBB * sauvePtr; // sauvegarde du pointeur de tenseur
PtTenseurBB (const PtTenseurBB * x1,const TenseurBB * x2)
// ici on ruse pour endormir le compilateur avec un cast
// par la suite on utilisera t1 et x2 pour effacer le tenseur !!
// ce qui n'est pas possible avec const !! normalement
{ t1 = (PtTenseurBB *) x1;sauvePtr = (TenseurBB *) x2;};
~PtTenseurBB ()
{ delete sauvePtr;};
};
// nouveau maillon
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsBB::NouveauMaillon(const TenseurBB * ptr)
{ PtTenseurBB * poi = maille;
maille = new PtTenseurBB ( poi,ptr);
#ifdef MISE_AU_POINT
nbmailBB++;
#endif
};
// liberation de l'espace
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsBB::Libere()
{ while (maille != NULL)
{ PtTenseurBB * poi = maille->t1; // recup de l'adresse precedente
delete maille; // liberation de l'espace de tenseur
maille = poi;
};
#ifdef MISE_AU_POINT
nbmailBB = 0;
#endif
};
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBB::Message(int dim, string mes) const
{ cout << " \n erreur de calcul sur les tenseurs BB ***** ";
cout << " dimension de l'argument = " << dim ;
cout << '\n' << mes << endl;
Sortie(1);
}
// retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
Mat_pleine& TenseurBB::Matrice_composante(Mat_pleine& res)const
{ int dim = Abs (this->dimension);
if ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
res.Initialise(dim,dim,0.);
if (this->dimension < 0)
// cas non symétrique
{for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
}
else // cas symétrique
{for (int i=1;i<=dim;i++) // le résultat peut ne pas être symétrique
for (int j=1;j<=i;j++)
res(j,i) = res(i,j) = (*this)(i,j);
};
return res;
};
// opération inverse: affectation en fonction d'une matrice
// donc sans vérif de variance !
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBB::Affectation(const Mat_pleine& a)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
{ cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
<< a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
<< " a une dimension differente du tenseur = "
<< dim << " !\n";
cout << "TenseurBB::Affectation(const Mat_pleine& a) \n";
Sortie(1);
};
#endif
if (this->dimension < 0)
// cas non symétrique
{for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
(*this).Coor(i,j) = a(i,j);
}
else // cas symétrique
{
#ifdef MISE_AU_POINT
if (!a.Symetrie())
{cout << "\n attention la matrice est non symetrique et on veut l'affecter a"
<< " un tenseur symetrique !";
a.Affiche();
cout << "TenseurBB::Affectation(const Mat_pleine& a) \n";
};
#endif
for (int i=1;i<=dim;i++) // le résultat peut ne pas être symétrique
for (int j=1;j<=i;j++) // à un chouat près
(*this).Coor(i,j) = 0.5* (a(j,i) + a(i,j));
};
};
// changement de base (cf. théorie) : la matrice beta est telle que:
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// et la matrice gamma telle que:
// gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
// c-a-d= gp^i = gamma^i_j * g^j
// rappel des différentes relations entre beta et gamma
// [beta]^{-1} = [gamma]^T ; [beta]^{-1T} = [gamma]
// [beta] = [gamma]^{-1T} ; [beta]^{T} = [gamma]^{-1}
// changement de base pour de deux fois covariants:
// [Ap_kl] = [beta] * [A_ij] * [beta]^T
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBB::ChBase(const Mat_pleine& beta)
{ // marche également si le tenseur n'est pas symétrique
int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBB::ChBase(Mat_pleine& beta) \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
// calcul de la matrice inverse
// on n'en a pas besoin mais comme cela le determinant est verifie
// non on vire car c'est débile de faire autant de calcul pour un déterminant
// Mat_pleine beta = beta.Inverse();
// calcul de la transformation
// res = ((beta.Transpose()) * res) * beta; ==> vieux calcul erreur
res = (beta * res) * beta.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
this->Coor(i,j) = res(i,j) ;
};
// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
// this : le tenseur
// var_tensBB : en entrée: la variation du tenseur dans la base initiale qu'on appelle g^i
// var_tensBB : en sortie: la variation du tenseur dans la base finale qu'on appelle gp^i
// beta : en entrée gpB(i) = beta(i,j) * gB(j)
// var_beta : en entrée : la variation de beta
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBB::Var_tenseur_dans_nouvelle_base
(const Mat_pleine& beta,TenseurBB& var_tensBB, const Mat_pleine& var_beta)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBB::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
if ((var_beta.Nb_ligne() != dim) || (var_beta.Nb_colonne() != dim))
{ cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBB::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
Mat_pleine res1(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = (*this)(i,j);
// [Ap_kl] = [beta] * [A_ij] * [beta]^T
// donc d[Ap_kl] = d[beta] * [A_ij] * [beta]^T
// + [beta] * d[A_ij] * [beta]^T
// + [beta] * [A_ij] * d[beta]^T
res = var_beta * res1 * beta.Transpose();
res += beta * res1 * var_beta.Transpose();
// récup des variations
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = var_tensBB(i,j);
// ajout de la ligne restante: res maintenant correspond à d[A_ij]
res += beta * res1 * beta.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
var_tensBB.Coor(i,j) = res(i,j) ;
};
// calcul des composantes du tenseur dans la base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant, retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & TenseurBB::BaseAbsolue(TenseurBB & A,const BaseH & gi) const
{ // on traite en fonction également de la symétrie de stockage ou pas
if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
else { return produit2(A,(*this),gi);};
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & TenseurBB::BaseAbsolue(TenseurHH & A,const BaseH & gi) const
{ // on traite en fonction également de la symétrie de stockage ou pas
if ((A.Dimension() > 0) && (dimension > 0)) { return produit1(A,(*this),gi);}
else { return produit2(A,(*this),gi);};
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & TenseurBB::BaseAbsolue(TenseurBH & A,const BaseH & gi) const
{ return produit2(A,(*this),gi);
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & TenseurBB::BaseAbsolue(TenseurHB & A,const BaseH & gi) const
{ return produit2(A,(*this),gi);
};
// calcul des composantes locales du tenseur considéré s'exprimé dans une base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes , retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & TenseurBB::Baselocale(TenseurBB & A,const BaseB & gi) const
{ // on traite en fonction également de la symétrie de stockage ou pas
if ((A.Dimension() > 0) && (dimension > 0)) { return produit4(A,(*this),gi);}
else { return produit6(A,(*this),gi);};
};
// ------------------tenseur BH ----------------------------------
//PtTenseurBH * LesMaillonsBH::maille = NULL; // initialisation de la première maille
#ifdef MISE_AU_POINT
int LesMaillonsBH::nbmailBH =0; // initialisation du nombre de tenseur intermediaire
#endif
// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurBH
{ public :
PtTenseurBH * t1; // adresse du maillon precedent
TenseurBH * sauvePtr; // sauvegarde du pointeur de tenseur
PtTenseurBH (const PtTenseurBH * x1,const TenseurBH * x2)
// ici on ruse pour endormir le compilateur avec un cast
// par la suite on utilisera t1 et x2 pour effacer le tenseur !!
// ce qui n'est pas possible avec const !! normalement
{ t1 = (PtTenseurBH *) x1;sauvePtr = (TenseurBH *) x2;};
~PtTenseurBH ()
{ delete sauvePtr;};
};
// nouveau maillon
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsBH::NouveauMaillon(const TenseurBH * ptr)
{ PtTenseurBH * poi = maille;
maille = new PtTenseurBH ( poi,ptr);
#ifdef MISE_AU_POINT
nbmailBH++;
#endif
};
// liberation de l'espace
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsBH::Libere()
{ while (maille != NULL)
{ PtTenseurBH * poi = maille->t1; // recup de l'adresse precedente
delete maille; // liberation de l'espace de tenseur
maille = poi;
};
#ifdef MISE_AU_POINT
nbmailBH = 0;
#endif
};
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBH::Message(int dim, string mes) const
{ cout << " \n erreur de calcul sur les tenseurs BH ***** ";
cout << " dimension de l'argument = " << dim ;
cout << '\n' << mes << endl;
Sortie(1);
}
// retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
Mat_pleine& TenseurBH::Matrice_composante(Mat_pleine& res)const
{ int dim = Abs (this->dimension);
if ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
res.Initialise(dim,dim,0.);
// toujours non symétrique
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
return res;
};
// opération inverse: affectation en fonction d'une matrice
// donc sans vérif de variance !
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBH::Affectation(const Mat_pleine& a)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
{ cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
<< a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
<< " a une dimension differente du tenseur = "
<< dim << " !\n";
cout << "TenseurBH::Affectation(const Mat_pleine& a) \n";
Sortie(1);
};
#endif
// on est toujours en non symétrique
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
(*this).Coor(i,j) = a(i,j);
};
/* on vire la méthode générale qui utilise NRC, car il y a une implantation spécifique bien plus
performante !!
// valeurs propre dans le vecteur de retour, classée par ordres "décroissants"
// cas indique le cas de valeur propre:
// quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
// dans ce cas les valeurs propres de retour sont nulles par défaut
// dim = 1, cas=1 pour une valeur propre;
// dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
// dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
// , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
// , cas = 2 si val(1)=val(2) ( val(1), et val(3) dans les 2 premières composantes du retour)
// , cas = 3 si val(2)=val(3) ( val(1), et val(2) dans les 2 premières composantes du retour)
#ifndef MISE_AU_POINT
inline
#endif
Coordonnee TenseurBH::ValPropre(int& cas) const
{ int dim = this->dimension;
Mat_pleine mat(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
mat(i,j) = (*this)(i,j);
int cass; Coordonnee vec = VectValPropre(mat,cass);
// classement des valeurs propres
// cas ou dim == 1 on ne fait rien
Tableau <int> ii(dim); // indicateur de l'ordre final
switch (dim)
{ case 1: ii(1)=1; break;
case 2:
{ if (vec(2) > vec(1)) { ii(1)=2;ii(2)=1;}
else { ii(1)=1;ii(2)=2;};
break;
}
case 3:
{ if (vec(1) >= vec(2))
{ if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
else { ii(1)=3; ii(2)=1; ii(3)=2;};
}
else
{ if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
else { ii(1)=3; ii(2)=1; ii(3)=2;};
}
break;
}
};
// maintenant on met dans l'ordre
Vecteur vec_bis(vec);
for (int i=1;i<=dim;i++) vec(i)=vec_bis(ii(i));
// définition du cas de valeurs propres
if ( difftrespetit(vec(1),vec(2)) && difftrespetit(vec(2),vec(3))) {cas = 1;}
else if ( !(difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 0;}
else if ( (difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 2;}
else if ( !(difftrespetit(vec(1),vec(2))) && difftrespetit(vec(2),vec(3))) {cas = 3;}
// s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
if (cass <= 0) cas=-1;
// retour de vec et de cas
return vec;
};
*/
mise à jour 7.010
2023-05-03 17:23:49 +02:00
//// idem met en retour la matrice mat contiend par colonne les vecteurs propre
//// elle doit avoir la dimension du tenseur
//// les vecteurs propre sont exprime dans le repere dual (contrairement au HB)
// // valeurs propre dans le vecteur de retour, classée par ordres "décroissants"
// // cas indique le cas de valeur propre:
// // quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
// // dans ce cas les valeurs propres de retour sont nulles par défaut
// // dim = 1, cas=1 pour une valeur propre;
// // dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
// // dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
// // , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
// // , cas = 2 si val(1)=val(2) ( val(1), et val(3) dans les 2 premières composantes du retour)
// // , cas = 3 si val(2)=val(3) ( val(1), et val(2) dans les 2 premières composantes du retour)
// // les vecteurs propre sont exprime dans le repere dual (contrairement au HB) pour les tenseurs dim 3
// // pour dim=2:le premier vecteur propre est exprime dans le repere dual
// // le second vecteur propre est exprimé dans le repère naturel
// // pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
//
//#ifndef MISE_AU_POINT
// inline
//#endif
//Coordonnee TenseurBH::ValPropre_int(int& cas, Mat_pleine& mat) const
// { int dim = this->dimension;
//
////---- a priori c'est méthode n'est plus utilisée par les tenseurs donc à virer par la suite ------
//
// #ifdef MISE_AU_POINT
// if ((mat.Nb_ligne() != dim) || (mat.Nb_colonne() != dim))
// { cout << "\nErreur : la matrice en parametre doit etre de dimension "
// << dim << " !\n";
// cout << "TenseurBH::ValPropre(Mat_pleine& mat) \n";
// Sortie(1);
// };
// #endif
// for (int i=1;i<=dim;i++)
// for (int j=1;j<=dim;j++)
// mat(i,j) = (*this)(i,j);
// int cass;
// Vecteur vec(MathUtil2::VectValPropre(mat,cass)); // calcul des vecteurs et valeurs propres
// // classement des valeurs propres
// // cas ou dim == 1 on ne fait rien
// Tableau <int> ii(dim); // indicateur de l'ordre final
// switch (dim)
// { case 1: ii(1)=1; break;
// case 2:
// { if (vec(2) > vec(1)) { ii(1)=2;ii(2)=1;}
// else { ii(1)=1;ii(2)=2;};
// break;
// }
// case 3:
// { if (vec(1) >= vec(2)) // d'abord le plus grand
// { if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
// else { ii(1)=3; ii(2)=1; ii(3)=2;};
// }
// else
// { if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
// else { ii(1)=3; ii(2)=1; ii(3)=2;};
// }
// // puis les deux derniers
// if (vec(ii(2)) < vec(ii(3))) {int ix=ii(2); ii(2)=ii(3);ii(3)=ix;};
// break;
// }
// };
// // maintenant on met dans l'ordre
// Mat_pleine mat_bis(mat);
// Coordonnee coor(vec.Taille());
// for (int i=1;i<=dim;i++)
// {coor(i)=vec(ii(i)); mat.RemplaceColonneSpe(i,mat_bis.Colonne (ii(i)));};
//
// // définition du cas de valeurs propres
// switch (dim)
// { case 1: break; // il n'y a rien à faire
// case 2:
// { if (coor(2) > coor(1))
// { double x = coor(1); coor(1) = coor(2); coor(2)=x;};
// // def du cas
// if ( difftrespetit(coor(1),coor(2))) { cas = 1;} else {cas = 0;};
// break;
// }
// case 3:
// { if ( difftrespetit(coor(1),coor(2)) && difftrespetit(coor(2),coor(3))) {cas = 1;}
// else if ( !(difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 0;}
// else if ( (difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 2;}
// else if ( !(difftrespetit(coor(1),coor(2))) && difftrespetit(coor(2),coor(3))) {cas = 3;}
// {cout << "\n classement des valeurs propres : non realise !! "
// << "\n valeurs propres: coor(1)= "<<coor(1) << ", coor(2)= "<<coor(2) << ", coor(3)= "<<coor(3)
// << "\n TenseurBH::ValPropre_int(... "<<endl;
// };
// break;
// };
// default: ; // sinon rien (ne devrait pas arriver !
// };
//
// // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
// if (cass <= 0) cas=-1;
// // retour de vec et de cas
// return coor;
// };
ajout des fichiers .h et .cc du répertoire tenseurs_mai99
2021-09-07 09:39:21 +02:00
// changement de base (cf. théorie) : la matrice beta est telle que:
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// et la matrice gamma telle que:
// gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
// c-a-d= gp^i = gamma^i_j * g^j
// rappel des différentes relations entre beta et gamma
// [beta]^{-1} = [gamma]^T ; [beta]^{-1T} = [gamma]
// [beta] = [gamma]^{-1T} ; [beta]^{T} = [gamma]^{-1}
// formule de changement de base
// [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
// beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
// gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBH::ChBase( const Mat_pleine& beta,const Mat_pleine& gamma)
{ // marche également si le tenseur n'est pas symétrique
int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBH::ChBase(Mat_pleine& beta... \n";
Sortie(1);
};
if ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBH::ChBase(Mat_pleine& beta, Mat_pleine& gamma) \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
// calcul de la matrice inverse transposée
// Mat_pleine gamma = (beta.Inverse()).Transpose();
// calcul de la transformation
// res = (gamma.Transpose() * res) * beta; ===> vieux calcul erreur voir théorie
// [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
res = (beta * res) * gamma.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
this->Coor(i,j) = res(i,j) ;
};
// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
// this : le tenseur
// var_tensBB : en entrée: la variation du tenseur dans les bases initiales qu'on appelle g_i et g^i
// var_tensBH : en sortie: la variation du tenseur dans les bases finales qu'on appelle gp_i et gp^j
// beta : en entrée gpB(i) = beta(i,j) * gB(j)
// var_beta : en entrée : la variation de beta
// gamma : en entrée gpH(i) = gamma(i,j) * gH(j)
// var_gamma : en entrée : la variation de gamma
// [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
#ifndef MISE_AU_POINT
inline
#endif
void TenseurBH::Var_tenseur_dans_nouvelle_base
(const Mat_pleine& beta,TenseurBH& var_tensBH, const Mat_pleine& var_beta
,const Mat_pleine& gamma,const Mat_pleine& var_gamma)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBH::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
if ((var_beta.Nb_ligne() != dim) || (var_beta.Nb_colonne() != dim))
{ cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurBH::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
Mat_pleine res1(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = (*this)(i,j);
// [Ap_k^l] = [beta] * [A_i^j] * [gamma]^T
// donc d[Ap_k^l] = d[beta] * [A_i^j] * [gamma]^T
// + [beta] * [A_i^j] * d[gamma]^T
// + [beta] * d[A_i^j] * [gamma]^T
res = var_beta * res1 * gamma.Transpose();
res += beta * res1 * var_gamma.Transpose();
// récup des variations
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = var_tensBH(i,j);
// ajout de la ligne restante: res maintenant correspond à d[A_i^j]
res += beta * res1 * gamma.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
var_tensBH.Coor(i,j) = res(i,j) ;
};
// calcul des composantes du tenseur dans la base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant, retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & TenseurBH::BaseAbsolue(TenseurBH & A,const BaseB & giB,const BaseH & giH) const
{ return produit3(A,(*this),giB,giH);
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & TenseurBH::BaseAbsolue(TenseurHB & A,const BaseB & giB,const BaseH & giH) const
{ return produit3(A,(*this),giB,giH);
};
// calcul des composantes locales du tenseur considéré s'exprimé dans une base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes , retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & TenseurBH::Baselocale(TenseurBH & A,const BaseH & giH,const BaseB & giB) const
{ return produit5(A,(*this),giB,giH);
};
// ------------------tenseur HB ----------------------------------
//PtTenseurHB * LesMaillonsHB::maille = NULL; // initialisation de la première maille
#ifdef MISE_AU_POINT
int LesMaillonsHB::nbmailHB =0; // initialisation du nombre de tenseur intermediaire
#endif
// def d'un maillon de liste chainee pour memoriser les differents tenseurs intermediaires
class PtTenseurHB
{ public :
PtTenseurHB * t1; // adresse du maillon precedent
TenseurHB * sauvePtr; // sauvegarde du pointeur de tenseur
PtTenseurHB (const PtTenseurHB * x1,const TenseurHB * x2)
// ici on ruse pour endormir le compilateur avec un cast
// par la suite on utilisera t1 et x2 pour effacer le tenseur !!
// ce qui n'est pas possible avec const !! normalement
{ t1 = (PtTenseurHB *) x1;sauvePtr = (TenseurHB *) x2;};
~PtTenseurHB ()
{ delete sauvePtr;};
};
// nouveau maillon
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsHB::NouveauMaillon(const TenseurHB * ptr)
{ PtTenseurHB * poi = maille;
maille = new PtTenseurHB ( poi,ptr);
#ifdef MISE_AU_POINT
nbmailHB++;
#endif
};
// liberation de l'espace
#ifndef MISE_AU_POINT
inline
#endif
void LesMaillonsHB::Libere()
{ while (maille != NULL)
{ PtTenseurHB * poi = maille->t1; // recup de l'adresse precedente
delete maille; // liberation de l'espace de tenseur
maille = poi;
};
#ifdef MISE_AU_POINT
nbmailHB = 0;
#endif
};
// sortie d'un message standard
// dim = dimension du tenseur argument
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHB::Message(int dim, string mes) const
{ cout << " \n erreur de calcul sur les tenseurs HB ***** ";
cout << " dimension de l'argument = " << dim ;
cout << '\n' << mes << endl;
Sortie(1);
}
// retourne la matrice contenant les composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
Mat_pleine& TenseurHB::Matrice_composante(Mat_pleine& res)const
{ int dim = Abs (this->dimension);
if ((res.Nb_ligne() != dim) || (res.Nb_colonne() != dim))
res.Initialise(dim,dim,0.);
// toujours non symétrique
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
return res;
};
// opération inverse: affectation en fonction d'une matrice
// donc sans vérif de variance !
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHB::Affectation(const Mat_pleine& a)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((a.Nb_ligne() != dim) || (a.Nb_colonne() != dim))
{ cout << "\nErreur : la matrice d'affectation: a.Nb_ligne()= "
<< a.Nb_ligne() << " a.Nb_colonne()= " << a.Nb_colonne()
<< " a une dimension differente du tenseur = "
<< dim << " !\n";
cout << "TenseurHB::Affectation(const Mat_pleine& a) \n";
Sortie(1);
};
#endif
// on est toujours en non symétrique
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
(*this).Coor(i,j) = a(i,j);
};
/* on vire la méthodes générale qui utilise NRC, car il y a une implantation spécifique bien plus
performante !!
// valeurs propre dans le vecteur de retour, classée par ordres "décroissants"
// cas indique le cas de valeur propre:
// quelque soit la dimension: cas = -1, si l'extraction des valeurs propres n'a pas pu se faire
// dans ce cas les valeurs propres de retour sont nulles par défaut
// dim = 1, cas=1 pour une valeur propre;
// dim = 2 , cas = 1 si deux valeurs propres distinctes, cas = 0 si deux val propres identiques
// dim = 3 , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
// , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
// , cas = 2 si val(1)=val(2) ( val(1), et val(3) dans les 2 premières composantes du retour)
// , cas = 3 si val(2)=val(3) ( val(1), et val(2) dans les 2 premières composantes du retour)
#ifndef MISE_AU_POINT
inline
#endif
Coordonnee TenseurHB::ValPropre(int& cas) const
{ int dim = this->dimension;
Mat_pleine mat(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
mat(i,j) = (*this)(i,j);
int caas;
Coordonnee vec = VectValPropre(mat,caas);
// classement des valeurs propres
// cas ou dim == 1 on ne fait rien
Tableau <int> ii(dim); // indicateur de l'ordre final
switch (dim)
{ case 1: ii(1)=1; break;
case 2:
{ if (vec(2) > vec(1)) { ii(1)=2;ii(2)=1;}
else { ii(1)=1;ii(2)=2;};
break;
}
case 3:
{ if (vec(1) >= vec(2))
{ if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
else { ii(1)=3; ii(2)=1; ii(3)=2;};
}
else
{ if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
else { ii(1)=3; ii(2)=1; ii(3)=2;};
}
break;
}
};
// maintenant on met dans l'ordre
Vecteur vec_bis(vec);
for (int i=1;i<=dim;i++) vec(i)=vec_bis(ii(i));
// définition du cas de valeurs propres
if ( difftrespetit(vec(1),vec(2)) && difftrespetit(vec(2),vec(3))) {cas = 1;}
else if ( !(difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 0;}
else if ( (difftrespetit(vec(1),vec(2))) && !difftrespetit(vec(2),vec(3))) {cas = 2;}
else if ( !(difftrespetit(vec(1),vec(2))) && difftrespetit(vec(2),vec(3))) {cas = 3;}
// s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
if (caas <= 0) cas=-1;
// retour de vec et de cas
return vec;
};
*/
mise à jour 7.010
2023-05-03 17:23:49 +02:00
//// idem met en retour la matrice mat contiend par colonne les vecteurs propre
//// elle doit avoir la dimension du tenseur
//// les vecteurs propre sont exprime dans le repere naturel
//#ifndef MISE_AU_POINT
// inline
//#endif
//Coordonnee TenseurHB::ValPropre_int(int& cas, Mat_pleine& mat) const
// { int dim = this->dimension;
//
////---- a priori c'est méthode n'est plus utilisée par les tenseurs donc à virer par la suite ------
//
// #ifdef MISE_AU_POINT
// if ((mat.Nb_ligne() != dim) || (mat.Nb_colonne() != dim))
// { cout << "\nErreur : la matrice en parametre doit etre de dimension "
// << dim << " !\n";
// cout << "TenseurBH::ValPropre(Mat_pleine& mat) \n";
// Sortie(1);
// };
// #endif
// for (int i=1;i<=dim;i++)
// for (int j=1;j<=dim;j++)
// mat(i,j) = (*this)(i,j);
// int caas;
// Vecteur vec(MathUtil2::VectValPropre(mat,caas)); // calcul des vecteurs et valeurs propres
// // classement des valeurs propres
// // cas ou dim == 1 on ne fait rien
// Tableau <int> ii(dim); // indicateur de l'ordre final
// switch (dim)
// { case 1: ii(1)=1; break;
// case 2:
// { if (vec(2) > vec(1)) { ii(1)=2;ii(2)=1;}
// else { ii(1)=1;ii(2)=2;};
// break;
// }
// case 3:
// { if (vec(1) >= vec(2)) // d'abord le plus grand
// { if (vec(1) >= vec(3)) { ii(1)=1; ii(2)=2; ii(3)=3;}
// else { ii(1)=3; ii(2)=1; ii(3)=2;};
// }
// else
// { if (vec(2) >= vec(3)) { ii(1)=2; ii(2)=1; ii(3)=3;}
// else { ii(1)=3; ii(2)=1; ii(3)=2;};
// }
// // puis les deux derniers
// if (vec(ii(2)) < vec(ii(3))) {int ix=ii(2); ii(2)=ii(3);ii(3)=ix;};
// break;
// }
// };
// // maintenant on met dans l'ordre
// Mat_pleine mat_bis(mat);
// Coordonnee coor(vec.Taille());
// for (int i=1;i<=dim;i++)
// {coor(i)=vec(ii(i)); mat.RemplaceColonneSpe(i,mat_bis.Colonne (ii(i)));};
// // définition du cas de valeurs propres
//
// // définition du cas de valeurs propres
// switch (dim)
// { case 1: break; // il n'y a rien à faire
// case 2:
// { if (coor(2) > coor(1))
// { double x = coor(1); coor(1) = coor(2); coor(2)=x;};
// // def du cas
// if ( difftrespetit(coor(1),coor(2))) { cas = 1;} else {cas = 0;};
// break;
// }
// case 3:
// { if ( difftrespetit(coor(1),coor(2)) && difftrespetit(coor(2),coor(3))) {cas = 1;}
// else if ( !(difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 0;}
// else if ( (difftrespetit(coor(1),coor(2))) && !difftrespetit(coor(2),coor(3))) {cas = 2;}
// else if ( !(difftrespetit(coor(1),coor(2))) && difftrespetit(coor(2),coor(3))) {cas = 3;}
// else
// {cout << "\n classement des valeurs propres : non realise !! "
// << "\n valeurs propres: coor(1)= "<<coor(1) << ", coor(2)= "<<coor(2) << ", coor(3)= "<<coor(3)
// << "\n TenseurHB::ValPropre_int(... "<<endl;
// };
// break;
// };
// default: ; // sinon rien (ne devrait pas arriver !
// };
//
// // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
// if (caas <= 0) cas=-1;
// // retour de coor et de cas
// return coor;
// };
//
ajout des fichiers .h et .cc du répertoire tenseurs_mai99
2021-09-07 09:39:21 +02:00
// changement de base (cf. théorie) : la matrice beta est telle que:
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// et la matrice gamma telle que:
// gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
// c-a-d= gp^i = gamma^i_j * g^j
// rappel des différentes relations entre beta et gamma
// [beta]^{-1} = [gamma]^T ; [beta]^{-1T} = [gamma]
// [beta] = [gamma]^{-1T} ; [beta]^{T} = [gamma]^{-1}
// formule de changement de base
// [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
// beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
// gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHB::ChBase( const Mat_pleine& beta,const Mat_pleine& gamma)
{ // marche également si le tenseur n'est pas symétrique
int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurHB::ChBase(Mat_pleine& beta) \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res(i,j) = (*this)(i,j);
// calcul de la matrice inverse transposée
// Mat_pleine gamma = (beta.Inverse()).Transpose(); // calcul de l'inverse de beta
// calcul de la transformation
// res = (beta.Transpose() * res) * gamma; ===> vieux calcul faux
// [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
res = (gamma * res) * beta.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
this->Coor(i,j) = res(i,j) ;
};
// il s'agit ici de calculer la variation d'un tenseur dans une nouvelle base
// connaissant sa variation dans la base actuelle
// this : le tenseur
// var_tensBB : en entrée: la variation du tenseur dans les bases initiales qu'on appelle g_i et g^i
// var_tensHB : en sortie: la variation du tenseur dans les bases finales qu'on appelle gp_i et gp^j
// beta : en entrée gpB(i) = beta(i,j) * gB(j)
// var_beta : en entrée : la variation de beta
// gamma : en entrée gpH(i) = gamma(i,j) * gH(j)
// var_gamma : en entrée : la variation de gamma
// [Ap^k_l] = [gamma] * [A^i^_j] * [beta]^T
#ifndef MISE_AU_POINT
inline
#endif
void TenseurHB::Var_tenseur_dans_nouvelle_base
(const Mat_pleine& beta,TenseurBH& var_tensHB, const Mat_pleine& var_beta
,const Mat_pleine& gamma,const Mat_pleine& var_gamma)
{ int dim = Abs (this->dimension);
#ifdef MISE_AU_POINT
if ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
{ cout << "\nErreur : la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurHB::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
if ((var_beta.Nb_ligne() != dim) || (var_beta.Nb_colonne() != dim))
{ cout << "\nErreur : la variation de la nouvelle base doit etre de dimension "
<< dim << " !\n";
cout << "TenseurHB::Var_tenseur_dans_nouvelle_base(.. \n";
Sortie(1);
};
#endif
Mat_pleine res(dim,dim); // matrice intermediaire
Mat_pleine res1(dim,dim); // matrice intermediaire
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = (*this)(i,j);
// [Ap^k_l] = [gamma] * [A^i_j] * [beta]^T
// donc d[Ap^k_l] = d[gamma] * [A^i_j] * [beta]^T
// + [gamma] * [A^i_j] * d[beta]^T
// + [gamma] * d[A^i_j] * [beta]^T
res = var_gamma * res1 * beta.Transpose();
res += gamma * res1 * var_beta.Transpose();
// récup des variations
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
res1(i,j) = var_tensHB(i,j);
// ajout de la ligne restante: res maintenant correspond à d[A^i_j]
res += gamma * res1 * beta.Transpose();
// retour au tenseur
for (int i=1;i<=dim;i++)
for (int j=1;j<=dim;j++)
var_tensHB.Coor(i,j) = res(i,j) ;
};
// calcul des composantes du tenseur dans la base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant, retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & TenseurHB::BaseAbsolue(TenseurHB & A,const BaseH & giH,const BaseB & giB) const
{ return produit3(A,(*this),giH,giB);
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & TenseurHB::BaseAbsolue(TenseurBH & A,const BaseH & giH,const BaseB & giB) const
{ return produit3(A,(*this),giH,giB);
};
// calcul des composantes locales du tenseur considéré s'exprimé dans une base absolue
// en argument : A -> une reference sur le tenseur résultat qui peut avoir une dimension
// différente du tenseur courant suivant que la dimension absolue et la dimension locale
// sont égales ou différentes , retour d'une reference sur A
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & TenseurHB::Baselocale(TenseurHB& A,const BaseB & giB,const BaseH & giH) const
{ return produit5(A,(*this),giH,giB);
};
// ----------------- fonction externe globale -----------------------
// liberation des tenseurs intermediaires
#ifndef MISE_AU_POINT
inline
#endif
void LibereTenseur()
{ // appel du membre static Libere des differents tenseurs
LesMaillonsHH::Libere();
LesMaillonsBB::Libere();
LesMaillonsHB::Libere();
LesMaillonsBH::Libere();
};
#endif // fin pour l'insertion dans le .h
Reference in a new issue Copy permalink