// 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 "MathUtil2.h"
//#include "nrutil.h"
//#include "jacobi.h"
#include "Util.h"
#include "ConstMath.h"
#include "MathUtil.h"




#ifndef  MATHUTIL2_H_deja_inclus

//// vecteurs propre et valeurs propres d'une matrice nxn c-a-d : en particulier 3 3 ou 2 2 ou 1 1
//// avec une limitation sur n (< 50)
//// *** dans le cas où les matrices sont symétriques !! ***
//// la méthode est itérative, --> méthode de Jacobi
//// la matrice A est ecrasee, a la sortie chaque colonne represente un vecteur propre
//// le vecteur de retour contiend les valeurs propres
//// s'il y a une erreur dans le calcul : cas = -1
//#ifndef MISE_AU_POINT
//  inline 
//#endif
//Vecteur MathUtil2::VectValPropre(Mat_pleine& A,int& cas)
//  { // appel d'une routine jacobi
//    int nbcol = A.Nb_colonne();
//    #ifndef MISE_AU_POINT
//     int nblign = A.Nb_ligne();
//     if (nbcol != nblign) 
//       { cout << "\n erreur dans VectValPropre(Mat_pleine& A), la matrice n'est pas carré"
//             << endl;
//        Sortie(1);
//       }      
//    #endif
//
//    // 1) def de la matrice et du vecteur contenant les vecteurs propre et les valeurs propres
//    Vecteur d(nbcol);
//    Mat_pleine V(nbcol,nbcol);
//    int nrot =0; // nombre de rotation dans jacobi
//    
//    // 3) appel du programme de la bibliothèque (fortement modifié)
//    Jacobi jaco; // une instance par défaut
//    
//    jaco.jacobi(A,d,V,nrot);
//    // traitement erreur
//    if (nrot == -1) {cas = -1;} else {cas = 1;};
//  
//    // 4) mise en ordre des valeurs et vecteurs propres
//    jaco.eigsrt(d,V); 
//    // 5) sauvegarde et retour
//    A = V;     
//    return d;
//  };
//

 
//calcul des vecteurs propres pour une matrice 3x3 réel uniquement, pas forcément symétrique
// dans le cas où on connait déjà les valeurs propres
// --> Utilisable que pour des matrices carrées
// en entrée : VP, de dimension 3, Vp contient les valeurs propres
//                 les 3 valeurs propres sont considérées classées: VP(1) >= VP(2) >= VP(3)
//             cas : qui indique le type de valeurs propres
//               , 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 Vp(1)=Vp(2) ( Vp(1), et Vp(3)  dans les 2 premières composantes du retour)
//               , cas = 3 si Vp(2)=Vp(3) ( Vp(1), et Vp(2)  dans les 2 premières composantes du retour)
// en sortie :  le tableau des vecteurs propres, rangés selon l'ordre d'entrée
//              cas = le cas de l'entrée, sauf s'il y a eu un pb, dans ce cas: cas = -1, et les vecteurs propres sont
//              dans ce cas quelconques
//
#ifndef MISE_AU_POINT
  inline 
