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

# include "Met_pout2D.h"   
# include "Util.h"

// constructeurs

Met_Pout2D::Met_Pout2D () : // constructeur par defaut mettant les pointeurs a NULL
  Met_PiPoCo()
		{ };  
	    // constructeur permettant de dimensionner  certaine variables
	    // dim = dimension de l'espace, nbvec = nb de vecteur des bases, tab = liste
	    // des variables a initialiser
Met_Pout2D::Met_Pout2D (int dim_base,int nbvec,const  DdlElement& tabddl,
	                   const  Tableau<Enum_variable_metrique>& tabb,int nomb_noeud):
	    Met_PiPoCo(dim_base,nbvec,tabddl,tabb,nomb_noeud)             
	    {  };
	    // Constructeur de copie :
Met_Pout2D::Met_Pout2D (const   Met_Pout2D& a) :
  Met_PiPoCo(a)
  { }
Met_Pout2D::~Met_Pout2D ()
		{  };
//============================ METHODES PUBLIQUES : ============================
    // Surcharge de l'operateur = : realise l'affectation
    // dérivant de virtuel, a ne pas employer -> message d'erreur
Met_abstraite& Met_Pout2D::operator= (const Met_abstraite& )
  {	// normalement ne devrait pas être utilisé
    cout << "\n erreur , l operateur d affectation a utiliser  doit etre celui explicite "
         << " de la classe Met_Pout2D \n"
         << " Met_abstraite& Met_Pout2D::operator= (const Met_abstraite& met) \n";
    Sortie(1);     
    return (*this);	
   };
Met_PiPoCo& Met_Pout2D::operator= (const Met_PiPoCo& )
  {	// normalement ne devrait pas être utilisé
    cout << "\n erreur , l operateur d affectation a utiliser  doit etre celui explicite "
         << " de la classe Met_Pout2D \n"
         << " Met_PiPoCo& Met_Pout2D::operator= (const Met_PiPoCo& met) \n";
    Sortie(1);     
    return (*this);	
   };
	// normale
Met_Pout2D& Met_Pout2D::operator= (const Met_Pout2D& met)
  {	(*this) = Met_PiPoCo::operator=(met);
    tabD2phi = met.tabD2phi;
    return (*this);	
  };
	  // passage de la dérivée seconde
void Met_Pout2D::DeriveeSeconde(Vecteur  const & taD2phi)
 { tabD2phi = &taD2phi;
  };  

// ============================ méthodes protegees ===========================
// calcul des normales a la facette
void Met_Pout2D::Calcul_N_0 ()
  { N_0(1) = -((*aiB_0)(1))(2); N_0(2) = ((*aiB_0)(1))(1);
    N_0.Normer (); 
  };
void Met_Pout2D::Calcul_N_t ()
  { N_t(1) = -((*aiB_t)(1))(2); N_t(2) = ((*aiB_t)(1))(1);
    N_t.Normer (); 
  };
void Met_Pout2D::Calcul_N_tdt ()
  { N_tdt(1) = -((*aiB_tdt)(1))(2); N_tdt(2) = ((*aiB_tdt)(1))(1);
    N_tdt.Normer (); 
  };
// calcul de la base  naturel a t0  
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de l'axe médian
void Met_Pout2D::Calcul_giB_0
 (  const   Tableau<Noeud *>& , const  Mat_pleine& , int ,const   Vecteur& phi)
{
  #ifdef MISE_AU_POINT
	  if  (giB_0 == NULL)
     { cout << "\nErreur : la base a t=0 n'est pas dimensionne !\n";
       cout << "void Met_Pout2D::Calcul_giB_0 \n";
       Sortie(1);
     };
  #endif
  //  derivee du vecteur normal
  // l'écriture est un peu particulière pour éviter de mettre en oeuvre plein de constructeur
  CoordonneeB dN1; dN1.ConstructionAPartirDe_H(curb_0(1) * (*aiH_0)(1)) ;  
  // vecteur de base gi
  // le premier element de phi est en fait est = a : " z  dans l'epaisseur " 
  giB_0->CoordoB(1) = (*aiB_0)(1) + phi(1) * dN1;
};

// calcul de la base  naturel a t
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de l'axe médian
void Met_Pout2D::Calcul_giB_t
 ( const   Tableau<Noeud *>& , const  Mat_pleine& , int,const   Vecteur& phi)
{
  #ifdef MISE_AU_POINT
	  if  (giB_t == NULL)
    { cout << "\nErreur : la base a t n'est pas dimensionne !\n";
      cout << "void Met_Pout2D::Calcul_giB_t \n";
      Sortie(1);
    };
  #endif
  //  derivee du vecteur normal
  // l'écriture est un peu particulière pour éviter de mettre en oeuvre plein de constructeur
  CoordonneeB dN1; dN1.ConstructionAPartirDe_H(curb_t(1) * (*aiH_t)(1)) ; 
  // vecteur de base gi
  // le premier element de phi est en fait est = a : " z  dans l'epaisseur " 
  giB_t->CoordoB(1) = (*aiB_t)(1) + phi(1) * dN1;
};

