182 lines
9 KiB
C
182 lines
9 KiB
C
|
// 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/>.
|
||
|
|
|
||
|
|
/************************************************************************
|
||
|
|
* DATE: 23/01/97 *
|
||
|
|
* $ *
|
||
|
|
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
|
||
|
|
* $ *
|
||
|
|
* PROJET: Herezh++ *
|
||
|
|
* $ *
|
||
|
|
************************************************************************
|
||
|
|
* BUT: Utilitaires divers mathematiques *
|
||
|
|
* $ *
|
||
|
|
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
|
||
|
|
* VERIFICATION: *
|
||
|
|
* *
|
||
|
|
* ! date ! auteur ! but ! *
|
||
|
|
* ------------------------------------------------------------ *
|
||
|
|
* ! ! ! ! *
|
||
|
|
* $ *
|
||
|
|
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
|
||
|
|
* MODIFICATIONS: *
|
||
|
|
* ! date ! auteur ! but ! *
|
||
|
|
* ------------------------------------------------------------ *
|
||
|
|
* $ *
|
||
|
|
************************************************************************/
|
||
|
|
#ifndef MATHUTIL2_H
|
||
|
|
#define MATHUTIL2_H
|
||
|
|
|
||
|
|
#include "Vecteur.h"
|
||
|
|
#include "Coordonnee.h"
|
||
|
|
|
||
|
|
#include "Mat_pleine.h"
|
||
|
|
|
||
|
|
/** @defgroup Classes_utilitaires_sur_vecteurs_et_matrices Classes_utilitaires_sur_vecteurs_et_matrices
|
||
|
|
*
|
||
|
|
* \author Gérard Rio
|
||
|
|
* \version 1.0
|
||
|
|
* \date 23/01/97
|
||
|
|
* \brief Def d'utilitaires divers mathematiques sur vecteurs et matrices
|
||
|
|
*
|
||
|
|
*/
|
||
|
|
|
||
|
|
/// @addtogroup Classes_utilitaires_sur_vecteurs_et_matrices
|
||
|
|
/// @{
|
||
|
|
///
|
||
|
|
|
||
|
|
/// cas d'une erreur survenue à cause d'une non convergence
|
||
|
|
|
||
|
|
class ErrMathUtil2
|
||
|
|
/// =0 cas courant, pas d'information particulière
|
||
|
|
{ public :
|
||
|
|
/// constructeur par défaut (pas d'erreur)
|
||
|
|
ErrMathUtil2 () : cas(0) {} ;
|
||
|
|
/// constructeur à partir d'un indice d'erreur ca donné
|
||
|
|
ErrMathUtil2 (int ca) : cas(ca) {} ;
|
||
|
|
|
||
|
|
/// contient l'erreur
|
||
|
|
int cas;
|
||
|
|
};
|
||
|
|
/// -------- fin gestion des erreurs ---------------
|
||
|
|
/// @} // end of group
|
||
|
|
|
||
|
|
/// @addtogroup Classes_utilitaires_sur_vecteurs_et_matrices
|
||
|
|
/// @{
|
||
|
|
///
|
||
|
|
|
||
|
|
class MathUtil2
|
||
|
|
{ public :
|
||
|
|
/// vecteurs propre et valeurs propres d'une matrice nxn c-a-d : en particulier 3 3 ou 2 2 ou 1 1
|
||
|
|
/// avec une limitation sur n (< 50)
|
||
|
|
/// *** dans le cas où les matrices sont symétriques !! ***
|
||
|
|
/// la méthode est itérative, --> méthode de Jacobi
|
||
|
|
/// la matrice A est ecrasee, a la sortie chaque colonne represente un vecteur propre
|
||
|
|
/// le vecteur de retour contiend les valeurs propres
|
||
|
|
/// s'il y a une erreur dans le calcul : cas = -1
|
||
|
|
static Vecteur VectValPropre(Mat_pleine& A,int& cas);
|
||
|
|
|
||
|
|
///calcul des vecteurs propres pour une matrice 3x3 réel uniquement, pas forcément symétrique
|
||
|
|
/// dans le cas où on connait déjà les valeurs propres
|
||
|
|
/// la méthode est directe: --> pas d'itération, donc très rapide
|
||
|
|
/// --> Utilisable que pour des matrices carrées
|
||
|
|
/// en entrée : VP, de dimension 3, Vp contient les valeurs propres
|
||
|
|
/// les 3 valeurs propres sont considérées classées: VP(1) >= VP(2) >= VP(3)
|
||
|
|
/// cas : qui indique le type de valeurs propres
|
||
|
|
/// , cas = 1 si 3 val propres différentes (= 3 composantes du vecteur de retour)
|
||
|
|
/// , cas = 0 si les 3 sont identiques (= la première composantes du vecteur de retour),
|
||
|
|
/// , cas = 2 si Vp(1)=Vp(2) ( Vp(1), et Vp(3) dans les 2 premières composantes du retour)
|
||
|
|
/// , cas = 3 si Vp(2)=Vp(3) ( Vp(1), et Vp(2) dans les 2 premières composantes du retour)
|
||
|
|
/// en sortie : le tableau des vecteurs propres, rangés selon l'ordre d'entrée
|
||
|
|
/// cas = le cas de l'entrée, sauf s'il y a eu un pb, dans ce cas: cas = -1, et les vecteurs propres sont
|
||
|
|
/// dans ce cas quelconques
|
||
|
|
static Tableau <Coordonnee> V_Propres3x3(const Mat_pleine& A,const Coordonnee & VP, int & cas );
|
||
|
|
|
||
|
|
|
||
|
|
/// calcul du rang d'une matrice 3x3 singulière, qui peut être 0, 1 ou 2
|
||
|
|
/// modifie la matrice en conséquence
|
||
|
|
/// en sortie:
|
||
|
|
/// rang : = le rang de la matrice
|
||
|
|
/// rang = 0 ==> veut dire que c'est la matrice nulle (à ConstMath::pasmalpetit près)
|
||
|
|
/// rang = 1 ==> la première ligne de la matrice contient la direction libre, génératrice
|
||
|
|
/// rang = 2 ==> les deux premières ligne de la matrice contiennent les deux directions libres
|
||
|
|
static int CalculRangMatrice3x3Singuliere(Mat_pleine& A);
|
||
|
|
|
||
|
|
/// --- cas d'une dimension d'espace == 3
|
||
|
|
/// calcul deux vecteurs perpendiculaires à un troisième, générateurs d'un plan, et orthonormées
|
||
|
|
/// Les deux vecteurs résultats sont quelconques a priori
|
||
|
|
static void Def_vecteurs_plan(const Coordonnee& V1, Coordonnee& V2, Coordonnee& V3);
|
||
|
|
static void Def_vecteurs_plan(const CoordonneeB& V1, CoordonneeB& V2, CoordonneeB& V3);
|
||
|
|
static void Def_vecteurs_plan(const CoordonneeH& V1, CoordonneeH& V2, CoordonneeH& V3);
|
||
|
|
|
||
|
|
/// --- cas d'une dimension d'espace == 2
|
||
|
|
/// calcul d'un vecteur perpendiculaire à un premier, et orthonormées
|
||
|
|
static void Def_vecteurs_plan(const Coordonnee& V1, Coordonnee& V2);
|
||
|
|
static void Def_vecteurs_plan(const CoordonneeB& V1, CoordonneeB& V2);
|
||
|
|
static void Def_vecteurs_plan(const CoordonneeH& V1, CoordonneeH& V2);
|
||
|
|
|
||
|
|
/// calcul deux vecteurs perpendiculaires à un troisième, générateurs d'un plan, et orthonormées
|
||
|
|
/// Les deux vecteurs résultats sont quelconques a priori
|
||
|
|
/// calcul également la variation des vecteurs résultats
|
||
|
|
/// en entrée: V1 et d_Vi(1)(ddl) : vecteur donné et sa variation
|
||
|
|
/// en sortie: V2 et V3, d_Vi(2)(ddl) et d_Vi(3)(ddl) : vecteurs résultats et leurs variations
|
||
|
|
//void Def_vecteurs_plan_et_variation(const CoordonneeB& V1, CoordonneeB& V2, CoordonneeB& V3
|
||
|
|
// ,Tableau <BaseB>& d_Vi);
|
||
|
|
|
||
|
|
/// calcul d'un changement de base: ceci n'est pas fait dans la classe Coordonnee car il faut y adjoindre
|
||
|
|
/// la classe Mat_pleine qui intègre beaucoup de chose, du coup la classe Coordonnee deviendrait une usine
|
||
|
|
/// changement de base (cf. théorie) : la matrice beta est telle que:
|
||
|
|
/// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
|
||
|
|
/// et la matrice gamma telle que:
|
||
|
|
/// gamma(i,j) represente les coordonnees de la nouvelle base duale gpH dans l'ancienne gH
|
||
|
|
/// gpH(i) = gamma(i,j) * gH(j), i indice de ligne, j indice de colonne
|
||
|
|
/// c-a-d= gp^i = gamma^i_j * g^j
|
||
|
|
/// rappel des différentes relations entre beta et gamma
|
||
|
|
/// [beta]^{-1} = [gamma]^T ; [beta]^{-1T} = [gamma]
|
||
|
|
/// [beta] = [gamma]^{-1T} ; [beta]^{T} = [gamma]^{-1}
|
||
|
|
/// changement de base pour de une fois covariant:
|
||
|
|
/// [Ap_k] = [beta] * [A_i]
|
||
|
|
|
||
|
|
static void ChBase(const CoordonneeB& A_B,const Mat_pleine& beta,CoordonneeB& Ap_B);
|
||
|
|
|
||
|
|
/// changement de base pour de une fois contravariant:
|
||
|
|
/// [Ap^k] = [gamma] * [A^i]
|
||
|
|
static void ChBase(const CoordonneeH& A_H,const Mat_pleine& gamma,CoordonneeH& Ap_H);
|
||
|
|
|
||
|
|
};
|
||
|
|
/// @} // end of group
|
||
|
|
|
||
|
|
|
||
|
|
#ifndef MISE_AU_POINT
|
||
|
|
#include "MathUtil2.cc"
|
||
|
|
#define MATHUTIL2_H_deja_inclus
|
||
|
|
#endif
|
||
|
|
|
||
|
|
|
||
|
|
#endif
|