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