// 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 "PiPoCo.h"
#include <iomanip>
#include "ConstMath.h"
#include "Loi_Umat.h"
#include "ExceptionsElemMeca.h"
#include "ExceptionsLoiComp.h"

// =========================== variables communes ===============================
// a priori les pointeurs de fonction pointent sur rien au debut
TenseurHH* PiPoCo::sig_bulk_pourPiPoCo_HH = NULL;  // variable de travail pour le bulk
                    
// =========================== fin variables communes ===========================

// =========================== methodes publiques ===============================

// calcul si un point est a l'interieur de l'element ou non
// il faut que M est la dimension globale
// les trois fonctions sont pour l'etude a t=0, t et tdt
// retour : =0 le point est externe, =1 le point est interne , 
//          = 2 le point est sur la frontière à la précision près 
// coor_locales : s'il est différent de NULL, est affecté des coordonnées locales calculées,
//                 uniquement précises si le point est interne
int PiPoCo::Interne_0(const Coordonnee& M,Coordonnee* )
  { cout << "\n fonction non implante, bool PiPoCo::Interne_0(Coordonnee& M)";
    Sortie(1);
    return Interne(TEMPS_0,M);    
   }; 
int PiPoCo::Interne_t(const Coordonnee& M,Coordonnee* ) 
  { cout << "\n fonction non implante, bool PiPoCo::Interne_t(Coordonnee& M)";
    Sortie(1);
    return Interne(TEMPS_t,M);    
   }; 
int PiPoCo::Interne_tdt(const Coordonnee& M,Coordonnee* ) 
  { cout << "\n fonction non implante, bool PiPoCo::Interne_tdt(Coordonnee& M)";
    Sortie(1);
    return Interne(TEMPS_tdt,M);    
   }; 
   
    // METHODES VIRTUELLES:
      // --------- calculs utils dans le cadre de la recherche du flambement linéaire		
// dans un premier temps uniquement virtuelles, ensuite se sera virtuelle pure pour éviter
// les oubli de définition ----> AMODIFIER !!!!!!!!    
		// Calcul de la matrice géométrique et initiale
ElemMeca::MatGeomInit PiPoCo::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & ) 
  { cout << "\n attention le calcul de la matrice g\'eom\'etrique  et de la matrice " 
         << " initiale   ne sont pas implant\'ees \n" ;
    this->Affiche(1);
    cout << "\nvoid PiPoCo::MatricesGeometrique_Et_Initiale (Mat_pleine &,Mat_pleine &)"
         << endl;
    Sortie(1);
    return MatGeomInit(NULL,NULL); // pour supprimer le warning
   };
    
// =========================== methodes protegees ===============================