#endif
Tableau <Coordonnee>   MathUtil2::V_Propres3x3(const Mat_pleine& A,const Coordonnee & VP, int & cas )
 {
  #ifdef MISE_AU_POINT 
  if ((A.Nb_colonne() != 3)|| (A.Nb_ligne() != 3))
      {cout << "\n erreur de dimension: Nb_colonne()= " << A.Nb_colonne() << ", Nb_ligne()= ";
       cout << A.Nb_ligne() << " alors que cette methode ne s'applique qu'aux matrices 3x3 ! "
            << "\n V_Propres3x3(..."  << endl;
       Sortie(1);
      };
  if (VP.Dimension() != 3)
      {cout << "\n erreur il n'y a pas 3 valeurs propres !  " << endl;
       Sortie(1);
      };
  
  if ( (VP(1) < VP(2)) || (VP(2) < VP(3)) || (VP(1) < VP(3)) )
      {cout << "\n erreur les trois valeurs propres ne sont pas dans le bon ordre (VP(1) >= VP(2) >= VP(3)) "
            << "\n VP(1)= " << VP(1) << ", VP(2)= " << VP(2) << ", VP(3)= " << VP(3)
            << "\n V_Propres3x3(..."  << endl;
       Sortie(1);
      };
  #endif
  // initialisation
  Coordonnee3 zero(3);
  Tableau <Coordonnee> Vpropre(3,zero); // les vecteurs propres
  
  // on utilise une copie de la matrice d'entrée
  Mat_pleine B(A);
  // dans tous les cas il y a au moins une valeur propre, on démarre avec celle là
  // on calcul la matrice B-lambda*Id qui normalement doit-être singulière
////---debug
//cout << "\n debug V_Propres3x3 ";
//cout <<"\n B: ";
//B.Affiche();
//// --- fin debug

  for (int i=1;i<=3;i++)
    B(i,i) -= VP(1);

////---debug
//cout << "\n debug V_Propres3x3 apres B(i,i) -= VP(1)";
//cout <<"\n B: ";
//B.Affiche();
//// --- fin debug
 
  // on calcul le rang de la matrice
  int rang_B = CalculRangMatrice3x3Singuliere(B);
//---debug
//cout << "\n debug V_Propres3x3 ";
//cout <<"\n rangB: "<<rang_B;
//// --- fin debug
  switch (rang_B)
   { case 0:  // si le rang = 0, la matrice B est uniformément nulle, et on peut choisir n'importe quel base
              // on prend la base identité
              // cela veut dire aussi que la matrice A initiale était sphérique !! donc que les trois valeurs propres sont identiques
      {Vpropre(1)(1)=1.;Vpropre(2)(2)=1.;Vpropre(3)(3)=1.;
       #ifdef MISE_AU_POINT
        {// en débug on vérifie que le rang est également nulle pour les deux autres valeurs
         Mat_pleine B2(A);
         for (int i=1;i<=3;i++)
           B2(i,i) -= VP(2);
         Mat_pleine B3(A);
         for (int i=1;i<=3;i++)
           B3(i,i) -= VP(3);
         if ((CalculRangMatrice3x3Singuliere(B2) != 0) || (CalculRangMatrice3x3Singuliere(B3) != 0))
           {cout << "\n ** erreur au niveau des rangs de matrice: rang_B == 0,"
                 << " mais rang_B2 =  " << CalculRangMatrice3x3Singuliere(B2)
                 << " et rang_B3= " << CalculRangMatrice3x3Singuliere(B3)
                 << "\n V_Propres3x3(..."  << endl;
            cas = -1; // on signale l'erreur
            break;    // et on arrête
           };
         if (cas != 0)
           {cout << "\n ** erreur au niveau des rangs de matrice: rang_B == 0,"
                 << " normalement on devrait avoir cas = 0, et en fait on a : cas=   " << cas
                 << " (pour voir les matrices utiliser un niveau de commentaire > 7) "
                 << "\n VP: ";
            if (ParaGlob::NiveauImpression() > 7)
               { VP.Affiche();
                 cout <<"\n A: ";
                 A.Affiche();
                 cout <<"\n B: ";
                 B.Affiche();
                 cout <<"\n rangB: "<<rang_B ;
                 cout << "\n cas= "<< cas << " VP= "<<VP;
                 cout << "\n les val propres: " << VP ;
               };
            cout  << "\n V_Propres3x3(..."  << endl;
            cas = -1; // on signale l'erreur
            break;    // et on arrête
           };
        };
       #endif
       // vérif pour le cas non débug
       if (cas != 0)
         {cas = -1; // on signale l'erreur
          break;    // et on arrête
         };
       break;
      }

     case 1:  // on n'a qu'une seule direction de libre
              // cela veut dire également que l'on a deux valeurs propres identiques
      { // on récupère la première ligne de la matrice qui est la direction libre
        Coordonnee3 V(B(1).Coordo()); // c-a-d la première direction libre
        // on définit les 2 vecteurs du plan normal
        Def_vecteurs_plan(V, Vpropre(2), Vpropre(3));
        // et le premier, qui correspond à la première valeur propre (celle qui a servi au calcul de B)
        Vpropre(1) = Util::ProdVec_coor(Vpropre(2),Vpropre(3));
        #ifdef MISE_AU_POINT
         // en débug que l'on a bien deux valeurs propres identiques
         if ((cas == 0) || (cas == 1))
           {cout << "\n erreur au niveau du rang de matrice: rang_B == 1,"
                 << " on devrait donc avoir deux valeurs propres identiques alors que cas =  " << cas
                 << " (pour voir les matrices utiliser un niveau de commentaire > 7) "
                 << "\n VP: ";
            if (ParaGlob::NiveauImpression() > 7)
               { VP.Affiche();
                 cout <<"\n A: ";
                 A.Affiche();
                 cout <<"\n B: ";
                 B.Affiche();
                 cout <<"\n rangB: "<<rang_B ;
                 cout << "\n cas= "<< cas << " VP= "<<VP;
                 cout << "\n les val propres: " << VP ;
               };
            cout  << "\n V_Propres3x3(..."  << endl;
            cas = -1; // on signale l'erreur
            break;    // et on arrête
           };
        #endif
        // vérif pour le cas non débug
        if ((cas == 0) || (cas == 1))
          {cas = -1; // on signale l'erreur
           break;    // et on arrête
          };
        break;
      }
     case 2:
      { // on récupère les deux première lignes de la matrice qui forme le plan normal à la direction principale
        Coordonnee3 V1 = B(1).Coordo(); // premier vecteur
        Coordonnee3 V2 = B(2).Coordo(); // second vecteur
        // on définit le vecteur du plan normal
        Vpropre(1) = Util::ProdVec_coor(V1,V2);
        Vpropre(1).Normer();
       
        // on s'occupe maintenant de la seconde valeur propre
        Mat_pleine B2(A);
        for (int i=1;i<=3;i++)
          B2(i,i) -= VP(2);
        int rang_B2 = CalculRangMatrice3x3Singuliere(B2);
        #ifdef MISE_AU_POINT
         // rang_B2 ne peut pas être nul
         if (rang_B2 == 0)
           {cout << "\n erreur au niveau du rang de matrice: rang_B == 2,"
                 << " et on a trouve rang_B2= 0 !! "
                 << " (pour voir les matrices utiliser un niveau de commentaire > 7) "
                 << "\n VP: ";
            if (ParaGlob::NiveauImpression() > 7)
               { VP.Affiche();
                 cout <<"\n A: ";
                 A.Affiche();
                 cout <<"\n B: ";
                 B.Affiche();
                 cout <<"\n rangB: "<<rang_B << ", rangB2= "<< rang_B2;
                 cout << "\n cas= "<< cas << " VP= "<<VP;
                 cout << "\n les val propres: " << VP ;
               };
            cout  << "\n V_Propres3x3(..."  << endl;
            cas = -1; // on signale l'erreur
            break;    // et on arrête
           };
        #endif
        // vérif pour le cas non débug
        if (rang_B2 == 0)
          {cas = -1; // on signale l'erreur
           break;    // et on arrête
          };
        if (rang_B2 == 1)
         {// cela veut dire que les deux dernières valeurs propres sont identiques
          // donc on peut choisir un plan normal au premier vecteur
          // produit 2 vecteurs normés
          Def_vecteurs_plan(Vpropre(1), Vpropre(2), Vpropre(3));
          #ifdef MISE_AU_POINT
           // en débug que l'on a bien deux valeurs propres identiques
           if ((cas == 0) || (cas == 1))
             {cout << "\n erreur au niveau du rang de matrice: rang_B1 == 2 et rang_B2 == 1,"
                   << " on devrait donc avoir deux valeurs propres identiques alors que cas =  " << cas
                 << " (pour voir les matrices utiliser un niveau de commentaire > 7) "
                 << "\n VP: ";
              if (ParaGlob::NiveauImpression() > 7)
                 { VP.Affiche();
                   cout <<"\n A: ";
                   A.Affiche();
                   cout <<"\n B: ";
                   B.Affiche();
                   cout <<"\n rangB: "<<rang_B << ", rangB2= "<< rang_B2;
                   cout << "\n cas= "<< cas << " VP= "<<VP;
                   cout << "\n les val propres: " << VP ;
                 };
              cout  << "\n V_Propres3x3(..."  << endl;
              cas = -1; // on signale l'erreur
              break;    // et on arrête
             };
          #endif
          // vérif pour le cas non débug
          if ((cas == 0) || (cas == 1))
            {cas = -1; // on signale l'erreur
             break;    // et on arrête
            };
         }
        else
         { // sinon rang_B2 == 2
           Coordonnee3 W1 = B2(1).Coordo(); // premier vecteur
           Coordonnee3 W2 = B2(2).Coordo(); // second vecteur
           // on définit le vecteur du plan normal
           Vpropre(2) = Util::ProdVec_coor(W1,W2);
           Vpropre(2).Normer();
           // puis le dernier vecteur propre: produit de deux vecteurs normés
           // donc a priori est normé également par le produit vectoriel
           Vpropre(3) = Util::ProdVec_coor(Vpropre(1),Vpropre(2));
           #ifdef MISE_AU_POINT
           if (cas == 0)
             {cout << "\n warning au niveau du rang de matrice: rang_B == 2 et rang_B2 == 2,"
                   << " bizarre les trois valeurs propres sont identiques cas =  " << cas
                   << "\n V_Propres3x3(..."  << endl;
              if (ParaGlob::NiveauImpression() > 7)
               {cout << "\n matrice initiale: ";
                A.Affiche();
                cout << "\n B: ";
                B.Affiche();
                cout << "\n B2: ";
                B2.Affiche();
                cout <<"\n rangB: "<<rang_B << ", rangB2= "<< rang_B2;
                cout << "\n cas= "<< cas << " VP= "<<VP;
                cout << "\n les val propres: " << VP ;
               };
             };
           #endif
           // vérif pour le cas non débug
//           if (cas == 0)
//             {cas = -1; // on signale l'erreur
//              break;    // et on arrête
//             };
         };
        break;
      }
     default:
      {cout << "\n erreur au niveau du rang de matrice: rang_B doit etre compris entre 0 et 2 "
            << " ici =  " << rang_B
            << "\n V_Propres3x3(..."  << endl;
       Sortie(1);
      };
      break;
   };
  
  #ifdef MISE_AU_POINT
             // en débug que les valeurs propres fonctionnent correctement
//               , 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 Vp(1)=Vp(2) ( Vp(1), et Vp(3)  dans les 2 premières composantes du retour)
//               , cas = 3 si Vp(2)=Vp(3) ( Vp(1), et Vp(2)  dans les 2 premières composantes du retour)
  { Vecteur Vp(3);
    switch (cas)
     { case 0: Vp(1)=Vp(2)=Vp(3)=VP(1); break;
       case 1: Vp = VP.Vect(); break;
       case 2: Vp(1)=Vp(2)=VP(1); Vp(3)=VP(3); break;
       case 3: Vp(2)=Vp(3)=VP(2); Vp(1)=VP(1); break;
     };
    for (int i=1;i<=3;i++)
      { Coordonnee3 W6 = A.Prod_mat_vec(Vpropre(i).Vect());
        // normalement W6 = Vp(i)*Vpropre(i)
        double diff = (W6 - (Vpropre(i) * Vp(i))).Norme();
        if (Dabs(diff) > ConstMath::unpeupetit)
          {cout << "\n erreur le vecteur (i= " << i <<") propre ";
           Vpropre(i).Affiche();
           cout << "\n n'est pas un vecteur propre pour la valeur propre: "
                << Vp(i) << " concernant la matrice: ";
           A.Affiche();
           cout << "\n A*V = "<<W6 << "(Vpropre(i) * Vp(i))= " << (Vpropre(i) * Vp(i))
                << " diff= "<<diff;
           
           cout << "\n cas= "<< cas << " VP= "<<VP;
           cout << "\n les val propres: " << Vp ;
           
           cout  << "\n V_Propres3x3(..."  << endl;
          };
      };
  };
  
  #endif
  // retour
  return Vpropre;
  
 };
  

