//#include "debug.h"
#include "Base.h"
#include "Util.h"


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

//------- contient que les méthodes qui ne sont pas inline ---------


//===============================================================
// Les base covariantes et absolus
//===============================================================

    

    // soient une base globale orthonormée : I_a et une base locale orthonormée I'_alpha
    // pour l'instant: I_a est en 3 dimensions et alpha varie de 1 à 2
    // soient une base naturelle g_alpha et duale associée g^beta
    // telles que : I'_alpha appartient à l'espace des g_alpha ou g^beta
    // on complète avec: I'_3 = I'_1 vectoriel I'_2
    // A) on a:
    // g^alpha = gamma(alpha,beta) * I'^beta  c-a-d : gamma^alpha_beta
    // NB: 1) comme I'est orthonormée: I'^beta = I'_beta
    //     2) chaque ligne de gamma représente un g^alpha
    // B) dans ce contexte on a:
    // inverse de gamma -> une matrice beta et chaque colonne de beta représente
    // les coordonnées de g_alpha dans I'^alpha
    // c-a-d g_alpha = beta (delta,alpha) * I'_delta
    // C) on peut également calculer les coordonnées de I'_beta dans le repère globale
    // I'_beta = gamma(alpha,beta) * g_alpha
    //
    // la méthode calcule à partir de this qui correspond à g_alpha dans I_a:
    //   apB :  g_alpha dans I'^alpha
    //   apH :  g^alpha dans I'^alpha
    //   IpB :  I'_beta dans I_a puis I'_3 = I'_1 vectoriel I'_2
void BaseB::ChangeBase_curviligne( const Mat_pleine& gamma, BaseB& apB , BaseH& apH, BaseB& IpB ) const
 {// vérif en cas de mise au point
  #ifdef MISE_AU_POINT
  if (dimension!= 3)
    { cout << "\n erreur, la methode n'est pour l'instant utilisable que pour la dimension 3 "
           << "\nBaseB::ChangeBase_curviligne(..";
      Sortie(1);
    };
  if ((gamma.Nb_ligne () != 2) || (gamma.Nb_colonne () != 2))
    { cout << "\n erreur, la matrice de passage gamma de nombre de ligne : " << gamma.Nb_ligne ()
           << ",  et de nombre de colonne  " << gamma.Nb_colonne ()
           << " doit avoir une dimension 2x2 : on ne traite pas les autres cas pour l'instant"
           << "\n BaseB::ChangeBase_curviligne(... ";
      Sortie(1);
    };
  #endif
  // 1) calcul de la base apB c-a-d de g_alpha dans la base I'
  //    calcul de alpha (delta,alpha)
  Mat_pleine beta =  gamma.Inverse();
  // g_alpha = beta (delta,alpha) * I'_delta
  int dim_i = gamma.Nb_ligne(); int dim_j = gamma.Nb_colonne(); // pour la généralisation éventuelle
  for (int al = 1; al <= dim_i; al++)
    for (int be = 1; be <= dim_j; be++)
      apB.CoordoB(al)(be) = beta(be,al);
  // 2) écriture de la base apH c-a-d de g^alpha dans la base I'
  // g^alpha = gamma(alpha,beta) * I'^beta
  for (int al = 1; al <= dim_i; al++)
    for (int be = 1; be <= dim_j; be++)
      apH.CoordoH(al)(be) = gamma(al,be);
  // 3) calcul des coordonnées de I' dans I_a
  for (int al = 1; al <= dim_i; al++)
   {IpB.CoordoB(al).Zero();
    for (int be = 1; be <= dim_j; be++)
      IpB.CoordoB(al) += gamma(be,al) * v(be);
   };
  // dans le cas en dimension 3 on complète
  if (ParaGlob::Dimension() == 3)
    { apB.CoordoB(1)(3) = apB.CoordoB(2)(3) = 0; // rien concernant la direction 3
      if (apB.NbVecteur() == 3)
         apB.CoordoB(3).Zero();
      apH.CoordoH(1)(3) = apH.CoordoH(2)(3) = 0;
      if (apH.NbVecteur() == 3)
         apH.CoordoH(3).Zero();
      IpB.CoordoB(3) = Util::ProdVec_coorB(IpB(1),IpB(2));
    };
 };


 
