// 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 "Debug.h"
# include "DeformationPP.h"
# include "Met_PiPoCo.h"
#include "ConstMath.h"
#include "ParaGlob.h"
#include "Util.h"

// il faut regarder et mettre à jour le calcul des différentes déformations ********
// il faut regarder dans déformation, en particulier il faut mettre à jour,           ***********
// si nécessaire différents cas de déformation              **********************
// il faut également mettre à jour l'influence de l'indicateur premier_calcul ????????????
// dans le cas du premier passage, pour l'instant pas fait !!!!!!!!!!!!!!!!!!!!!!!!!!
// -------------------------- variables static -------------------

int DeformationPP::numInteg_ep = 0; 
int DeformationPP::nbtotalep = 0; 
int DeformationPP::nbtotalaxpl = 0; 
Vecteur DeformationPP::theta_z(1); 
// on utilise deux types d'interpolation
Tableau <Mat_pleine const *>  DeformationPP::taDphi(2); 
Tableau <Vecteur const *>  DeformationPP::taPhi(2);  

// ---------- constructeur----------------------------------------

DeformationPP::DeformationPP ()   // constructeur ne devant pas etre utilise
   {
     #ifdef MISE_AU_POINT	 	 
		{ cout << "\nErreur : le constructeur par defaut ne doit pa etre utilise !\n";
			  cout << "DeformationPP::DeformationPP () \n";
			  Sortie(1);
			};
     #endif
   };
// constructeur normal dans le cas d'un ou de plusieurs pt d'integration             	  
DeformationPP::DeformationPP (Met_abstraite & a,Tableau<Noeud *>& tabnoeud    
           ,Tableau <Mat_pleine>  const & tDphH,Tableau <Vecteur>  const & tPhH
           ,Tableau <Mat_pleine>  const & tDphS,Tableau <Vecteur>  const & tPhS ): 
   Deformation (a,tabnoeud,tDphS,tPhS )
   {  type_deformation = DEFORMATION_POUTRE_PLAQUE_STANDART;
      tabDphiH = &(tDphH);
      tabPhiH = &(tPhH);
      // a priori pas  de numero d'integration affectes
      numInteg = 0;
      numInteg_ep = 0;
    };
   
// constructeur de copie
DeformationPP::DeformationPP (const  DeformationPP& a) :
   Deformation(a) 
   { numInteg_ep = a.numInteg_ep;
     tabDphiH = a.tabDphiH;
     tabPhiH = a.tabPhiH;
    };
     
DeformationPP::~DeformationPP () 
   { }; 
    
//  ============      METHODES PUBLIQUES :  ==================
    // Surcharge de l'operateur = : realise l'affectation
    // fonction virtuelle
Deformation& DeformationPP::operator= (const Deformation& )
  {	// normalement ne devrait pas être utilisé
    cout << "\n erreur , l operateur d affectation a utiliser  doit etre celui explicite "
         << " de la classe DeformationPP \n"
         << " Deformation& DeformationPP::operator= (const Deformation& def) \n";
    Sortie(1);     
    return (*this);	
   };
    // fonction normale
DeformationPP& DeformationPP::operator= (const DeformationPP& def)
  {	tabDphiH = def.tabDphiH;
    tabPhiH = def.tabPhiH;
    (*this) = Deformation::operator=(def);
    return (*this);	
   };
    // change les numeros d'integration de surface et d'epaisseur courant    
    void DeformationPP::ChangeNumIntegSH
            (int nisurf, int niepaiss)
      {  //PiPoCo & elem = *((PiPoCo*) element);
         numInteg_ep = niepaiss; numInteg = nisurf;
         // on fait les changements dus aux nouveaux pts integ
         AppliquePtInteg();
       } ;

     // calcul explicit à t :tous les parametres sont de resultats