// calcul du rang d'une matrice 3x3 singulière, qui peut être 0, 1 ou 2
// modifie la matrice en conséquence
// en sortie:
//  rang : = le rang de la matrice
//       rang = 0 ==> veut dire que c'est la matrice nulle (à ConstMath::pasmalpetit près)
//       rang = 1 ==> la première ligne de la matrice contient la direction libre, génératrice
//       rang = 2 ==> les deux premières ligne de la matrice contiennent les deux directions libres
#ifndef MISE_AU_POINT
  inline 
#endif
int MathUtil2::CalculRangMatrice3x3Singuliere( Mat_pleine& A)
  {
    #ifdef MISE_AU_POINT
    if ((A.Nb_colonne() != 3)|| (A.Nb_ligne() != 3))
        {cout << "erreur de dimension: Nb_colonne()= " << A.Nb_colonne() << ", Nb_ligne()= ";
         cout << A.Nb_ligne() << " alors que cette methode ne s'applique qu'aux matrices 3x3 ! "
              << "\n CalculRangMatrice(..."  << endl;
         Sortie(1);
        }
    #endif
    int rang = 0; // init par défaut
    // on récupère la norme maxi de la matrice
    int imax,jmax;
    double maxMat = A.MaxiValAbs(imax,jmax);
    // si maxMat == 0, la matrice est nulle -> rang = 0
    if (maxMat <= ConstMath::pasmalpetit)
     { rang = 0; }
    else
     {// sinon le rang sera au moins égal à 1 et maxMat n'est pas trop petit
      double unSurmaxMat = 1./maxMat;
      // on permute la première ligne, pour qu'elle contienne le maxi
      if (imax != 1)
       {// sauvegarde de la ligne du maxi, avec normalisation à 1
        Vecteur ligne_du_maxi = unSurmaxMat * A(imax);
        A.Ligne_set(imax) = A(1);  // permutation
        A.Ligne_set(1) = ligne_du_maxi;       // """"
       };
      // maintenant on fait une réduction, c-a-d on transforme la matrice pour qu'il n'y ait que (1,0,0) dans la colonne du maxi
      Vecteur col_du_maxi = A(jmax); // récup de la colonne du maxi
      A.Ligne_set(2) = A(2) - col_du_maxi(2)*A(1);
      A.Ligne_set(3) = A(3) - col_du_maxi(3)*A(1);
      // arrivé ici normalement on a A(imax,2) = A(imax,3) = 0.
      // on calcul le maxi des composantes de la deuxième et troisième ligne de la matrice
      // -- tout d'abord la deuxième ligne
      int imax2=2; // init
      int jmax2; double maxMat2 = A(2).Max_val_abs(jmax2);
      // -- idem troisième ligne
      int imax3=3; // init
      int jmax3; double maxMat3 = A(3).Max_val_abs(jmax3);
      // les maxi finaux sont avec les indices 2 (imax2, jmax2, maxMat2)
      if (maxMat3 > maxMat2)
        {imax2=imax3;jmax2=jmax3;maxMat2 = maxMat3;
        };
      // si le maxi est très petit, on n'a qu'une seule direction de libre
      if (maxMat2 <= ConstMath::pasmalpetit)
       { rang = 1; }
      else
       { // sinon on refait la même manip entre les lignes 2 et 3
         double unSurmaxMat2 = 1./maxMat2;
         // on permute la dernière ligne, si besoin, pour que la deuxième contienne le maxi
         if (imax2 != 2)
          {// sauvegarde de la ligne du maxi, avec normalisation à 1
           Vecteur ligne_du_maxi = unSurmaxMat2 * A(imax2);
           A.Ligne_set(imax2) = A(2);  // permutation
           A.Ligne_set(2) = ligne_du_maxi;       // """"
          };
         // le rang est 2 et les 2 premières lignes contiennent deux vecteurs indépendants
         rang = 2;
       }
     };
    // retour
    return rang;
  };


    // --- cas d'une dimension d'espace == 3

