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

/************************************************************************
 *     DATE:        28/07/2006                                          *
 *                                                                $     *
 *     AUTEUR:      G RIO   (mailto:gerardrio56@free.fr)                *
 *                                                                $     *
 *     PROJET:      Herezh++                                            *
 *                                                                $     *
 ************************************************************************
 *     BUT:    Algorithme de calcul dynamique explicite, pour de la     *
 *             mecanique du solide déformable en coordonnees materielles*
 *             entrainees.                                              *
 *            Correspond à l'implantation de l'algo de Bonelli et Bursi.*
 *            Le modèle est de type galerkin discontinu en tempsP1-P1,  *
 *            pour la discrétisation. L'algo est de type prédiction     *
 *            suivi de 1 ou plusieurs cycle de correction.              *
 *            Il est explicite, dans le sens où on n'utilise pas de     *
 *            comportement tangent du matériau et que l'on contrôle     *
 *            exactement le nombre d'itération du processus.            *
 *                                                                      *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *                                                                      *
 *     VERIFICATION:                                                    *
 *                                                                      *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *     !        !            !                                    !     *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *
 *     MODIFICATIONS:                                                   *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *                                                                $     *
 ************************************************************************/
#ifndef AGORI_BONELLI_H
#define AGORI_BONELLI_H


#include "Algori.h"
#include "Assemblage.h"
#include "GeomSeg.h"

/// @addtogroup Les_algorithmes_de_resolutions_globales
///  @{
///

///     BUT:    Algorithme de calcul dynamique explicite, pour de la
///             mecanique du solide déformable en coordonnees materielles
///             entrainees.
///            Correspond à l'implantation de l'algo de Bonelli et Bursi.
///            Le modèle est de type galerkin discontinu en tempsP1-P1,
///            pour la discrétisation. L'algo est de type prédiction
///            suivi de 1 ou plusieurs cycle de correction.
///            Il est explicite, dans le sens où on n'utilise pas de
///            comportement tangent du matériau et que l'on contrôle
///            exactement le nombre d'itération du processus.

class AlgoBonelli : public Algori
{
  public :
    // CONSTRUCTEURS :
    AlgoBonelli () ; // par defaut
    
    // constructeur en fonction du type de calcul
    // du sous type (pour les erreurs, remaillage etc...)
    // il y a ici lecture des parametres attaches au type
    AlgoBonelli (const bool avec_typeDeCal
           ,const list <EnumSousTypeCalcul>& soustype
           ,const list <bool>& avec_soustypeDeCal
           ,UtilLecture& entreePrinc);
    // constructeur de copie
    AlgoBonelli (const AlgoBonelli& algo);
    
    // constructeur de copie à partie d'une instance indifférenciée
    Algori * New_idem(const Algori* algo) const
     {// on vérifie qu'il s'agit bien d'une instance
      if (algo->TypeDeCalcul() != DYNA_EXP_BONELLI)
        { cout << "\n *** erreur lors de la creation par copie d'un algo DYNA_EXP_BONELLI "
               << " l'algo passe en parametre est en fait : " << Nom_TypeCalcul(algo->TypeDeCalcul())
               << " arret !! " << flush;
          Sortie(1);
          return NULL; 
        }
      else
        { AlgoBonelli* inter = (AlgoBonelli*) algo;
          return ((Algori *) new AlgoBonelli(*inter));
        };
     };

    // DESTRUCTEUR :
     ~AlgoBonelli () ;  
    
    // METHODES PUBLIQUES :
    // execution de l'algorithme explicite dans le cas  dynamique sans contact
    void Execution(ParaGlob * ,LesMaillages *,LesReferences*,LesCourbes1D* ,LesFonctions_nD*
                  ,VariablesExporter* varExpor,LesLoisDeComp*
                  ,DiversStockage*,Charge*,LesCondLim*,LesContacts*,Resultats* );
        
    //------- décomposition en 3 du calcul d'équilibre -------------
    // a priori   : InitAlgorithme  et FinCalcul ne s'appellent qu'une fois,
    // par contre : CalEquilibre peut s'appeler plusieurs fois, le résultat sera différent si entre deux calcul
    //              certaines variables ont-été changés
    
