// 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-2021 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:        12/02/2006                                          *
 *                                                                $     *
 *     AUTEUR:      G RIO   (mailto:gerardrio56@free.fr)                *
 *                                                                $     *
 *     PROJET:      Herezh++                                            *
 *                                                                $     *
 ************************************************************************
 *     BUT:   Algorithmes de base pour la résolution d'équations        *
 *            différentielles.                                          *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *
 *                                                                $     *
 ************************************************************************/
#ifndef ALGO_EDP_H
#define ALGO_EDP_H

#include "Vecteur.h"
#include "Mat_abstraite.h"
#include "Racine.h"

/** @defgroup Les_classes_algo
*
*     BUT:   Algorithmes de base pour la résolution d'équations
*            différentielles, recherche de zéros, intégration.
*
*
* \author    Gérard Rio
* \version   1.0
* \date      12/02/2006
* \brief       Algorithmes de base pour la résolution d'équations différentielles, recherche de zéros, intégration.
*
*/

/// @addtogroup Les_classes_algo
///  @{
///

///     BUT:   Algorithmes de base pour la résolution d'équations
///            différentielles.
class Algo_edp
{
  public :
    // CONSTRUCTEURS :
    // constructeur par défaut
    Algo_edp();
    // constructeur de copie 
    Algo_edp(const Algo_edp& a);    
    // DESTRUCTEUR :
    ~Algo_edp();
    
    // METHODES PUBLIQUES :
    

    /// résolution par runge kutta imbriqué --> estimation d'erreur
    /// il y a 3 calcul de la fonction dérivée
    /// 2 et 3 ième ordre imbriquée
    /// calcul de la solution à t+dt en fonction de la solution à t,
    /// en entrée:
    ///  *Ptder_fonc  : pointeur de la fonction de calcul  de la dérivées
    ///           en entrée: t et f : temps et la valeur de la fonction a t
    ///           en sortie: df     : dérivée de la fonction a t
    ///                      erreur : si diff de 0, indique qu'il y a eu une erreur
    ///  *Ptverif_fonc : pointeur de fonction, pour tester val_finale après la prédiction explicite finale
    ///
    ///  val_initiale : valeur des fonctions en t0
    ///  der_initiale  : valeur de la dérivée pour t0
    ///  t0           : temps initiale
    ///  deltat       : incrément de temps demandé
    ///  en sortie :
    ///  val_finale    : valeur des fonctions à tdt
    ///  estime_erreur : estimation d'erreur pour chaque fonction
    ///   --> si estime_erreur(i) est >= à ConstMath::tresgrand, la valeur finale = la valeur initiale (c-a-d pas de calcul valide)
    template <class T> void  Runge_Kutta_step23(T& instance
                         ,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
                         ,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
                         ,const Vecteur& val_initiale,const Vecteur& der_initiale
                         ,const double& t0,const double& deltat
                         ,Vecteur& val_finale,Vecteur& estime_erreur) ;


    /// résolution par runge kutta imbriqué --> estimation d'erreur
    /// il y a 4 calcul de la fonction dérivée
    /// 3 et 4 ième ordre imbriquée (Fehlberg)
    /// calcul de la solution à t+dt en fonction de la solution à t,
    /// en entrée:
    ///  *Ptder_fonc  : pointeur de la fonction de calcul  de la dérivées
    ///           en entrée: t et f : temps et la valeur de la fonction a t
    ///           en sortie: df     : dérivée de la fonction a t
    ///                      erreur : si diff de 0, indique qu'il y a eu une erreur
    ///  *Ptverif_fonc : pointeur de fonction, pour tester val_finale après la prédiction explicite finale
    ///
    ///  val_initiale : valeur des fonctions en t0
    ///  der_initiale  : valeur de la dérivée pour t0
    ///  t0           : temps initiale
    ///  deltat       : incrément de temps demandé
    ///  en sortie :
    ///  val_finale    : valeur des fonctions à tdt
    ///  estime_erreur : estimation d'erreur pour chaque fonction
    ///   --> si estime_erreur(i) est >= à ConstMath::tresgrand, la valeur finale = la valeur initiale (c-a-d pas de calcul valide)
    template <class T> void  Runge_Kutta_step34(T& instance
                         ,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
                         ,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
                         ,const Vecteur& val_initiale,const Vecteur& der_initiale
                         ,const double& t0,const double& deltat
                         ,Vecteur& val_finale,Vecteur& estime_erreur) ;

