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

1316 lines
50 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 "Debug.h"
#include "Util.h"
#include "Coordonnee3.h"
#include "Coordonnee2.h"
#include "Coordonnee1.h"
#include "ConstMath.h"
#include "MathUtil.h"
#ifndef Util_H_deja_inclus
//#ifndef MISE_AU_POINT
// inline
//#endif
//// PRODUIT VECTORIEL DE DEUX VECTEURS EN COORDONNEES ASBOLUS
//Vecteur Util::ProdVec( const Vecteur & v1, const Vecteur & v2)
// {
// #ifdef MISE_AU_POINT
// if (v1.Taille() != v2.Taille())
// { cout << "\nErreur : la taille des deux vecteur n'est pas identiques !\n";
// cout << "Util::ProdVec(Vecteur & v1,Vecteur & v2) \n";
// Sortie(1);
// };
// if (v1.Taille() != 3)
// { cout << "\nErreur : la taille des deux vecteur est differente de trois !\n";
// cout << "Util::ProdVec(Vecteur & v1,Vecteur & v2) \n";
// Sortie(1);
// };
// #endif
// Vecteur v(3);
// v(1) = v1(2) * v2(3) - v1(3) * v2(2);
// v(2) = v1(3) * v2(1) - v1(1) * v2(3);
// v(3) = v1(1) * v2(2) - v1(2) * v2(1);
// return v;
// };
//
#ifndef MISE_AU_POINT
inline
#endif
// idem en coordonnées absolues avec le type Coordonnee
Coordonnee Util::ProdVec_coor( const Coordonnee & v1, const Coordonnee & v2)
{
#ifdef MISE_AU_POINT
if (v1.Dimension() != v2.Dimension())
{ cout << "\nErreur : la taille des deux Coordonnees n'est pas identiques !\n";
cout << "Util::ProdVec(const Coordonnee & v1, const Coordonnee & v2) \n";
Sortie(1);
};
if (v1.Dimension() != 3)
{ cout << "\nErreur : la taille des deux vecteur est differente de trois !\n";
cout << "Util::ProdVec(Vecteur & v1,Vecteur & v2) \n";
Sortie(1);
};
#endif
Coordonnee v(3);
v(1) = v1(2) * v2(3) - v1(3) * v2(2);
v(2) = v1(3) * v2(1) - v1(1) * v2(3);
v(3) = v1(1) * v2(2) - v1(2) * v2(1);
return v;
};
// idem en coordonnées contravariantes avec le type CoordonneeH
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeH Util::ProdVec_coorH( const CoordonneeH & v1, const CoordonneeH & v2)
{
#ifdef MISE_AU_POINT
if (v1.Dimension() != v2.Dimension())
{ cout << "\nErreur : la taille des deux CoordonneesH n'est pas identiques !\n";
cout << "Util::ProdVec_coorH( const CoordonneeH & v1, const CoordonneeH & v2) \n";
Sortie(1);
};
if (v1.Dimension() != 3)
{ cout << "\nErreur : la taille des deux CoordonneesH est differente de trois !\n";
cout << "Util::ProdVec_coorH( const CoordonneeH & v1, const CoordonneeH & v2) \n";
Sortie(1);
};
#endif
CoordonneeH v(3);
v(1) = v1(2) * v2(3) - v1(3) * v2(2);
v(2) = v1(3) * v2(1) - v1(1) * v2(3);
v(3) = v1(1) * v2(2) - v1(2) * v2(1);
return v;
};
// idem en coordonnées covariantes avec le type CoordonneeB
#ifndef MISE_AU_POINT
inline
#endif
CoordonneeB Util::ProdVec_coorB( const CoordonneeB & v1, const CoordonneeB & v2)
{
#ifdef MISE_AU_POINT
if (v1.Dimension() != v2.Dimension())
{ cout << "\nErreur : la taille des deux CoordonneesB n'est pas identiques !\n";
cout << "Util::ProdVec_coorB( const CoordonneeB & v1, const CoordonneeB & v2) \n";
Sortie(1);
};
if (v1.Dimension() != 3)
{ cout << "\nErreur : la taille des deux CoordonneesB est differente de trois !\n";
cout << "Util::ProdVec_coorB( const CoordonneeB & v1, const CoordonneeB & v2) \n";
Sortie(1);
};
#endif
CoordonneeB v(3);
v(1) = v1(2) * v2(3) - v1(3) * v2(2);
v(2) = v1(3) * v2(1) - v1(1) * v2(3);
v(3) = v1(1) * v2(2) - v1(2) * v2(1);
return v;
};
// idem en coordonnées covariantes avec le type CoordonneeB, et en retour un coordonnee normal
#ifndef MISE_AU_POINT
inline
#endif
Coordonnee Util::ProdVec_coorBN( const CoordonneeB & v1, const CoordonneeB & v2)
{
#ifdef MISE_AU_POINT
if (v1.Dimension() != v2.Dimension())
{ cout << "\nErreur : la taille des deux CoordonneesB n'est pas identiques !\n";
cout << "Util::ProdVec_coorBN( const CoordonneeB & v1, const CoordonneeB & v2) \n";
Sortie(1);
};
if (v1.Dimension() != 3)
{ cout << "\nErreur : la taille des deux CoordonneesB est differente de trois !\n";
cout << "Util::ProdVec_coorBN( const CoordonneeB & v1, const CoordonneeB & v2) \n";
Sortie(1);
};
#endif
Coordonnee v(3);
v(1) = v1(2) * v2(3) - v1(3) * v2(2);
v(2) = v1(3) * v2(1) - v1(1) * v2(3);
v(3) = v1(1) * v2(2) - v1(2) * v2(1);
return v;
};
//#ifndef MISE_AU_POINT
// inline
//#endif
//// CALCUL DU DETERMINANT DE TROIS VECTEURS
//double Util::Determinant( const Vecteur & v1, const Vecteur & v2, const Vecteur & v3)
// {
// #ifdef MISE_AU_POINT
// if ((v1.Taille() != 3) || (v2.Taille() != 3) || (v3.Taille() != 3))
// { cout << "\nErreur : la taille des trois vecteur est differente de trois !\n";
// cout << "Util::Determinant(Vecteur & v1,Vecteur & v2,Vecteur & v3)\n";
// Sortie(1);
// }
// #endif
// return (ProdVec(v1,v2) * v3);
// };
// CALCUL DU DETERMINANT DE TROIS VECTEURS en coordonnées covariantes
#ifndef MISE_AU_POINT
inline
#endif
double Util::DeterminantB( const CoordonneeB & v1, const CoordonneeB & v2, const CoordonneeB & v3)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != 3) || (v2.Dimension() != 3) || (v3.Dimension() != 3))
{ cout << "\nErreur : la taille des trois Coordonee est differente de trois !\n";
cout << "Util::DeterminantB(CoordonneeB & v1,CoordonneeB & v2,CoordonneeB & v3)\n";
Sortie(1);
}
#endif
return ((ProdVec_coorB(v1,v2)).ScalBB(v3));
};
// CALCUL DU DETERMINANT DE TROIS VECTEURS en coordonnées absolues avec le type Coordonnee
#ifndef MISE_AU_POINT
inline
#endif
double Util::Determinant( const Coordonnee & v1, const Coordonnee & v2, const Coordonnee & v3)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != 3) || (v2.Dimension() != 3) || (v3.Dimension() != 3))
{ cout << "\nErreur : la taille des trois Coordonee est differente de trois !\n";
cout << "Util::Determinant(Coordonnee & v1,Coordonnee & v2,Coordonnee & v3)\n";
Sortie(1);
}
#endif
return (ProdVec_coor(v1,v2) * v3);
};
//#ifndef MISE_AU_POINT
// inline
//#endif
//// CALCUL DU DETERMINANT DE DEUX VECTEURS
//double Util::Determinant( const Vecteur & v1, const Vecteur & v2)
// {
// #ifdef MISE_AU_POINT
// if (((v1.Taille() != 3) && (v1.Taille() != 2))
// || ((v2.Taille() != 3) && (v2.Taille() != 2)) )
// { cout << "\nErreur : la dimension des deux vecteur est differente de trois et deux !\n";
// cout << "Util::Determinant(Vecteur & v1,Vecteur & v2)\n";
// Sortie(1);
// }
// #endif
// double res;
// switch (v1.Taille())
// { case 2 : res = v1(1) * v2(2) - v1(2) * v2(1); break;
// case 3 : res = !(ProdVec(v1,v2)); break;
// default : cout << "erreur Determinant(Vecteur & v1,Vecteur & v2)\n";
// Sortie(1);
// }
// return res;
// };
//
// CALCUL DU DETERMINANT DE DEUX VECTEURS en coordonnées covariantes
#ifndef MISE_AU_POINT
inline
#endif
double Util::DeterminantB( const CoordonneeB & v1, const CoordonneeB & v2)
{
#ifdef MISE_AU_POINT
if (((v1.Dimension() != 3) && (v1.Dimension() != 2))
|| ((v2.Dimension() != 3) && (v2.Dimension() != 2)) )
{ cout << "\nErreur : la dimension des deux Coordonnees est differente de trois et deux !\n";
cout << "Util::DeterminantB(CoordonneeB & v1,CoordonneeB & v2)\n";
Sortie(1);
}
#endif
double res;
switch (v1.Dimension())
{ case 2 : res = v1(1) * v2(2) - v1(2) * v2(1); break;
case 3 : res = (ProdVec_coorB(v1,v2)).Norme(); break;
default : cout << "erreur DeterminantB(CoordonneeB & v1,CoordonneeB & v2)\n";
Sortie(1);
}
return res;
};
// CALCUL DU DETERMINANT DE DEUX VECTEURS en coordonnées absolues avec le type Coordonnee
#ifndef MISE_AU_POINT
inline
#endif
double Util::Determinant( const Coordonnee & v1, const Coordonnee & v2)
{
#ifdef MISE_AU_POINT
if (((v1.Dimension() != 3) && (v1.Dimension() != 2))
|| ((v2.Dimension() != 3) && (v2.Dimension() != 2)) )
{ cout << "\nErreur : la dimension des deux Coordonnees est differente de trois et deux !\n";
cout << "Util::Determinant(Coordonnee & v1,Coordonnee & v2)\n";
Sortie(1);
}
#endif
double res;
switch (v1.Dimension())
{ case 2 : res = v1(1) * v2(2) - v1(2) * v2(1); break;
case 3 : res = (ProdVec_coor(v1,v2)).Norme(); break;
default : cout << "erreur Determinant(Coordonnee & v1,Coordonnee & v2)\n";
Sortie(1);
}
return res;
};
//#ifndef MISE_AU_POINT
// inline
//#endif
//// CALCUL DU DETERMINANT DE UN VECTEUR
//double Util::Determinant( const Vecteur & v1)
// {
// return v1.Norme(); //v1 * v1;
// };
//
// CALCUL DU DETERMINANT DE UN VECTEUR en coordonnée covariante
#ifndef MISE_AU_POINT
inline
#endif
double Util::DeterminantB( const CoordonneeB & v1)
{
return v1.Norme(); //v1 * v1;
};
// CALCUL DU DETERMINANT DE UN VECTEUR en coordonnées absolues avec le type Coordonnee
#ifndef MISE_AU_POINT
inline
#endif
double Util::Determinant( const Coordonnee & v1)
{
return v1.Norme(); //v1 * v1;
};
#ifndef MISE_AU_POINT
inline
#endif
// calcul de la variation d'un vecteur unitaire connaissant la variation du
// vecteur non norme
// v : le vecteur non norme, Dv : la variation de v, nor : la norme de v
// en retour : la variation de vecteur
Vecteur Util::VarUnVect( const Vecteur & v, const Vecteur & Dv, double nor)
{ double nor3 = nor * nor * nor;
Vecteur De = Dv / nor - v * ( (v * Dv) / nor3);
return De;
};
#ifndef MISE_AU_POINT
inline
#endif
// idem avec coordonnee en retour
Coordonnee Util::VarUnVect_coor( const Coordonnee & v, const Coordonnee & Dv, double nor)
{ double nor3 = nor * nor * nor;
Coordonnee De = Dv / nor - v * ( (v * Dv) / nor3);
return De;
};
#ifndef MISE_AU_POINT
inline
#endif
// idem avec coordonneeB en in-out
//!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
// ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
// d'une métrique associée à un repère non orthonormé
// ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
CoordonneeB Util::VarUnVect_coorB( const CoordonneeB & v, const CoordonneeB & Dv, double nor)
{ double nor3 = nor * nor * nor;
CoordonneeB De = Dv / nor - v * ( (v.ScalBB(Dv)) / nor3);
return De;
};
#ifndef MISE_AU_POINT
inline
#endif
// idem avec coordonneeH en in-out
//!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
// ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
// d'une métrique associée à un repère non orthonormé
// ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
CoordonneeH Util::VarUnVect_coorH( const CoordonneeH & v, const CoordonneeH & Dv, double nor)
{ double nor3 = nor * nor * nor;
CoordonneeH De = Dv / nor - v * ( (v.ScalHH(Dv)) / nor3);
return De;
};
#ifndef MISE_AU_POINT
inline
#endif
// calcul du tableau de variation d'un vecteur unitaire connaissant le tableau de variation du
// vecteur non norme
// v : le vecteur non norme, Dv : la variation de v, nor : la norme de v
// en retour : le tableau de variation
Tableau <Vecteur> Util::VarUnVect( const Vecteur & v, const Tableau <Vecteur>& Dv, double nor)
{ Tableau <Vecteur> De(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ De(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
}
return De;
};
// idem avec des coordonnees
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee> Util::VarUnVect_coor( const Coordonnee & v, const Tableau <Coordonnee>& Dv, double nor)
{ Tableau <Coordonnee> De(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ De(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
}
return De;
};
// idem avec des coordonneesB
#ifndef MISE_AU_POINT
inline
#endif
//!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
// ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
// d'une métrique associée à un repère non orthonormé
// ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
Tableau <CoordonneeB> Util::VarUnVect_coorB( const CoordonneeB & v, const Tableau <CoordonneeB>& Dv, double nor)
{ Tableau <CoordonneeB> De(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ De(i) = Dv(i) / nor - v * ( ( v.ScalBB(Dv(i)) ) / nor3);
}
return De;
};
#ifndef MISE_AU_POINT
inline
#endif
// idem et le tableau de retour passé en paramètre
// calcul du tableau de variation d'un vecteur unitaire connaissant le tableau de variation du
// vecteur non norme
// v : le vecteur non norme, Dv : la variation de v, nor : la norme de v
// en retour : le tableau de variation
Tableau <Vecteur>& Util::VarUnVect( const Vecteur & v, const Tableau <Vecteur>& Dv, double nor, Tableau <Vecteur>& retour)
{ retour.Change_taille(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ retour(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
}
return retour;
};
// idem avec des coordonnees
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee>& Util::VarUnVect_coor( const Coordonnee & v, const Tableau <Coordonnee>& Dv, double nor, Tableau <Coordonnee>& retour )
{ retour.Change_taille(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ retour(i) = Dv(i) / nor - v * ( (v * Dv(i)) / nor3);
}
return retour;
};
// idem avec des coordonneesB
#ifndef MISE_AU_POINT
inline
#endif
//!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
// ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
// d'une métrique associée à un repère non orthonormé
// ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
Tableau <CoordonneeB>& Util::VarUnVect_coorB( const CoordonneeB & v, const Tableau <CoordonneeB>& Dv, double nor,Tableau <CoordonneeB>& retour)
{ retour.Change_taille(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ retour(i) = Dv(i) / nor - v * ( ( v.ScalBB(Dv(i)) ) / nor3);
}
return retour;
};
// idem avec des coordonnees
// la variation du vecteur est supposé se trouver dans le premier vecteur de la base
#ifndef MISE_AU_POINT
inline
#endif
//!!!!!!!!!!!!! très important: il doit s'agir de vecteur exprimé dans un repère orthonormé
// ceci est vrai quelque soit la variance affichée: car ici on ne prend pas en compte la variation
// d'une métrique associée à un repère non orthonormé
// ex: les variations des gi sont ok car en fait les gi représentent les coorddonées dans un repère absolu
Tableau <Coordonnee>& Util::VarUnVect_coorBN( const CoordonneeB & v, const Tableau <BaseB>& Dv, double nor,Tableau <Coordonnee>& retour)
{ retour.Change_taille(Dv.Taille());
double nor3 = nor * nor * nor;
int DvTaille=Dv.Taille();
for (int i=1;i<= DvTaille;i++)
{ //retour(i) = (Dv(i) / nor - v * ( ( v.ScalBB(Dv(i)(1)) ) / nor3)).Coor_const();
const Coordonnee & dv = Dv(i).Coordo(1); // pour simplifier
retour(i) = (dv / nor - v.Coor_const() * ( v.Coor_const() * dv) / nor3);
}
return retour;
};
//#ifndef MISE_AU_POINT
// inline
//#endif
//// calcul de la variation d'un produit vectoriel
//// vi et Dvi les vecteurs du produit vectoriel et leurs variations
//Tableau <Vecteur> Util::VarProdVect( const Vecteur & v1, const Vecteur & v2,
// const Tableau <Vecteur >& Dv1, const Tableau <Vecteur >& Dv2)
// {
// #ifdef MISE_AU_POINT
// if ((v1.Taille() != v2.Taille()) || (Dv1.Taille() != Dv2.Taille())
// || (Dv1.Taille() == 0)
// || (v1.Taille() != Dv1(1).Taille()) || (Dv1(1).Taille() != Dv1(2).Taille()) )
// { cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
// cout << "Util::VarProdVect(Vecteur & v1,Vecteur & v2," <<
// "Tableau <Vecteur >& Dv1,Tableau <Vecteur >& Dv2)\n";
// Sortie(1);
// }
// #endif
// int Dv1Taille = Dv1.Taille();
// Tableau <Vecteur> De(Dv1Taille);
// for (int i=1;i<= Dv1Taille;i++)
// { De(i) = ProdVec(Dv1(i),v2) + ProdVec(v1,Dv2(i));
// }
// return De;
// };
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee> Util::VarProdVect_coor( const Coordonnee & v1, const Coordonnee & v2,
const Tableau <Coordonnee >& Dv1, const Tableau <Coordonnee >& Dv2)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv1.Taille() != Dv2.Taille())
|| (Dv1.Taille() == 0)
|| (v1.Dimension() != Dv1(1).Dimension()) || (Dv1(1).Dimension() != Dv1(2).Dimension()) )
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect_coor(Coordonnee & v1,Coordonnee & v2," <<
"Tableau <Coordonnee >& Dv1,Tableau <Coordonnee >& Dv2)\n";
Sortie(1);
}
#endif
int Dv1taille=Dv1.Taille();
Tableau <Coordonnee> De(Dv1taille);
for (int i=1;i<= Dv1taille;i++)
{ De(i) = ProdVec_coor(Dv1(i),v2) + ProdVec_coor(v1,Dv2(i));
}
return De;
};
// idem que le précédent module mais avec un retour en coordonnée, et un tableau de BaseB
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee> Util::VarProdVect_coor(const Coordonnee & v1, const Coordonnee & v2,
const Tableau <BaseB >& Dv)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
|| (Dv.Taille() == 0))
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
<< " const Tableau <BaseB >& Dv) \n";
Sortie(1);
}
#endif
int DvTaille = Dv.Taille();
Tableau <Coordonnee> De(DvTaille);
for (int i=1;i<= DvTaille;i++)
{ De(i) = ProdVec_coor(Dv(i).Coordo(1),v2) + ProdVec_coor(v1,Dv(i).Coordo(2));
}
return De;
};
// idem que le précédent module mais avec un in-out en coordonnéeB, et un tableau de BaseB
#ifndef MISE_AU_POINT
inline
#endif
Tableau <CoordonneeB> Util::VarProdVect_coorB(const CoordonneeB & v1, const CoordonneeB & v2,
const Tableau <BaseB >& Dv)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
|| (Dv.Taille() == 0))
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
<< " const Tableau <BaseB >& Dv) \n";
Sortie(1);
}
#endif
int DvTaille = Dv.Taille();
Tableau <CoordonneeB> De(DvTaille);
for (int i=1;i<= DvTaille;i++)
{ De(i) = ProdVec_coorB(Dv(i)(1),v2) + ProdVec_coorB(v1,Dv(i)(2));
}
return De;
};
// idem mais avec le retour passé en paramètre
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee>& Util::VarProdVect_coor( const Coordonnee & v1, const Coordonnee & v2,
const Tableau <Coordonnee >& Dv1, const Tableau <Coordonnee >& Dv2,Tableau <Coordonnee>& retour)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv1.Taille() != Dv2.Taille())
|| (Dv1.Taille() == 0)
|| (v1.Dimension() != Dv1(1).Dimension()) || (Dv1(1).Dimension() != Dv1(2).Dimension()) )
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect_coor(Coordonnee & v1,Coordonnee & v2," <<
"Tableau <Coordonnee >& Dv1,Tableau <Coordonnee >& Dv2)\n";
Sortie(1);
}
#endif
int Dv1taille=Dv1.Taille();
retour.Change_taille(Dv1taille);
for (int i=1;i<= Dv1taille;i++)
{ retour(i) = ProdVec_coor(Dv1(i),v2) + ProdVec_coor(v1,Dv2(i));
}
return retour;
};
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee>& Util::VarProdVect_coor(const Coordonnee & v1, const Coordonnee & v2,
const Tableau <BaseB >& Dv,Tableau <Coordonnee>& retour)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
|| (Dv.Taille() == 0))
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
<< " const Tableau <BaseB >& Dv) \n";
Sortie(1);
}
#endif
int DvTaille = Dv.Taille();
retour.Change_taille(DvTaille);
for (int i=1;i<= DvTaille;i++)
{ retour(i) = ProdVec_coor(Dv(i).Coordo(1),v2) + ProdVec_coor(v1,Dv(i).Coordo(2));
}
return retour;
};
#ifndef MISE_AU_POINT
inline
#endif
Tableau <CoordonneeB>& Util::VarProdVect_coorB(const CoordonneeB & v1, const CoordonneeB & v2,
const Tableau <BaseB >& Dv,Tableau <CoordonneeB>& retour)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
|| (Dv.Taille() == 0))
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
<< " const Tableau <BaseB >& Dv) \n";
Sortie(1);
}
#endif
int DvTaille = Dv.Taille();
retour.Change_taille(DvTaille);
for (int i=1;i<= DvTaille;i++)
{ retour(i) = ProdVec_coorB(Dv(i)(1),v2) + ProdVec_coorB(v1,Dv(i)(2));
}
return retour;
};
#ifndef MISE_AU_POINT
inline
#endif
Tableau <Coordonnee>& Util::VarProdVect_coorBN(const CoordonneeB & v1, const CoordonneeB & v2,
const Tableau <BaseB >& Dv,Tableau <Coordonnee>& retour)
{
#ifdef MISE_AU_POINT
if ((v1.Dimension() != v2.Dimension()) || (Dv(1)(1).Dimension() != v1.Dimension())
|| (Dv.Taille() == 0))
{ cout << "\nErreur : dans le calcul de la variation du produit vectoriel de 2 vecteurs";
cout << "Util::VarProdVect( int toto, const Vecteur & v1, const Vecteur & v2,"
<< " const Tableau <BaseB >& Dv) \n";
Sortie(1);
}
#endif
int DvTaille = Dv.Taille();
retour.Change_taille(DvTaille);
for (int i=1;i<= DvTaille;i++)
{ retour(i) = ProdVec_coorBN(Dv(i)(1),v2) + ProdVec_coorBN(v1,Dv(i)(2));
}
return retour;
};
#ifndef MISE_AU_POINT
inline
#endif
// retourne le numero du ddl recherche identifie par en, dans le tableau passé en paramètre
// s'il existe sinon 0
int Util::Existe(const Tableau<Ddl >& tab_ddl,Enum_ddl en)
{ static int posiddl; // initialisée par le système à zéro
int tab_ddl_Taille = tab_ddl.Taille();
// on commence par tester a partir de la derniere position sauvegardee
for (int i=posiddl+1;i<= tab_ddl_Taille;i++)
if (en == tab_ddl(i).Id_nom())
{ posiddl=i; return i; }
// on continue la recherche a partir du debut
for (int i=1;i<= posiddl;i++)
if (en == tab_ddl(i).Id_nom())
{ posiddl=i; return i; }
// sinon retour 0
return 0;
};
#ifndef MISE_AU_POINT
inline
#endif
// calcul du produit mixte des vecteurs d'une base naturelle en coordonnée covariante
double Util::ProduitMixte(const BaseB & tab_v )
{ int nbvecteur= tab_v.NbVecteur(); // le nombre de vecteur
// on utilise le fait que les vecteurs de la base sont uniquement 1D 2D ou 3D pour optimiser
switch (nbvecteur)
{ case 1: // on calcul la norme du vecteur qui est donc toujours positive
return DeterminantB(tab_v(1));
break;
case 2: // cas de deux vecteurs, on ramène la norme du produit vectoriel
return DeterminantB(tab_v(1),tab_v(2));
break;
case 3: // produit mixte normal qui peut-être négatif
return DeterminantB(tab_v(1),tab_v(2),tab_v(3));
break;
default:
cout << "\n nombre de vecteur de la base non correcte : " << nbvecteur
<< "\n double Util::ProduitMixte(const BaseB & tab_v )";
Sortie(1);
}
// pour le compilo
Sortie(1);
return 0.; // ne doit jamais passer ici
};
// calcul de la variation de la norme d'un vecteur, connaissant la variation du vecteur,
// le vecteur, et sa norme
#ifndef MISE_AU_POINT
inline
#endif
Tableau <double> Util::VarNorme( const Tableau <Coordonnee >& Dv,const Coordonnee& V,const double& norme)
{ int nbddl = Dv.Taille();
Tableau <double> tab_Nv(nbddl);
if (norme < ConstMath::trespetit) {return tab_Nv;}; // retour dans le cas où la vitesse est nulle
double unsurnorme = 1./norme;
for (int i=1;i<=nbddl;i++)
{ tab_Nv(i) = (V*Dv(i))*unsurnorme;};
return tab_Nv;
};
// idem mais avec le type vecteur
#ifndef MISE_AU_POINT
inline
#endif
Tableau <double> Util::VarNorme( const Tableau <Vecteur >& Dv,const Vecteur& V,const double& norme)
{ int nbddl = Dv.Taille();
Tableau <double> tab_Nv(nbddl);
if (norme < ConstMath::trespetit) {return tab_Nv;}; // retour dans le cas où la vitesse est nulle
double unsurnorme = 1./norme;
for (int i=1;i<=nbddl;i++)
{ tab_Nv(i) = (V*Dv(i))*unsurnorme;};
return tab_Nv;
};
// calcul de l'inverse d'une matrice 3x3 donnée par ces coordonnées, en retour la matrice
// d'entrée est remplacée par la matrice inverse
// 1) cas d'une matrice non symétrique (quelconque), rangement des valeurs:
// (1,1) ; (1,2) ; (1,3) ; (2,1) ; (2,2) ; (2,3) ; (3,1) ; (3,2) ; (3,3)
#ifndef MISE_AU_POINT
inline
#endif
void Util::Inverse_mat3x3(listdouble9Iter & i9Iter)
// le codage est très basique, pour essayer d'avoir un calcul très efficace !!
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
double t[9];
double * tex = (i9Iter)->donnees; // idem pour la matrice passée en paramètre
for ( int i=0; i< 9; i++)
t[i] = tex[i]; // recopie des valeurs
// pour plus de clarté on va définir des identifiants explicites
double* M11 =&t[0];
double* M12 =&t[1];
double* M13 =&t[2];
double* M21 =&t[3];
double* M22 =&t[4];
double* M23 =&t[5];
double* M31 =&t[6];
double* M32 =&t[7];
double* M33 =&t[8];
// on calcul un vecteur d'équilibrage V
double V[3];
V[0] = DabsMaX(*M11,DabsMaX(*M12, *M13));
V[1] = DabsMaX(*M21,DabsMaX(*M22, *M23));
V[2] = DabsMaX(*M31,DabsMaX(*M32, *M33));
if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) || (V[2] < ConstMath::trespetit ))
{ cout << "\n erreur *** la matrice 3x3 est singuliere et ne peut pas etre inversee"
<< "\n m11= "<< *M11 <<" m12= "<< *M12 <<" m13= "<< *M13
<< "\n m21= "<< *M21 <<" m22= "<< *M22 <<" m23= "<< *M23
<< "\n m31= "<< *M31 <<" m32= "<< *M32 <<" m33= "<< *M33 ;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n Util::Inverse_mat3x3(listdouble9Iter & i9Iter)";
Sortie(1);
};
V[0] = 1./V[0]; V[1] = 1./V[1];V[2] = 1./V[2];
// def d'un vecteur de permutation
short int indi[3]; // par exemple indi[1]+1 = indique pour la 2ième ligne de la nouvelle matrice alph_beta,
// à quelle ligne de l'ancienne matrice cela correspond
short int jdi[3]; // par exemple jdi[1]+1 = indique pour la 2ième ligne de l'ancienne matrice,
// à quelle ligne de la nouvelle matrice alpha_beta cela correspond
// permet de faire le retour inverse de la permutation
//-----------------------------------------------------------------------
// --- A triangulation de la matrice avec équilibrage suivant les lignes
//-----------------------------------------------------------------------
// Beta11 beta12 beta13
// alpha21 beta22 beta23
// alpha31 alpha32 beta33
// 1) == cas de la première nouvelle colonne
// on commence par chercher la ligne telle que |Mi1|*Vi est maxi
double mvi1 = Dabs(*M11) * V[0];
double mvi2 = Dabs(*M21) * V[1]; double mvi3 = Dabs(*M31) * V[2];
if (mvi1 >= mvi2) { if (mvi1 >= mvi3) { indi[0]=0;jdi[0]=0;} // mvi1 est le plus grand
else { indi[0]=2;jdi[2]=0;}; // mvi3 est le plus grand
}
else { if (mvi2 >= mvi3) { indi[0]=1;jdi[1]=0;} // mvi2 est le plus grand
else { indi[0]=2;jdi[2]=0;}; // mvi3 est le plus grand
};
// on permutte les lignes 1 et indi[0]+1 si indi[0] n'est pas égal à 0
double inter = 0.;
switch (indi[0])
{ case 1: // permutation des lignes 1 et 2
{inter = *M11; *M11 = *M21; *M21 = inter;
inter = *M12; *M12 = *M22; *M22 = inter;
inter = *M13; *M13 = *M23; *M23 = inter;
inter = V[0]; V[0] = V[1]; V[1] = inter;
}
break;
case 2: // permutation des lignes 1 et 3
{inter = *M11; *M11 = *M31; *M31 = inter;
inter = *M12; *M12 = *M32; *M32 = inter;
inter = *M13; *M13 = *M33; *M33 = inter;
inter = V[0]; V[0] = V[2]; V[2] = inter;
}
break;
};
// calcule de la première colonne
// beta11 = M11 donc rien n'a faire pour le premier terme
*M21 /= *M11; // alpha21 = M21 / beta11
*M31 /= *M11; // alpha31 = M31 / beta11
// 2) == cas de la deuxième nouvelle colonne
//Beta12 = M12 : donc rien n'a faire pour Beta12
// choix entre la 2ième et 3ième ligne
double p22 = *M22 - *M21 * *M12;
double p32 = *M32 - *M31 * *M12;
if(V[1]*Dabs(p22) < V[2]*Dabs(p32))
{ // cas ou l'on inverse les lignes 2 et 3
indi[1] = 2;jdi[2]=1;indi[2] = 1;jdi[1]=2;
inter = *M21; *M21 = *M31; *M31 = inter;
inter = *M22; *M22 = *M32; *M32 = inter;
inter = *M23; *M23 = *M33; *M33 = inter;
inter = p22; p22 = p32; p32 = inter;
}
else // sinon c'est ok
{ indi[1] = 1;jdi[1]=1;indi[2] = 2;jdi[2]=2;};
// calcul de la deuxième colonne
// beta12 = M12 donc rien n'a faire pour ce premier terme
*M22 = p22; // beta22 = M22 - alpha21 * beta12
*M32 = p32 / *M22; // alpha32 = 1/beta22 * (M32 - alpha31*beta12)
// 3) == cas de la troisième colonne
// beta13 = M13 donc rien n'a faire avec ce terme
*M23 -= *M21 * *M13 ; // beta23 = M23 - alpha21*beta13
*M33 -= *M31 * *M13 - *M32 * *M23 ; // beta33 = M33 - alpha31*beta13 - alpha32*beta23
//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 3 problèmes
//-----------------------------------------------------------------------
// on a les formules suivantes:
// [alpha](Z)=(Y) sachant que alphaii=1
// -> Z1=Y1; Z2=Y2-alpha21*Z1; Z3=Y3-alpha31*Z1-alpha32*Z2
//
// ensuite [beta](X)=(Z) d'où
// -> X3=Z3/beta33; X2=1/beta22 * (Z2-beta23*X3);
// X1=1./beta11 * (Z1-beta12*X2-beta13*X3)
//
// il faut tenir compte des permutations éventuelles Yi -> Y_i
// du coup on ne sait pas exactement quelles sont les grandeurs nulles ou pas
//1) première colonne
// Y1=1.,Y2=0.,Y3=0.; et sans permutation on a X1=tex[0], X2 = tex[3] et X3 = tex[6]
double Y[3];double * pt[3];
{Y[0]=1.;Y[1]=0.;Y[2]=0.;
double Z1 = Y[indi[0]]; // Z1 = Y_1
double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]];
double Z2= Y_2- *M21 * Z1; double Z3=Y_3- *M31 * Z1- *M32 * Z2;
pt[0] = &tex[0]; pt[1] = &tex[3]; pt[2] = &tex[6];
*(pt[jdi[2]]) = Z3/ *M33; *(pt[jdi[1]]) = (Z2 - *M23 * (*(pt[jdi[2]])) )/ *M22;
*(pt[jdi[0]]) = (Z1 - *M12 * (*(pt[jdi[1]])) - *M13 * (*(pt[jdi[2]])) )/ *M11;
}
//2) deuxième colonne
// Y1=0.,Y2=1.,Y3=0.; et sans permutation on a X1=tex[1], X2 = tex[4] et X3 = tex[7]
{Y[0]=0.;Y[1]=1.;Y[2]=0.;
double Z1 = Y[indi[0]]; // Z1 = Y_1
double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]];
double Z2= Y_2- *M21 * Z1; double Z3=Y_3- *M31 * Z1- *M32 * Z2;
pt[0] = &tex[1]; pt[1] = &tex[4]; pt[2] = &tex[7];
*(pt[jdi[2]]) = Z3/ *M33; *(pt[jdi[1]]) = (Z2 - *M23 * (*(pt[jdi[2]])) )/ *M22;
*(pt[jdi[0]]) = (Z1 - *M12 * (*(pt[jdi[1]])) - *M13 * (*(pt[jdi[2]])) )/ *M11;
}
//2) troisième colonne
// Y1=0.,Y2=0.,Y3=1.; et sans permutation on a X1=tex[2], X2 = tex[5] et X3 = tex[8]
{Y[0]=0.;Y[1]=0.;Y[2]=1.;
double Z1 = Y[indi[0]]; // Z1 = Y_1
double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]];
double Z2= Y_2- *M21 * Z1; double Z3=Y_3- *M31 * Z1- *M32 * Z2;
pt[0] = &tex[2]; pt[1] = &tex[5]; pt[2] = &tex[8];
*(pt[jdi[2]]) = Z3/ *M33; *(pt[jdi[1]]) = (Z2 - *M23 * (*(pt[jdi[2]])) )/ *M22;
*(pt[jdi[0]]) = (Z1 - *M12 * (*(pt[jdi[1]])) - *M13 * (*(pt[jdi[2]])) )/ *M11;
}
};
// calcul de l'inverse d'une matrice 3x3 donnée par ces coordonnées, en retour la matrice
// d'entrée est remplacée par la matrice inverse
// 2) cas d'une matrice symétrique , rangement des valeurs:
// (1,1) ; (2,2) ; (3,3) ; (2,1) = (1,2) ; (3,2) = (2,3) ; (3,1) = (1,3) ;
#ifndef MISE_AU_POINT
inline
#endif
void Util::Inverse_mat3x3(listdouble6Iter & i6Iter)
// le codage est très basique, pour essayer d'avoir un calcul très efficace !!
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
double t[9];
double * tex = (i6Iter)->donnees; // idem pour la matrice passée en paramètre
// recopie des valeurs
t[0] = tex[0]; t[4] = tex[1]; t[8] = tex[2];
t[1] = t[3] = tex[3]; t[5] = t[7]= tex[4]; t[6] = t[2] = tex[5];
// pour plus de clarté on va définir des identifiants explicites
double* A11 = &t[0];
double* A12 = &t[1];
double* A13 = &t[2];
double* A21 = &t[3];
double* A22 = &t[4];
double* A23 = &t[5];
double* A31 = &t[6];
double* A32 = &t[7];
double* A33 = &t[8];
//cout << "\n " << *A11 <<" "<< *A12 <<" "<< *A13 <<" "<< *A21 <<" "<< *A22 <<" "<< *A23 <<" "<< *A31 <<" "<< *A32 <<" "<< *A33 <<endl;
// on calcul un vecteur d'équilibrage V
double V[3];
V[0] = DabsMaX(*A11,DabsMaX(*A12, *A13));
V[1] = DabsMaX(*A21,DabsMaX(*A22, *A23));
V[2] = DabsMaX(*A31,DabsMaX(*A32, *A33));
if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) || (V[2] < ConstMath::trespetit ))
{ cout << "\n erreur *** la matrice 3x3 est singuliere et ne peut pas etre inversee"
<< "\n m11= "<< *A11 <<" m12= "<< *A12 <<" m13= "<< *A13
<< "\n m21= "<< *A21 <<" m22= "<< *A22 <<" m23= "<< *A23
<< "\n m31= "<< *A31 <<" m32= "<< *A32 <<" m33= "<< *A33 ;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n Util::Inverse_mat3x3(listdouble9Iter & i9Iter)";
Sortie(1);
};
V[0] = 1./V[0]; V[1] = 1./V[1];V[2] = 1./V[2];
// def d'un vecteur de permutation
short int indi[3]; // par exemple indi[1]+1 = indique pour la 2ième ligne de la nouvelle matrice alph_beta,
// à quelle ligne de l'ancienne matrice cela correspond
short int jdi[3]; // par exemple jdi[1]+1 = indique pour la 2ième ligne de l'ancienne matrice,
// à quelle ligne de la nouvelle matrice alpha_beta cela correspond
// permet de faire le retour inverse de la permutation
//-----------------------------------------------------------------------
// --- A triangulation de la matrice avec équilibrage suivant les lignes
//-----------------------------------------------------------------------
// Beta11 beta12 beta13
// alpha21 beta22 beta23
// alpha31 alpha32 beta33
// 1) == cas de la première nouvelle colonne
// on commence par chercher la ligne telle que |Mi1|*Vi est maxi
double mvi1 = Dabs(*A11) * V[0];
double mvi2 = Dabs(*A21) * V[1]; double mvi3 = Dabs(*A31) * V[2];
if (mvi1 >= mvi2) { if (mvi1 >= mvi3) { indi[0]=0;jdi[0]=0;} // mvi1 est le plus grand
else { indi[0]=2;jdi[2]=0;}; // mvi3 est le plus grand
}
else {if (mvi2 >= mvi3) { indi[0]=1;jdi[1]=0;} // mvi2 est le plus grand
else { indi[0]=2;jdi[2]=0;}; // mvi3 est le plus grand
};
// on permutte les lignes 1 et indi[0] si indi[0] n'est pas égal à 1
double inter = 0.;
switch (indi[0])
{ case 1: // permutation des lignes 1 et 2
{inter = *A11; *A11 = *A21; *A21 = inter;
inter = *A12; *A12 = *A22; *A22 = inter;
inter = *A13; *A13 = *A23; *A23 = inter;
inter = V[0]; V[0] = V[1]; V[1] = inter;
}
break;
case 2: // permutation des lignes 1 et 3
{inter = *A11; *A11 = *A31; *A31 = inter;
inter = *A12; *A12 = *A32; *A32 = inter;
inter = *A13; *A13 = *A33; *A33 = inter;
inter = V[0]; V[0] = V[2]; V[2] = inter;
}
break;
};
// calcule de la première colonne
// beta11 = A11 donc rien n'a faire pour le premier terme
*A21 /= *A11; // alpha21 = A21 / beta11
*A31 /= *A11; // alpha31 = A31 / beta11
// 2) == cas de la deuxième nouvelle colonne
//Beta12 = A12 : donc rien n'a faire pour Beta12
// choix entre la 2ième et 3ième ligne
double p22 = *A22 - *A21 * *A12;
double p32 = *A32 - *A31 * *A12;
if(V[1]*Dabs(p22) < V[2]*Dabs(p32))
{ // cas ou l'on inverse les lignes 2 et 3
indi[1] = 2;jdi[2]=1;indi[2] = 1;jdi[1]=2;
inter = *A21; *A21 = *A31; *A31 = inter;
inter = *A22; *A22 = *A32; *A32 = inter;
inter = *A23; *A23 = *A33; *A33 = inter;
inter = p22; p22 = p32; p32 = inter;
}
else // sinon c'est ok
{ indi[1] = 1;jdi[1]=1;indi[2] = 2;jdi[2]=2;};
// calcul de la deuxième colonne
// beta12 = A12 donc rien n'a faire pour ce premier terme
*A22 = p22; // beta22 = A22 - alpha21 * beta12
*A32 = p32 / *A22; // alpha32 = 1/beta22 * (A32 - alpha31*beta12)
// 3) == cas de la troisième colonne
// beta13 = A13 donc rien n'a faire avec ce terme
*A23 -= *A21 * *A13 ; // beta23 = A23 - alpha21*beta13
*A33 -= *A31 * *A13 - *A32 * *A23 ; // beta33 = A33 - alpha31*beta13 - alpha32*beta23
//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 3 problèmes
//-----------------------------------------------------------------------
// on a les formules suivantes:
// [alpha](Z)=(Y) sachant que alphaii=1
// -> Z1=Y1; Z2=Y2-alpha21*Z1; Z3=Y3-alpha31*Z1-alpha32*Z2
//
// ensuite [beta](X)=(Z) d'où
// -> X3=Z3/beta33; X2=1/beta22 * (Z2-beta23*X3);
// X1=1./beta11 * (Z1-beta12*X2-beta13*X3)
//
// il faut tenir compte des permutations éventuelles Yi -> Y_i
// du coup on ne sait pas exactement quelles sont les grandeurs nulles ou pas
//1) première colonne
// Y1=1.,Y2=0.,Y3=0.; et sans permutation on a X1=tex[0], X2 = tex[3] et X3 = tex[5]
double Y[3];double * pt[3];
{Y[0]=1.;Y[1]=0.;Y[2]=0.;
double Z1 = Y[indi[0]]; // Z1 = Y_1
double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]];
double Z2= Y_2- *A21 * Z1; double Z3=Y_3- *A31 * Z1- *A32 * Z2;
pt[0] = &tex[0]; pt[1] = &tex[3]; pt[2] = &tex[5];
*(pt[jdi[2]]) = Z3/ *A33; *(pt[jdi[1]]) = (Z2 - *A23 * (*(pt[jdi[2]])) )/ *A22;
*(pt[jdi[0]]) = (Z1 - *A12 * (*(pt[jdi[1]])) - *A13 * (*(pt[jdi[2]])) )/ *A11;
}
//2) deuxième colonne
// Y1=0.,Y2=1.,Y3=0.; et sans permutation on a X1=tex[3] déjà calculé, X2 = tex[1] et X3 = tex[4]
// on pourrait sans doute faire des éconnomies due à la symétrie, cependant due aux permutations
// ce n'est pas évident d'utiliser correctement les symétries, donc on recalcul tout
{Y[0]=0.;Y[1]=1.;Y[2]=0.;
double Z1 = Y[indi[0]]; // Z1 = Y_1
double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]];
double Z2= Y_2- *A21 * Z1; double Z3=Y_3- *A31 * Z1- *A32 * Z2;
pt[0] = &tex[3]; pt[1] = &tex[1]; pt[2] = &tex[4];
*(pt[jdi[2]]) = Z3/ *A33; *(pt[jdi[1]]) = (Z2 - *A23 * (*(pt[jdi[2]])) )/ *A22;
*(pt[jdi[0]]) = (Z1 - *A12 * (*(pt[jdi[1]])) - *A13 * (*(pt[jdi[2]])) )/ *A11;
}
//2) troisième colonne
// Y1=0.,Y2=0.,Y3=1.; et sans permutation on a X1=tex[5] déjà calculé, X2 = tex[4] déjà calculé et X3 = tex[2]
{Y[0]=0.;Y[1]=0.;Y[2]=1.;
double Z1 = Y[indi[0]]; // Z1 = Y_1
double Y_2 = Y[indi[1]];double Y_3 = Y[indi[2]];
double Z2= Y_2- *A21 * Z1; double Z3=Y_3- *A31 * Z1- *A32 * Z2;
pt[0] = &tex[5]; pt[1] = &tex[4]; pt[2] = &tex[2];
*(pt[jdi[2]]) = Z3/ *A33; *(pt[jdi[1]]) = (Z2 - *A23 * (*(pt[jdi[2]])) )/ *A22;
*(pt[jdi[0]]) = (Z1 - *A12 * (*(pt[jdi[1]])) - *A13 * (*(pt[jdi[2]])) )/ *A11;
}
};
// calcul de l'inverse d'une matrice 2x2 donnée par ces coordonnées, en retour la matrice
// d'entrée est remplacée par la matrice inverse
// 1) cas d'une matrice non symétrique (quelconque), rangement des valeurs:
// (1,1) ; (2,2) ; (2,1) ; (1,2) ;
#ifndef MISE_AU_POINT
inline
#endif
void Util::Inverse_mat2x2(listdouble4Iter & i4Iter)
// le codage est très basique, pour essayer d'avoir un calcul très efficace !!
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
double t[4];
double * tex = (i4Iter)->donnees; // idem pour la matrice passée en paramètre
// recopie des valeurs
t[0] = tex[0]; t[3] = tex[1]; t[2] = tex[2]; t[1] = tex[4];
// pour plus de clarté on va définir des identifiants explicites
// *** finalement on n'utilise pas de define car cela donne une grandeur définie pour tous Herezh
// et dans certaine méthodes on utilise des variables qui ont ce nom là !!
// #define B11 t[0]
// #define B12 t[1]
// #define B21 t[2]
// #define B22 t[3]
double* B11 = &t[0];
double* B12 = &t[1];
double* B21 = &t[2];
double* B22 = &t[3];
//cout << "\n " << *B11 <<" "<< *B12 <<" "<< *B21 <<" "<< *B22 << endl;
// on calcul un vecteur d'équilibrage V
double V[2];
V[0] = DabsMaX(*B11,*B12);
V[1] = DabsMaX(*B21,*B22);
if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) )
{ cout << "\n erreur *** la matrice 2x2 est singuliere et ne peut pas etre inversee"
<< "\n m11= "<< *B11 <<" m12= "<< *B12
<< "\n m21= "<< *B21 <<" m22= "<< *B22 ;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n Util::Inverse_mat2x2(listdouble4Iter & i4Iter)";
Sortie(1);
};
V[0] = 1./V[0]; V[1] = 1./V[1];
// def d'un indicateur de permutation
bool permutation = false; // indique si oui ou non on permutte la première ligne avec la seconde
// 1) == cas de la première nouvelle colonne
// on commence par chercher la ligne telle que |Mi1|*Vi est maxi
double mvi1 = Dabs(*B11) * V[0];
double mvi2 = Dabs(*B21) * V[1];
double inter = 0.;
if (mvi1 < mvi2)
{ permutation = true; // mvi2 est le plus grand
// permutation des lignes 1 et 2
inter = *B11; *B11 = *B21; *B21 = inter;
inter = *B12; *B12 = *B22; *B22 = inter;
inter = V[0]; V[0] = V[1]; V[1] = inter;
};
// calcule de la première colonne
// beta11 = A11 donc rien n'a faire pour le premier terme
*B21 /= *B11; // alpha21 = B21 / beta11
// 2) == cas de la deuxième nouvelle colonne
//Beta12 = B12 : donc rien n'a faire pour Beta12
// calcul du deuxième terme de la deuxième colonne
*B22 = *B22 - *B21 * *B12; // beta22 = B22 - alpha21 * beta12
//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 2 problèmes
//-----------------------------------------------------------------------
// on a les formules suivantes:
// [alpha](Z)=(Y) sachant que alphaii=1
// -> Z1=Y1; Z2=Y2-alpha21*Z1;
//
// ensuite [beta](X)=(Z) d'où
// X2=Z2/beta22 ;
// X1=1./beta11 * (Z1-beta12*X2)
//
// il faut tenir compte d'une permutation éventuelle
if (!permutation)
{ //1) première colonne
// Y1=1.,Y2=0., et sans permutation on a X1=tex[0], X2 = tex[2]
double Y2=0.;
double Z1 = 1.;
double Z2 = Y2 - *B21 * Z1;
tex[2] = Z2 / *B22; tex[0] = (Z1 - *B12 * tex[2] )/ *B11;
//2) deuxième colonne
// Y1=0.,Y2=1.; et sans permutation on a X1=tex[3], X2 = tex[1]
// on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
Y2=1.;
Z1 = 0.; // Z1 = Y_1
Z2 = Y2- *B21 * Z1;
tex[1] = Z2 / *B22; tex[3] = (Z1 - *B12 * tex[1] )/ *B11;
}
else
{ //avec permutation: 1) première colonne
// Y1=0.,Y2=1., et avec permutation on a X1=tex[2], X2 = tex[0]
double Y2=1.;
double Z1 = 0.;
double Z2 = Y2 - *B21 * Z1;
tex[0] = Z2 / *B22; tex[2] = (Z1 - *B12 * tex[0] )/ *B11;
//2) deuxième colonne
// Y1=1.,Y2=0.; et sans permutation on a X1=tex[1] , X2 = tex[3]
// on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
Y2=0.;
Z1 = 1.; // Z1 = Y_1
Z2 = Y2- *B21 * Z1;
tex[3] = Z2 / *B22; tex[1] = (Z1 - *B12 * tex[3] )/ *B11;
};
};
// calcul de l'inverse d'une matrice 2x2 donnée par ces coordonnées, en retour la matrice
// d'entrée est remplacée par la matrice inverse
// 2) cas d'une matrice symétrique , rangement des valeurs:
// (1,1) ; (2,2) ; (2,1) = (1,2) ;
#ifndef MISE_AU_POINT
inline
#endif
void Util::Inverse_mat2x2(listdouble3Iter & i3Iter)
// le codage est très basique, pour essayer d'avoir un calcul très efficace !!
{ // on commence par créer une matrice qui contient une copie de la matrice a inverser
double t[4];
double * tex = (i3Iter)->donnees; // idem pour la matrice passée en paramètre
// recopie des valeurs
t[0] = tex[0]; t[3] = tex[1]; t[1] = t[2] = tex[2];
// pour plus de clarté on va définir des identifiants explicites
double* C11 = &t[0];
double* C12 = &t[1];
double* C21 = &t[2];
double* C22 = &t[3];
//cout << "\n " << *C11 <<" "<< *C12 <<" "<< *C21 <<" "<< *C22 << endl;
// on calcul un vecteur d'équilibrage V
double V[2];
V[0] = DabsMaX(*C11,*C12);
V[1] = DabsMaX(*C21,*C22);
if ( (V[0] < ConstMath::trespetit ) || (V[1] < ConstMath::trespetit ) )
{ cout << "\n erreur *** la matrice 2x2 est singuliere et ne peut pas etre inversee"
<< "\n m11= "<< *C11 <<" m12= "<< *C12
<< "\n m21= "<< *C21 <<" m22= "<< *C22 ;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n Util::Inverse_mat2x2(listdouble3Iter & i3Iter)";
Sortie(1);
};
V[0] = 1./V[0]; V[1] = 1./V[1];
// def d'un indicateur de permutation
bool permutation = false; // indique si oui ou non on permutte la première ligne avec la seconde
//-----------------------------------------------------------------------
// --- A triangulation de la matrice avec équilibrage suivant les lignes
//-----------------------------------------------------------------------
// Beta11 beta12
// alpha21 beta22
// 1) == cas de la première nouvelle colonne
// on commence par chercher la ligne telle que |Mi1|*Vi est maxi
double mvi1 = Dabs(*C11) * V[0];
double mvi2 = Dabs(*C21) * V[1];
double inter = 0.;
if (mvi1 < mvi2)
{ permutation = true; // mvi2 est le plus grand
// permutation des lignes 1 et 2
inter = *C11; *C11 = *C21; *C21 = inter;
inter = *C12; *C12 = *C22; *C22 = inter;
inter = V[0]; V[0] = V[1]; V[1] = inter;
};
// calcule de la première colonne
// beta11 = A11 donc rien n'a faire pour le premier terme
*C21 /= *C11; // alpha21 = C21 / beta11
// 2) == cas de la deuxième nouvelle colonne
//Beta12 = C12 : donc rien n'a faire pour Beta12
// calcul du deuxième terme de la deuxième colonne
*C22 = *C22 - *C21 * *C12; // beta22 = C22 - alpha21 * beta12
//-----------------------------------------------------------------------
// --- B calcul de l'inverse de la matrice en résolvant 2 problèmes
//-----------------------------------------------------------------------
// on a les formules suivantes:
// [alpha](Z)=(Y) sachant que alphaii=1
// -> Z1=Y1; Z2=Y2-alpha21*Z1;
//
// ensuite [beta](X)=(Z) d'où
// X2=Z2/beta22 ;
// X1=1./beta11 * (Z1-beta12*X2)
//
// il faut tenir compte d'une permutation éventuelle
if (!permutation)
{ //1) première colonne
// Y1=1.,Y2=0., et sans permutation on a X1=tex[0], X2 = tex[2]
double Y2=0.;
double Z1 = 1.;
double Z2 = Y2 - *C21 * Z1;
tex[2] = Z2 / *C22; tex[0] = (Z1 - *C12 * tex[2] )/ *C11;
//2) deuxième colonne
// Y1=0.,Y2=1.; et sans permutation on a X1=tex[2] déjà calculé, X2 = tex[1]
// on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
Y2=1.;
Z1 = 0.; // Z1 = Y_1
Z2 = Y2- *C21 * Z1;
tex[1] = Z2 / *C22; tex[2] = (Z1 - *C12 * tex[1] )/ *C11;
}
else
{ //avec permutation: 1) première colonne
// Y1=0.,Y2=1., et avec permutation on a X1=tex[2], X2 = tex[0]
double Y2=1.;
double Z1 = 0.;
double Z2 = Y2 - *C21 * Z1;
tex[0] = Z2 / *C22; tex[2] = (Z1 - *C12 * tex[0] )/ *C11;
//2) deuxième colonne
// Y1=1.,Y2=0.; et sans permutation on a X1=tex[1] déjà calculé, X2 = tex[2]
// on pourrait faire des éconnomies due à la symétrie, cependant due aux permutations on recalcul tout
Y2=0.;
Z1 = 1.; // Z1 = Y_1
Z2 = Y2- *C21 * Z1;
tex[2] = Z2 / *C22; tex[1] = (Z1 - *C12 * tex[2] )/ *C11;
};
};
#endif
Reference in a new issue Copy permalink