// calcul deux vecteurs perpendiculaires à un troisième, générateurs d'un plan, et orthonormée
// Les deux vecteurs résultats sont quelconques a priori
#ifndef MISE_AU_POINT
  inline 
#endif
void  MathUtil2::Def_vecteurs_plan(const Coordonnee& V1, Coordonnee& V2, Coordonnee& V3)
  { double norme_V1= V1.Norme();
    #ifdef MISE_AU_POINT
    if (norme_V1 < ConstMath::trespetit)
        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
         cout << " le calcul du plan normal est impossible ! "
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        };
    if ((V1.Dimension() != 3)|| (V2.Dimension() != 3) || (V3.Dimension() != 3))
        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 3: V1.Taille()= " << V1.Dimension()
              << ", V2.Dimension()= " << V2.Dimension() << ", V3.Dimension()= " << V3.Dimension()
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        }
    #endif
    Coordonnee V(V1/norme_V1); // on normalise

    // on va chercher l'ordre des composantes de V
    int un,deux,trois;
    if (Dabs(V(2)) > Dabs(V(1)))
     { if (Dabs(V(2)) > Dabs(V(3)))
         {un = 2;
          if (Dabs(V(1)) > Dabs(V(3)))
           {deux = 1; trois = 3;}
          else
           {deux = 3; trois = 1;}
         }
       else
         {un = 3; deux= 2; trois = 1;}
     }
    else
     { if ((Dabs(V(3)) > Dabs(V(1))))
        {un = 3; deux = 1; trois = 2; }
       else
        {if ((Dabs(V(3)) > Dabs(V(2))))
          {un = 1; deux = 3; trois = 2;}
         else
          {un = 1; deux = 2; trois = 3;};
        }
     };
    // on calcule les vecteurs en fonction
    double unsurlongue = 1./sqrt(V(un)*V(un) + V(trois)*V(trois));
    V2(un) = -V(trois) * unsurlongue;
    V2(deux) = 0.;
    V2(trois) = V(un) * unsurlongue;
    
    V3(un) = V(deux) * V2(trois);
    V3(deux) = V(trois) * V2(un) - V(un) * V2(trois);
    V3(trois) = -V(deux) * V2(un);
    
    #ifdef MISE_AU_POINT
    if (ParaGlob::NiveauImpression() > 8)
     {// on vérifie l'algo de tri
      bool erreur = false;
      if ((Dabs(V(deux)) > Dabs(V(un))) || (Dabs(V(trois)) > Dabs(V(un)))
          || (Dabs(V(trois)) > Dabs(V(deux)))
         )
       {cout << "\n *** erreur dans le tri des composantes de V: "
             << "un= "<< un << "deux= "<< deux <<" trois= "<< trois
             << " V= "<< V ;
        erreur = true;
       };
      // on vérifie les normes et les produits scalaires
      if ((V.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V a une norme " << V.Norme() << " diff de 1 ";
        erreur = true;
       }
      if ((V2.Norme()-1) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V2 a une norme " << V2.Norme() << " diff de 1 ";
        erreur = true;
       }
      if ((V3.Norme()-1) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 a une norme " << V3.Norme() << " diff de 1 ";
        erreur = true;
       }
      if (Dabs(V * V2) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V2 n'est pas normal a V1 "
            << "\n V1= " << V1 << " V2= " << V2 ;
        erreur = true;
       }
      if (Dabs(V * V3) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 n'est pas normal a V1 "
            << "\n V1= " << V << " V3= " << V3 ;
        erreur = true;
       }
      if (Dabs(V2 * V3) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 n'est pas normal a V2 "
            << "\n V3= " << V3 << " V2= " << V2 ;
        erreur = true;
       }
      if (erreur)
       { cout << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << flush;
         Sortie(1);
       };
     };
    #endif

/* // calcul original, qui ne fonctionne pas pour le cas V(1) et V(3) nulle ensemble !!


    //
    if (Dabs(V(1)) >= Dabs(V(3)))
     { double unsurlongue = 1./sqrt(V(1)*V(1) + V(3)*V(3));
       // si V(3) est nulle, alors V(1) est forcément non nulle ==>longue est non nul
       // si (V(3) est non nulle, alors ==>longue est non nul
       V2(1) = -V(3) * unsurlongue;
       V2(2) = 0.;
       V2(3) = V(1) * unsurlongue;
      
       V3(1) = V(2) * V2(3);
       V3(2) = V(3) * V2(1) - V(1) * V2(3);
       V3(3) = -V(2) * V2(1);
     }
    else
     { double unsurlongue = 1./sqrt(V(2)*V(2) + V(3)*V(3));
       // ici V(2) prend le rôle de V(1)
       V2(1) = 0.;
       V2(2) = V(3) * unsurlongue;
       V2(3) = -V(2) * unsurlongue;
      
       V3(1) = V(2) * V2(3) - V(3) * V2(2);
       V3(2) = -V(1) * V2(3);
       V3(3) = V(1) * V2(2);
     };
 */
  };
