Herezh_dev/Util/Algo_edp.cc

// 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 "Algo_edp.h"

// contiend que les méthodes qui ne sont pas des templates

  #include "MathUtil.h"
  #include "ConstMath.h"

// ----- def des variables statiques ---------
    // ---- variables intermédiaires utilisées par la méthode Runge_Kutta_step23
double Algo_edp::AB1=1./6.,Algo_edp::AB2=2./3.,Algo_edp::AB3=1./6.
      ,Algo_edp::A2=1.,Algo_edp::aa2=0.5,Algo_edp::aa3=1.
      ,Algo_edp::bb21=0.5,Algo_edp::bb31=-1.,Algo_edp::bb32=2.
      ,Algo_edp::dA1=1./6,Algo_edp::dA2=-1./3.,Algo_edp::dA3=1./6.;
    // ---- variables intermédiaires utilisées par la méthode Runge_Kutta_step34
double Algo_edp::A_B1=229./1470.,Algo_edp::A_B3=1125./1813.,Algo_edp::A_B4=13718./81585.
      ,Algo_edp::A_B5=1./18.
      ,Algo_edp::A_1=79./490.,Algo_edp::A_3=2175./3626.,Algo_edp::A_4=2166./9065.
      ,Algo_edp::a_a2=2./7.,Algo_edp::a_a3=7./15.,Algo_edp::a_a4=35./38.,Algo_edp::a_a5=1.
      ,Algo_edp::b_b21=2./7.
      ,Algo_edp::b_b31=77./900.,Algo_edp::b_b32=343./900.
      ,Algo_edp::b_b41=805./1444.,Algo_edp::b_b42=-77175./54872.,Algo_edp::b_b43=97125./54872.
      ,Algo_edp::b_b51=79./490.,Algo_edp::b_b53=2175./3626.,Algo_edp::b_b54=2166./9065.
      ,Algo_edp::d_A1=0.,Algo_edp::d_A3=0.,Algo_edp::d_A4=0.,Algo_edp::d_A5=1./18.;
    // --- variables intermédiaires utilisées par la méthode Runge_Kutta_step45
double // d'abord les ai
       Algo_edp::a2=0.2,Algo_edp::a3=0.3,Algo_edp::a4=0.6,Algo_edp::a5=1.0,Algo_edp::a6=7./8.;
       // puis les bij
double       
Algo_edp::b21=0.2
,Algo_edp::b31=3.0/40.0,Algo_edp::b32=9.0/40.0
,Algo_edp::b41=0.3,Algo_edp::b42 = -0.9,Algo_edp::b43=1.2
,Algo_edp::b51 = -11.0/54.0, Algo_edp::b52=2.5,Algo_edp::b53 = -70.0/27.0,Algo_edp::b54=35.0/27.0
,Algo_edp::b61=1631.0/55296.0,Algo_edp::b62=175.0/512.0,Algo_edp::b63=575.0/13824.0
		      ,Algo_edp::b64=44275.0/110592.0,Algo_edp::b65=253.0/4096.0;
	   // puis les ci	      
double
Algo_edp::c1=37.0/378.0,Algo_edp::c3=250.0/621.0,Algo_edp::c4=125.0/594.0,Algo_edp::c6=512.0/1771.0;
       // et  c5-c_étoile_5
double Algo_edp::ce5 = -277.00/14336.0;
// initialisés dans le premier appel du construteur
double Algo_edp::ce1=0.,Algo_edp::ce3=0.,Algo_edp::ce4=0.,Algo_edp::ce6=0.;
		      
// pour la méthode Pilotage_kutta (cf.rkf45)
double Algo_edp::remin = 1.0E-12;
	      
		      
// -------- fin des variables statiques --------

    // CONSTRUCTEURS :
    // constructeur par défaut