// calcul de la base  naturel a tdt
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de l'axe médian
void Met_Pout2D::Calcul_giB_tdt
 ( const   Tableau<Noeud *>& ,const   Mat_pleine& , int,const   Vecteur& phi)
{
  #ifdef MISE_AU_POINT
	  if  (giB_tdt == NULL)
     { cout << "\nErreur : la base a t+dt n'est pas dimensionne !\n";
       cout << "void Met_Pout2D::Calcul_giB_tdt \n";
       Sortie(1);
     };
  #endif
  // derivee du vecteur normal
  // l'écriture est un peu particulière pour éviter de mettre en oeuvre plein de constructeur
  CoordonneeB dN1; dN1.ConstructionAPartirDe_H(curb_tdt(1) * (*aiH_tdt)(1)); 
  // vecteur de base gi
  // le premier element de phi est en fait est = a : " z  dans l'epaisseur " 
  giB_tdt->CoordoB(1) = (*aiB_tdt)(1) + phi(1) * dN1;  
};

//------------// variation des vecteurs de base
 void Met_Pout2D::D_giB_t(const   Mat_pleine& , int ,const   Vecteur & phi)
  { for (int iddl=1;iddl<= 6;iddl++)
     { // variation de la derivee du vecteur normal
	      // l'écriture est un peu particulière pour éviter de mettre en oeuvre plein de constructeur
	      CoordonneeB ddN1; ddN1.ConstructionAPartirDe_H(dcurb_t(iddl)(1) * (*aiH_t)(1)
	                       + curb_t(1) * (*d_aiH_t)(iddl)(1)); 
	      // vecteur de base gi
	      // le premier element de phi est en fait est = a : " z * phi dans l'epaisseur "
	      (*d_giB_t)(iddl).CoordoB(1) = (*d_aiB_t)(iddl)(1) + phi(1) * ddN1;
     };
  };

 void Met_Pout2D::D_giB_tdt(const   Mat_pleine& , int ,const   Vecteur & phi)
  { for (int iddl=1;iddl<= 6;iddl++)
     { // variation de la derivee du vecteur normal
       // l'écriture est un peu particulière pour éviter de mettre en oeuvre plein de constructeur
	      CoordonneeB ddN1; ddN1.ConstructionAPartirDe_H(dcurb_tdt(iddl)(1) * (*aiH_tdt)(1)
	                                       + curb_tdt(1) * (*d_aiH_tdt)(iddl)(1)); 
	      // vecteur de base gi
	      // le premier element de phi est en fait est = a : " z * phi dans l'epaisseur "
	      (*d_giB_tdt)(iddl).CoordoB(1) = (*d_aiB_tdt)(iddl)(1) + phi(1) * ddN1;
     };
  };

	
    //-----------// calcul du tenseur de courbure  dans la base naturelle
    // plusieurs cas sont etudies suivant l'instant considere
    // a l'instant t = 0
Vecteur& Met_Pout2D::courbure_0 (const  Tableau<Noeud *>& tab_noeud)
  { const CoordonneeB & aiB1 = (*aiB_0)(1); // pour simplifier
    const CoordonneeH & aiH1 = (*aiH_0)(1); // pour simplifier
    Tableau<Coordonnee> tab_coor(nomb_noeud);
    for (int i=1;i<=nomb_noeud;i++)
      tab_coor(i) = tab_noeud(i)->Coord0();
    curb_0(1) =  courbure(tab_coor,aiB1,aiH1);
    return curb_0;
   };      

    // a l'instant t
Vecteur& Met_Pout2D::courbure_t (const  Tableau<Noeud *>& tab_noeud)
  { const CoordonneeB & aiB1 = (*aiB_t)(1); // pour simplifier
    const CoordonneeH & aiH1 = (*aiH_t)(1); // pour simplifier
    Tableau<Coordonnee> tab_coor(nomb_noeud);
    for (int i=1;i<=nomb_noeud;i++)
      tab_coor(i) = tab_noeud(i)->Coord1();
    curb_t(1) =  courbure(tab_coor,aiB1,aiH1);
    return curb_t;
   };      

    // a l'instant t+dt
Vecteur& Met_Pout2D::courbure_tdt (const  Tableau<Noeud *>& tab_noeud)
  { const CoordonneeB & aiB1 = (*aiB_tdt)(1); // pour simplifier
    const CoordonneeH & aiH1 = (*aiH_tdt)(1); // pour simplifier
    Tableau<Coordonnee> tab_coor(nomb_noeud);
    for (int i=1;i<=nomb_noeud;i++)
      tab_coor(i) = tab_noeud(i)->Coord2();
    curb_tdt(1) =  courbure(tab_coor,aiB1,aiH1);
    return curb_tdt;
   };      

    // routine generale de calcul de la courbure