const Met_abstraite::Expli&  DeformationPP::Cal_explicit_t
                       ( const Tableau <double>& def_equi_t,TenseurBB & epsBB_t,Tableau <TenseurBB *> & d_epsBB
							   ,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& DeltaEpsBB,bool premier_calcul)
  {  bool gradV_instantane = false; // ************ pour l'instant figé
     // appel de la metrique
     Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
     const Met_abstraite::Expli& ex = met_pipoco->Cal_explicit_t
       (*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
     
 /// voir Deformation car ça a changé !!
/*    // epsBB_t est dimensionne dans l'element              
     epsBB_t = 0.5 * ( *(ex.gijBB_t) - *(ex.gijBB_0)); 
     // ici on considère que le delta_epsBB=epsBB_t ! c-a-d le delta de 0 à t
     DeltaEpsBB = epsBB_t;
     // calcul simplifié éventuelle de la vitesse de déformation si les ddl vitesses ne sont pas présent
     if (pas_de_gradV)
      { // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
	    const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
	    // dans le cas où l'incrément de temps est nul on garde la vitesse actuelle
	    double deltat=vartemps.IncreTempsCourant();
	    if (deltat >= ConstMath::trespetit)
    	    DepsBB = DeltaEpsBB/deltat;
	   }
	 else // cas où le gradient a été calculé
	   DepsBB = 0.5 * ((*ex.gradVBB_t) + ex.gradVBB_t->Transpose());
	 // variation de la déformation / au ddl  
     int  d_epsBBTaille = d_epsBB.Taille();                                          
     for (int i=1; i<= d_epsBBTaille; i++)
         *(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_t))(i)));*/
     return ex;           
  };

     // calcul explicit à tdt :tous les parametres sont de resultats
const Met_abstraite::Expli_t_tdt&  DeformationPP::Cal_explicit_tdt
                  ( const Tableau <double>& def_equi_t,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB
						  ,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& delta_epsBB_tdt,bool premier_calcul)

  {  bool gradV_instantane = false; // ************ pour l'instant figé
     // appel de la metrique
     Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
     const Met_abstraite::Expli_t_tdt& ex= met_pipoco->Cal_explicit_tdt
                (*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);

 /// voir Deformation car ça a changé !!
/*     // epsBB_tdt est dimensionne auparavant              
     epsBB_tdt = 0.5 * ( *(ex.gijBB_t) - *(ex.gijBB_0)); 
     // epsBB_tdt est dimensionne auparavant             
     delta_epsBB_tdt = 0.5 * ( *(ex.gijBB_tdt) - *(ex.gijBB_t)); 
     if (pas_de_gradV)
      { // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
	    const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
	    // dans le cas où l'incrément de temps est nul on garde la vitesse actuelle
	    double deltat=vartemps.IncreTempsCourant();
	    if (deltat >= ConstMath::trespetit)
    	    DepsBB = delta_epsBB_tdt/deltat;
	   } 
	 else // cas où le gradient a été calculé
	   DepsBB = 0.5 * ((*ex.gradVBB_tdt) + ex.gradVBB_tdt->Transpose());
	 // variation de la déformation / au ddl  
     int  d_epsBBTaille = d_epsBB.Taille();                                          
     for (int i=1; i<= d_epsBBTaille; i++)
         *(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));*/
     return ex;           
  };

  	  // cas implicite : tous les parametres sont de resultats	
const Met_abstraite::Impli& DeformationPP::Cal_implicit
         ( const Tableau <double>& def_equi_t,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
			  ,Tableau <double>& def_equi,Tableau2 <TenseurBB *>& d2_epsBB_tdt
			  ,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
  {  bool gradV_instantane = false; // ************ pour l'instant figé
     // recup d'un pointeur sur la metrique sfe1
     Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
     // toutes les variables de passage a metrique apres l'appel
     // pointeront sur des variables  deja dimensionnees
     // pour les Tableau <> il y a dimensionnement auto a l'affectation
     const Met_abstraite::Impli& ex =met_pipoco->Cal_implicit ( *tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul); 
		  	
 /// voir Deformation car ça a changé !!
 
 
 
/*      epsBB_tdt = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_0));
     delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
     if (pas_de_gradV)
      { // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
	    const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
	    // dans le cas où l'incrément de temps est nul on garde la vitesse actuelle
	    double deltat=vartemps.IncreTempsCourant();
	    if (deltat >= ConstMath::trespetit)
	        DepsBB = delta_epsBB/deltat;
	   } 
	 else // cas où le gradient a été calculé
	   DepsBB = 0.5 * ((*ex.gradVBB_tdt) + ex.gradVBB_tdt->Transpose());
	 // variation de la déformation / au ddl  
     int d_epsBB_tdtTaille = d_epsBB_tdt.Taille();		      
     for (int i=1; i<= d_epsBB_tdtTaille; i++)
         *(d_epsBB_tdt(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
     // calcul éventuelle de la dérivée seconde de la déformation
     if (pa.Var_D())    
	  {// recup des variations secondes de la déformation
	   int d2_epsBB_tdtTaille1 = d2_epsBB_tdt.Taille1();
       for (int i=1; i<= d2_epsBB_tdtTaille1; i++)
        for (int j=1; j<= i; j++)  // symétrie du  tableau et du tenseur
	      { *(d2_epsBB_tdt(i,j)) = 0.5 * (*((*(ex.d2_gijBB_tdt))(i,j))) ;
	        *(d2_epsBB_tdt(j,i)) = *(d2_epsBB_tdt(i,j)) ; 
	       }
	   }  */  
     return ex;
  };    
// ---------------- calcul des variables primaires autre que pour la mécanique -------- 
// ------------     donc par de retour relatif aux déformations   
     // calcul explicit à t :tous les parametres sont de resultats
const Met_abstraite::Expli&  DeformationPP::Cal_explicit_t(bool premier_calcul )

  {  bool gradV_instantane = false; // ************ pour l'instant figé
     // appel de la metrique
     Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
     // appel de la metrique
     const Met_abstraite::Expli& ex = met_pipoco->Cal_explicit_t(*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
     return ex;           
  };

     // calcul explicit à tdt :tous les parametres sont de resultats
const Met_abstraite::Expli_t_tdt&  DeformationPP::Cal_explicit_tdt(bool premier_calcul )

  {  bool gradV_instantane = false; // ************ pour l'instant figé
     // appel de la metrique
     Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
     // appel de la metrique
     const Met_abstraite::Expli_t_tdt& ex = met_pipoco->Cal_explicit_tdt(*tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul);
     return ex;           
  };

  	  // cas implicite : tous les parametres sont de resultats	
const Met_abstraite::Impli& DeformationPP::Cal_implicit( bool premier_calcul)
  {  bool gradV_instantane = false; // ************ pour l'instant figé
     // recup d'un pointeur sur la metrique sfe1
     Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
     // toutes les variables de passage a metrique apres l'appel
     // pointeront sur des variables  deja dimensionnees
     // pour les Tableau <> il y a dimensionnement auto a l'affectation
     // appel de la metrique
     const Met_abstraite::Impli& ex =met_pipoco->Cal_implicit ( *tabnoeud,gradV_instantane,taDphi,nbNoeud,taPhi,premier_calcul); 
		  	
     return ex;
  };    
  
// ========== remontee aux informations =========================
		
// cas sortie d'un calcul implicit
// Aa(i,a) = Aa^i_{.a}, avec  g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
const Met_abstraite::InfoImp DeformationPP::RemontImp(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
  { // recup d'un pointeur sur la metrique sfe1
    Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
    Met_abstraite::InfoImp ex = met_pipoco->Cal_InfoImp
        ( *tabnoeud,taDphi,taPhi,nbNoeud);
    // determination des matrices de transformation de base
    const BaseB & giB0 = *(ex.giB_0);
    const BaseB & giB = *(ex.giB_tdt);
    const BaseH & giH0 = *(ex.giH_0);
    const BaseH & giH = *(ex.giH_tdt);
    // on différencie le cas des poutres et celui des plaques
    int dim = giB0.NbVecteur();
    if (dim == 1) // cas de la dimension 1 c-a-dire poutre
        Deformation::BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin); // générale
    else         
        BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);  // particulier plaques
    return ex;        
   };

// cas sortie d'un calcul implicit
// sans les matrices de passage
const Met_abstraite::InfoImp DeformationPP::RemontImp()
  { // recup d'un pointeur sur la metrique sfe1
    Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
    Met_abstraite::InfoImp ex = met_pipoco->Cal_InfoImp
        ( *tabnoeud,taDphi,taPhi,nbNoeud);
    // determination des matrices de transformation de base
    const BaseB & giB0 = *(ex.giB_0);
    const BaseB & giB = *(ex.giB_tdt);
    const BaseH & giH0 = *(ex.giH_0);
    const BaseH & giH = *(ex.giH_tdt);
    return ex;
   };

// cas sortie d'un calcul explicit à t
// Aa(i,a) = Aa^i_{.a}, avec  g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
const Met_abstraite::InfoExp_t DeformationPP::RemontExp_t(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
  { // recup d'un pointeur sur la metrique sfe1
    Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
    Met_abstraite::InfoExp_t ex = met_pipoco->Cal_InfoExp_t
                 ( *tabnoeud,taDphi,taPhi,nbNoeud);
    // determination des matrices de transformation de base
    const BaseB & giB0 = *(ex.giB_0);
    const BaseB & giB = *(ex.giB_t);
    const BaseH & giH0 = *(ex.giH_0);
    const BaseH & giH = *(ex.giH_t);
    BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
    return ex;      
  };

// cas sortie d'un calcul explicit à t
// sans les matrices de passage
const Met_abstraite::InfoExp_t DeformationPP::RemontExp_t()
  { // recup d'un pointeur sur la metrique sfe1
    Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
    Met_abstraite::InfoExp_t ex = met_pipoco->Cal_InfoExp_t
                 ( *tabnoeud,taDphi,taPhi,nbNoeud);
    // determination des matrices de transformation de base
    const BaseB & giB0 = *(ex.giB_0);
    const BaseB & giB = *(ex.giB_t);
    const BaseH & giH0 = *(ex.giH_0);
    const BaseH & giH = *(ex.giH_t);
    return ex;
  };
   
// cas sortie d'un calcul explicit à tdt
// Aa(i,a) = Aa^i_{.a}, avec  g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
const Met_abstraite::InfoExp_tdt DeformationPP::RemontExp_tdt(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
  { // recup d'un pointeur sur la metrique sfe1
    Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
    Met_abstraite::InfoExp_tdt ex = met_pipoco->Cal_InfoExp_tdt
                 ( *tabnoeud,taDphi,taPhi,nbNoeud);
    // determination des matrices de transformation de base
    const BaseB & giB0 = *(ex.giB_0);
    const BaseB & giB = *(ex.giB_tdt);
    const BaseH & giH0 = *(ex.giH_0);
    const BaseH & giH = *(ex.giH_tdt);
    BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
    return ex;      
  };

// cas sortie d'un calcul explicit à tdt
// sans les matrices de passage
const Met_abstraite::InfoExp_tdt DeformationPP::RemontExp_tdt()
  { // recup d'un pointeur sur la metrique sfe1
    Met_PiPoCo* met_pipoco = (Met_PiPoCo*) metrique;
    Met_abstraite::InfoExp_tdt ex = met_pipoco->Cal_InfoExp_tdt
                 ( *tabnoeud,taDphi,taPhi,nbNoeud);
    // determination des matrices de transformation de base
    const BaseB & giB0 = *(ex.giB_0);
    const BaseB & giB = *(ex.giB_tdt);
    const BaseH & giH0 = *(ex.giH_0);
    const BaseH & giH = *(ex.giH_tdt);
    return ex;
  };
  
// gestion du parcours de tous les points d'integration
// bien voir que numInteg est relatif aux pt d'intégration de l'axe
void DeformationPP::PremierPtInteg()
   {  nbtotalaxpl = tabPhi->Taille();
      nbtotalep = tabPhiH->Taille();
      numInteg = 1 ; numInteg_ep = 1; 
     // on fait les changements dus aux nouveaux pts integ
     AppliquePtInteg();
    };
bool DeformationPP::DernierPtInteg()
   { if (numInteg <= nbtotalaxpl) 
       return true;
     else
       return false;  
    };
void DeformationPP::NevezPtInteg()
   { if (numInteg_ep < nbtotalep)
       numInteg_ep++;
     else
      { numInteg++;
        numInteg_ep = 1;
       }
     // on fait les changements dus aux nouveaux pts integ
     // dans le cas ou ils sont valides
     if (numInteg <= nbtotalaxpl)
        AppliquePtInteg();
    };
    // ------------------------ METHODES PROTEGEES : -------------------------------
  
// Aa(i,a) = Aa^i_{.a}, avec  g^i = Aa^i_{.a} * I^a
// tout ce passe comme si I^a est la nouvelle base vers laquelle on veut évoluer
void DeformationPP::BasePassage
        (bool absolue,const BaseB & giB0,const BaseB & giB,const BaseH & giH0,const BaseH & giH,
         Mat_pleine& Aa0,Mat_pleine& Aafin)
 {  
    // determination des matrices de transformation de base
    // l'objectif est de determiner un repere pertinant dans le cas de plaque , coque .
    // Dans le cas d'une dimension globale 3D avec des plaques ou coque ->
    // choix : un repere qui appartiend a la facette, obtenu apres projection
    // du repere global
    //------ cas de la config initiale, on regarde si la projection de I1 n'est pas trop petite
    // def de la normale a la facette
    Coordonnee N = (Util::ProdVec_coorBN(giB0(1),giB0(2))).Normer();
    Coordonnee ve(0.,-N(3),-N(2)); // = ProdVec(N,I1)
    double norve = ve.Norme();
    Tableau <Coordonnee> vi(2); // les vecteurs plans orthonormes de la facette
    if (norve >= ConstMath::petit)
     { vi(2) = ve.Normer();
       vi(1) = Util::ProdVec_coor(vi(2),N);
      }
    else
      { ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
        vi(1) = ve.Normer();
        vi(2) = Util::ProdVec_coor(N,vi(1));
       } 
    for (int al=1 ;al<=2;al++)
      {Aa0(1,al) = giH0.Coordo(1) * vi(al) ;
       Aa0(2,al) = giH0.Coordo(2) * vi(al) ;
       } 
    //------ cas de la config finale, 
    N = (Util::ProdVec_coorBN(giB(1),giB(2))).Normer();
    ve.Change_Coordonnee(3,0.,N(3),-N(2)); // = ProdVec(N,I1)
    norve = ve.Norme();
    Tableau <Coordonnee> Aa(2);
    if (norve >= ConstMath::petit)
     { Aa(2) = ve.Normer();
       Aa(1) = Util::ProdVec_coor(Aa(2),N);
      }
    else
      { ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
        Aa(1) = ve.Normer();
        Aa(2) = Util::ProdVec_coor(N,Aa(1));
       } 
    for (int be=1;be<=2;be++)
      { Aafin(1,be) = giH.Coordo(1) * Aa(be);
        Aafin(2,be) = giH.Coordo(2) * Aa(be);
       } 
    return;      
  };           

// applique les conséquences d'un changement de numéro de point d'intégration
// par exemple effectué via NevezPtInteg() ou ChangeNumIntegSH
    void DeformationPP::AppliquePtInteg()
      {  //PiPoCo & elem = *((PiPoCo*) element);
         // on est obligé de préciser le caractère réel du tableau ?
//------- pour l'instant c'est inopérationnelle: les 3 lignes qui suivent donne une erreur à la compilation
          
//         taPhi(1) = & ((Vecteur const)((*tabPhi)(numInteg)));
//         taDphi(1) = &((Mat_pleine const)((*tabDphi)(numInteg)));
//         taDphi(2) = &((Mat_pleine const)((*tabDphiH)(numInteg_ep)));
         // calcul de la cote en epaisseur
// les trois lignes suivantes ont été mise en commentaire car on ne peut plus
// accèder à l'élément directement, a moins par exemple de le passer en 
// paramètre ou de trouver une autre technique pour avoir les infos a doc
cout << "\n pour l'instant la méthode n'est pas disponible cf le listing";
cout << " void DeformationPP::AppliquePtInteg()";
Sortie(1);        
 //        double epaisseur = elem.H();  
 //        theta_z(1) = 0.5 * epaisseur * (elem.KSIepais(numInteg_ep));
 //        taPhi(2) = & theta_z;
       } ;