Algo_edp::Algo_edp() :
  k2(1),k3(1),k4(1),k5(1),k6(1)//,f_inter(1)
 ,val_inter(1),der_inter(1),estime_erreur(1),coef(0.9)
 ,nb_cas_test(1),permet_affichage(0)
    { if (ce1 == 0.) // au premier appel on initialise
       {ce1=c1-2825.0/27648.0;
        ce3=c3-18575.0/48384.0;
		      ce4=c4-13525.0/55296.0;
		      ce6=c6-0.25;
       };
      if (d_A1==0.) 
      	{d_A1=A_B1-A_1;
        d_A3=A_B3-A_3;
        d_A4=A_B4-A_4;
        d_A5=A_B5;
      	};
      // initialisation des paramètres par défaut
      Init_param_val_defaut(); 
      // --- test des routines ---- (en fonctionnement normal: commenté)
//      Test_template(); 
    };
    // constructeur de copie 
Algo_edp::Algo_edp(const Algo_edp& a) :
  k2(a.k2),k3(a.k3),k4(a.k4),k5(a.k5),k6(a.k6)//,f_inter(a.f_inter)
  ,val_inter(a.val_inter),der_inter(a.der_inter),estime_erreur(a.estime_erreur),coef(a.coef)
  ,estime_prec_erreur(a.estime_prec_erreur)
  ,nb_cas_test(1),permet_affichage(0)
    { if (ce1 == 0.) // au premier appel on initialise
       {ce1=c1-2825.0/27648.0;
        ce3=c3-18575.0/48384.0;
		      ce4=c4-13525.0/55296.0;
		      ce6=c6-0.25;
       };
      // initialisation des paramètres par défaut
      Init_param_val_defaut(); 
    };
       
    // DESTRUCTEUR :
Algo_edp::~Algo_edp()
   {}
    
    // METHODES PUBLIQUES :
        // initialisation de tous les paramètres à leurs valeurs par défaut
void Algo_edp::Init_param_val_defaut()        
 {  nbMaxiAppel = 3000; }; 

// explication erreur
// affiche à l'écran l'explication de l'erreur de Pilotage_kutta
void Algo_edp::Affiche_Explication_erreur_RG(int type_erreur)
  { switch (type_erreur)
     {case 2  : cout << " calcul correct "; break;
      case 3  : cout << " la precision relative demandee est trop petite "; break;
      case 4  : cout << " nombre maxi d'appels de la fonction, atteints "; break;
      case 6  : cout << " integration impossible, due aux precisions demandees, on doit augmenter ces precisions ";break;
      case 8  : cout << " cas_kutta ne correspond pas a une routine de base existante ";break;
      case 9  : cout << " calcul correct, mais la derniere valeur a t, n'a pu etre calculee, "
                     << " c'est la valeur a t - quelque chose qui est renvoyee "; break;
      case 0  : cout << " erreur inconnue ";break;
      default : cout << " type d'erreur inconnue ";break;
     }
  };


//==================================== méthodes privées ==============================

    // on défini une fonction a intégrer
Vecteur& Algo_edp::FoncDeriv(const double & t,const Vecteur& f,Vecteur& df,int& erreur)
	{ switch (nb_cas_test)
		{case 1: 
		  	{df = f; // juste pour dimensionner
		  	 df(1)=cos(t);
      erreur = 0;
		  	}break;
		 case 2: 
		  	{df = f; // juste pour dimensionner
		  	 df(1)=cos(t);
      erreur = 0;
		  	}break; 	
		 case 3: // y" + y=0 -> y'=f, f'=-y,  solution cos 
		  	{df(1) = f(2); // juste pour dimensionner
		  	 df(2) = -f(1); 
      erreur = 0;
		  	}break; 	
		 case 4: // deux fonctions sin cos indépendantes
		  	{df(1) = cos(t); // juste pour dimensionner
		  	 df(2) = -sin(t); 
      erreur = 0;
		  	}break;
   default:
     erreur = 1;
     cout << "\n cas de test non implemente !! nb_cas_test="<<nb_cas_test
          << "\n lgo_edp::FoncDeriv(..." << endl;
     Sortie(1);
		};
	  return df;	
	};