double Met_Pout2D::courbure(const  Tableau<Coordonnee>& tab_coor,const  CoordonneeB & aiB1,const  CoordonneeH & )
   { // dans le cas de la poutre 2D, il n'y a qu'une seule courbure
     double  curb;
     // calcul de la normale qui est normale à aiB1
     Coordonnee N(2); N(1) = -aiB1(2);N(2) = aiB1(1);  N.Normer ();
     // calcul de la dérivée du vecteur aiB1
     Coordonnee asiB1(2);
     for (int a=1;a<= 2;a++)
	     {asiB1(a) = 0.;
		     for (int r=1;r<=tabD2phi->Taille();r++)
			      asiB1(a) += tab_coor(r)(a) * (*tabD2phi)(r);
      };
     // calcul de la courbure
     curb = asiB1 * N;
     return curb;
   };
     			
    //--------------// calcul du tenseur de courbure et de sa variation 
    // plusieurs cas sont etudies suivant l'instant considere
    // a l'instant t 
void Met_Pout2D::Dcourbure_t(const  Tableau<Noeud *>& tab_noeud,
               Vecteur& curb,TabOper<Vecteur>& dcurb)
  { const CoordonneeB & aiB1 = (*aiB_t)(1); // pour simplifier
    const CoordonneeH & aiH1 = (*aiH_t)(1);  // pour simplifier
    Tableau<Coordonnee> tab_coor(nomb_noeud);
    for (int i=1;i<=nomb_noeud;i++)
      tab_coor(i) = tab_noeud(i)->Coord1();
    Dcourbure (tab_coor,aiB1,aiH1,*d_aiB_t,curb,dcurb);
   };      

    // a l'instant tdt
void Met_Pout2D::Dcourbure_tdt(const  Tableau<Noeud *>& tab_noeud,
                Vecteur& curb,TabOper<Vecteur>& dcurb)
  { const CoordonneeB & aiB1 = (*aiB_tdt)(1); // pour simplifier
    const CoordonneeH & aiH1 = (*aiH_tdt)(1);  // pour simplifier
    Tableau<Coordonnee> tab_coor(nomb_noeud);
    for (int i=1;i<=nomb_noeud;i++)
      tab_coor(i) = tab_noeud(i)->Coord2();
    Dcourbure (tab_coor,aiB1,aiH1,*d_aiB_tdt,curb,dcurb);
   };      

    // routine generale de calcul de la courbure et de sa variation
    // en sortie : curb , la courbure  b11
    //            dcurb , la variation de courbure.
void Met_Pout2D::Dcourbure (const  Tableau<Coordonnee>& tab_coor,const CoordonneeB & aiB1
                    ,const CoordonneeH & ,Tableau <BaseB>& DaiB
                    ,Vecteur& curb,TabOper <Vecteur>& dcurb)
   { // dans le cas de la poutre 2D, il n'y a qu'une seule courbure
     // calcul de la normale qui est normale à aiB1  
     Coordonnee N(2); N(1) = -aiB1(2);N(2) = aiB1(1);  N.Normer ();
     // calcul de la dérivée du vecteur aiB1
     Coordonnee asiB1(2);
     for (int a=1;a<= 2;a++)
	      {asiB1(a) = 0.;
		      for (int r=1;r<=tabD2phi->Taille();r++)
			        asiB1(a) += tab_coor(r)(a) * (*tabD2phi)(r);
	      };
     // calcul de la courbure
     curb(1) = asiB1 * N;
     // maintenant le cas des variations
     int nbddl = dim_base * nomb_noeud;
     Coordonnee bidon(2) ; // pour le dimensionnement de dN et de dasiB1
     TabOper <Coordonnee> dN(nbddl,bidon); // variation de la normale
     TabOper <Coordonnee> dasiB1(nbddl,bidon) ;// " " de la dérivée du vecteur a1
     // on fait une boucle sur les degrés de liberté
     int indice;
     // normalement dim_base = 2
     for (int b = 1; b<= dim_base; b++)
      for (int r = 1; r <= nomb_noeud; r++)
         { indice = (r-1)*dim_base + b;
           dN(indice)(1) = -DaiB(indice)(1)(2);
           dN(indice)(2) =  DaiB(indice)(1)(1);
           dasiB1(indice)(b) = (*tabD2phi)(r);
         };
     for (indice =1;indice <= nbddl;indice++)
        dcurb(indice)(1) = dasiB1(indice) * N + asiB1 * dN(indice);      
   };