// 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"
#ifndef Tenseur3_H_deja_inclus
// variables globales
// initialisation dans EnteteTenseur.h , utilisé dans le progr principal
//------------------------------------------------------------------
// cas des composantes deux fois contravariantes
//------------------------------------------------------------------
// --- gestion de changement d'index ----
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::ChangementIndex::ChangementIndex() :
idx_i(6),idx_j(6),odVect(3)
// ici pour un stockage triangle inférieur
{ idx_i(1)=1;idx_i(2)=2;idx_i(3)=3;idx_i(4)=2;idx_i(5)=3;idx_i(6)=3;
idx_j(1)=1;idx_j(2)=2;idx_j(3)=3;idx_j(4)=1;idx_j(5)=2;idx_j(6)=1;
odVect(1,1)=1;odVect(1,2)=4;odVect(1,3)=6;
odVect(2,1)=4;odVect(2,2)=2;odVect(2,3)=5;
odVect(3,1)=6;odVect(3,2)=5;odVect(3,3)=3;
};
// Constructeur
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::Tenseur3HH(): // par défaut
ipointe()
{ dimension = 3;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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] = 0.; t[1] = 0.; t[2] = 0.;
t[3] = 0.; t[4] = 0.; t[5] = 0.;
};
// initialisation de toutes les composantes a une meme valeur
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::Tenseur3HH( const double val) :
ipointe()
{ dimension = 3;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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] =val; t[1] =val; t[2] =val;
t[3] =val; t[4] =val; t[5] =val;
};
// initialisation avec 6 valeurs différentes
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::Tenseur3HH
( const double val1, const double val2, const double val3,
const double val4, const double val5, const double val6) :
ipointe()
{ dimension = 3;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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;
};
// DESTRUCTEUR :
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::~Tenseur3HH()
{ //if(listdouble6.end() != listdouble6.begin()) // si la liste n'est pas vide
listdouble6.erase(ipointe);} ; // suppression de l'élément de la liste
// constructeur a partir d'une instance non differenciee
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::Tenseur3HH ( const TenseurHH & B):
ipointe()
{ dimension = B.dimension;
#ifdef MISE_AU_POINT
if (Dabs(dimension) != 3)
{ cout << "\n erreur de dimension, elle devrait etre = 3 ";
cout << "\n Tenseur3HH::Tenseur3HH ( TenseurHH &) " << endl;
Sortie(1);
}
#endif
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
if (dimension == 3) // cas symetrique
for (int i=0;i< 6;i++)
t[i] = B.t[i];
else // cas non symetrique de stockage
// pour les 3D cela ne signifie pas : non symetrique car du aux couts de la verif
// de symetrie, on transforme systematiquement en non symetrie, la verif est donc
// a faire ici
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if ((!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
|| (Abs(Z) < ConstMath::petit) )
#endif
// on a retrouve un tenseur symetrique
{ dimension = 3;
t[0] = B.t[0];
t[1] = B.t[4];
t[2] = B.t[8];
t[3] = 0.5*(B.t[3]+B.t[1]) ;
t[4] = 0.5*(B.t[7]+B.t[5]) ;
t[5] = 0.5*(B.t[6]+B.t[2]) ;
}
#ifdef MISE_AU_POINT
else // erreur d'affectation
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3HH&) B);
Message(3,
"erreur d\'affectation, tenseur non symetrique, Tenseur3HH::Tenseur3HH ( TenseurHH &)");
}
#endif
}
};
// constructeur de copie
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3HH::Tenseur3HH ( const Tenseur3HH & B):
ipointe()
{ this->dimension = B.dimension;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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< 6;i++)
this->t[i] = B.t[i];
};
// METHODES PUBLIQUES :
// initialise toutes les composantes à val
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HH::Inita(double val)
{ t[0] =val; t[1] =val; t[2] =val;
t[3] =val; t[4] =val; t[5] =val;
};
// operations
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator + ( const TenseurHH & B) const
{ TenseurHH * res;
if (B.Dimension() == -3) return (B + *this); // cas B non symetrique
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HH::operator + ( etc..");
#endif
res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] + B.t[i]; //somme des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HH::operator += ( const TenseurHH & B)
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HH::operator += ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
for (int i = 0; i< 6; i++)
this->t[i] += B.t[i];
else // cas non symetrique
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if ((!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
|| (Abs(Z) < ConstMath::petit) )
#endif
{ this->t[0] += B.t[0];
this->t[1] += B.t[4];
this->t[2] += B.t[8];
this->t[3] += 0.5*(B.t[3]+B.t[1]) ;
this->t[4] += 0.5*(B.t[7]+B.t[5]) ;
this->t[5] += 0.5*(B.t[6]+B.t[2]) ;
}
#ifdef MISE_AU_POINT
else
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3HH&) B);
Message(3,
"erreur d\'affectation, tenseur non symetrique, Tenseur3HH::operator += ( etc..");
}
#endif
}
LesMaillonsHH::Libere(); // destruction des tenseurs intermediaires
}; //somme des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator - () const
{ TenseurHH * res;
res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = - this->t[i]; //oppose
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator - ( const TenseurHH & B) const
{ TenseurHH * res;
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HH::operator - ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
{ res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] - B.t[i]; //soustraction des données
}
else // cas non symetrique
{ res = new Tenseur_ns3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = this->t[0] - B.t[0]; //soustraction des données
res->t[1] = this->t[3] - B.t[1]; //soustraction des données
res->t[2] = this->t[5] - B.t[2]; //soustraction des données
res->t[3] = this->t[3] - B.t[3]; //soustraction des données
res->t[4] = this->t[1] - B.t[4]; //soustraction des données
res->t[5] = this->t[4] - B.t[5]; //soustraction des données
res->t[6] = this->t[5] - B.t[6]; //soustraction des données
res->t[7] = this->t[4] - B.t[7]; //soustraction des données
res->t[8] = this->t[2] - B.t[8]; //soustraction des données
}
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HH::operator -= ( const TenseurHH & B)
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HH::operator -= ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
for (int i = 0; i< 6; i++)
this->t[i] -= B.t[i];
else // cas non symetrique
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if ((!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
|| (Abs(Z) < ConstMath::petit) )
#endif
{ this->t[0] -= B.t[0];
this->t[1] -= B.t[4];
this->t[2] -= B.t[8];
this->t[3] -= 0.5*(B.t[3]+B.t[1]) ;
this->t[4] -= 0.5*(B.t[7]+B.t[5]) ;
this->t[5] -= 0.5*(B.t[6]+B.t[2]) ;
}
#ifdef MISE_AU_POINT
else
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3HH&) B);
Message(3,
"erreur d\'affectation, tenseur non symetrique, Tenseur3HH::operator -= ( etc..");
}
#endif
}
LesMaillonsHH::Libere(); // destruction des tenseurs intermediaires
}; //soustraction des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator = ( const TenseurHH & B)
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(2,"Tenseur3HH::operator = ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
{ for (int i = 0; i< 6; i++)
this->t[i] = B.t[i];
}
else // cas non symetrique de stockage
// pour les 3D cela ne signifie pas : non symetrique car du aux couts de la verif
// de symetrie, on transforme systematiquement en non symetrie, la verif est donc
// a faire a l'affectation c-a-d ici
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if (!(!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
&& (Abs(Z) < ConstMath::petit) )
// erreur d'affectation
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3HH&) B);
cout << "WARNING ** erreur d\'affectation, tenseur non symetrique, Tenseur3HH::operator = ( etc..";
}
#endif
// on a retrouve un tenseur symetrique (ou il y a eu un message)
this->t[0] = B.t[0];
this->t[1] = B.t[4];
this->t[2] = B.t[8];
this->t[3] = 0.5*(B.t[3]+B.t[1]) ;
this->t[4] = 0.5*(B.t[7]+B.t[5]) ;
this->t[5] = 0.5*(B.t[6]+B.t[2]) ;
}
LesMaillonsHH::Libere(); // destruction des tenseurs intermediaires
return *this;
}; //affectation des données;
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator * ( const double & b) const
{ TenseurHH * res;
res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] * b; //multiplication des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HH::operator *= ( const double & b)
{for (int i = 0; i< 6; i++)
this->t[i] *= b ;}; //multiplication des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator / ( const double & b) const
{ TenseurHH * res;
res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] / b; //division des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3HH::operator /= ( const double & b)
{ for (int i = 0; i< 6; 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 Tenseur3HH::Affectation_2D_a_3D(const Tenseur2HH & B,bool plusZero)
{ this->t[0] = B.t[0];this->t[1] = B.t[1];this->t[3] = B.t[2];
if (plusZero)
{ this->t[2] = this->t[4] = this->t[5] = 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 Tenseur3HH::Affectation_trans_dimension(const TenseurHH & B,bool plusZero)
{ switch (B.Dimension())
{case 3: case -3: *this = B; break; // affectation normale
case 2:
{ const Tenseur2HH & bn = *((Tenseur2HH *) &B);
this->Affectation_2D_a_3D(bn,plusZero);
break;
}
case -2:
{ const Tenseur2HH bn = B; // on crée un nouveau tenseur transitoire
this->Affectation_2D_a_3D(bn,plusZero); // qui est tout de suite supprimé
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, Tenseur3HH::Affectation_trans_dimension( const TenseurHH & B, ..");
}
};
// produit contracte avec un vecteur
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeH Tenseur3HH::operator * ( const CoordonneeB & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != dimension)
{ cout << "\nErreur : dimensions vecteur tenseur non egales !\n";
cout << " Tenseur3HH::operator *\n";
Sortie(1);
};
#endif
CoordonneeH v(dimension);
v(1) = this->t[0] * B(1) + this->t[3] * B(2) + this->t[5] * B(3);
v(2) = this->t[3] * B(1) + this->t[1] * B(2) + this->t[4] * B(3);
v(3) = this->t[5] * B(1) + this->t[4] * B(2) + this->t[2] * B(3);
return v;
};
// produit contracte contracté une fois puis deux fois
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::operator * ( const TenseurBH & B) const
{ TenseurHH * res;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HH::operator * ( etc..");
#endif
res = new Tenseur_ns3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
// on cree systematiquement un tenseur non symetrique
res->t[0] = this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[1] * B.t[3] + this->t[4] * B.t[6];
res->t[6] = this->t[5] * B.t[0] + this->t[4] * B.t[3] + this->t[2] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[3] * B.t[4] + this->t[5] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[1] * B.t[4] + this->t[4] * B.t[7];
res->t[7] = this->t[5] * B.t[1] + this->t[4] * B.t[4] + this->t[2] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[3] * B.t[5] + this->t[5] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[1] * B.t[5] + this->t[4] * B.t[8];
res->t[8] = this->t[5] * B.t[2] + this->t[4] * B.t[5] + this->t[2] * B.t[8];
return *res; };
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB & Tenseur3HH::operator * ( const TenseurBB & B) const
{ TenseurHB * res;
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HH::operator * ( etc..");
#endif
res = new Tenseur3HB;
LesMaillonsHB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
if (B.Dimension() == 3) // cas symetrique
{ res->t[0] = this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[5];
res->t[3] = this->t[3] * B.t[0] + this->t[1] * B.t[3] + this->t[4] * B.t[5];
res->t[6] = this->t[5] * B.t[0] + this->t[4] * B.t[3] + this->t[2] * B.t[5];
res->t[1] = this->t[0] * B.t[3] + this->t[3] * B.t[1] + this->t[5] * B.t[4];
res->t[4] = this->t[3] * B.t[3] + this->t[1] * B.t[1] + this->t[4] * B.t[4];
res->t[7] = this->t[5] * B.t[3] + this->t[4] * B.t[1] + this->t[2] * B.t[4];
res->t[2] = this->t[0] * B.t[5] + this->t[3] * B.t[4] + this->t[5] * B.t[2];
res->t[5] = this->t[3] * B.t[5] + this->t[1] * B.t[4] + this->t[4] * B.t[2];
res->t[8] = this->t[5] * B.t[5] + this->t[4] * B.t[4] + this->t[2] * B.t[2];
}
else // cas ou B est en stockage non symetrique
{ res->t[0] = this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[1] * B.t[3] + this->t[4] * B.t[6];
res->t[6] = this->t[5] * B.t[0] + this->t[4] * B.t[3] + this->t[2] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[3] * B.t[4] + this->t[5] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[1] * B.t[4] + this->t[4] * B.t[7];
res->t[7] = this->t[5] * B.t[1] + this->t[4] * B.t[4] + this->t[2] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[3] * B.t[5] + this->t[5] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[1] * B.t[5] + this->t[4] * B.t[8];
res->t[8] = this->t[5] * B.t[2] + this->t[4] * B.t[5] + this->t[2] * B.t[8];
}
return *res; };
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HH::operator && ( const TenseurBB & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3HH::operator && ( etc..");
#endif
double b = 0;
if (B.Dimension() == 3) // cas ou B est en stockage symetrique
{ b += this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[5];
b += this->t[3] * B.t[3] + this->t[1] * B.t[1] + this->t[4] * B.t[4];
b += this->t[5] * B.t[5] + this->t[4] * B.t[4] + this->t[2] * B.t[2];
}
else // cas ou B est en stockage non symetrique
{ b += this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[6];
b += this->t[3] * B.t[1] + this->t[1] * B.t[4] + this->t[4] * B.t[7];
b += this->t[5] * B.t[2] + this->t[4] * B.t[5] + this->t[2] * B.t[8];
}
return b; };
// test
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3HH::operator == ( const TenseurHH & B) const
{ int res = 1;
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HH::operator == ( etc..");
#endif
if (B.Dimension() == 3)
{ for (int i = 0; i< 6; i++)
if (this->t[i] != B.t[i]) res = 0 ;
return res;
}
else
return (B == *this);
};
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3HH::operator != ( const TenseurHH & B) const
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3HH::operator != ( etc..");
#endif
if ((*this) == B)
return 0;
else
return 1; };
// determinant de la matrice des coordonnees
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HH::Det() const
{ double b = 0;
b += this->t[0] * (this->t[1] * this->t[2] - this->t[4] * this->t[4]) ;
b -= this->t[3] * (this->t[3] * this->t[2] - this->t[5] * this->t[4]) ;
b += this->t[5] * (this->t[3] * this->t[4] - this->t[5] * this->t[1]) ;
return b;
};
// calcul du tenseur inverse par rapport au produit contracte
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3HH::Inverse() const
{TenseurBB * res;
res = new Tenseur3BB;
LesMaillonsBB::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< 6; i++)
res->t[i] = t[i];
// appel de l'inversion
Util::Inverse_mat3x3(((Tenseur3BB*) res)->ipointe);
}
break;
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 << "Tenseur3HH::Inverse() \n";
Sortie(1);
};
#endif
det =1./det;
res->t[0] = (this->t[1]*this->t[2] - this->t[4]*this->t[4])*det;
res->t[3] = (this->t[5]*this->t[4] - this->t[3]*this->t[2])*det;
res->t[5] = (this->t[3]*this->t[4] - this->t[5]*this->t[1])*det;
res->t[1] = (this->t[0]*this->t[2] - this->t[5]*this->t[5])*det;
res->t[4] = (this->t[3]*this->t[5] - this->t[4]*this->t[0])*det;
res->t[2] = (this->t[0]*this->t[1] - this->t[3]*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 << "Tenseur3HH::Inverse() \n";
Sortie(1);
};
break;
};
return *res;
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::Transpose() const
{ TenseurHH * res;
res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon(res); // ajout d'un tenseur intermediaire
// *(res->t) = *(this->t); // transposé d'un tenseur symétrique = idem
res->t[0] = this->t[0]; // transposé d'un tenseur symétrique = idem
res->t[1] = this->t[1];
res->t[2] = this->t[2];
res->t[3] = this->t[3];
res->t[4] = this->t[4];
res->t[5] = this->t[5];
return *res;};
// ---- manipulation d'indice ---- -> création de nouveaux tenseurs
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB& Tenseur3HH::Baisse2Indices() const
{ TenseurBB * res;
res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon(res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i];
return *res;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH& Tenseur3HH::BaissePremierIndice() 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[5];
res->t[3] = this->t[3];
res->t[4] = this->t[1];
res->t[5] = this->t[4];
res->t[6] = this->t[5];
res->t[7] = this->t[4];
res->t[8] = this->t[2];
return *res;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB& Tenseur3HH::BaisseDernierIndice() 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[5];
res->t[3] = this->t[3];
res->t[4] = this->t[1];
res->t[5] = this->t[4];
res->t[6] = this->t[5];
res->t[7] = this->t[4];
res->t[8] = this->t[2];
return *res;};
// calcul du maximum en valeur absolu des composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3HH::MaxiComposante() const
{ return DabsMaxiTab(t, 6) ;
};
// retourne la composante i,j en lecture et écriture
#ifndef MISE_AU_POINT
inline
#endif
double& Tenseur3HH::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 << "Tenseur3HH::Coor(int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[3]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[1]; break;
case 3 : return t[4]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[5]; break;
case 2 : return t[4]; break;
case 3 : return t[2]; 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 Tenseur3HH::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 << "Tenseur3HH::OPERATOR() (int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[3]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[1]; break;
case 3 : return t[4]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[5]; break;
case 2 : return t[4]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
//fonctions static définissant le produit tensoriel de deux vecteurs
// si les vecteurs sont égaux le tenseur est symétrique sinon il est non symétrique
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3HH::Prod_tensoriel(const CoordonneeH & aH, const CoordonneeH & bH)
{ TenseurHH * res;
#ifdef MISE_AU_POINT
if ((aH.Dimension() != 3) || (bH.Dimension() != 3))
{ cout << "\n erreur de dimension dans les coordonnees d'entree, dim1 et dim2 ="
<< aH.Dimension() << " " << bH.Dimension()
<< "\n Tenseur2HH::Prod_tensoriel( etc.." << endl;
Sortie(1);
}
#endif
if (aH == bH) // cas d'un résultat symetrique
{res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = aH(1) * aH(1);
res->t[3] = aH(2) * aH(1); res->t[1] = aH(2) * aH(2);
res->t[5] = aH(3) * aH(1); res->t[4] = aH(3) * aH(2); res->t[2] = aH(3) * aH(3);
}
else // cas d'un résultat non symétrique
{res = new Tenseur_ns3HH;
LesMaillonsHH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = aH(1) * bH(1); res->t[1] = aH(1) * bH(2); res->t[2] = aH(1) * bH(3);
res->t[3] = aH(2) * bH(1); res->t[4] = aH(2) * bH(2); res->t[5] = aH(2) * bH(3);
res->t[6] = aH(3) * bH(1); res->t[7] = aH(3) * bH(2); res->t[8] = aH(3) * bH(3);
}
return *res ;};
// lecture et écriture de données
#ifndef MISE_AU_POINT
inline
#endif
istream & Tenseur3HH::Lecture(istream & entree)
{ // lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3HH")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i< 6; i++)
entree >> this->t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
ostream & Tenseur3HH::Ecriture(ostream & sort) const
{ // écriture du type
sort << "Tenseur3HH ";
int nb_diggit = ParaGlob::NbdigdoCA();
// puis les datas
for (int i = 0; i < 6; i++)
sort << setprecision(nb_diggit) << this->t[i] << " ";
return sort;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator de lecture
istream & operator >> (istream & entree, Tenseur3HH & A)
{ int dim = A.Dimension();
#ifdef MISE_AU_POINT
if (dim != 3) A.Message(3,"operator >> (istream & entree, Tenseur3HH & A)");
#endif
// lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3HH")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i< 6; i++)
entree >> A.t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator d'ecriture
ostream & operator << (ostream & sort , const Tenseur3HH & A)
{ //int dim = A.Dimension();
// écriture du type
sort << "Tenseur3HH ";
int nb_diggit = ParaGlob::NbdigdoCA();
// puis les datas
for (int i = 0; i<= 5; i++)
sort << setprecision(nb_diggit) << A.t[i] << " ";
return sort;
};
//------------------------------------------------------------------
// cas des composantes deux fois covariantes
//------------------------------------------------------------------
// --- gestion de changement d'index ----
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::ChangementIndex::ChangementIndex() :
idx_i(6),idx_j(6),odVect(3)
// ici pour un stockage triangle inférieur
{ idx_i(1)=1;idx_i(2)=2;idx_i(3)=3;idx_i(4)=2;idx_i(5)=3;idx_i(6)=3;
idx_j(1)=1;idx_j(2)=2;idx_j(3)=3;idx_j(4)=1;idx_j(5)=2;idx_j(6)=1;
odVect(1,1)=1;odVect(1,2)=4;odVect(1,3)=6;
odVect(2,1)=4;odVect(2,2)=2;odVect(2,3)=5;
odVect(3,1)=6;odVect(3,2)=5;odVect(3,3)=3;
};
// Constructeur
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::Tenseur3BB() : // par défaut
ipointe()
{ dimension = 3;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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] = 0.; t[1] = 0.; t[2] = 0.;
t[3] = 0.; t[4] = 0.; t[5] = 0.;
};
// initialisation de toutes les composantes a une meme valeur
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::Tenseur3BB( const double val):
ipointe()
{ dimension = 3;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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] =val; t[1] =val; t[2] =val;
t[3] =val; t[4] =val; t[5] =val;
};
// initialisation avec 6 valeurs différentes
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::Tenseur3BB
( const double val1, const double val2, const double val3,
const double val4, const double val5, const double val6) :
ipointe()
{ dimension = 3;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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;
};
// DESTRUCTEUR :
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::~Tenseur3BB()
{//if(listdouble6.end() != listdouble6.begin()) // si la liste n'est pas vide
listdouble6.erase(ipointe);} ; // suppression de l'élément de la liste
// constructeur a partir d'une instance non differenciee
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::Tenseur3BB ( const TenseurBB & B) :
ipointe()
{ this->dimension = B.dimension;
#ifdef MISE_AU_POINT
if (Dabs(dimension) != 3)
{ cout << "\n erreur de dimension, elle devrait etre = 3 ";
cout << "\n Tenseur3BB::Tenseur3BB ( TenseurBB &) " << endl;
Sortie(1);
}
#endif
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.begin(); // recup de la position de la maille dans la liste
t = (ipointe)->donnees; // recup de la position des datas dans la maille
if (dimension == 3) // cas symetrique
for (int i=0;i< 6;i++)
t[i] = B.t[i];
else // cas non symetrique de stockage
// pour les 3D cela ne signifie pas : non symetrique car du aux couts de la verif
// de symetrie, on transforme systematiquement en non symetrie, la verif est donc
// a faire ici
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if ((!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
|| (Abs(Z) < ConstMath::petit) )
#endif
// on a retrouve un tenseur symetrique
{ dimension = 3;
t[0] = B.t[0];
t[1] = B.t[4];
t[2] = B.t[8];
t[3] = 0.5*(B.t[3]+B.t[1]) ;
t[4] = 0.5*(B.t[7]+B.t[5]) ;
t[5] = 0.5*(B.t[6]+B.t[2]) ;
}
#ifdef MISE_AU_POINT
else // erreur d'affectation
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3BB&) B);
Message(3,
"erreur d\'affectation, tenseur non symetrique, Tenseur3BB::Tenseur3BB ( TenseurBB &)");
}
#endif
}
};
// constructeur de copie
#ifndef MISE_AU_POINT
inline
#endif
Tenseur3BB::Tenseur3BB ( const Tenseur3BB & B):
ipointe()
{ this->dimension = B.dimension;
listdouble6.push_front(Reels6()); // allocation
ipointe = listdouble6.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< 6;i++)
this->t[i] = B.t[i];
};
// METHODES PUBLIQUES :
// initialise toutes les composantes à val
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BB::Inita(double val)
{ t[0] =val; t[1] =val; t[2] =val;
t[3] =val; t[4] =val; t[5] =val;
};
// operations
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator + ( const TenseurBB & B) const
{ TenseurBB * res;
if (B.Dimension() == -3) return (B + *this); // cas B non symetrique
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BB::operator + ( etc..");
#endif
res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] + B.t[i]; //somme des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BB::operator += ( const TenseurBB & B)
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BB::operator += ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
for (int i = 0; i< 6; i++)
this->t[i] += B.t[i];
else // cas non symetrique
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if ((!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
|| (Abs(Z) < ConstMath::petit) )
#endif
{ this->t[0] += B.t[0];
this->t[1] += B.t[4];
this->t[2] += B.t[8];
this->t[3] += 0.5*(B.t[3]+B.t[1]) ;
this->t[4] += 0.5*(B.t[7]+B.t[5]) ;
this->t[5] += 0.5*(B.t[6]+B.t[2]) ;
}
#ifdef MISE_AU_POINT
else
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3BB&) B);
Message(3,
"erreur d\'affectation, tenseur non symetrique, Tenseur3BB::operator += ( etc..");
}
#endif
}
LesMaillonsHH::Libere(); // destruction des tenseurs intermediaires
}; //somme des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator - () const
{ TenseurBB * res;
res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = - this->t[i]; //oppose des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator - ( const TenseurBB & B) const
{ TenseurBB * res;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BB::operator - ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
{ res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] - B.t[i]; //soustraction des données
}
else // cas non symetrique
{ res = new Tenseur_ns3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = this->t[0] - B.t[0]; //soustraction des données
res->t[1] = this->t[3] - B.t[1]; //soustraction des données
res->t[2] = this->t[5] - B.t[2]; //soustraction des données
res->t[3] = this->t[3] - B.t[3]; //soustraction des données
res->t[4] = this->t[1] - B.t[4]; //soustraction des données
res->t[5] = this->t[4] - B.t[5]; //soustraction des données
res->t[6] = this->t[5] - B.t[6]; //soustraction des données
res->t[7] = this->t[4] - B.t[7]; //soustraction des données
res->t[8] = this->t[2] - B.t[8]; //soustraction des données
}
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BB::operator -= ( const TenseurBB & B)
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BB::operator -= ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
for (int i = 0; i< 6; i++)
this->t[i] -= B.t[i];
else // cas non symetrique
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if ((!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
|| (Abs(Z) < ConstMath::petit) )
#endif
{ this->t[0] -= B.t[0];
this->t[1] -= B.t[4];
this->t[2] -= B.t[8];
this->t[3] -= 0.5*(B.t[3]+B.t[1]) ;
this->t[4] -= 0.5*(B.t[7]+B.t[5]) ;
this->t[5] -= 0.5*(B.t[6]+B.t[2]) ;
}
#ifdef MISE_AU_POINT
else
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3BB&) B);
Message(3,
"erreur d\'affectation, tenseur non symetrique, Tenseur3BB::operator -= ( etc..");
}
#endif
}
LesMaillonsHH::Libere(); // destruction des tenseurs intermediaires
}; //soustraction des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator = ( const TenseurBB & B)
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BB::operator = ( etc..");
#endif
if (B.Dimension() == 3) // cas symetrique
{ for (int i = 0; i< 6; i++)
this->t[i] = B.t[i];
}
else // cas non symetrique de stockage
// pour les 3D cela ne signifie pas : non symetrique car du aux couts de la verif
// de symetrie, on transforme systematiquement en non symetrie, la verif est donc
// a faire a l'affectation c-a-d ici
{
#ifdef MISE_AU_POINT
double Z = B.MaxiComposante();
if (!(!diffpourcent(B(2,1),B(1,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(2,3),B(3,2),Z,ConstMath::unpeupetit)
&& !diffpourcent(B(3,1),B(1,3),Z,ConstMath::unpeupetit))
&& (Abs(Z) < ConstMath::petit) )
{ if (ParaGlob::NiveauImpression() > 5)
cout << "\n tenseur = " << ((Tenseur_ns3BB&) B);
cout << "WARNING ** erreur d\'affectation, tenseur non symetrique, Tenseur3BB::operator = ( etc..";
}
#endif
// on a retrouve un tenseur symetrique (ou il y a eu un message)
this->t[0] = B.t[0];
this->t[1] = B.t[4];
this->t[2] = B.t[8];
this->t[3] = 0.5*(B.t[3]+B.t[1]) ;
this->t[4] = 0.5*(B.t[7]+B.t[5]) ;
this->t[5] = 0.5*(B.t[6]+B.t[2]) ;
}
LesMaillonsBB::Libere(); // destruction des tenseurs intermediaires
return *this; }; //affectation des données;
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator * ( const double & b) const
{ TenseurBB * res;
res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] * b; //multiplication des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BB::operator *= ( const double & b)
{ for (int i = 0; i< 6; i++)
this->t[i] *= b ;}; //multiplication des données
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator / ( const double & b) const
{ TenseurBB * res;
res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i] / b; //division des données
return *res ;};
#ifndef MISE_AU_POINT
inline
#endif
void Tenseur3BB::operator /= ( const double & b)
{ for (int i = 0; i< 6; 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 Tenseur3BB::Affectation_2D_a_3D(const Tenseur2BB & B,bool plusZero)
{ this->t[0] = B.t[0];this->t[1] = B.t[1];this->t[3] = B.t[2];
if (plusZero)
{ this->t[2] = this->t[4] = this->t[5] = 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 Tenseur3BB::Affectation_trans_dimension(const TenseurBB & B,bool plusZero)
{ switch (B.Dimension())
{case 3: case -3: *this = B; break; // affectation normale
case 2:
{ const Tenseur2BB & bn = *((Tenseur2BB *) &B);
this->Affectation_2D_a_3D(bn,plusZero);
break;
}
case -2:
{ const Tenseur2BB bn = B; // on crée un nouveau tenseur transitoire
this->Affectation_2D_a_3D(bn,plusZero); // qui est tout de suite supprimé
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, Tenseur3BB::Affectation_trans_dimension( const TenseurBB & B, ..");
}
};
// produit contracte avec un vecteur
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeB Tenseur3BB::operator * ( const CoordonneeH & B) const
{
#ifdef MISE_AU_POINT
if (B.Dimension() != dimension)
{ cout << "\nErreur : dimensions vecteur tenseur non egales !\n";
cout << " Tenseur3BB::operator *\n";
Sortie(1);
};
#endif
CoordonneeB v(dimension);
v(1) = this->t[0] * B(1) + this->t[3] * B(2) + this->t[5] * B(3);
v(2) = this->t[3] * B(1) + this->t[1] * B(2) + this->t[4] * B(3);
v(3) = this->t[5] * B(1) + this->t[4] * B(2) + this->t[2] * B(3);
return v;
};
// produit contracte contracté une fois puis deux fois
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::operator * ( const TenseurHB & B) const
{ TenseurBB * res;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BB::operator * ( etc..");
#endif
res = new Tenseur_ns3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
// on cree systematiquement un tenseur non symetrique
res->t[0] = this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[1] * B.t[3] + this->t[4] * B.t[6];
res->t[6] = this->t[5] * B.t[0] + this->t[4] * B.t[3] + this->t[2] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[3] * B.t[4] + this->t[5] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[1] * B.t[4] + this->t[4] * B.t[7];
res->t[7] = this->t[5] * B.t[1] + this->t[4] * B.t[4] + this->t[2] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[3] * B.t[5] + this->t[5] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[1] * B.t[5] + this->t[4] * B.t[8];
res->t[8] = this->t[5] * B.t[2] + this->t[4] * B.t[5] + this->t[2] * B.t[8];
return *res; };
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH & Tenseur3BB::operator * ( const TenseurHH & B) const
{ TenseurBH * res;
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BB::operator * ( etc..");
#endif
res = new Tenseur3BH;
LesMaillonsBH::NouveauMaillon( res); // ajout d'un tenseur intermediaire
if (B.Dimension() == 3) // cas symetrique
{ res->t[0] = this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[5];
res->t[3] = this->t[3] * B.t[0] + this->t[1] * B.t[3] + this->t[4] * B.t[5];
res->t[6] = this->t[5] * B.t[0] + this->t[4] * B.t[3] + this->t[2] * B.t[5];
res->t[1] = this->t[0] * B.t[3] + this->t[3] * B.t[1] + this->t[5] * B.t[4];
res->t[4] = this->t[3] * B.t[3] + this->t[1] * B.t[1] + this->t[4] * B.t[4];
res->t[7] = this->t[5] * B.t[3] + this->t[4] * B.t[1] + this->t[2] * B.t[4];
res->t[2] = this->t[0] * B.t[5] + this->t[3] * B.t[4] + this->t[5] * B.t[2];
res->t[5] = this->t[3] * B.t[5] + this->t[1] * B.t[4] + this->t[4] * B.t[2];
res->t[8] = this->t[5] * B.t[5] + this->t[4] * B.t[4] + this->t[2] * B.t[2];
}
else // cas ou B est en stockage non symetrique
{ res->t[0] = this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[6];
res->t[3] = this->t[3] * B.t[0] + this->t[1] * B.t[3] + this->t[4] * B.t[6];
res->t[6] = this->t[5] * B.t[0] + this->t[4] * B.t[3] + this->t[2] * B.t[6];
res->t[1] = this->t[0] * B.t[1] + this->t[3] * B.t[4] + this->t[5] * B.t[7];
res->t[4] = this->t[3] * B.t[1] + this->t[1] * B.t[4] + this->t[4] * B.t[7];
res->t[7] = this->t[5] * B.t[1] + this->t[4] * B.t[4] + this->t[2] * B.t[7];
res->t[2] = this->t[0] * B.t[2] + this->t[3] * B.t[5] + this->t[5] * B.t[8];
res->t[5] = this->t[3] * B.t[2] + this->t[1] * B.t[5] + this->t[4] * B.t[8];
res->t[8] = this->t[5] * B.t[2] + this->t[4] * B.t[5] + this->t[2] * B.t[8];
}
return *res; };
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BB::operator && ( const TenseurHH & B) const
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BB::operator && ( etc..");
#endif
double b = 0;
if (B.Dimension() == 3) // cas ou B est en stockage symetrique
{ b += this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[5];
b += this->t[3] * B.t[3] + this->t[1] * B.t[1] + this->t[4] * B.t[4];
b += this->t[5] * B.t[5] + this->t[4] * B.t[4] + this->t[2] * B.t[2];
}
else // cas ou B est en stockage non symetrique
{ b += this->t[0] * B.t[0] + this->t[3] * B.t[3] + this->t[5] * B.t[6];
b += this->t[3] * B.t[1] + this->t[1] * B.t[4] + this->t[4] * B.t[7];
b += this->t[5] * B.t[2] + this->t[4] * B.t[5] + this->t[2] * B.t[8];
}
return b; };
// test
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3BB::operator == ( const TenseurBB & B) const
{ int res = 1;
#ifdef MISE_AU_POINT
if (B.Dimension() != 3) Message(3,"Tenseur3BB::operator == ( etc..");
#endif
if (B.Dimension() == 3)
{ for (int i = 0; i< 6; i++)
if (this->t[i] != B.t[i]) res = 0 ;
return res;
}
else
return (B == *this);
};
#ifndef MISE_AU_POINT
inline
#endif
int Tenseur3BB::operator != ( const TenseurBB & B) const
{
#ifdef MISE_AU_POINT
if (Dabs(B.Dimension()) != 3) Message(3,"Tenseur3BB::operator != ( etc..");
#endif
if ((*this) == B)
return 0;
else
return 1; };
// calcul du determinant de la matrice des coordonnees
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BB::Det() const
{ double b = 0;
b += this->t[0] * (this->t[1] * this->t[2] - this->t[4] * this->t[4]) ;
b -= this->t[3] * (this->t[3] * this->t[2] - this->t[5] * this->t[4]) ;
b += this->t[5] * (this->t[3] * this->t[4] - this->t[5] * this->t[1]) ;
return b;
};
// calcul du tenseur inverse par rapport au produit contracte
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH & Tenseur3BB::Inverse() const
{TenseurHH * res;
res = new Tenseur3HH;
LesMaillonsHH::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< 6; i++)
res->t[i] = t[i];
// pour le débug
//res->t[0]=3.; res->t[1]=2.;res->t[2]=1.;res->t[3]=1.;res->t[4]=0.;res->t[5]=0.;
// appel de l'inversion
Util::Inverse_mat3x3(((Tenseur3HH*) 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 << "Tenseur3BB::Inverse() \n";
Sortie(1);
};
#endif
det =1./det;
res->t[0] = (this->t[1]*this->t[2] - this->t[4]*this->t[4])*det;
res->t[3] = (this->t[5]*this->t[4] - this->t[3]*this->t[2])*det;
res->t[5] = (this->t[3]*this->t[4] - this->t[5]*this->t[1])*det;
res->t[1] = (this->t[0]*this->t[2] - this->t[5]*this->t[5])*det;
res->t[4] = (this->t[3]*this->t[5] - this->t[4]*this->t[0])*det;
res->t[2] = (this->t[0]*this->t[1] - this->t[3]*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 << "Tenseur3BB::Inverse() \n";
Sortie(1);
};
break;
};
// res->Ecriture(cout); // pour le debug
return *res;
};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::Transpose() const
{ TenseurBB * res;
res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon(res); // ajout d'un tenseur intermediaire
// *(res->t) = *(this->t); // transposé d'un tenseur symétrique = idem
res->t[0] = this->t[0]; // transposé d'un tenseur symétrique = idem
res->t[1] = this->t[1];
res->t[2] = this->t[2];
res->t[3] = this->t[3];
res->t[4] = this->t[4];
res->t[5] = this->t[5];
return *res;};
// ---- manipulation d'indice ---- -> création de nouveaux tenseurs
#ifndef MISE_AU_POINT
inline
#endif
TenseurHH& Tenseur3BB::Monte2Indices() const
{ TenseurHH * res;
res = new Tenseur3HH;
LesMaillonsHH::NouveauMaillon(res); // ajout d'un tenseur intermediaire
for (int i = 0; i< 6; i++)
res->t[i] = this->t[i];
return *res;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurBH& Tenseur3BB::MonteDernierIndice() 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[5];
res->t[3] = this->t[3];
res->t[4] = this->t[1];
res->t[5] = this->t[4];
res->t[6] = this->t[5];
res->t[7] = this->t[4];
res->t[8] = this->t[2];
return *res;};
#ifndef MISE_AU_POINT
inline
#endif
TenseurHB& Tenseur3BB::MontePremierIndice() 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[5];
res->t[3] = this->t[3];
res->t[4] = this->t[1];
res->t[5] = this->t[4];
res->t[6] = this->t[5];
res->t[7] = this->t[4];
res->t[8] = this->t[2];
return *res;};
// calcul du maximum en valeur absolu des composantes du tenseur
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BB::MaxiComposante() const
{ return DabsMaxiTab(t, 6) ;
};
// retourne la composante i,j en lecture et écriture
#ifndef MISE_AU_POINT
inline
#endif
double& Tenseur3BB::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 << "Tenseur3BB::Coor(int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[3]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[1]; break;
case 3 : return t[4]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[5]; break;
case 2 : return t[4]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
// retourne la composante i,j en lecture uniquement
#ifndef MISE_AU_POINT
inline
#endif
double Tenseur3BB::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 << "Tenseur3BB::OPERATOR() (int,int ) \n";
Sortie(1);
};
#endif
switch (i)
{ case 1 : { switch (j)
{ case 1 : return t[0]; break;
case 2 : return t[3]; break;
case 3 : return t[5]; break;
default : return t[0]; }
break;}
case 2 : { switch (j)
{ case 1 : return t[3]; break;
case 2 : return t[1]; break;
case 3 : return t[4]; break;
default : return t[0]; }
break;}
case 3 : { switch (j)
{ case 1 : return t[5]; break;
case 2 : return t[4]; break;
case 3 : return t[2]; break;
default : return t[0]; }
break;}
default : return t[0];
}
};
//fonctions static définissant le produit tensoriel de deux vecteurs
// si les vecteurs sont égaux le tenseur est symétrique sinon il est non symétrique
#ifndef MISE_AU_POINT
inline
#endif
TenseurBB & Tenseur3BB::Prod_tensoriel(const CoordonneeB & aB, const CoordonneeB & bB)
{ TenseurBB * res;
#ifdef MISE_AU_POINT
if ((aB.Dimension() != 3) || (bB.Dimension() != 3))
{ cout << "\n erreur de dimension dans les coordonnees d'entree, dim1 et dim2 ="
<< aB.Dimension() << " " << bB.Dimension()
<< "\n Tenseur2BB::Prod_tensoriel( etc.." << endl;
Sortie(1);
}
#endif
if (aB == bB) // cas d'un résultat symetrique
{res = new Tenseur3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = aB(1) * aB(1);
res->t[3] = aB(2) * aB(1); res->t[1] = aB(2) * aB(2);
res->t[5] = aB(3) * aB(1); res->t[4] = aB(3) * aB(2); res->t[2] = aB(3) * aB(3);
}
else // cas d'un résultat non symétrique
{res = new Tenseur_ns3BB;
LesMaillonsBB::NouveauMaillon( res); // ajout d'un tenseur intermediaire
res->t[0] = aB(1) * bB(1); res->t[1] = aB(1) * bB(2); res->t[2] = aB(1) * bB(3);
res->t[3] = aB(2) * bB(1); res->t[4] = aB(2) * bB(2); res->t[5] = aB(2) * bB(3);
res->t[6] = aB(3) * bB(1); res->t[7] = aB(3) * bB(2); res->t[8] = aB(3) * bB(3);
}
return *res ;};
// lecture et écriture de données
#ifndef MISE_AU_POINT
inline
#endif
istream & Tenseur3BB::Lecture(istream & entree)
{ // lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3BB")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i<= 5; i++)
entree >> this->t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
ostream & Tenseur3BB::Ecriture(ostream & sort) const
{ // écriture du type
sort << "Tenseur3BB ";
// puis les datas
for (int i = 0; i<= 5; 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, Tenseur3BB & A)
{ int dim = A.Dimension();
#ifdef MISE_AU_POINT
if (dim != 3) A.Message(3,"operator >> (istream & entree, Tenseur3BB & A)");
#endif
// lecture et vérification du type
string nom_type;
entree >> nom_type;
if (nom_type != "Tenseur3BB")
{ Sortie(1);
return entree;
}
// lecture des coordonnées
for (int i = 0; i<= 5; i++)
entree >> A.t[i];
return entree;
};
#ifndef MISE_AU_POINT
inline
#endif
// surcharge de l'operator d'ecriture
ostream & operator << (ostream & sort , const Tenseur3BB & A)
{ //int dim = A.Dimension();
// écriture du type
sort << "Tenseur3BB ";
// puis les datas
for (int i = 0; i<= 5; i++)
sort << setprecision(ParaGlob::NbdigdoCA()) << A.t[i] << " ";
return sort;
};
#endif