// calcul deux vecteurs perpendiculaires à un troisième, générateurs d'un plan, et orthonormée
// Les deux vecteurs résultats sont quelconques a priori
#ifndef MISE_AU_POINT
  inline 
#endif
void  MathUtil2::Def_vecteurs_plan(const CoordonneeB& V1, CoordonneeB& V2, CoordonneeB& V3)
  { double norme_V1= V1.Norme();
    #ifdef MISE_AU_POINT
    if (norme_V1 < ConstMath::trespetit)
        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
         cout << " le calcul du plan normal est impossible ! "
              << "\n Def_vecteurs_plan(const CoordonneeB& V1,..."  << endl;
         Sortie(1);
        };
    if ((V1.Dimension() != 3)|| (V2.Dimension() != 3) || (V3.Dimension() != 3))
        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 3: V1.Taille()= " << V1.Dimension()
              << ", V2.Dimension()= " << V2.Dimension() << ", V3.Dimension()= " << V3.Dimension()
              << "\n Def_vecteurs_plan(const CoordonneeB& V1,..."  << endl;
         Sortie(1);
        }
    #endif
    CoordonneeB V(V1/norme_V1); // on normalise
   
    // on va chercher l'ordre des composantes de V
    int un,deux,trois;
    if (Dabs(V(2)) > Dabs(V(1)))
     { if (Dabs(V(2)) > Dabs(V(3)))
         {un = 2;
          if (Dabs(V(1)) > Dabs(V(3)))
           {deux = 1; trois = 3;}
          else
           {deux = 3; trois = 1;}
         }
       else
         {un = 3; deux= 2; trois = 1;}
     }
    else
     { if ((Dabs(V(3)) > Dabs(V(1))))
        {un = 3; deux = 1; trois = 2; }
       else
        {if ((Dabs(V(3)) > Dabs(V(2))))
          {un = 1; deux = 3; trois = 2;}
         else
          {un = 1; deux = 2; trois = 3;};
        }
     };
    // on calcule les vecteurs en fonction
    double unsurlongue = 1./sqrt(V(un)*V(un) + V(trois)*V(trois));
    V2(un) = -V(trois) * unsurlongue;
    V2(deux) = 0.;
    V2(trois) = V(un) * unsurlongue;
    
    V3(un) = V(deux) * V2(trois);
    V3(deux) = V(trois) * V2(un) - V(un) * V2(trois);
    V3(trois) = -V(deux) * V2(un);
    
    #ifdef MISE_AU_POINT
    if (ParaGlob::NiveauImpression() > 8)
     {// on vérifie l'algo de tri
      bool erreur = false;
      if ((Dabs(V(deux)) > Dabs(V(un))) || (Dabs(V(trois)) > Dabs(V(un)))
          || (Dabs(V(trois)) > Dabs(V(deux)))
         )
       {cout << "\n *** erreur dans le tri des composantes de V: "
             << "un= "<< un << "deux= "<< deux <<" trois= "<< trois
             << " V= "<< V ;
        erreur = true;
       };
      // on vérifie les normes et les produits scalaires
      if ((V.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V a une norme " << V.Norme() << " diff de 1 ";
        erreur = true;
       }
      if ((V2.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V2 a une norme " << V2.Norme() << " diff de 1 ";
        erreur = true;
       }
      if ((V3.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 a une norme " << V3.Norme() << " diff de 1 ";
        erreur = true;
       }
      if (Dabs(V.Coor() * V2.Coor()) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V2 n'est pas normal a V1 "
            << "\n V1= " << V1 << " V2= " << V2 ;
        erreur = true;
       }
      if (Dabs(V.Coor() * V3.Coor()) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 n'est pas normal a V1 "
            << "\n V1= " << V << " V3= " << V3 ;
        erreur = true;
       }
      if (Dabs(V2.Coor() * V3.Coor()) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 n'est pas normal a V2 "
            << "\n V3= " << V3 << " V2= " << V2 ;
        erreur = true;
       }
      if (erreur)
       { cout << "\n Def_vecteurs_plan(const CoordonneeB& V1,..."  << flush;
         Sortie(1);
       };
     }
    #endif

    
/* // calcul original, qui ne fonctionne pas pour le cas V(1) et V(3) nulle ensemble !!
    if (Dabs(V(1)) >= Dabs(V(3)))
       { double unsurlongue = 1./sqrt(V(1)*V(1) + V(3)*V(3));
         // si V(3) est nulle, alors V(1) est forcément non nulle ==>longue est non nul
         // si (V(3) est non nulle, alors ==>longue est non nul
         V2(1) = -V(3) * unsurlongue;
         V2(2) = 0.;
         V2(3) = V(1) * unsurlongue;
        
         V3(1) = V(2) * V2(3);
         V3(2) = V(3) * V2(1) - V(1) * V2(3);
         V3(3) = -V(2) * V2(1);
       }
      else
       { double unsurlongue = 1./sqrt(V(2)*V(2) + V(3)*V(3));
         // ici V(2) prend le rôle de V(1)
         V2(1) = 0.;
         V2(2) = V(3) * unsurlongue;
         V2(3) = -V(2) * unsurlongue;
        
         V3(1) = V(2) * V2(3) - V(3) * V2(2);
         V3(2) = -V(1) * V2(3);
         V3(3) = V(1) * V2(2);
       };
*/
       
  };
// calcul deux vecteurs perpendiculaires à un troisième, générateurs d'un plan, et orthonormée
// Les deux vecteurs résultats sont quelconques a priori
#ifndef MISE_AU_POINT
  inline
#endif
void  MathUtil2::Def_vecteurs_plan(const CoordonneeH& V1, CoordonneeH& V2, CoordonneeH& V3)
  { double norme_V1= V1.Norme();
    #ifdef MISE_AU_POINT
    if (norme_V1 < ConstMath::trespetit)
        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
         cout << " le calcul du plan normal est impossible ! "
              << "\n Def_vecteurs_plan(const CoordonneeH& V1,..."  << endl;
         Sortie(1);
        };
    if ((V1.Dimension() != 3)|| (V2.Dimension() != 3) || (V3.Dimension() != 3))
        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 3: V1.Taille()= " << V1.Dimension()
              << ", V2.Dimension()= " << V2.Dimension() << ", V3.Dimension()= " << V3.Dimension()
              << "\n Def_vecteurs_plan(const CoordonneeH& V1,..."  << endl;
         Sortie(1);
        }
    #endif
    CoordonneeH V(V1/norme_V1); // on normalise
   
    // on va chercher l'ordre des composantes de V
    int un,deux,trois;
    if (Dabs(V(2)) > Dabs(V(1)))
     { if (Dabs(V(2)) > Dabs(V(3)))
         {un = 2;
          if (Dabs(V(1)) > Dabs(V(3)))
           {deux = 1; trois = 3;}
          else
           {deux = 3; trois = 1;}
         }
       else
         {un = 3; deux= 2; trois = 1;}
     }
    else
     { if ((Dabs(V(3)) > Dabs(V(1))))
        {un = 3; deux = 1; trois = 2; }
       else
        {if ((Dabs(V(3)) > Dabs(V(2))))
          {un = 1; deux = 3; trois = 2;}
         else
          {un = 1; deux = 2; trois = 3;};
        }
     };
    // on calcule les vecteurs en fonction
    double unsurlongue = 1./sqrt(V(un)*V(un) + V(trois)*V(trois));
    V2(un) = -V(trois) * unsurlongue;
    V2(deux) = 0.;
    V2(trois) = V(un) * unsurlongue;
   
    V3(un) = V(deux) * V2(trois);
    V3(deux) = V(trois) * V2(un) - V(un) * V2(trois);
    V3(trois) = -V(deux) * V2(un);
   
    #ifdef MISE_AU_POINT
    if (ParaGlob::NiveauImpression() > 8)
     {// on vérifie l'algo de tri
      bool erreur = false;
      if ((Dabs(V(deux)) > Dabs(V(un))) || (Dabs(V(trois)) > Dabs(V(un)))
          || (Dabs(V(trois)) > Dabs(V(deux)))
         )
       {cout << "\n *** erreur dans le tri des composantes de V: "
             << "un= "<< un << "deux= "<< deux <<" trois= "<< trois
             << " V= "<< V ;
        erreur = true;
       };
      // on vérifie les normes et les produits scalaires
      if ((V.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V a une norme " << V.Norme() << " diff de 1 ";
        erreur = true;
       }
      if ((V2.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V2 a une norme " << V2.Norme() << " diff de 1 ";
        erreur = true;
       }
      if ((V3.Norme()-1.) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 a une norme " << V3.Norme() << " diff de 1 ";
        erreur = true;
       }
      if (Dabs(V.Coor() * V2.Coor()) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V2 n'est pas normal a V1 "
            << "\n V1= " << V1 << " V2= " << V2 ;
        erreur = true;
       }
      if (Dabs(V.Coor() * V3.Coor()) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 n'est pas normal a V1 "
            << "\n V1= " << V << " V3= " << V3 ;
        erreur = true;
       }
      if (Dabs(V2.Coor() * V3.Coor()) > ConstMath::petit)
       {cout << "erreur le nouveau vecteur V3 n'est pas normal a V2 "
            << "\n V3= " << V3 << " V2= " << V2 ;
        erreur = true;
       }
      if (erreur)
       { cout << "\n Def_vecteurs_plan(const CoordonneeH& V1,..."  << flush;
         Sortie(1);
       };
     }
    #endif

   
/* // calcul original, qui ne fonctionne pas pour le cas V(1) et V(3) nulle ensemble !!
    if (Dabs(V(1)) >= Dabs(V(3)))
       { double unsurlongue = 1./sqrt(V(1)*V(1) + V(3)*V(3));
         // si V(3) est nulle, alors V(1) est forcément non nulle ==>longue est non nul
         // si (V(3) est non nulle, alors ==>longue est non nul
         V2(1) = -V(3) * unsurlongue;
         V2(2) = 0.;
         V2(3) = V(1) * unsurlongue;
 
         V3(1) = V(2) * V2(3);
         V3(2) = V(3) * V2(1) - V(1) * V2(3);
         V3(3) = -V(2) * V2(1);
       }
      else
       { double unsurlongue = 1./sqrt(V(2)*V(2) + V(3)*V(3));
         // ici V(2) prend le rôle de V(1)
         V2(1) = 0.;
         V2(2) = V(3) * unsurlongue;
         V2(3) = -V(2) * unsurlongue;
 
         V3(1) = V(2) * V2(3) - V(3) * V2(2);
         V3(2) = -V(1) * V2(3);
         V3(3) = V(1) * V2(2);
       };
*/
   
  };

    // --- cas d'une dimension d'espace == 2

// calcul d'un vecteur perpendiculaire à un premier, et orthonormée
#ifndef MISE_AU_POINT
  inline
#endif
void  MathUtil2::Def_vecteurs_plan(const Coordonnee& V1, Coordonnee& V2)
  { double norme_V1= V1.Norme();
    #ifdef MISE_AU_POINT
    if (norme_V1 < ConstMath::trespetit)
        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
         cout << " le calcul du plan normal est impossible ! "
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        };
    if ((V1.Dimension() != 2)|| (V2.Dimension() != 2) )
        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 2: V1.Taille()= " << V1.Dimension()
              << ", V2.Dimension()= " << V2.Dimension()
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        }
    #endif
    Coordonnee V(V1/norme_V1); // on normalise
   
    // le vecteur V2 est normale à V
    // et V,V2 -> sens direct (inverse de l'horloge)
    V2(1) = - V1(2);
    V2(2) = V1(1);
  };

