1316 lines
50 KiB
C++
Executable file
1316 lines
50 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 "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
|