rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/tenseurs_mai99/Tenseur/Tenseur3-2.cc
Gérard Rio 0460624e94 mise à jour 7.010
2023-05-03 17:23:49 +02:00

1846 lines
71 KiB
C++
Raw Blame History

// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
//#include "Debug.h"
#include "Tenseur3.h"
#include "ConstMath.h"
#include "MathUtil.h"
#include "Util.h"
#include "MathUtil2.h"
#ifndef Tenseur3_H_deja_inclus
// variables globales
// initialisation dans EnteteTenseur.h , utilisé dans le progr principal
//------------------------------------------------------------------
// cas des composantes mixtes BH
//------------------------------------------------------------------
// --- gestion de changement d'index ----
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::ChangementIndex::ChangementIndex() :
idx_i(9),idx_j(9),odVect(3)
{ idx_i(1)=1;idx_i(2)=1;idx_i(3)=1; idx_j(1)=1;idx_j(2)=2;idx_j(3)=3;
idx_i(4)=2;idx_i(5)=2;idx_i(6)=2; idx_j(4)=1;idx_j(5)=2;idx_j(6)=3;
idx_i(7)=3;idx_i(8)=3;idx_i(9)=3; idx_j(7)=1;idx_j(8)=2;idx_j(9)=3;
odVect(1,1)=1;odVect(1,2)=2;odVect(1,3)=3;
odVect(2,1)=4;odVect(2,2)=5;odVect(2,3)=6;
odVect(3,1)=7;odVect(3,2)=8;odVect(3,3)=9;
};
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::Tenseur3BH() :
ipointe() // par defaut
{ dimension = 3;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = 0.;
};
// initialisation de toutes les composantes a une meme valeur
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::Tenseur3BH( const double val):
ipointe()
{ dimension = 3;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = val;
};
// initialisation avec 9 valeurs différentes
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::Tenseur3BH
( const double val1, const double val2, const double val3, // 1ere ligne
const double val4, const double val5, const double val6, // 2ieme ligne
const double val7, const double val8, const double val9) // 3ieme ligne
:
ipointe()
{ dimension = 3;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
t[0] = val1;t[1] = val2;t[2] = val3;
t[3] = val4;t[4] = val5;t[5] = val6;
t[6] = val7;t[7] = val8;t[8] = val9;
};
// constructeur a partir d'une instance non differenciee
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::Tenseur3BH ( const TenseurBH & B) :
ipointe()
{ this->dimension = B.dimension ;
#ifdef MISE_AU_POINT
if (dimension != 3)
{ cout << "\n erreur de dimension, elle devrait etre = 3 ";
cout << "\n Tenseur3BH::Tenseur3BH ( TenseurBH &) " << endl;
Sortie(1);
}
#endif
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = B.t[i];
};
// constructeur de copie
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::Tenseur3BH ( const Tenseur3BH & B):
ipointe()
{ this->dimension = B.dimension ;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = B.t[i];
};
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BH::~Tenseur3BH()
{//if(listdouble9.end() != listdouble9.begin()) // si la liste n'est pas vide
listdouble9.erase(ipointe);} ; // suppression de l'élément de la liste
// initialise toutes les composantes à val
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BH::Inita(double val)
{ t[0] =val; t[1] =val; t[2] =val;
t[3] =val; t[4] =val; t[5] =val;
t[6] =val; t[7] =val; t[8] =val;
};
// operations
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator + ( const TenseurBH & B) const
{ TenseurBH * ptr;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator + ( etc..");
#endif
ptr = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( ptr); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
(*ptr).t[i] = this->t[i] + B.t[i]; //somme des données
return (*ptr) ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BH::operator += ( const TenseurBH & B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator + ( etc..");
#endif
for (int i = 0; i<=8; i++)
this->t[i] += B.t[i];}; //somme des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator - () const
{ TenseurBH * ptr;
ptr = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( ptr); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
(*ptr).t[i] = - this->t[i]; // soustraction des données
return (*ptr) ;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator - ( const TenseurBH & B) const
{ TenseurBH * ptr;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator - ( etc..");
#endif
ptr = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( ptr); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
(*ptr).t[i] = this->t[i] - B.t[i]; // soustraction des données
return (*ptr) ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BH::operator -= ( const TenseurBH & B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator -= ( etc..");
#endif
for (int i = 0; i<=8; i++)
this->t[i] -= B.t[i];}; //soustraction des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator = ( const TenseurBH & B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator = ( etc..");
#endif
for (int i = 0; i<=8; i++)
this->t[i] = B.t[i];
LesMaillonsBH::Libere(); // destruction des tenseurs intermediaires
return *this; }; //affectation des données;
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator * ( const double & b) const
{ TenseurBH * res;
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
res->t[i] = this->t[i] * b; //multiplication des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BH::operator *= ( const double & b)
{ for (int i = 0; i<=8; i++)
this->t[i] *= b;}; //multiplication des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator / ( const double & b) const
{ TenseurBH * res;
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
res->t[i] = this->t[i] / b; //division des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BH::operator /= ( const double & b)
{ for (int i = 0; i<=8; i++)
this->t[i] /= b;}; //division des données
#ifndef MISE_AU_POINT
inline
#endif
// affectation de B dans this, plusZero = false: les données manquantes sont inchangées,
// plusZero = true: les données manquantes sont mises à 0
void Tenseur3BH::Affectation_2D_a_3D(const Tenseur2BH & B,bool plusZero)
{ this->t[0] = B.t[0];this->t[4] = B.t[1];this->t[3] = B.t[2];this->t[1] = B.t[3];
if (plusZero)
{ this->t[2] = this->t[5] = this->t[8] = this->t[7] = this->t[6] = 0.;};
};
#ifndef MISE_AU_POINT
inline
#endif
// affectation de B dans this, plusZero = false: les données manquantes sont inchangées,
// plusZero = true: les données manquantes sont mises à 0
// si au contraire la dimension de B est plus grande que *this, il y a uniquement affectation
// des données possibles
void Tenseur3BH::Affectation_trans_dimension(const TenseurBH & B,bool plusZero)
{ switch (B.Dimension())
{case 3: *this = B; break; // affectation normale
case 2:
{ const Tenseur2BH & bn = *((Tenseur2BH *) &B);
this->Affectation_2D_a_3D(bn,plusZero);
break;
}
case 1:
{ if (plusZero)
this->Inita(0.);
this->t[0] = B.t[0]; //on affecte le seul terme
break;
}
default:
cout << "\n this= " << *this << " B= "; B.Ecriture(cout);
Message(3,
"erreur d\'affectation, Tenseur3BH::Affectation_trans_dimension( const TenseurBH & B, ..");
};
};
// produit contracte avec un vecteur
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeB Tenseur3BH::operator * ( const CoordonneeB & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != dimension)
{ cout << "\nErreur : dimensions vecteur tenseur non egales !\n";
cout << " Tenseur3BH::operator *\n";
Sortie(1);
};
#endif
CoordonneeB v(dimension);
v(1) = this->t[0] * B(1) + this->t[1] * B(2) + this->t[2] * B(3);
v(2) = this->t[3] * B(1) + this->t[4] * B(2) + this->t[5] * B(3);
v(3) = this->t[6] * B(1) + this->t[7] * B(2) + this->t[8] * B(3);
return v;
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BH::operator * ( const TenseurBB & B) const // produit une fois contracte
{ TenseurBB * res;
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BH::operator * ( etc..");
#endif
res = new Tenseur_ns3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
// on cree systematiquement un tenseur non symetrique
if (B.Dimension() == 3) // cas symetrique
{
res->t[0] = this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[5];
res->t[3] = this->t[3] * B.t[0] + this->t[4] * B.t[3] + this->t[5] * B.t[5];
res->t[6] = this->t[6] * B.t[0] + this->t[7] * B.t[3] + this->t[8] * B.t[5];
res->t[1] = this->t[0] * B.t[3] + this->t[1] * B.t[1] + this->t[2] * B.t[4];
res->t[4] = this->t[3] * B.t[3] + this->t[4] * B.t[1] + this->t[5] * B.t[4];
res->t[7] = this->t[6] * B.t[3] + this->t[7] * B.t[1] + this->t[8] * B.t[4];
res->t[2] = this->t[0] * B.t[5] + this->t[1] * B.t[4] + this->t[2] * B.t[2];
res->t[5] = this->t[3] * B.t[5] + this->t[4] * B.t[4] + this->t[5] * B.t[2];
res->t[8] = this->t[6] * B.t[5] + this->t[7] * B.t[4] + this->t[8] * B.t[2];
}
else // cas ou B n'est pas symetrique
{
res->t[0] = this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[4] * B.t[3] + this->t[5] * B.t[6];
res->t[6] = this->t[6] * B.t[0] + this->t[7] * B.t[3] + this->t[8] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[1] * B.t[4] + this->t[2] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[4] * B.t[4] + this->t[5] * B.t[7];
res->t[7] = this->t[6] * B.t[1] + this->t[7] * B.t[4] + this->t[8] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[1] * B.t[5] + this->t[2] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[4] * B.t[5] + this->t[5] * B.t[8];
res->t[8] = this->t[6] * B.t[2] + this->t[7] * B.t[5] + this->t[8] * B.t[8];
}
return (*res); };
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::operator * ( const TenseurBH & B) const // produit une fois contracte
{ TenseurBH * res;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator * ( etc..");
#endif
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[4] * B.t[3] + this->t[5] * B.t[6];
res->t[6] = this->t[6] * B.t[0] + this->t[7] * B.t[3] + this->t[8] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[1] * B.t[4] + this->t[2] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[4] * B.t[4] + this->t[5] * B.t[7];
res->t[7] = this->t[6] * B.t[1] + this->t[7] * B.t[4] + this->t[8] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[1] * B.t[5] + this->t[2] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[4] * B.t[5] + this->t[5] * B.t[8];
res->t[8] = this->t[6] * B.t[2] + this->t[7] * B.t[5] + this->t[8] * B.t[8];
return (*res); };
#ifndef MISE_AU_POINT
inline
#endif
// produit contracte contracté deux fois A(i,j)*B(j,i)=A..B c-c-d : A_i^{.k} * B_k^{.i}
// -> on contracte d'abord l'indice du milieu puis l'indice externe
// cela permet d'utiliser les mêmes variances pour les deux tenseurs
double Tenseur3BH::operator && ( const TenseurBH & B) const // produit deux fois contracte
{ double b = 0;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator && ( etc..");
#endif
b += this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[6];
b += this->t[3] * B.t[1] + this->t[4] * B.t[4] + this->t[5] * B.t[7];
b += this->t[6] * B.t[2] + this->t[7] * B.t[5] + this->t[8] * B.t[8];
return b;};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BH::Trace() const // trace du tenseur
{ double b;
b = this->t[0] + this->t[4] + this->t[8];
return b;};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BH::II() const // second invariant
{ //return Dabs((*this && *this)); // on met Dabs car dans les petits nombres on peut avoir des nombres
double ret = (*this && *this);
// if (Dabs(ret) < 0.1*ConstMath::pasmalpetit)
// on part du principe qu'il s'agit du produit doublement contracté d'un tenseur symétrique
// donc normalement il devrait donner un nombre positif, donc s'il est très petit, on le met arbitrairement à 0
if (ret < 0.)
{ret = 0.;}
return ret;
};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BH::III() const // troisieme invariant
{ Tenseur3BH a; a = ((*this) * (*this));
return (a && (*this)) ;};
// 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 = 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 Tenseur3BH::ValPropre(int& cas) const
{ Coordonnee ret(3);
// si le tenseur est sphérique, les trois valeurs propres sont identique au 1/3 de la trace
double b = this->Trace(); // la trace
double sqrt_II_b = sqrt(Dabs(this->II() - Sqr(b) / 3.)) ; // l'intensité du déviateur
if (sqrt_II_b <= ConstMath::unpeupetit * Dabs(b) )
{ cas = 0; ret(1)=ret(2)=ret(3)=b/3.; return ret;};
// sinon ..
b=-b; double c=0.5*(b*b-((*this)*(*this)).Trace()) ; double d=-this->Det();
// maintenant on va essayer de normaliser les termes de l'équation du troisième degré
// au début on a une équation x^3+B.x^2+C.x+D=0=(x-x1)(x-x2)(x-x3)
// on veut tenir compte de l'ordre de grandeur des xi avec une transformation: xp=a.x
// a étant un facteur de normalisation. L'équation en xpi devient:
// x^3+(a.B).x^2+(a^2.C).x+(a^3.D)=0=x^3+Bp.x^2+Cp.x+Dp, Bp Cp Dp étant les nouveaux coefficients
// on choisit pour a un coef fonction des valeurs des composantes du tenseur
double a_N = 1./ this->MaxiComposante(); // coeff de normalisation normalement non nulle
double bp = b * a_N; double cp = c * a_N*a_N; double dp = d * a_N*a_N*a_N;
// int caas;
// alg_zero.TroisiemeDegre dans certain ne fonctionne pas parfaitement, because , sans doute des arrondis
// j'essaie donc une résolution directe (voire après) qui a l'air de toujours fonctionner (a voir dans le temps)
// alg_zero.TroisiemeDegre(1.,b,c,d,ret(1),ret(2),ret(3),caas);
// #ifdef MISE_AU_POINT
// if (caas <= 0)
// { cout << "\n warning **** recherche de valeurs propre du tenseur: pas de racine correcte ! "
// << "\n tenseur= " << *this << "racines= " << ret
// << "\n on laisse le calcul se poursuivre cependant"
// << "\n Coordonnee Tenseur3BH::ValPropre()";
// };
// #endif
//
// ret /= a_N;
// cout << "\n \n " << ret(1) << " " << ret(2) << " " << ret(3);
//
// on va essayer d'utiliser les formules directes analytiques:
// 1) les invariants: Yavuz Basar and Dieter Weichert, "Nonlinear continuum Mechanis of solids"
// Sringer, p 33, ISBN 3-540-66601-X, 2000
double IA=-bp; double IIA=cp; double IIIA=-dp;
double tp1=IA*IA-3.*IIA; // grandeur intermédiaire qui doit-être positive
double t1 = 2.*sqrt(MaX(0.,tp1)); // on la met systématiquement à positif et on prend la racine
// si t1 est nul cela signifie qu'il y a trois valeurs propres égales
double untiers = 1./3.;
if (Abs(t1) <= ConstMath::trespetit)
{ret(1)=ret(2)=ret(3)=untiers*IA;}
else // sinon cas où l'on peut calculer t2
{// maintenant une seconde grandeur intermédiaire qui est un cosinus
double t2 = (2.*IA*IA*IA-9.*IA*IIA+27.*IIIA)/(t1*tp1);
if (Abs(t2) > 1.) t2 = Signe(t2); // limite à 1 au cas où
double theta = acos(t2);
// calcul des valeurs propres
ret(1)= untiers*(IA+t1*cos(untiers*(theta + 2.* ConstMath::Pi)));
ret(2)= untiers*(IA+t1*cos(untiers*(theta + 4.* ConstMath::Pi)));
ret(3)= untiers*(IA+t1*cos(untiers*(theta + 6.* ConstMath::Pi)));
};
ret /= a_N;
// caas = 5; // qui est le cas 3 valeurs propres, mais en fait elles peuvent être égales par exemple
/// cout << "\n " << ret(1) << " " << ret(2) << " " << ret(3);
// classement des valeurs propres
double vv1,vv2,vv3;
// if (caas == 5) // cas où il y a 3 racines réelles
// il n'y a plus qu'un seul cas
{// d'abord on regarde le plus grand que l'on met dans ret(1)
if (ret(1) >= ret(2))
{ if (ret(1) >= ret(3)) { vv1=ret(1); vv2=ret(2); vv3=ret(3);}
else { vv2=ret(1); ret(1)=vv1=ret(3); vv3=ret(2);};
}
else
{ if (ret(2) >= ret(3)) { vv2=ret(1);ret(1)=vv1=ret(2); vv3=ret(3);}
else { vv2=ret(1);ret(1)=vv1=ret(3); vv3=ret(2);};
};
// puis les deux derniers
if (vv2 < vv3) { ret(2) = vv3; ret(3) = vv2;vv2=vv3;vv3=ret(2);}
else { ret(2) = vv2; ret(3) = vv3;};
// définition du cas de valeurs propres
// 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)
if ( difftrespetit(vv1,vv2) && difftrespetit(vv2,vv3)) {cas = 1;}
else if ( (EgaleApeuPres(vv1,vv2)) && EgaleApeuPres(vv2,vv3)) {cas = 0;}
else if ( (EgaleApeuPres(vv1,vv2)) && !EgaleApeuPres(vv2,vv3)) {cas = 2;}
else if ( !(EgaleApeuPres(vv1,vv2)) && EgaleApeuPres(vv2,vv3)) {cas = 3;}
else
{cout << "\n classement des valeurs propres : non realise !! "
<< "\n vv1= "<<vv1 <<", vv2= "<< vv2 <<", vv3= "<< vv3
<< "\n r(1)= "<<ret(1) << ", r(2)= "<<ret(2) << ", r(3)= "<<ret(3)
<< "\n Tenseur3BH::ValPropre(... "<<endl;
};
}
// else if ((caas == 4) || (caas == 1)) // cas où il y a 2 racines
// {//on regarde le plus grand
// if (ret(1) < ret(2))
// { vv2=ret(1); ret(1)=ret(2);ret(2)=vv2; cas=2;};
// // sinon l'ordre est déjà bon
// cas=3;
// }
// else if ((caas == 2) || (caas == 3))
// { // cas d'une seule racine
// cas=0;
// };
// // s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
// if (caas <= 0) cas=-1;
// if (cas <= 0)
// { cout << "\n warning **** recherche de valeurs propre du tenseur: pas de racine correcte ! "
// << "\n tenseur= " << *this << "racines= " << ret
// << "\n on laisse le calcul se poursuivre cependant"
// << "\n Coordonnee Tenseur3BH::ValPropre()";
// };
// retour
return ret;
};
// 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 à HB)
#ifndef MISE_AU_POINT
inline
#endif
Coordonnee Tenseur3BH::ValPropre(int& cas, Mat_pleine& mat) const
{ Coordonnee ret(3);
// si le tenseur est sphérique, les trois valeurs propres sont identique au 1/3 de la trace
// on ramène le tenseur identité en absolue comme base de vecteurs propres car toute base est ok
double b = this->Trace(); // la trace
double sqrt_II_b = sqrt(Dabs(this->II() - Sqr(b) / 3.)) ; // l'intensité du déviateur
if (sqrt_II_b <= ConstMath::unpeupetit * Dabs(b) )
{ mat.Initialise(0.); mat(1,1)=1.;mat(2,2)=1.;mat(3,3)=1.;
cas=0;ret(1)=ret(2)=ret(3)=b/3.;return ret;
};
// sinon on commence par calculer les valeurs propres
Coordonnee valPropre = Tenseur3BH::ValPropre(cas);
// --- gestion des erreurs éventuelles
if (cas == -1) // cas d'une erreur
{mat.Zero(); // mise à 0 de la matrice
ret.Zero(); // mise à 0 des valeurs propres
return ret;
};
// sinon ok on continue
// on crée une matrice pour utiliser une routine générale
#ifdef MISE_AU_POINT
if ((mat.Nb_ligne() != 3) || (mat.Nb_colonne() != 3))
{ cout << "\nErreur : la matrice en parametre doit etre de dimension 3x3 ";
cout << "\n Tenseur3BH::ValPropre(... \n";
Sortie(1);
};
#endif
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
mat(i,j) = (*this)(i,j);
// int cass;
Tableau <Coordonnee> V_P = MathUtil2::V_Propres3x3(mat,valPropre,cas);
// --- gestion des erreurs éventuelles
if (cas == -1) // cas d'une erreur
{mat.Zero(); // mise à 0 de la matrice
ret.Zero(); // mise à 0 des valeurs propres
return ret;
};
// sinon c'est ok, on continue
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
mat(j,i) = V_P(i)(j); // le vecteur propre i, sa coordonnée j
//retour
return ret;
};
// ici il s'agit uniquement de calculer les vecteurs propres, les valeurs propres
// étant déjà connues
// en retour VP les vecteurs propre : doivent avoir la dimension du tenseur
// les vecteurs propre sont exprime dans le repere naturel (pour les tenseurs dim 3
// pour dim=2:le premier vecteur propre est exprime dans le repere naturel
// le second vecteur propre est exprimé dans le repère dual
// pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
// en sortie cas = -1 s'il y a eu un problème, dans ce cas, V_P est quelconque
// sinon si tout est ok, cas est identique en sortie avec l'entrée
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BH::VecteursPropres(const Coordonnee& valPropre,int& cas, Tableau <Coordonnee>& V_P) const
{ // on dimensionne en conséquence
Coordonnee toto(3); // un vecteur nul
V_P.Change_taille(3,toto); // tableau initialisé à 0
// si le tenseur est sphérique, les trois valeurs propres sont identique au 1/3 de la trace
// on ramène le tenseur identité en absolue comme base de vecteurs propres car toute base est ok
double b = this->Trace(); // la trace
double sqrt_II_b = sqrt(Dabs(this->II() - Sqr(b) / 3.)) ; // l'intensité du déviateur
if (sqrt_II_b <= ConstMath::unpeupetit * Dabs(b) )
{
#ifdef MISE_AU_POINT
if (cas != 0)
{ cout << "\nErreur : les valeurs propres devraient etre identiques alors qu'en entree: cas= " << cas ;
cout << "\n Tenseur3BH::VecteursPropres(... \n";
cas=-1;
return;
};
#endif
V_P(1).Zero(); V_P(1)(1)=1.;
V_P(2).Zero(); V_P(2)(2)=1.;
V_P(3).Zero(); V_P(3)(3)=1.;
cas=0;return;
};
// sinon ok on continue
// on crée une matrice pour utiliser une routine générale
Mat_pleine mat(3,3,0.);
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
mat(i,j) = (*this)(i,j);
// int cass;
V_P = MathUtil2::V_Propres3x3(mat,valPropre,cas);
};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BH::Det() const // determinant de la matrice des coordonnees
{ double b = 0;
b += this->t[0] *(this->t[4] * this->t[8] - this->t[5] * this->t[7]);
b -= this->t[3] *(this->t[1] * this->t[8] - this->t[2] * this->t[7]);
b += this->t[6] *(this->t[1] * this->t[5] - this->t[2] * this->t[4]);
return b;};
// test
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3BH::operator == ( const TenseurBH & B) const
{ int res = 1;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator == ( etc..");
#endif
for (int i = 0; i<=8; i++)
if (this->t[i] != B.t[i]) res = 0 ;
return res; };
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3BH::operator != ( const TenseurBH & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BH::operator != ( etc..");
#endif
if ((*this) == B)
return 0;
else
return 1; };
// calcul du tenseur inverse par rapport au produit contracte
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::Inverse() const
{TenseurBH * res;
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
// choix sur la méthode d'inversion
switch (ParaGlob::param->ParaAlgoControleActifs().Type_calnum_inversion_metrique())
{ case LU_EQUILIBRE:
{ // on recopie this dans le nouveau tenseur
for (int i = 0; i< 9; i++)
res->t[i] = t[i];
// pour le débug
//res->t[0]=3.; res->t[1]=2.;res->t[2]=1.;
// appel de l'inversion
Util::Inverse_mat3x3(((Tenseur3BH *) res)->ipointe);
}
break;
//cout << "\n comp \n ";
// res->Ecriture(cout); cout << "\n";
case CRAMER : // méthode historique
{ // calcul du determinant
double det = Det();
#ifdef MISE_AU_POINT
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nErreur : le determinant du tenseur est nul !\n";
cout << "Tenseur3BH::Inverse() \n";
Sortie(1);
};
#endif
det =1./det;
res->t[0] = (this->t[4]*this->t[8] - this->t[5]*this->t[7])*det;
res->t[3] = (this->t[5]*this->t[6] - this->t[3]*this->t[8])*det;
res->t[6] = (this->t[3]*this->t[7] - this->t[4]*this->t[6])*det;
res->t[1] = (this->t[2]*this->t[7] - this->t[1]*this->t[8])*det;
res->t[4] = (this->t[0]*this->t[8] - this->t[2]*this->t[6])*det;
res->t[7] = (this->t[1]*this->t[6] - this->t[0]*this->t[7])*det;
res->t[2] = (this->t[1]*this->t[5] - this->t[2]*this->t[4])*det;
res->t[5] = (this->t[2]*this->t[3] - this->t[0]*this->t[5])*det;
res->t[8] = (this->t[0]*this->t[4] - this->t[1]*this->t[3])*det;
}
break;
default:
{ cout << "\nErreur **** : la methode de resolution de l'inversion de tenseur "
<< ParaGlob::param->ParaAlgoControleActifs().Type_calnum_inversion_metrique() << " n'est pas implante \n";
cout << "Tenseur3BH::Inverse() \n";
Sortie(1);
};
break;
};
// res->Ecriture(cout); // pour le debug
return *res;
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3BH::Transpose() const
{ TenseurHB * res;
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon(res); // ajout d'un tenseur intermediaire
res->t[0] = this->t[0];
res->t[1] = this->t[3];
res->t[2] = this->t[6];
res->t[3] = this->t[1];
res->t[4] = this->t[4];
res->t[5] = this->t[7];
res->t[6] = this->t[2];
res->t[7] = this->t[5];
res->t[8] = this->t[8];
return *res;};
#ifndef MISE_AU_POINT
inline
#endif
// permute Bas Haut, mais reste dans le même tenseur
void Tenseur3BH::PermuteHautBas()
{ // t[0] = t[0];
double t1=t[1]; t[1] = t[3];
double t2=t[2]; t[2] = t[6];
t[3] = t1;
//[4] = [4];
double t5=t[5]; t[5] = t[7];
t[6] = t2;
t[7] = t5;
//t[8] = t[8];
};
// calcul du maximum en valeur absolu des composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BH::MaxiComposante() const
{ return DabsMaxiTab(t, 9) ;
};
// retourne la composante i,j en lecture et écriture
#ifndef MISE_AU_POINT
inline
#endif
double& Tenseur3BH::Coor( const int i, const int j)
{
#ifdef MISE_AU_POINT
if ( ((i!=1)&&(i!=2)&&(i!=3)) || ((j!=1)&&(j!=2)&&(j!=3)) )
{ cout << "\nErreur : composante inexistante !\n";
cout << " i = " << i << ", j = " << j << '\n';
cout << "Tenseur3BH::OPERATOR() (int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[1]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[4]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[6]; break;
case 2 : return t[7]; break;
case 3 : return t[8]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
// retourne la composante i,j en lecture seulement
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BH::operator () ( const int i, const int j) const
{
#ifdef MISE_AU_POINT
if ( ((i!=1)&&(i!=2)&&(i!=3)) || ((j!=1)&&(j!=2)&&(j!=3)) )
{ cout << "\nErreur : composante inexistante !\n";
cout << " i = " << i << ", j = " << j << '\n';
cout << "Tenseur3BH::OPERATOR() (int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[1]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[4]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[6]; break;
case 2 : return t[7]; break;
case 3 : return t[8]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
//fonctions static définissant le produit tensoriel de deux vecteurs
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BH::Prod_tensoriel(const CoordonneeB & aB, const CoordonneeH & bH)
{ TenseurBH * res;
#ifdef MISE_AU_POINT
if ((aB.Dimension() != 3) || (bH.Dimension() != 3))
{ cout << "\n erreur de dimension dans les coordonnées d'entrée, dim1 et dim2 ="
<< aB.Dimension() << " " << bH.Dimension()
<< "\n Tenseur3BH::Prod_tensoriel( etc.." << endl;
Sortie(1);
}
#endif
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = aB(1) * bH(1); res->t[1] = aB(1) * bH(2); res->t[2] = aB(1) * bH(3);
res->t[3] = aB(2) * bH(1); res->t[4] = aB(2) * bH(2); res->t[5] = aB(2) * bH(3);
res->t[6] = aB(3) * bH(1); res->t[7] = aB(3) * bH(2); res->t[8] = aB(3) * bH(3);
return *res ;};
// lecture et écriture de données
#ifndef MISE_AU_POINT
inline
#endif
istream & Tenseur3BH::Lecture(istream & entree)
{ // lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3BH")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i< 9; i++)
entree >> this->t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
ostream & Tenseur3BH::Ecriture(ostream & sort) const
{ // écriture du type
sort << "Tenseur3BH ";
// puis les datas
for (int i = 0; i< 9; i++)
sort << setprecision(ParaGlob::NbdigdoCA()) << this->t[i] << " ";
return sort;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator de lecture
istream & operator >> (istream & entree, Tenseur3BH & A)
{ int dim = A.Dimension();
#ifdef MISE_AU_POINT
if (dim != 3) A.Message(3,"operator >> (istream & entree, Tenseur3BH & A)");
#endif
// lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3BH")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i< 9 ; i++)
entree >> A.t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator d'ecriture
ostream & operator << (ostream & sort , const Tenseur3BH & A)
{ //int dim = A.Dimension();
// écriture du type
sort << "Tenseur3BH ";
// puis les datas
for (int i = 0; i< 9; i++)
sort << setprecision(ParaGlob::NbdigdoCA()) << A.t[i] << " ";
return sort;
};
//------------------------------------------------------------------
// cas des composantes mixtes HB
//------------------------------------------------------------------
// --- gestion de changement d'index ----
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::ChangementIndex::ChangementIndex() :
idx_i(9),idx_j(9),odVect(3)
{ idx_i(1)=1;idx_i(2)=1;idx_i(3)=1; idx_j(1)=1;idx_j(2)=2;idx_j(3)=3;
idx_i(4)=2;idx_i(5)=2;idx_i(6)=2; idx_j(4)=1;idx_j(5)=2;idx_j(6)=3;
idx_i(7)=3;idx_i(8)=3;idx_i(9)=3; idx_j(7)=1;idx_j(8)=2;idx_j(9)=3;
odVect(1,1)=1;odVect(1,2)=2;odVect(1,3)=3;
odVect(2,1)=4;odVect(2,2)=5;odVect(2,3)=6;
odVect(3,1)=7;odVect(3,2)=8;odVect(3,3)=9;
};
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::Tenseur3HB() :
ipointe() // par defaut
{ dimension = 3;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = 0.;
};
// initialisation de toutes les composantes a une meme valeur
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::Tenseur3HB( const double val) :
ipointe()
{ dimension = 3;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = val;
};
// initialisation avec 9 valeurs différentes
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::Tenseur3HB
( const double val1, const double val2, const double val3, // 1ere ligne
const double val4, const double val5, const double val6, // 2ieme ligne
const double val7, const double val8, const double val9) // 3ieme ligne
:
ipointe()
{ dimension = 3;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
t[0] = val1;t[1] = val2;t[2] = val3;
t[3] = val4;t[4] = val5;t[5] = val6;
t[6] = val7;t[7] = val8;t[8] = val9;
};
// constructeur a partir d'une instance non differenciee
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::Tenseur3HB ( const TenseurHB & B):
ipointe()
{ this->dimension = B.dimension ;
#ifdef MISE_AU_POINT
if (dimension != 3)
{ cout << "\n erreur de dimension, elle devrait etre = 3 ";
cout << "\n Tenseur3HB::Tenseur3HB ( TenseurHB &) " << endl;
Sortie(1);
}
#endif
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = B.t[i];
};
// constructeur de copie
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::Tenseur3HB ( const Tenseur3HB & B):
ipointe()
{ this->dimension = B.dimension ;
listdouble9.push_front(Reels9()); // allocation
ipointe = listdouble9.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
for ( int i=0; i< 9; i++)
t[i] = B.t[i];
};
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HB::~Tenseur3HB()
{//if(listdouble9.end() != listdouble9.begin()) // si la liste n'est pas vide
listdouble9.erase(ipointe);} ; // suppression de l'élément de la liste
// initialise toutes les composantes à val
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HB::Inita(double val)
{ t[0] =val; t[1] =val; t[2] =val;
t[3] =val; t[4] =val; t[5] =val;
t[6] =val; t[7] =val; t[8] =val;
};
// operations
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator + (const TenseurHB & B) const
{ TenseurHB * ptr;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator + ( etc..");
#endif
ptr = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( ptr); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
(*ptr).t[i] = this->t[i] + B.t[i]; //somme des données
return (*ptr) ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HB::operator += ( const TenseurHB & B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator += ( etc..");
#endif
for (int i = 0; i<=8; i++)
this->t[i] += B.t[i];}; //somme des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator - () const
{ TenseurHB * ptr;
ptr = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( ptr); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
(*ptr).t[i] = - this->t[i]; // soustraction des données
return (*ptr) ;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator - ( const TenseurHB & B) const
{ TenseurHB * ptr;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator - ( etc..");
#endif
ptr = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( ptr); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
(*ptr).t[i] = this->t[i] - B.t[i]; // soustraction des données
return (*ptr) ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HB::operator -= ( const TenseurHB & B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator -= ( etc..");
#endif
for (int i = 0; i<=8; i++)
this->t[i] -= B.t[i];}; //soustraction des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator = ( const TenseurHB & B)
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator = ( etc..");
#endif
for (int i = 0; i<=8; i++)
this->t[i] = B.t[i];
LesMaillonsHB::Libere(); // destruction des tenseurs intermediaires
return *this; }; //affectation des données;
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator * ( const double & b) const
{ TenseurHB * res;
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
res->t[i] = this->t[i] * b; //multiplication des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HB::operator *= ( const double & b)
{ for (int i = 0; i<=8; i++)
this->t[i] *= b;}; //multiplication des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator / ( const double & b) const
{ TenseurHB * res;
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i<=8; i++)
res->t[i] = this->t[i] / b; //division des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HB::operator /= ( const double & b)
{ for (int i = 0; i<=8; i++)
this->t[i] /= b;}; //division des données
#ifndef MISE_AU_POINT
inline
#endif
// affectation de B dans this, plusZero = false: les données manquantes sont inchangées,
// plusZero = true: les données manquantes sont mises à 0
void Tenseur3HB::Affectation_2D_a_3D(const Tenseur2HB & B,bool plusZero)
{ this->t[0] = B.t[0];this->t[4] = B.t[1];this->t[3] = B.t[2];this->t[1] = B.t[3];
if (plusZero)
{ this->t[2] = this->t[5] = this->t[8] = this->t[7] = this->t[6] = 0.;};
};
#ifndef MISE_AU_POINT
inline
#endif
// affectation de B dans this, plusZero = false: les données manquantes sont inchangées,
// plusZero = true: les données manquantes sont mises à 0
// si au contraire la dimension de B est plus grande que *this, il y a uniquement affectation
// des données possibles
void Tenseur3HB::Affectation_trans_dimension(const TenseurHB & B,bool plusZero)
{ switch (B.Dimension())
{case 3: *this = B; break; // affectation normale
case 2:
{ const Tenseur2HB & bn = *((Tenseur2HB *) &B);
this->Affectation_2D_a_3D(bn,plusZero);
break;
}
case 1:
{ if (plusZero)
this->Inita(0.);
this->t[0] = B.t[0]; //on affecte le seul terme
break;
}
default:
cout << "\n this= " << *this << " B= "; B.Ecriture(cout);
Message(3,
"erreur d\'affectation, Tenseur3HB::Affectation_trans_dimension( const TenseurHB & B, ..");
};
};
// produit contracte avec un vecteur
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeH Tenseur3HB::operator * ( const CoordonneeH & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != dimension)
{ cout << "\nErreur : dimensions vecteur tenseur non egales !\n";
cout << " Tenseur3HB::operator *\n";
Sortie(1);
};
#endif
CoordonneeH v(dimension);
v(1) = this->t[0] * B(1) + this->t[1] * B(2) + this->t[2] * B(3);
v(2) = this->t[3] * B(1) + this->t[4] * B(2) + this->t[5] * B(3);
v(3) = this->t[6] * B(1) + this->t[7] * B(2) + this->t[8] * B(3);
return v;
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HB::operator * ( const TenseurHH & B) const // produit une fois contracte
{ TenseurHH * res;
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HB::operator * ( etc..");
#endif
res = new Tenseur_ns3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
// on cree systematiquement un tenseur non symetrique
if (B.Dimension() == 3) // cas symetrique
{
res->t[0] = this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[5];
res->t[3] = this->t[3] * B.t[0] + this->t[4] * B.t[3] + this->t[5] * B.t[5];
res->t[6] = this->t[6] * B.t[0] + this->t[7] * B.t[3] + this->t[8] * B.t[5];
res->t[1] = this->t[0] * B.t[3] + this->t[1] * B.t[1] + this->t[2] * B.t[4];
res->t[4] = this->t[3] * B.t[3] + this->t[4] * B.t[1] + this->t[5] * B.t[4];
res->t[7] = this->t[6] * B.t[3] + this->t[7] * B.t[1] + this->t[8] * B.t[4];
res->t[2] = this->t[0] * B.t[5] + this->t[1] * B.t[4] + this->t[2] * B.t[2];
res->t[5] = this->t[3] * B.t[5] + this->t[4] * B.t[4] + this->t[5] * B.t[2];
res->t[8] = this->t[6] * B.t[5] + this->t[7] * B.t[4] + this->t[8] * B.t[2];
}
else // cas ou B n'est pas symetrique
{
res->t[0] = this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[4] * B.t[3] + this->t[5] * B.t[6];
res->t[6] = this->t[6] * B.t[0] + this->t[7] * B.t[3] + this->t[8] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[1] * B.t[4] + this->t[2] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[4] * B.t[4] + this->t[5] * B.t[7];
res->t[7] = this->t[6] * B.t[1] + this->t[7] * B.t[4] + this->t[8] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[1] * B.t[5] + this->t[2] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[4] * B.t[5] + this->t[5] * B.t[8];
res->t[8] = this->t[6] * B.t[2] + this->t[7] * B.t[5] + this->t[8] * B.t[8];
}
return (*res); };
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::operator * ( const TenseurHB & B) const // produit une fois contracte
{ TenseurHB * res;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator * ( etc..");
#endif
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[4] * B.t[3] + this->t[5] * B.t[6];
res->t[6] = this->t[6] * B.t[0] + this->t[7] * B.t[3] + this->t[8] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[1] * B.t[4] + this->t[2] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[4] * B.t[4] + this->t[5] * B.t[7];
res->t[7] = this->t[6] * B.t[1] + this->t[7] * B.t[4] + this->t[8] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[1] * B.t[5] + this->t[2] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[4] * B.t[5] + this->t[5] * B.t[8];
res->t[8] = this->t[6] * B.t[2] + this->t[7] * B.t[5] + this->t[8] * B.t[8];
return (*res); };
#ifndef MISE_AU_POINT
inline
#endif
// produit contracte contracté deux fois A(i,j)*B(j,i)=A..B c-c-d : A^i_{.k} * B^k_{.i}
// -> on contracte d'abord l'indice du milieu puis l'indice externe
// cela permet d'utiliser les mêmes variances pour les deux tenseurs
double Tenseur3HB::operator && ( const TenseurHB & B) const // produit deux fois contracte
{ double b = 0;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator && ( etc..");
#endif
b += this->t[0] * B.t[0] + this->t[1] * B.t[3] + this->t[2] * B.t[6];
b += this->t[3] * B.t[1] + this->t[4] * B.t[4] + this->t[5] * B.t[7];
b += this->t[6] * B.t[2] + this->t[7] * B.t[5] + this->t[8] * B.t[8];
return b;};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HB::Trace() const // trace du tenseur
{ double b;
b = this->t[0] + this->t[4] + this->t[8];
return b;};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HB::II() const // second invariant
{ //return Dabs((*this && *this)); // on met Dabs car dans les petits nombres on peut avoir des nombres
double ret = (*this && *this);
// if (Dabs(ret) < 0.1*ConstMath::pasmalpetit)
// on part du principe qu'il s'agit du produit doublement contracté d'un tenseur symétrique
// donc normalement il devrait donner un nombre positif, donc s'il est très petit, on le met arbitrairement à 0
if (ret < 0.)
{ret = 0.;}
return ret;
};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HB::III() const // troisieme invariant
{ Tenseur3HB a; a = ((*this) * (*this));
return (a && (*this)) ;};
// valeurs propre dans le vecteur de retour, classée par ordres "décroissants"
// cas indique le cas de valeur propre:
// 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 = 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 Tenseur3HB::ValPropre(int& cas) const
{ Coordonnee ret(3);
// si le tenseur est sphérique, les trois valeurs propres sont identique au 1/3 de la trace
double b = this->Trace(); // la trace
double sqrt_II_b = sqrt(Dabs(this->II() - Sqr(b) / 3.)) ; // l'intensité du déviateur
if (sqrt_II_b <= ConstMath::unpeupetit * Dabs(b) )
{ cas = 0; ret(1)=ret(2)=ret(3)=b/3.; return ret;};
// sinon ..
b=-b; double c=0.5*(b*b-((*this)*(*this)).Trace()) ; double d=-this->Det();
// maintenant on va essayer de normaliser les termes de l'équation du troisième degré
// au début on a une équation x^3+B.x^2+C.x+D=0=(x-x1)(x-x2)(x-x3)
// on veut tenir compte de l'ordre de grandeur des xi avec une transformation: xp=a.x
// a étant un facteur de normalisation. L'équation en xpi devient:
// x^3+(a.B).x^2+(a^2.C).x+(a^3.D)=0=x^3+Bp.x^2+Cp.x+Dp, Bp Cp Dp étant les nouveaux coefficients
// on choisit pour a un coef fonction des valeurs des composantes du tenseur
double a_N = 1./ this->MaxiComposante(); // coeff de normalisation normalement non nulle
double bp = b * a_N; double cp = c * a_N*a_N; double dp = d * a_N*a_N*a_N;
// int caas; // pour la gestion des différents cas de retour de alg_zero.TroisiemeDegre
// alg_zero.TroisiemeDegre dans certain ne fonctionne pas parfaitement, because , sans doute des arrondis
// j'essaie donc une résolution directe (voire après) qui a l'air de toujours fonctionner (a voir dans le temps)
// alg_zero.TroisiemeDegre(1.,bp,cp,dp,ret(1),ret(2),ret(3),caas);
// #ifdef MISE_AU_POINT
// if (caas <= 0)
// { cout << "\n warning **** recherche de valeurs propre du tenseur: pas de racine correcte ! "
// << "\n tenseur= " << *this << "racines= " << ret
// << "\n on laisse le calcul se poursuivre cependant"
// << "\n Coordonnee Tenseur3HB::ValPropre()";
// cas=-1;
// ret.Zero();
// // retour
// return ret;
// };
// #endif
// ret /= a_N;
// cout << "\n \n " << ret(1) << " " << ret(2) << " " << ret(3);
//
// on va essayer d'utiliser les formules directes analytiques:
// 1) les invariants: Yavuz Basar and Dieter Weichert, "Nonlinear continuum Mechanis of solids"
// Sringer, p 33, ISBN 3-540-66601-X, 2000
double IA=-bp; double IIA=cp; double IIIA=-dp;
double tp1=IA*IA-3.*IIA; // grandeur intermédiaire qui doit-être positive
double t1 = 2.*sqrt(MaX(0.,tp1)); // on la met systématiquement à positif et on prend la racine
// si t1 est nul cela signifie qu'il y a trois valeurs propres égales
double untiers = 1./3.;
if (Abs(t1) <= ConstMath::trespetit)
{ret(1)=ret(2)=ret(3)=untiers*IA;}
else // sinon cas où l'on peut calculer t2
{// maintenant une seconde grandeur intermédiaire qui est un cosinus
double t2 = (2.*IA*IA*IA-9.*IA*IIA+27.*IIIA)/(t1*tp1);
if (Abs(t2) > 1.) t2 = Signe(t2); // limite à 1 au cas où
double theta = acos(t2);
// calcul des valeurs propres
ret(1)= untiers*(IA+t1*cos(untiers*(theta + 2.* ConstMath::Pi)));
ret(2)= untiers*(IA+t1*cos(untiers*(theta + 4.* ConstMath::Pi)));
ret(3)= untiers*(IA+t1*cos(untiers*(theta + 6.* ConstMath::Pi)));
};
ret /= a_N;
/// caas = 5; // qui est le cas 3 valeurs propres, mais en fait elles peuvent être égales par exemple
/// cout << "\n " << ret(1) << " " << ret(2) << " " << ret(3);
// // dans certains cas meme avec caas ok, il y a des pb ??
// // dans ce cas on utilise lapack
// bool pb=false;
// switch(caas)
// {
//// !!!!!!!!! a revoir car ne fonctionne pas en fait !!!!!!!!!
// case 5: if(ret(3)== ConstMath::tresgrand) pb=true;
// case 4 : case 1 : if(ret(2)== ConstMath::tresgrand) pb=true;
// case 2 : case 3 : if(ret(1)== ConstMath::tresgrand) pb=true;
// default: break;
// };
// if (pb)
// { // pour l'instant on brame !!
// #ifdef MISE_AU_POINT
// if (caas <= 0)
// { cout << "\n warning **** recherche de valeurs propre du tenseur: pas de racine correcte ! "
// << "\n tenseur= " << *this << "racines= " << ret
// << "\n on laisse le calcul se poursuivre cependant"
// << "\n Coordonnee Tenseur3HB::ValPropre()";
// };
// #endif
//
// }
// classement des valeurs propres
double vv1,vv2,vv3;
// if (caas == 5) // cas où il y a 3 racines réelles
// il n'y a plus qu'un seul cas
{// d'abord on regarde le plus grand que l'on met dans ret(1)
if (ret(1) >= ret(2))
{ if (ret(1) >= ret(3)) { vv1=ret(1); vv2=ret(2); vv3=ret(3);}
else { vv2=ret(1); ret(1)=vv1=ret(3); vv3=ret(2);};
}
else
{ if (ret(2) >= ret(3)) { vv2=ret(1);ret(1)=vv1=ret(2); vv3=ret(3);}
else { vv2=ret(1);ret(1)=vv1=ret(3); vv3=ret(2);};
};
// puis les deux derniers
if (vv2 < vv3) { ret(2) = vv3; ret(3) = vv2;vv2=vv3;vv3=ret(2);}
else { ret(2) = vv2; ret(3) = vv3;};
// définition du cas de valeurs propres
// 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)
if ( difftrespetit(vv1,vv2) && difftrespetit(vv2,vv3)) {cas = 1;}
else if ( (EgaleApeuPres(vv1,vv2)) && EgaleApeuPres(vv2,vv3)) {cas = 0;}
else if ( (EgaleApeuPres(vv1,vv2)) && !EgaleApeuPres(vv2,vv3)) {cas = 2;}
else if ( !(EgaleApeuPres(vv1,vv2)) && EgaleApeuPres(vv2,vv3)) {cas = 3;}
else
{cout << "\n classement des valeurs propres : non realise !! "
<< "\n vv1= "<<vv1 <<", vv2= "<< vv2 <<", vv3= "<< vv3
<< "\n r(1)= "<<ret(1) << ", r(2)= "<<ret(2) << ", r(3)= "<<ret(3)
<< "\n Tenseur3HB::ValPropre(... "<<endl;
};
}
// else if ((caas == 4) || (caas == 1)) // cas où il y a 2 racines
// {//on regarde le plus grand
// if (ret(1) < ret(2))
// { vv2=ret(1); ret(1)=ret(2);ret(2)=vv2; cas=2;};
// // sinon l'ordre est déjà bon
// cas=3;
// }
// else if ((caas == 2) || (caas == 3))
// {// cas d'une seule racine
// cas=0;
// };
// s'il y a eu erreur dans la recherche de valeurs propre on modifie le cas
// if (caas <= 0) {cas=-1;ret.Zero();};
// if (cas <= 0)
// { cout << "\n warning **** recherche de valeurs propre du tenseur: pas de racine correcte ! "
// << "\n tenseur= " << *this << "racines= " << ret
// << "\n on laisse le calcul se poursuivre cependant"
// << "\n Coordonnee Tenseur3HB::ValPropre()";
// };
// retour
return ret;
};
// 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 Tenseur3HB::ValPropre(int& cas, Mat_pleine& mat) const
{ Coordonnee ret(3);
// si le tenseur est sphérique, les trois valeurs propres sont identique au 1/3 de la trace
// on ramène le tenseur identité en absolue comme base de vecteurs propres car toute base est ok
double b = this->Trace(); // la trace
double sqrt_II_b = sqrt(Dabs(this->II() - Sqr(b) / 3.)) ; // l'intensité du déviateur
if (sqrt_II_b <= ConstMath::unpeupetit * Dabs(b) )
{ mat.Initialise(0.); mat(1,1)=1.;mat(2,2)=1.;mat(3,3)=1.;
cas=0;ret(1)=ret(2)=ret(3)=b/3.;return ret;
};
// sinon on commence par rechercher les valeurs propres
// sinon on commence par calculer les valeurs propres
Coordonnee valPropre = Tenseur3HB::ValPropre(cas);
// --- gestion des erreurs éventuelles
if (cas == -1) // cas d'une erreur
{mat.Zero(); // mise à 0 de la matrice
ret.Zero(); // mise à 0 des valeurs propres
return ret;
};
// sinon ok on continue
// on crée une matrice pour utiliser une routine générale
#ifdef MISE_AU_POINT
if ((mat.Nb_ligne() != 3) || (mat.Nb_colonne() != 3))
{ cout << "\nErreur : la matrice en parametre doit etre de dimension 3x3 ";
cout << "\n Tenseur3BH::ValPropre(... \n";
Sortie(1);
};
#endif
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
mat(i,j) = (*this)(i,j);
// int cass;
Tableau <Coordonnee> V_P = MathUtil2::V_Propres3x3(mat,valPropre,cas);
// --- gestion des erreurs éventuelles
if (cas == -1) // cas d'une erreur
{mat.Zero(); // mise à 0 de la matrice
ret.Zero(); // mise à 0 des valeurs propres
return ret;
};
// sinon c'est ok, on continue
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
mat(j,i) = V_P(i)(j); // le vecteur propre i, sa coordonnée j
//retour
return ret;
};
// ici il s'agit uniquement de calculer les vecteurs propres, les valeurs propres
// étant déjà connues
// en retour VP les vecteurs propre: doivent avoir la dimension du tenseur
// les vecteurs propre sont exprime dans le repere naturel (pour les tenseurs dim 3
// pour dim=2:le premier vecteur propre est exprime dans le repere naturel
// le second vecteur propre est exprimé dans le repère dual
// pour dim=1 le vecteur est dans la base naturelle (mais ça importe pas)
// en sortie cas = -1 s'il y a eu un problème, dans ce cas, V_P est quelconque
// sinon si tout est ok, cas est identique en sortie avec l'entrée
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HB::VecteursPropres(const Coordonnee& valPropre,int& cas, Tableau <Coordonnee>& V_P) const
{ // on dimensionne en conséquence
Coordonnee toto(3); // un vecteur nul
V_P.Change_taille(3,toto); // tableau initialisé à 0
// si le tenseur est sphérique, les trois valeurs propres sont identique au 1/3 de la trace
// on ramène le tenseur identité en absolue comme base de vecteurs propres car toute base est ok
double b = this->Trace(); // la trace
double sqrt_II_b = sqrt(Dabs(this->II() - Sqr(b) / 3.)) ; // l'intensité du déviateur
if (sqrt_II_b <= ConstMath::unpeupetit * Dabs(b) )
{
#ifdef MISE_AU_POINT
if (cas != 0)
{ cout << "\nErreur : les valeurs propres devraient etre identiques alors qu'en entree: cas= " << cas ;
cout << "\n Tenseur3HB::VecteursPropres(... \n";
cas=-1;
return;
};
#endif
V_P(1).Zero(); V_P(1)(1)=1.;
V_P(2).Zero(); V_P(2)(2)=1.;
V_P(3).Zero(); V_P(3)(3)=1.;
cas=0;return;
};
// sinon ok on continue
// on crée une matrice pour utiliser une routine générale
Mat_pleine mat(3,3,0.);
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
mat(i,j) = (*this)(i,j);
// int cass;
V_P = MathUtil2::V_Propres3x3(mat,valPropre,cas);
};
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HB::Det() const // determinant de la matrice des coordonnees
{ double b = 0;
b += this->t[0] *(this->t[4] * this->t[8] - this->t[5] * this->t[7]);
b -= this->t[3] *(this->t[1] * this->t[8] - this->t[2] * this->t[7]);
b += this->t[6] *(this->t[1] * this->t[5] - this->t[2] * this->t[4]);
return b;};
// test
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3HB::operator == ( const TenseurHB & B) const
{ int res = 1;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator == ( etc..");
#endif
for (int i = 0; i<=8; i++)
if (this->t[i] != B.t[i]) res = 0 ;
return res; };
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3HB::operator != ( const TenseurHB & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HB::operator != ( etc..");
#endif
if ((*this) == B)
return 0;
else
return 1; };
// calcul du tenseur inverse par rapport au produit contracte
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::Inverse() const
{TenseurHB * res;
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
// choix sur la méthode d'inversion
switch (ParaGlob::param->ParaAlgoControleActifs().Type_calnum_inversion_metrique())
{ case LU_EQUILIBRE:
{ // on recopie this dans le nouveau tenseur
for (int i = 0; i< 9; i++)
res->t[i] = t[i];
// pour le débug
//res->t[0]=3.; res->t[1]=2.;res->t[2]=1.;
// appel de l'inversion
Util::Inverse_mat3x3(((Tenseur3HB *) res)->ipointe);
}
break;
//cout << "\n comp \n ";
// res->Ecriture(cout); cout << "\n";
case CRAMER : // méthode historique
{ // calcul du determinant
double det = Det();
#ifdef MISE_AU_POINT
if (Dabs(det) <= ConstMath::trespetit)
{ cout << "\nErreur : le determinant du tenseur est nul !\n";
cout << "Tenseur3HB::Inverse() \n";
Sortie(1);
};
#endif
det =1./det;
res->t[0] = (this->t[4]*this->t[8] - this->t[5]*this->t[7])*det;
res->t[3] = (this->t[5]*this->t[6] - this->t[3]*this->t[8])*det;
res->t[6] = (this->t[3]*this->t[7] - this->t[4]*this->t[6])*det;
res->t[1] = (this->t[2]*this->t[7] - this->t[1]*this->t[8])*det;
res->t[4] = (this->t[0]*this->t[8] - this->t[2]*this->t[6])*det;
res->t[7] = (this->t[1]*this->t[6] - this->t[0]*this->t[7])*det;
res->t[2] = (this->t[1]*this->t[5] - this->t[2]*this->t[4])*det;
res->t[5] = (this->t[2]*this->t[3] - this->t[0]*this->t[5])*det;
res->t[8] = (this->t[0]*this->t[4] - this->t[1]*this->t[3])*det;
}
break;
default:
{ cout << "\nErreur **** : la methode de resolution de l'inversion de tenseur "
<< ParaGlob::param->ParaAlgoControleActifs().Type_calnum_inversion_metrique() << " n'est pas implante \n";
cout << "Tenseur3BH::Inverse() \n";
Sortie(1);
};
break;
};
// res->Ecriture(cout); // pour le debug
return *res;
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3HB::Transpose() const
{ TenseurBH * res;
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon(res); // ajout d'un tenseur intermediaire
res->t[0] = this->t[0];
res->t[1] = this->t[3];
res->t[2] = this->t[6];
res->t[3] = this->t[1];
res->t[4] = this->t[4];
res->t[5] = this->t[7];
res->t[6] = this->t[2];
res->t[7] = this->t[5];
res->t[8] = this->t[8];
return *res;};
#ifndef MISE_AU_POINT
inline
#endif
// permute Bas Haut, mais reste dans le même tenseur
void Tenseur3HB::PermuteHautBas()
{ // t[0] = t[0];
double t1=t[1]; t[1] = t[3];
double t2=t[2]; t[2] = t[6];
t[3] = t1;
//[4] = [4];
double t5=t[5]; t[5] = t[7];
t[6] = t2;
t[7] = t5;
//t[8] = t[8];
};
// calcul du maximum en valeur absolu des composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HB::MaxiComposante() const
{ return DabsMaxiTab(t, 9) ;
};
// retourne la composante i,j en lecture et écriture
#ifndef MISE_AU_POINT
inline
#endif
double& Tenseur3HB::Coor( const int i, const int j)
{
#ifdef MISE_AU_POINT
if ( ((i!=1)&&(i!=2)&&(i!=3)) || ((j!=1)&&(j!=2)&&(j!=3)) )
{ cout << "\nErreur : composante inexistante !\n";
cout << " i = " << i << ", j = " << j << '\n';
cout << "TenseurHH::Coor(int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[1]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[4]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[6]; break;
case 2 : return t[7]; break;
case 3 : return t[8]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
// retourne la composante i,j en lecture seulement
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HB::operator () ( const int i, const int j) const
{
#ifdef MISE_AU_POINT
if ( ((i!=1)&&(i!=2)&&(i!=3)) || ((j!=1)&&(j!=2)&&(j!=3)) )
{ cout << "\nErreur : composante inexistante !\n";
cout << " i = " << i << ", j = " << j << '\n';
cout << "TenseurHH::OPERATOR() (int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[1]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[4]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[6]; break;
case 2 : return t[7]; break;
case 3 : return t[8]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
//fonctions static définissant le produit tensoriel de deux vecteurs
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HB::Prod_tensoriel(const CoordonneeH & aH, const CoordonneeB & bB)
{ TenseurHB * res;
#ifdef MISE_AU_POINT
if ((aH.Dimension() != 3) || (bB.Dimension() != 3))
{ cout << "\n erreur de dimension dans les coordonnées d'entrée, dim1 et dim2 ="
<< aH.Dimension() << " " << bB.Dimension()
<< "\n Tenseur2HB::Prod_tensoriel( etc.." << endl;
Sortie(1);
}
#endif
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = aH(1) * bB(1); res->t[1] = aH(1) * bB(2); res->t[2] = aH(1) * bB(3);
res->t[3] = aH(2) * bB(1); res->t[4] = aH(2) * bB(2); res->t[5] = aH(2) * bB(3);
res->t[6] = aH(3) * bB(1); res->t[7] = aH(3) * bB(2); res->t[8] = aH(3) * bB(3);
return *res ;};
// lecture et écriture de données
#ifndef MISE_AU_POINT
inline
#endif
istream & Tenseur3HB::Lecture(istream & entree)
{ // lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3HB")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i< 9; i++)
entree >> this->t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
ostream & Tenseur3HB::Ecriture(ostream & sort) const
{ // écriture du type
sort << "Tenseur3HB ";
// puis les datas
for (int i = 0; i< 9; i++)
sort << setprecision(ParaGlob::NbdigdoCA()) << this->t[i] << " ";
return sort;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator de lecture
istream & operator >> (istream & entree, Tenseur3HB & A)
{ int dim = A.Dimension();
#ifdef MISE_AU_POINT
if (dim != 3) A.Message(3,"operator >> (istream & entree, Tenseur3HB & A)");
#endif
// lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3HB")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i< 9; i++)
entree >> A.t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator d'ecriture
ostream & operator << (ostream & sort , const Tenseur3HB & A)
{ //int dim = A.Dimension();
// écriture du type
sort << "Tenseur3HB ";
// puis les datas
for (int i = 0; i< 9; i++)
sort << setprecision(ParaGlob::NbdigdoCA()) << A.t[i] << " ";
return sort;
};
#endif
Reference in a new issue View git blame Copy permalink