// 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  