// calcul d'un vecteur perpendiculaire à un premier, et orthonormée
#ifndef MISE_AU_POINT
  inline
#endif
void  MathUtil2::Def_vecteurs_plan(const CoordonneeB& V1, CoordonneeB& V2)
  { double norme_V1= V1.Norme();
    #ifdef MISE_AU_POINT
    if (norme_V1 < ConstMath::trespetit)
        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
         cout << " le calcul du plan normal est impossible ! "
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        };
    if ((V1.Dimension() != 2)|| (V2.Dimension() != 2) )
        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 2: V1.Taille()= " << V1.Dimension()
              << ", V2.Dimension()= " << V2.Dimension()
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        }
    #endif
    CoordonneeB V(V1/norme_V1); // on normalise
   
    // le vecteur V2 est normale à V
    // et V,V2 -> sens direct (inverse de l'horloge)
    V2(1) = - V1(2);
    V2(2) = V1(1);
  };

// calcul d'un vecteur perpendiculaire à un premier, et orthonormée
#ifndef MISE_AU_POINT
  inline
#endif
void  MathUtil2::Def_vecteurs_plan(const CoordonneeH& V1, CoordonneeH& V2)
  { double norme_V1= V1.Norme();
    #ifdef MISE_AU_POINT
    if (norme_V1 < ConstMath::trespetit)
        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
         cout << " le calcul du plan normal est impossible ! "
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        };
    if ((V1.Dimension() != 2)|| (V2.Dimension() != 2) )
        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 2: V1.Taille()= " << V1.Dimension()
              << ", V2.Dimension()= " << V2.Dimension()
              << "\n Def_vecteurs_plan(const Coordonnee& V1,..."  << endl;
         Sortie(1);
        }
    #endif
    CoordonneeH V(V1/norme_V1); // on normalise
   
    // le vecteur V2 est normale à V
    // et V,V2 -> sens direct (inverse de l'horloge)
    V2(1) = - V1(2);
    V2(2) = V1(1);
  };