//===============================================================
// Les base duales
//===============================================================
    
    


    // soient une base globale orthonormée : I_a et une base locale orthonormée I'_alpha
    // pour l'instant: I_a est en 3 dimensions et alpha varie de 1 à 2
    // soient une base naturelle g_alpha et duale associée g^beta
    // telles que : I'_alpha appartient à l'espace des g_alpha ou g^beta
    // A) on a:
    // g^alpha = gamma(alpha,beta) * I'^beta  c-a-d : gamma^alpha_beta
    // NB: 1) comme I'est orthonormée: I'^beta = I'_beta
    //     2) chaque ligne de gamma représente un g^alpha
    // B) dans ce contexte on a:
    // inverse de gamma -> une matrice beta et chaque ligne de beta représente
    // les coordonnées de g_alpha dans I'^alpha
    // c-a-d g_alpha = beta (delta,alpha) * I'_delta
    // C) on peut également calculer les coordonnées de I'_beta dans le repère globale
    // I'^beta = beta(beta,al) * g^alpha
    //
    // la méthode calcule à partir de this qui correspond à g^alpha dans I_a:
    //   apB :  g_alpha dans I'^alpha
    //   apH :  g^alpha dans I'^alpha
    //   IpBH :  I'^beta dans I_a  puis I'_3 = I'_1 vectoriel I'_2
void BaseH::ChangeBase_curviligne( const Mat_pleine& gamma, BaseB& apB , BaseH& apH, BaseH& IpH ) const
 {// vérif en cas de mise au point
  #ifdef MISE_AU_POINT
  if (dimension!= 3)
    { cout << "\n erreur, la methode n'est pour l'instant utilisable que pour la dimension 3 "
           << "\nBaseH::ChangeBase_curviligne(..";
      Sortie(1);
    };
  if ((gamma.Nb_ligne () != 2) || (gamma.Nb_colonne () != 2))
    { cout << "\n erreur, la matrice de passage gamma de nombre de ligne : " << gamma.Nb_ligne ()
           << ",  et de nombre de colonne  " << gamma.Nb_colonne ()
           << " doit avoir une dimension 2x2 : on ne traite pas les autres cas pour l'instant"
           << "\n BaseH::ChangeBase_curviligne(... ";
      Sortie(1);
    };
  #endif
  // 1) calcul de la base apB c-a-d de g_alpha dans la base I'
  //    calcul de alpha (delta,alpha)
  Mat_pleine beta =  gamma.Inverse();
  // g_alpha = beta (delta,alpha) * I'_delta
  int dim_i = gamma.Nb_ligne(); int dim_j = gamma.Nb_colonne(); // pour la généralisation éventuelle
  for (int al = 1; al <= dim_i; al++)
    for (int be = 1; be <= dim_j; be++)
      apB.CoordoB(al)(be) = beta(be,al);
  // 2) écriture de la base apH c-a-d de g^alpha dans la base I'
  // g^alpha = gamma(alpha,beta) * I'^beta
  for (int al = 1; al <= dim_i; al++)
    for (int be = 1; be <= dim_j; be++)
      apH.CoordoH(al)(be) = gamma(al,be);
  // 3) calcul des coordonnées de I' dans I_a
  for (int al = 1; al <= dim_i; al++)
   { IpH.CoordoH(al).Zero();
    for (int be = 1; be <= dim_j; be++)
      IpH.CoordoH(al) += beta(al,be) * v(be);
   };
  // dans le cas en dimension 3 on complète
  if (ParaGlob::Dimension() == 3)
    { apB.CoordoB(1)(3) = apB.CoordoB(2)(3) = 0; // rien concernant la direction 3
      if (apB.NbVecteur() == 3)
         apB.CoordoB(3).Zero();
      apH.CoordoH(1)(3) = apH.CoordoH(2)(3) = 0;
      if (apH.NbVecteur() == 3)
         apH.CoordoH(3).Zero();
     IpH.CoordoH(3) = Util::ProdVec_coorH(IpH(1),IpH(2));
    };
 };