rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/Util/MathUtil2.cc

1091 lines
42 KiB
C++
Raw Normal View History

ajout des fichiers .h et .cc des répertoire Util et TypeBase
2021-09-07 09:51:43 +02:00
// 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
Reference in a new issue Copy permalink