// Calcul du residu local et de la raideur locale,
//  pour le schema implicite
void PiPoCo::Cal_implicitPiPoCo (DdlElement & tab_ddl,
                Tableau <TenseurBB *> & d_epsBB,Tableau < Tableau2 <TenseurBB *> > d2_epsBB_,
                Tableau <TenseurHH *> & d_sigHH,int nbintS,Vecteur& poidsS,int nbintH
                ,Vecteur& poidsH,const ParaAlgoControle & pa,bool cald_Dvirtuelle)
                
   {  int nbddl = tab_ddl.NbDdl();
      double jacobien; // def du jacobien
      Vecteur d_jacobien_tdt; 
      energie_totale.Inita(0.); // init de l'énergie totale sur l'élément
      TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
      if (ElemMeca::Bulk_visco_actif())  
          { if (sig_bulk_pourPiPoCo_HH == NULL)  {sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);}
            else if (sig_bulk_pourPiPoCo_HH->Dimension() != sigHH_t_1->Dimension())
               {delete sig_bulk_pourPiPoCo_HH; sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);};           
           };               
      
      int indice = 1;
      for (int nisur =1; nisur<= nbintS; nisur++)  // boucle sur les pt d'integ de surface
        for (int niepais =1; niepais<= nbintH; niepais++, indice++) // pt d'integ d'epaisseur
        { ((DeformationPP*) def)->ChangeNumIntegSH(nisur,niepais);
          ((DeformationPP*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
          PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
          TenseurHH & sigHH = *(ptIntegMeca.SigHH());
          TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
          EnergieMeca& energ =  tab_energ(indice);
          EnergieMeca& energ_t =  tab_energ_t(indice);
          double poids = poidsS(nisur) * poidsH(niepais)  ;
          CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
          if (loiTP != NULL) {sTP = tabSaveTP(indice);};         // au pt d'integ si elles existes 
          if ((loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
              || (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);

          // -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
          const Met_abstraite::Impli*  ex_pointe = NULL; // c'est pour construire ex
          try // ce qui permet de donner les numéros d'éléments et de pti
           { ex_pointe = &loiComp->Cal_implicit(tabSaveDon(indice), *def,tab_ddl,ptIntegMeca
                                                   ,d_epsBB,jacobien,d_jacobien_tdt,d_sigHH
                                                   ,pa,sTP,loiTP,dilatation,energ,energ_t
                                                   ,premier_calcul_meca_impli_expli);
           }
          catch (ErrNonConvergence_loiDeComportement excep)
           // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
           // d'une loi de comportement incrémentale
            { cout << "\n erreur de loi de comportement element= " << this->Num_elt()
                     << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw excep;
            }
          catch (ErrSortieFinale)
               // cas d'une direction voulue vers la sortie
               // on relance l'interuption pour le niveau supérieur
             { ErrSortieFinale toto;
               throw (toto);
             }
          catch ( ... )
            { cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
                     << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw (Err_inconnue_ElemMeca());
            };
          const Met_abstraite::Impli& ex=(*ex_pointe);

          // on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
          Tableau2 <TenseurBB *>& d2_epsBB =  d2_epsBB_(indice); // pour simplifier
          if (cald_Dvirtuelle) 
             def->Cal_var_def_virtuelle(false,d2_epsBB);
          // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
          if (Bulk_visco_actif()) 
             ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_tdt,(*sig_bulk_pourPiPoCo_HH) );

          if (jacobien <= 0.) // vérif du jacobien
            { cout << "\n **********  attention on a trouve un jacobien  negatif !! ********* " << endl; 
              cout << "\n on annulle sa contribution. Element nb: " <<  Num_elt() << " nb d'integ : " << indice;
              }
          else 
           { // cas avec toutes les variations
              {energie_totale += energ * (poids * jacobien);
               for (int j =1; j<= nbddl; j++)  // 1ere boucle sur les ddl
                {(*residu)(j) -= ((*(d_epsBB(j))) && sigHH ) * poids * jacobien;
                 for (int i =1; i<= nbddl; i++)  // 2ere boucle sur les ddl
                   (*raideur)(j,i) +=  ( ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien 
                                       +  (sigHH && (*(d2_epsBB(j,i)))) * jacobien
                                       +  (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i)
                                        ) * poids;
                 }
              }
           }//-- fin des différents cas
        if (Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
           { sigHH -= (*sig_bulk_pourPiPoCo_HH); }  // ceci pour la sauvegarde  ou autre utilisation       
 
        }; // fin boucle sur les points d'intégration d'épaisseur
     // liberation des tenseurs intermediaires
     LibereTenseur();   
   };

// Calcul uniquement du residu local a l'instant t ou pas suivant la variable atdt
// pour le schema explicit par exemple 
void PiPoCo::Cal_explicitPiPoCo (DdlElement & tab_ddl,Tableau <TenseurBB *> & d_epsBB,int nbintS,
                Vecteur& poidsS,int nbintH,Vecteur& poidsH,const ParaAlgoControle & ,bool atdt)

   {  double jacobien; // def du jacobien
      int nbddl = tab_ddl.NbDdl();
      int indice = 1;
      energie_totale.Inita(0.); // init de l'énergie totale sur l'élément
      TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
      if (ElemMeca::Bulk_visco_actif())  
          { if (sig_bulk_pourPiPoCo_HH == NULL)  {sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);}
            else if (sig_bulk_pourPiPoCo_HH->Dimension() != sigHH_t_1->Dimension())
               {delete sig_bulk_pourPiPoCo_HH; sig_bulk_pourPiPoCo_HH = NevezTenseurHH(*sigHH_t_1);};           
           };               
      for (int nisur =1; nisur<= nbintS; nisur++)  // boucle sur les pt d'integ de surface
        for (int niepais =1; niepais<= nbintH; niepais++, indice++) // pt d'integ d'epaisseur
        { ((DeformationPP*) def)->ChangeNumIntegSH(nisur,niepais);
          ((DeformationPP*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
          PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
          TenseurHH & sigHH = *(ptIntegMeca.SigHH());
          TenseurBB & DepsBB_ = *(ptIntegMeca.DepsBB());
          EnergieMeca& energ =  tab_energ(indice);
          EnergieMeca& energ_t =  tab_energ_t(indice);
          double poids = poidsS(nisur) * poidsH(niepais)  ;
          CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
          if (loiTP != NULL) {sTP = tabSaveTP(indice);};         // au pt d'integ si elles existes 
          if ((loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
              || (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);

          // -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
          try // ce qui permet de donner les numéros d'éléments et de pti
           {if (atdt)
               { const Met_abstraite::Expli_t_tdt& ex = loiComp->Cal_explicit_tdt 
                               (tabSaveDon(indice), *def,tab_ddl,ptIntegMeca,d_epsBB,jacobien
                               ,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
                 // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
                 if (Bulk_visco_actif()) 
                    ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_ ,(*sig_bulk_pourPiPoCo_HH));
               }
            else   
               { const Met_abstraite::Expli& ex = loiComp->Cal_explicit_t 
                               (tabSaveDon(indice), *def,tab_ddl,ptIntegMeca,d_epsBB,jacobien
                               ,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
                 // examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
                 if (Bulk_visco_actif()) 
                      ModifContrainteBulk(TEMPS_t,*ex.gijHH_t,sigHH,DepsBB_ ,(*sig_bulk_pourPiPoCo_HH));
               };
           }
          catch (ErrNonConvergence_loiDeComportement excep)
           // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
           // d'une loi de comportement incrémentale
            { cout << "\n erreur de loi de comportement element= " << this->Num_elt()
                   << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw excep;
            }
          catch (ErrSortieFinale)
               // cas d'une direction voulue vers la sortie
               // on relance l'interuption pour le niveau supérieur
             { ErrSortieFinale toto;
               throw (toto);
             }
          catch ( ... )
            { cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
                   << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw (Err_inconnue_ElemMeca());
            };

          if (jacobien <= 0.) // vérif du jacobien
            { cout << "\n **********  attention on a trouve un jacobien  negatif !! ********* " << endl; 
              cout << "\n on annulle sa contribution. Element nb: " <<  Num_elt() << " nb d'integ : " << indice;
              }                    
          else
           { energie_totale += energ * (poids * jacobien);
             for (int j =1; j<= nbddl; j++)
               (*residu)(j) -=  (sigHH && (*(d_epsBB(j)))) * jacobien * poids;
            } 
        if (Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
           { sigHH -= (*sig_bulk_pourPiPoCo_HH); }  // ceci pour la sauvegarde  ou autre utilisation       
 
        }; // fin boucle sur les points d'intégration
     // liberation des tenseurs intermediaires
     LibereTenseur();   
   };

// Calcul de la matrice géométrique et de la matrice initiale
// cette fonction est éventuellement appelée par les classes dérivées 
// ddl represente les degres de liberte specifiques a l'element
// epsBB = deformation, sigHH = contrainte, d_epsbb = variation des def
// nbint = nb maxi de pt d'integration , poids = poids d'integration
// S pour surface et H pour épaisseur -> nbint et poids  
void PiPoCo::Cal_matGeom_InitPiPoCo (Mat_pleine & matGeom, Mat_pleine & matInit,
                 DdlElement & tab_ddl,Tableau <TenseurBB *>& d_epsBB, Tableau < Tableau2 <TenseurBB *> > d2_epsBB_,
                Tableau <TenseurHH *> & d_sigHH,int nbintS
                ,Vecteur& poidsS,int nbintH,Vecteur& poidsH,const ParaAlgoControle & pa,bool cald_Dvirtuelle)
   {  // le début du programme est semblable au calcul implicit
      // seule la constitution des matrices finales est différente
      int nbddl = tab_ddl.NbDdl();
      double jacobien; // def du jacobien
      Vecteur d_jacobien_tdt; 

      int indice = 1;
      for (int nisur =1; nisur<= nbintS; nisur++)  // boucle sur les pt d'integ de surface
        for (int niepais =1; niepais<= nbintH; niepais++, indice++) // pt d'integ d'epaisseur
        { ((DeformationPP*) def)->ChangeNumIntegSH(nisur,niepais);
          ((DeformationPP*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
          PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
          TenseurHH & sigHH = *(ptIntegMeca.SigHH());
          TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
          EnergieMeca& energ =  tab_energ(indice);
          EnergieMeca& energ_t =  tab_energ_t(indice);
          double poids = poidsS(nisur) * poidsH(niepais)  ;
          CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
          if (loiTP != NULL) {sTP = tabSaveTP(indice);};         // au pt d'integ si elles existes 
          if ((loiComp->Id_comport()==LOI_VIA_UMAT)   // cas de l'utilisation d'une loi umat
             || (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
            ((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);

          // -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
          try // ce qui permet de donner les numéros d'éléments et de pti
           {loiComp->Cal_flamb_lin(tabSaveDon(indice), *def,tab_ddl,ptIntegMeca,d_epsBB
                             ,jacobien,d_jacobien_tdt, d_sigHH,pa,sTP,loiTP,dilatation
                             ,energ,energ_t,premier_calcul_meca_impli_expli);
           }
          catch (ErrNonConvergence_loiDeComportement excep)
           // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
           // d'une loi de comportement incrémentale
            { cout << "\n erreur de loi de comportement element= " << this->Num_elt()
                   << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw excep;
            }
          catch (ErrSortieFinale)
               // cas d'une direction voulue vers la sortie
               // on relance l'interuption pour le niveau supérieur
             { ErrSortieFinale toto;
               throw (toto);
             }
          catch ( ... )
            { cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
                   << " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
              throw (Err_inconnue_ElemMeca());
            };
         
          // on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
          Tableau2 <TenseurBB *>& d2_epsBB =  d2_epsBB_(indice); // pour simplifier
          if (cald_Dvirtuelle) 
             def->Cal_var_def_virtuelle(false,d2_epsBB);
          for (int j =1; j<= nbddl; j++)  // 1ere boucle sur les ddl
           for (int i =1; i<= nbddl; i++)  // 2ere boucle sur les ddl
              {matInit(j,i) += ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien * poids;
               matGeom(j,i) += (sigHH && (*(d2_epsBB(j,i)))) * jacobien * poids; 
              }
        }
     // liberation des tenseurs intermediaires
     LibereTenseur();   
   };