1091 lines
42 KiB
C++
Executable file
1091 lines
42 KiB
C++
Executable file
// 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-2021 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
|