// on défini une fonction de vérification d'intégrité
void Algo_edp::FoncVerif_fonc(const double & ,const Vecteur& f,int& erreur)
	{// ici la fonction ne vérifie quasiment rien
  switch (nb_cas_test)
		{case 1: 
		  	{erreur = 0;
		  	}break;
		 case 2: 
		  	{erreur = 0;
		  	}break; 	
		 case 3: // y" + y=0 -> y'=f, f'=-y,  solution cos 
		  	{erreur = 0;
		  	}break; 	
		 case 4: // deux fonctions sin cos indépendantes
		  	{erreur = 0;
		  	}break;
   default:
     erreur = 1;
     cout << "\n cas de test non implemente !! nb_cas_test="<<nb_cas_test
          << "\n lgo_edp::FoncVerif_fonc(..." << endl;
     Sortie(1);
		};
	};

    // méthode test
void Algo_edp::Test_template()
	{ int nb_max_test = 4; // nombre de test actuellement existant
	  // def de la taille du problème et des différentes constantes
	  int taille=0;
	  Vecteur val_initiale(0),der_initiale(0),val_finale(0),der_finale(0);
   // affichage de l'entête
   cout << "\n ----------- test de Pilotage kutta ------------- ";
   //boucle sur les différentes routines de base
   for (int cas_kutta=3;cas_kutta<=5;cas_kutta++)
	  // boucle sur les différents tests
	  for (nb_cas_test=1; nb_cas_test<= nb_max_test;nb_cas_test++)
    { cout << "\n \n    test " << nb_cas_test;
      // def de la taille
      switch (nb_cas_test)
       {case 1: case 2 : {taille = 1; break;}
        case 3: {taille = 2;break;}
        case 4: {taille = 2;break;}
       };
      // def des différentes variables
      double tdebut=0.2,tfin=3.3,erreurAbsolue=1.e-5,erreurRelative=1.e-6;
      double dernierTemps,dernierdeltat; // infos de retour
      int nombreAppelF,nb_step;          // infos de retour
      double erreur_maxi_global;         // infos de retour
      // dimensionement des vecteurs
      val_initiale.Change_taille(taille);
      der_initiale.Change_taille(taille);
      val_finale.Change_taille(taille);
      der_finale.Change_taille(taille);
      // init des valeurs initiales de la fonction
      int nb_boucle = 1; // combien de boucle on veut faire
      switch (nb_cas_test)
       {case 1:{val_initiale(1)=sin(tdebut);nb_boucle = 3;break;}
        case 2:{tfin = - tfin;
                 val_initiale(1)=sin(tdebut);break;}
        case 3:{tdebut = 0.; tfin = 3.* ConstMath::Pi;
                 val_initiale(1)=1.;val_initiale(2)=0.;
                 nb_boucle = 9;
                 break;}
        case 4:{ tdebut = 0.; tfin = 3.* ConstMath::Pi;
                 val_initiale(1)=sin(tdebut);val_initiale(2)=cos(tdebut);
                 nb_boucle = 9;
                 break;}
       };
      int erreur = 0; // init par défaut
      Algo_edp::FoncDeriv(tdebut,val_initiale,der_initiale,erreur);
      if (erreur) // gestion d'une erreur éventuelle
       { cout << "\n erreur 1 dans le calcul de FoncDeriv, le test n'est pas correct "
              << "\n Algo_edp::Test_template() " << endl;
         Sortie(1);
       };
      // premier appel directe de kutta
      double deltat = tfin-tdebut;
      switch (cas_kutta)
        {case 3: Runge_Kutta_step23(*this,& Algo_edp::FoncDeriv,& Algo_edp::FoncVerif_fonc
                                     ,val_initiale,der_initiale,tdebut,deltat
                                     ,val_finale,estime_erreur); break;
         case 4: Runge_Kutta_step34(*this,& Algo_edp::FoncDeriv,& Algo_edp::FoncVerif_fonc
                                     ,val_initiale,der_initiale,tdebut,deltat
                                     ,val_finale,estime_erreur); break;
         case 5: Runge_Kutta_step45(*this,& Algo_edp::FoncDeriv,& Algo_edp::FoncVerif_fonc
                                     ,val_initiale,der_initiale,tdebut,deltat
                                     ,val_finale,estime_erreur); break;            	
        };
      switch (nb_cas_test)
       {case 1:
          {cout << "\n --- appel directe de kutta test 1 --- "
                << "\n valeur finale "<< val_finale(1) << " exacte " << sin(tfin)
                << " erreur " << (val_finale(1)-sin(tfin)) << " estime " << estime_erreur(1);
           break;
           }
        case 2:
          {cout << "\n --- appel directe de kutta test 2 --- "
                << "\n valeur finale "<< val_finale(1) << " exacte " << sin(tfin)
                << " erreur " << (val_finale(1)-sin(tfin)) << " estime " << estime_erreur(1);
           break;
           }
       };
      // appel de la fonction maitre
      double tdeb=tdebut; double intertemps=(tfin-tdebut)/nb_boucle;
      double tfi=tdeb;
     
      for (int iboucle = 1; iboucle <= nb_boucle; iboucle++)		
       { tdeb=tfi; tfi+=intertemps;
         Algo_edp::FoncDeriv(tdeb,val_initiale,der_initiale,erreur);
         if (erreur) // gestion d'une erreur éventuelle
          { cout << "\n erreur dans le calcul de FoncDeriv, le test n'est pas correct "
                 << "\n Algo_edp::Test_template() " << endl;
            Sortie(1);
          };
//         int conver=this->Pilotage_kutta
         this->Pilotage_kutta
                        (cas_kutta,*this,& Algo_edp::FoncDeriv,& Algo_edp::FoncVerif_fonc
                                  ,val_initiale,der_initiale
                                  ,tdeb,tfi,erreurAbsolue,erreurRelative
                                  ,val_finale,der_finale,dernierTemps,dernierdeltat
                                  ,nombreAppelF,nb_step,erreur_maxi_global);
         // affichage des informations
         switch (nb_cas_test)
          {case 1:
             {cout << "\n iter= " << iboucle << " valeur finale "<< val_finale(1) << " exacte " << sin(tfi)
                   << " erreur " << (val_finale(1)-sin(tfi)) << " estime " << erreur_maxi_global
                   << " nb step " << nb_step << " nombre d'appel " << nombreAppelF
                   << " dernier temps " << dernierTemps << " dernier deltat " << dernierdeltat;
              break;
              }
           case 2:
             {cout << "\n iter= " << iboucle << " valeur finale "<< val_finale(1) << " exacte " << sin(tfi)
                   << " erreur " << (val_finale(1)-sin(tfi)) << " estime " << erreur_maxi_global
                   << " nb step " << nb_step << " nombre d'appel " << nombreAppelF
                   << " dernier temps " << dernierTemps << " dernier deltat " << dernierdeltat;
              break;
              }
           case 3 :
             {cout << "\n iter= " << iboucle << " valfin1 "<< val_finale(1) << " exact " << cos(tfi)
                   << " valfin2 "<< val_finale(2) << " exact " << -sin(tfi)
                   << " err1 " << (val_finale(1)-cos(tfi)) 
                   << " err2 " << (val_finale(2)+sin(tfi)) << " estime " << erreur_maxi_global
                   << " nb step " << nb_step << " nbappel " << nombreAppelF
                   << " dertemps " << dernierTemps << " derdeltat " << dernierdeltat;
              break;
              }
           case 4 :  
             {cout << "\n iter= " << iboucle << " valfin1 "<< val_finale(1) << " exact " << sin(tfi)
                   << " valfin2 "<< val_finale(2) << " exact " << cos(tfi)
                   << " err1 " << (val_finale(1)-sin(tfi)) 
                   << " err2 " << (val_finale(2)-cos(tfi)) << " estime " << erreur_maxi_global
                   << " nb step " << nb_step << " nbappel " << nombreAppelF
                   << " dertemps " << dernierTemps << " derdeltat " << dernierdeltat;
              break;
              }
          };
         // prépa de la boucle suivante
         val_initiale=val_finale; der_initiale=der_finale;		
       };
		};
   // fin programmé
   Sortie(1);                    
    
	};


//**********************************************************************
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Modified:
//
//    24 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
void Algo_edp::timestamp ( void )

{
#define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
#undef TIME_SIZE
};