    // initialisation
    void InitAlgorithme(ParaGlob * ,LesMaillages *,LesReferences*,LesCourbes1D*
                        ,LesFonctions_nD* ,VariablesExporter*,LesLoisDeComp*
                        ,DiversStockage*,Charge*,LesCondLim*,LesContacts*,Resultats* );
    // mise à jour
    void MiseAJourAlgo(ParaGlob * ,LesMaillages *,LesReferences*,LesCourbes1D*
                       ,LesFonctions_nD* ,VariablesExporter* ,LesLoisDeComp*
                       ,DiversStockage*,Charge*,LesCondLim*,LesContacts*,Resultats* );
    // calcul de l'équilibre
    // si tb_combiner est non null -> un tableau de 2 fonctions
    //  - la première fct dit si on doit valider ou non le calcul à convergence ok,
    //  - la seconde dit si on doit sortir de la boucle ou non à convergence ok
    //
    // si la validation est effectuée, la sauvegarde pour le post-traitement est également effectuée
    //                                 en fonction de la demande de sauvegard,
    // sinon pas de sauvegarde pour le post-traitement à moins que l'on a demandé un mode debug
    //   qui lui fonctionne indépendamment
    void CalEquilibre(ParaGlob * ,LesMaillages *,LesReferences*,LesCourbes1D*
                      ,LesFonctions_nD* ,VariablesExporter*,LesLoisDeComp*
                      ,DiversStockage*,Charge*,LesCondLim*,LesContacts*,Resultats*
                      ,Tableau < Fonction_nD* > * tb_combiner);
    // dernière passe
    void FinCalcul(ParaGlob * ,LesMaillages *,LesReferences*,LesCourbes1D*
                   ,LesFonctions_nD* ,VariablesExporter*,LesLoisDeComp*
                   ,DiversStockage*,Charge*,LesCondLim*,LesContacts*,Resultats* );

    // sortie du schemaXML: en fonction de enu
    void SchemaXML_Algori(ofstream& sort,const Enum_IO_XML enu) const ;
 
  protected :  
    // VARIABLES PROTEGEES :
    // paramètres de l'algorithme
    double a_a,b_b;
    int k_max;
    double omega_b; // omega de bifurcation < omega critique
    double rho_b; // rayon spectral à la bifurcation
    GeomSeg* pt_seg_temps; // un élement segment pour l'intégration temporelle avec gauss
    int nb_pt_int_t; // le nombre de points d'intégration temporelle 
    
    // liste de variables de travail déclarées ici pour éviter le passage de paramètre entre les 
    // méthodes internes à la classe
    double delta_t; // pas de temps instantané
    double deltat2,unsurdeltat,deltatSurDeux; // variable intermédiaire temps instantané
    double delta_t_total; // pas de temps total 
    double deltat2_total,unsurdeltat_total,deltatSurDeux_total; // variable intermédiaire temps totales
    
    
    // pour clarifier la lecture des différentes phases pour la mise en place des CL sur ddl
    enum enuTypePhase {PREDICTION_bonelli, INTEGRATION_bonelli
                       , CORRECTION_bonelli, FIN_DELTA_T_bonelli }; 

    // ---------------- variables de transferts internes  ------------------ 
    //  === pointeurs d'instance et classe particulières
    LesMaillages * lesMail_;  
    LesReferences* lesRef_;
    LesCourbes1D* lesCourbes1D_;
    LesFonctions_nD* lesFonctionsnD_ ;
    Charge* charge_;
    LesCondLim* lesCondLim_;
    LesContacts* lesContacts_;
    Assemblage * Ass1, * Ass2, * Ass3; // pointeurs d'assemblages
    
    //  === variables scalaires
    double maxPuissExt;   // maxi de la puissance des efforts externes
    double maxPuissInt;   // maxi de la puissance des efforts internes 
    double maxReaction;   // maxi des reactions 
    int inReaction;       // pointeur d'assemblage pour le maxi de reaction
    int inSol;            // pointeur d'assemblage du maxi de variation de ddl      
    double maxDeltaDdl;   // maxi de variation de ddl
    int cas_combi_ddl;   // def combinaison des ddl 
    int icas;        // idem cas_combi_ddl mais pour lesCondlim
    bool erreurSecondMembre; // pour la gestion des erreurs de calcul au second membre
    bool prepa_avec_remont; // comme son nom l'indique
    bool brestart;   // booleen qui indique si l'on est en restart ou pas
    OrdreVisu::EnumTypeIncre type_incre; // pour la visualisation au fil du calcul
    
    //  === vecteurs
    Vecteur q_0,p_0; // position et quantité de mouvement  à ti^{-}
    Vecteur q_1,p_1; // position et quantité de mouvement  à ti^{+}
    Vecteur q_2,p_2; // position et quantité de mouvement  à ti+1^{-}
    
    Vecteur pPoint_i; // quantité de mouvement à ti^{-}
    Vecteur P_i_1,P_i_2,P_q1_k,P_q2_k; // inter
    Vecteur r_1_k,r_2_k; // les résidus
    Vecteur vec_travail1,vec_travail2; // vecteurs de travail
    Vecteur dernier_phi; // fonction d'interpolation temporelle au dernier pt d'integ en temps
    
    
    Vecteur vglobin; // puissance interne : pour ddl accélération
    Vecteur vglobex; // puissance externe 
    Vecteur vglobaal; // puissance totale qui ecrase vglobin
    Vecteur vcontact; // puissance des forces de contact
    Vecteur vitesse_tplus;  // vitesses
    Vecteur X_Bl,V_Bl,G_Bl; // stockage transitoirement des X V GAMMA <-> CL
    Vecteur forces_vis_num; // forces visqueuses d'origines numériques