    /// résolution par runge kutta imbriqué --> estimation d'erreur
    /// il y a 5 calcul de la fonction dérivée
    /// l'algorithme s'appuis sur la présentation de numerical recipes
    /// fifth-order Cash-Karp Runge-Kutta
    /// calcul de la solution à t+dt en fonction de la solution à t,
    /// en entrée:
    ///  *Ptder_fonc   : pointeur de la fonction de calcul  de la dérivées
    ///           en entrée: t et f : temps et la valeur de la fonction a t
    ///           en sortie: df     : dérivée de la fonction a t
    ///                      erreur : si diff de 0, indique qu'il y a eu une erreur
    ///  *Ptverif_fonc : pointeur de fonction, pour tester val_finale après la prédiction explicite finale
    ///
    ///  val_initiale : valeur des fonctions en t0
    ///  der_initiale  : valeur de la dérivée pour t0
    ///  t0           : temps initiale
    ///  deltat       : incrément de temps demandé
    ///  en sortie :
    ///  val_finale    : valeur des fonctions à tdt
    ///  estime_erreur : estimation d'erreur pour chaque fonction
    ///   --> si estime_erreur(i) est >= à ConstMath::tresgrand, la valeur finale = la valeur initiale (c-a-d pas de calcul valide)
    template <class T> void  Runge_Kutta_step45(T& instance
                         ,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
                         ,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
                         ,const Vecteur& val_initiale,const Vecteur& der_initiale
                         ,const double& t0,const double& deltat
                         ,Vecteur& val_finale,Vecteur& estime_erreur) ;
 
    
    /// un programme de pilotage de l'intégration d'un pas (on s'inspire du programme libre rkf45)
    /// on distingue une erreur globale et une erreur relative
    /// l'erreur réelle utilisée sur chaque step est: erreurRelative * |f| + erreurAbsolue
    ///   ==== en entrée: =====
    ///   cas_kutta      : donne le type de routine kutta imbriqué a employer:
    ///                    =3 pour kutta23, =4 pour kutta34, =5 pour kutta45
    ///   tdebut et tfin : temps de début et de fin du calcul, demandé
    ///   val_initiale   : valeur initiale de la fonction (donc à tdebut)
    ///   der_initiale   : dérivée initiale de la fonction (donc à tdebut)
    ///   erreurAbsolue  : erreur absolue demandée
    ///   erreurRelative : erreur relative (à la fonction) demandée
    ///   Ptder_fonc     : pointeur de fonction
    ///           en entrée: t et f : temps et la valeur de la fonction a t
    ///           en sortie: df     : dérivée de la fonction a t
    ///           en entrée: t et f : temps et la valeur de la fonction a t
    ///           en sortie: df     : dérivée de la fonction a t
    ///                      erreur : si diff de 0, indique qu'il y a eu une erreur
    ///   Ptverif_fonc   : pointeur de fonction, pour tester l'intégrité du résultat
    ///           cela signifie que val_finale est testée avec cette fonction, systématiquement
    ///           - à la fin de chaque prédiction des kuttas imbriqués
    ///           - au retour du pilotage
    ///   ---- en sortie: -----
    /// dernierTemps    : temps final utilisé = tfin si calcul ok, sinon = le dernier temps calculé
    /// val_finale      : valeur de la fonction à "dernierTemps"
    /// der_finale      :  dérivée finale
    /// dernierdeltat   : le dernier deltat utilisé
    /// nombreAppelF    : nombre d'appel de la fonction réellement utilisé dans le programme
    /// erreurRelative  : si retour = 3 ==> erreur relative (à la fonction) minimum possible,
    /// nb_step         : nombre de step utilisé pour le calcul (non compté les steps avec trop d'erreur)
    /// erreur_maxi_global : une approximation de l'erreur maxi global cumulant les erreurs sur tous les steps
    ///   === en retour un entier qui donne le résultat du calcul : ====
    /// =2  : tout est ok : seul cas où toutes les valeurs de retour sont valides
    /// =3  : erreur: la précision relative demandée est trop petite, en retour la précision mini possible
    /// =4  : erreur: on a dépassé le nombre maxi d'appel fixé, d'où a peu près à (nombreAppelF/6) step de calcul.
    /// =6  : erreur: l'intégration pas possible, due aux précisions demandées, on doit augmenter ces précisions
    ///               souvent du à une solution qui varie très rapidemment --> problème potentiel
    /// =8  : erreur: cas_kutta ne correspond pas à une routine de base existante
    /// = 0 : erreur inconnue