//// calcul deux vecteurs perpendiculaires à un troisième, générateurs d'un plan, et orthonormées
//// Les deux vecteurs résultats sont quelconques a priori
//// calcul également la variation des vecteurs résultats
////  en entrée: V1 et d_Vi(1)(ddl) : vecteur donné et sa variation
////  en sortie: V2 et V3, d_Vi(2)(ddl) et d_Vi(3)(ddl) : vecteurs résultats et leurs variations
//#ifndef MISE_AU_POINT
//  inline 
//#endif
//void Def_vecteurs_plan_et_variation(const CoordonneeB& V1, CoordonneeB& V2, CoordonneeB& V3
//                                    ,Tableau <BaseB>&  d_Vi)
//  { double norme_V1= V1.Norme();
//    #ifdef MISE_AU_POINT
//    if (norme_V1 < ConstMath::trespetit)
//        {cout << "erreur le vecteur V1 a une norme nulle (= " << norme_V1 << ") ";
//         cout << " le calcul du plan normal est impossible ! "
//              << "\n Def_vecteurs_plan_et_variation(const CoordonneeB& V1,..."  << endl;
//         Sortie(1);
//        };
//    if ((V1.Dimension() != 3)|| (V2.Dimension() != 3) || (V3.Dimension() != 3))
//        {cout << "erreur de dimension dans les vecteurs, qui devraient etre de dimension 3: V1.Taille()= " << V1.Dimension()
//              << ", V2.Dimension()= " << V2.Dimension() << ", V3.Dimension()= " << V3.Dimension()
//              << "\n Def_vecteurs_plan_et_variation(const CoordonneeB& V1,..."  << endl;
//         Sortie(1);
//        }
//    #endif
//    CoordonneeB V(V1/norme_V1); // on normalise
//    //
//    if (Dabs(V(1)) >= Dabs(V(3)))
//     { double unsurlongue = 1./sqrt(V(1)*V(1) + V(3)*V(3));
//       // si V(3) est nulle, alors V(1) est forcément non nulle ==>longue est non nul
//       // si (V(3) est non nulle, alors ==>longue est non nul
//       V2(1) = -V(3) * unsurlongue;
//       V2(2) = 0.;
//       V2(3) = V(1) * unsurlongue;
//      
//       V3(1) = V(2) * V2(3);
//       V3(2) = V(3) * V2(1) - V(1) * V2(3);
//       V3(3) = -V(2) * V2(1);
//      
//       // et les variations
//       int taille = d_Vi.Taille();
//       for (int i=1; i<= taille; i++)
//        { BaseB & d_Vi_i = d_Vi(i); // pour simplifier
//          CoordonneeB & d_V1 = d_Vi_i(1);  // "
//          double d_denom = 2. * (d_V1(1) * V(1) + d_V1(3) * V(3));
//          double d_unsurlongue =
//           1./sqrt(V(1)*V(1) + V(3)*V(3));
//         
//          CoordonneeB & d_V2 = d_Vi_i(2);  // "
//       
//         
//        };
//     }
//    else
//     { double unsurlongue = 1./sqrt(V(2)*V(2) + V(3)*V(3));
//       // ici V(2) prend le rôle de V(1)
//       V2(1) = 0.;
//       V2(2) = V(3) * unsurlongue;
//       V2(3) = -V(2) * unsurlongue;
//      
//       V3(1) = V(2) * V2(3) - V(3) * V2(2);
//       V3(2) = -V(1) * V2(3);
//       V3(3) = V(1) * V2(2);
//     };
//  };
//