    //  === les listes
    list <LesCondLim::Gene_asso> li_gene_asso; // tableaux d'indices généraux des ddl bloqués
    //  === les tableaux
    Tableau <Nb_assemb> t_assemb; // tableau globalisant les numéros d'assemblage de X V gamma
    Tableau <Enum_ddl> tenuXVG; // les enum des inconnues
    //  === les matrices
    Mat_abstraite* mat_masse,* mat_masse_sauve; // choix de la matrice de masse
    Mat_abstraite* mat_C_pt; // matrice visqueuse numérique
    
    
    // METHODES PROTEGEES :

    // lecture des paramètres du calcul
    void lecture_Parametres(UtilLecture& entreePrinc);
    // écriture des paramètres dans la base info
    // = 1 : on écrit tout
    // = 2 : on écrot uniquement les données variables (supposées comme telles)
    void Ecrit_Base_info_Parametre(ostream& sort,const int& cas);
    // lecture des paramètres dans la base info
    // = 1 : on récupère tout
    // = 2 : on récupère uniquement les données variables (supposées comme telles)
    // choix = true  : fonctionnememt normal
    // choix = false : la méthode ne doit pas lire mais initialiser les données à leurs valeurs par défaut
    //                 car la lecture est impossible
    void Lecture_Base_info_Parametre(istream& ent,const int& cas,bool choix);
    // création d'un fichier de commande: cas des paramètres spécifiques
    void Info_commande_parametres(UtilLecture& entreePrinc);
    
    // gestion et vérification du pas de temps et modif en conséquence si nécessaire
    // cas = 0: premier passage à blanc, delta t = 0
    // cas = 1: initialisation du pas de temps et vérif / au pas de temps critique
    //          ceci pour le temps t=0
    // cas = 2: initialisation du pas de temps et vérif / au pas de temps critique
    //          ceci pour le temps t. et on divise par nbstep
    //          ceci pour garantir que l'on fait le calcul avec 1 step
    void Gestion_pas_de_temps(LesMaillages * lesMail,int cas,int nbstep);
    
    // modification transitoire du pas de temps et modif en conséquence si nécessaire
    // utilisée pour les différentes intégrales temporelles
    // delta_tau : nouveau pas de temps transitoire imposé
    void Modif_transi_pas_de_temps(double  delta_tau);
    
    // calcul des ddl avec prise en comptes des conditions limites suivant une technique 
    // s'appuyant sur une différence finie décentrée à droite
    // utilise les vecteurs globaux : q_2, q_1, q_0
    // --> mise à jour des ddl au niveau globaux sur les X, V, gamma puis au niveau des noeuds
    // --> ou mise à jour des composantes adoc des vecteurs globaux q_2, q_1, q_0, selon "phasage"
    // phii : interpolation en fonction du temps, ne sert que pour phasage = INTEGRATION_bonelli
    // phasage =
    //  PREDICTION_bonelli: phase de prédiction, on n'utilise "pas" l'interpolation en temps (phii)
    //                      --> calcul de q_1 et q_2 avec les conditions limites
    //  INTEGRATION_bonelli: phase d'intégration, "on utilise l'interpolation en temps" 
    //                      --> calcul de X_tdt et vitesse_tdt
    //  CORRECTION_bonelli: phase de correction, on n'utilise "pas" l'interpolation en temps (phii)
    //                      --> calcul de q_1 et q_2 avec les conditions limites
    //  FIN_DELTA_T_bonelli: phase final, on n'utilise "pas" l'interpolation en temps (phii)
    //                      --> calcul de X_tdt et vitesse_tdt pour le temps final
    void AvanceDDL_avec_CL(const Vecteur & phii,enuTypePhase phasage);
    
    // calcul  des différences énergie en jeux dans le cas Bonelli
    // affichage éventuelle des ces énergies et du bilan
    // V         : ddl de vitesse
    // coef_mass : en fait il faut utiliser 1/coef_mass * mat_mass pour la matrice masse, ceci due à la spécificité de l'algo
    // les différentes énergie et efforts généralisées sont des variables internes à l'algo
    // en sortie: énergie cinétique : E_cin_tdt,  énergie interne : E_int_tdt, énergie externe = E_ext_tdt, bilan : bilan_E
    //            idem avec les puissances
    //            énergie visqueuse due aux forces visqueuses numériques
    // icharge : compteur d'increment
    void CalEnergieAffichageBonelli(const double& coef_mass,int icharge);

    //---- gestion des commndes interactives --------------
    // écoute et prise en compte d'une commande interactive
    // ramène true tant qu'il y a des commandes en cours
    bool ActionInteractiveAlgo();
    
 };
 /// @}  // end of group

#endif