    template <class T> 
    int  Pilotage_kutta(int cas_kutta, T& instance
                        ,Vecteur& (T::*Ptder_fonc) (const double & t,const Vecteur& f,Vecteur& df,int& erreur)
                        ,void (T::*Ptverif_fonc) (const double & t,const Vecteur& f,int& erreur_final)
                        ,const Vecteur& val_initiale,const Vecteur& der_initiale
                        ,const double& tdebut,const double& tfin
                        ,const double& erreurAbsolue, double& erreurRelative
                        ,Vecteur& val_finale,Vecteur& der_finale
                        ,double& dernierTemps, double& dernierdeltat,int& nombreAppelF,int& nb_step
                        ,double& erreur_maxi_global);
    /// explication erreur
    /// affiche à l'écran l'explication de l'erreur de Pilotage_kutta
    void Affiche_Explication_erreur_RG(int type_erreur);
 
    /// --méthodes pour modifier les différents paramètres
        /// nombre maxi d'appel de fonction permis
        /// dans le cas du runge 4-5 -> 6 appels par step (à la louche)
    void Modif_nbMaxiAppel (int nb) {nbMaxiAppel = nb;};
        /// initialisation de tous les paramètres à leurs valeurs par défaut
    void Init_param_val_defaut();
        /// modifie le niveau d'affichage par exemple pour le débug
    void Change_niveau_affichage(int niveau) {permet_affichage=niveau;};
        
  private :  
    // VARIABLES PROTEGEES :
        
    // ----- controle de la sortie des informations
    int permet_affichage; // pour permettre un affichage spécifique dans les méthodes, pour les erreurs et warning
    int nbMaxiAppel; // nombre maxi d'Appel autorisé de la fonction
    // ---- variables intermédiaires utilisées par la méthode Runge_Kutta_step23
    static double AB1,AB2,AB3,A2,aa2,aa3
                 ,bb21,bb31,bb32
                 ,dA1,dA2,dA3;
    // ---- variables intermédiaires utilisées par la méthode Runge_Kutta_step34
    static double A_B1,A_B3,A_B4,A_B5,A_1,A_3,A_4
                 ,a_a2,a_a3,a_a4,a_a5
                 ,b_b21,b_b31,b_b32,b_b41,b_b42,b_b43,b_b51,b_b53,b_b54
                 ,d_A1,d_A3,d_A4,d_A5;
    // --- variables intermédiaires utilisées par la méthode Runge_Kutta_step45
    Vecteur k2,k3,k4,k5,k6;
//    Vecteur f_inter;
    static double a2,a3,a4,a5,a6    // les ai de Cash-karp paramètres
                  ,b21,b31,b32,b41,b42,b43               //
               ,b51,b52,b53,b54,b61,b62,b63,b64,b65   //  les bij 
               ,c1,c3,c4,c6;     // les ci
               // c2 =0, c5=0
    static double ce1,ce3,ce4,ce5,ce6; // ce2=0;  : les c-c_étoile_i
    // --- variables intermédiaires utilisées par la méthode Pilotage_kutta
    Vecteur val_inter; // valeur de la fonction, intermédiaire
    Vecteur der_inter; // dérivée de la fonction, intermédiaire
    Vecteur estime_erreur; // estimation d'erreur calculée
       // puis des sauvegardes pour pouvoir revenir en arrière
    Vecteur val_prec_inter; // sauve avant appel: valeur de la fonction, intermédiaire
    Vecteur der_prec_inter; // sauve avant appel: dérivée de la fonction, intermédiaire
    Vecteur estime_prec_erreur; // sauve avant appel: estimation d'erreur calculée

    double coef; // coeff de sécurité préconnisé pour les variations de temps
    static double remin; // pour calculer l'erreur mini acceptable dans le runge
    
    // METHODES PROTEGEES :
    // impression du temps
    void timestamp ();
    // --  méthode test ---
    // on défini une fonction a intégrer
    Vecteur& FoncDeriv(const double & t,const Vecteur& f,Vecteur& df,int& erreur);
    // on défini une fonction de vérification d'intégrité
    void FoncVerif_fonc(const double & t,const Vecteur& f,int& erreur);
    int nb_cas_test; // nb pour choisir entre différents cas test
    // la méthode qui fait les tests 
    void Test_template();

 };
 
  // pour faire de l'inline: nécessaire avec les templates
  // on n'inclut que les méthodes templates
  #include "Algo_edp_2.cc"
  #define  Algo_edp_deja_inclus
 /// @}  // end of group

#endif  