// calcul d'un changement de base: ceci n'est pas fait dans la classe Coordonnee car il faut y adjoindre
// la classe Mat_pleine qui intègre beaucoup de chose, du coup la classe Coordonnee deviendrait une usine
    // changement de base (cf. théorie) : la matrice beta est telle que:
    // gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
    // et la matrice gamma telle que:
    // gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
    // gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
    //   c-a-d= gp^i = gamma^i_j * g^j
    // rappel des différentes relations entre beta et gamma
    // [beta]^{-1} = [gamma]^T      ; [beta]^{-1T} = [gamma]
    // [beta] = [gamma]^{-1T}      ; [beta]^{T} = [gamma]^{-1}
    // changement de base pour de une fois covariant:
    // [Ap_k] = [beta] * [A_i]

#ifndef MISE_AU_POINT
  inline
#endif
void MathUtil2::ChBase(const CoordonneeB& A_B,const Mat_pleine& beta,CoordonneeB& Ap_B)
  { int dim = A_B.Dimension();
    #ifdef MISE_AU_POINT
    if  ((beta.Nb_ligne() != dim) || (beta.Nb_colonne() != dim))
     { cout << "\nErreur : matrice beta n'a pas les memes dimensions beta.Nb_ligne():  "
            << beta.Nb_ligne() <<" et beta.Nb_colonne(): " << beta.Nb_colonne()
            << " que le vecteur de dim: "
            << dim << " ! le changement de base n'est pas possible ! \n";
       cout << "ChBase(const CoordonneeB& A_B,const Mat_pleine& beta,CoordonneeB& Ap_B)  \n";
       Sortie(1);
     };
    #endif
    Vecteur inter = beta.Prod_mat_vec(A_B.Vect());
    for (int i=1; i<= dim;i++)
       Ap_B(i) = inter(i);
  };


    // changement de base pour de une fois contravariant:
    // [Ap^k] = [gamma] * [A^i]
#ifndef MISE_AU_POINT
  inline
#endif
void MathUtil2::ChBase(const CoordonneeH& A_H,const Mat_pleine& gamma,CoordonneeH& Ap_H)
  { int dim = A_H.Dimension();
    #ifdef MISE_AU_POINT
    if  ((gamma.Nb_ligne() != dim) || (gamma.Nb_colonne() != dim))
     { cout << "\nErreur : matrice gamma n'a pas les memes dimensions gamma.Nb_ligne():  "
            << gamma.Nb_ligne() <<" et gamma.Nb_colonne(): " << gamma.Nb_colonne()
            << " que le vecteur de dim: "
            << dim << " ! le changement de base n'est pas possible ! \n";
       cout << "ChBase(const CoordonneeH& A_H,const Mat_pleine& gamma,CoordonneeH& Ap_H)  \n";
       Sortie(1);
     };
    #endif
   
    Vecteur inter = gamma.Prod_mat_vec(A_H.Vect());
////-------- debug -------
//cout << "\n debug MathUtil2::ChBase( ...";
//cout << "\n A_H initiale: "<< A_H;
//cout << "\n gamma: "; gamma.Affiche();
//cout << "\n gamma * A_H " << inter << flush;
////-------- fin debug ---------
    for (int i=1; i<= dim;i++)
       Ap_H(i) = inter(i);
  };



#endif