// 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 "Algori.h"
#include "string"
#include "MathUtil.h"
#include <iostream>   //  pour utiliser la classe istrstream
#include <strstream>  //      nouveau dans CW5.3


//====================================================================================================== 
//                 troisième fichier pour l'implantation des méthodes d'Algori       
//======================================================================================================


// lecture des paramètres du calcul généraux pour tous les algos (dans le cas où il y en a)
// qui peuvent éventuellement ne pas servir !!
void Algori::lecture_Parametres(UtilLecture& entreePrinc)
{ bool a_lire = false;
  MotCle motCle; // ref aux mots cle
  // dans le cas où l'entête des paramètres n'a pas été lue par l'algorithme spécifique
  // on regarde s'il existe
  if (deja_lue_entete_parametre == 2)
   // on peut directement lire
   {a_lire = true;}
  else  // sinon
   { // deux cas suivant que l'algo spécifique a une méthode de lecture de paramètre ou non
     if (deja_lue_entete_parametre == 0)
       { // cas ou l'algo spécifique n'a pas de méthode de lecture de paramètre donc il n'y a pas
         // eu d'essai de lecture de l'entête 
         // on se positionne sur le prochain mot clé
         do { entreePrinc.NouvelleDonnee();}
         while ( !motCle.SimotCle(entreePrinc.tablcar)) ;
         // si le mot clé est "PARA_TYPE_DE_CALCUL" cela signifie 
         // qu'il y a une zone de paramètre éventuellement  à lire
         if (strstr(entreePrinc.tablcar,"PARA_TYPE_DE_CALCUL")!=NULL) 
          {entreePrinc.NouvelleDonnee();
           a_lire = true;
           }
         else // sinon le mot clé ne concerne pas les paramètres il n'y a rien n'a lire
          {a_lire = false;
          };
       }
     else if (deja_lue_entete_parametre == 1) 
       { // cas ou l'algo spécifique a une méthode de lecture de paramètre donc il a
           // déjà essayer de lire  l'entête et n'a rien trouvé: il n'y a rien n'a lire
           a_lire = false;
         }
     else
       { cout << "\n erreur deja_lue_entete_parametre = " << deja_lue_entete_parametre
              << " alors qu'il devrait etre soit 1 ou 0 ou 2 !! "
              << "\n Algori::lecture_Parametres(....";
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
         if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
         Sortie(1);
       };
   };
  // maintenant on lit en fonction du booléen a_lire     
  if (a_lire)
  {// on s'arrête lorsqu'il n'y a plus de mot clé spécifique à lire
		 while (    (strstr(entreePrinc.tablcar,"parametre_calcul_de_la_viscosite_")!=NULL)
				  || (strstr(entreePrinc.tablcar,"avec_amortissement_cinetique_")!=NULL)
				  || (strstr(entreePrinc.tablcar,"avec_au_noeud_amortissement_cinetique_")!=NULL)
				  || (strstr(entreePrinc.tablcar,"ARRET_A_EQUILIBRE_STATIQUE_")!=NULL)
				  || (strstr(entreePrinc.tablcar,"mode_debug_=")!=NULL)
      || (strstr(entreePrinc.tablcar,"permet_affichage_")!=NULL)
      || (strstr(entreePrinc.tablcar,"modulation_precision_equilibre_=")!=NULL)
				 ) 
						 
		 {bool lecture_effective = false;
    string nom_class_methode("Algori::lecture_Parametres((..."); // pour la lecture
    int min=1,max=0;string mot_cle;int int_defaut=0; double double_defaut=0;// variables inter qui servent pour la lecture
			 //=== parametres de controle du calcul de la matrice de viscosite =====
    if (strstr(entreePrinc.tablcar,"parametre_calcul_de_la_viscosite_")!=NULL)
     { mot_cle="parametre_calcul_de_la_viscosite_";
       lecture_effective = entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
       // si on arrive ici c'est que l'on a bien lue le mot clé "parametre_calcul_de_la_viscosite_"
       // donc il faudra passer une ligne de toute manière, donc ce n'est plus la peine de tester les booleen
       
       // --- lecture éventuelle de type_calcul_visqu_critique
       int_defaut=0; mot_cle="type_calcul_visqu_critique="; min=1;max=3; 
       entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,type_calcul_visqu_critique);
       // on change le paramètre de viscosité numérique général en conséquence
       if (type_calcul_visqu_critique)
         pa.Change_amort_visco_artificielle(4); // on indique qu'il s'agit d'une gestion par un amortissement visqueux critique 
        // utilisé par de la relaxation dynamique
       // --- lecture éventuelle de opt_cal_C_critique
       int_defaut=-1; mot_cle="opt_cal_C_critique="; min=-1;max=2; 
       entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,opt_cal_C_critique);
       // --- lecture éventuelle de f_
       double_defaut=1.; mot_cle="f_="; min=1;max=0; 
       entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,f_);
       // --- lecture éventuelle de ampli_visco
       double_defaut=1.; mot_cle="ampli_visco_="; min=1;max=0; 
       entreePrinc.Lecture_un_parametre_double(double_defaut,nom_class_methode,min,max,mot_cle,ampli_visco);
     };
//			else // sinon on met les valeurs par défaut
//				{ type_calcul_visqu_critique= 0;
//				  opt_cal_C_critique= 1; f_= 1.; ampli_visco= 1.;
//				};
    //===  cas où l'on veut de l'amortissement cinétique global avec éventuellement des paramètres particuliers de contrôle
    if (strstr(entreePrinc.tablcar,"avec_amortissement_cinetique_")!=NULL)
     { mot_cle="avec_amortissement_cinetique_";
       entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
       amortissement_cinetique=true; // on active l'amortissement cinétique
       string st1;
       // lecture des paramètres de l'amortissement cinétique
       while (st1 != "fi_parametre_amortissement_cinetique_")
       { entreePrinc.NouvelleDonnee();
       *(entreePrinc.entree) >> st1; 
       if (st1 == "max_nb_decroit_pourRelaxDyn_")
        {*(entreePrinc.entree) >> max_nb_decroit_pourRelaxDyn; }
       else if (st1 == "inter_nb_entre_relax_")
          {// il faut que l'on regarde s'il y a éventuellement un remplacement par une
           // fonction nD
           if (strstr(entreePrinc.tablcar,"nom_fct_nD_inter_nb_entre_relax_")!=NULL)
             {mot_cle="nom_fct_nD_inter_nb_entre_relax_";
              entreePrinc.Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_fct_nD_inter_nb_entre_relax);
             }
           else // cas sans fonction nD
             {nom_fct_nD_inter_nb_entre_relax="";
              *(entreePrinc.entree) >> inter_nb_entre_relax;
             };
          }
       else if (st1 == "taille_moyenne_glissante_")
          {*(entreePrinc.entree) >> taille_moyenne_glissante;}
       else if (st1 == "test_fraction_energie_")
          {*(entreePrinc.entree) >> test_fraction_energie;}		
       else if  (st1 == "coef_arret_pourRelaxDyn_")
        {*(entreePrinc.entree) >> coef_arret_pourRelaxDyn; }
       else if  (st1 == "coef_redemarrage_pourRelaxDyn_")
        {*(entreePrinc.entree) >> coef_redemarrage_pourRelaxDyn; }
       else if  (st1 == "max_deltaX_pourRelaxDyn_")
        {*(entreePrinc.entree) >> max_deltaX_pourRelaxDyn; }
       else if  (st1 == "nb_max_dX_OK_pourRelaxDyn_")
        {*(entreePrinc.entree) >> nb_max_dX_OK_pourRelaxDyn; }
       else if  (st1 == "nb_deb_testfin_pourRelaxDyn_")
        {*(entreePrinc.entree) >> nb_deb_testfin_pourRelaxDyn; }
       else if  (st1 == "nb_deb_test_amort_cinetique_")
        {*(entreePrinc.entree) >> nb_deb_test_amort_cinetique;};
                
       if ((st1 != "max_nb_decroit_pourRelaxDyn_")
        && (st1 != "inter_nb_entre_relax_")
        && (st1 != "taille_moyenne_glissante_")
        && (st1 != "test_fraction_energie_")
        && (st1 != "coef_arret_pourRelaxDyn_")
        && (st1 != "coef_redemarrage_pourRelaxDyn_")
        && (st1 != "max_deltaX_pourRelaxDyn_")
        && (st1 != "nb_max_dX_OK_pourRelaxDyn_")
        && (st1 != "nb_deb_testfin_pourRelaxDyn_") 
        && (st1 != "nb_deb_test_amort_cinetique_")
        && (st1 != "fi_parametre_amortissement_cinetique_"))
        {   cout << "\n erreur en lecture du coef d'arret pour l'amortissement cinetique"
                 << "\n on attendait un mot cle , au lieu de " << st1
                 << "\n Algori::lecture_Parametres( ... ";
            // dans le cas où un comptage du calcul est en cours on l'arrête
            if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
            if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
            Sortie(1);
        };
                 };
       // préparation pour une nouvelle donnée
       lecture_effective= true;
     }; 
    //===  idem mais pour de l'amortissement cinétique à chaque noeud
    if (strstr(entreePrinc.tablcar,"avec_au_noeud_amortissement_cinetique_")!=NULL)
     { mot_cle="avec_au_noeud_amortissement_cinetique_";
       entreePrinc.Lecture_et_verif_mot_cle(nom_class_methode,mot_cle);
       amortissement_cinetique_au_noeud=true; // on active l'amortissement cinétique
       string st1;
       // lecture des paramètres de l'amortissement cinétique
       while (st1 != "fi_parametre_au_noeud_amortissement_cinetique_")
       { entreePrinc.NouvelleDonnee();
         *(entreePrinc.entree) >> st1;
         if (st1 == "max_nb_decroit_pourRelaxDyn_noe_")
          {*(entreePrinc.entree) >> max_nb_decroit_pourRelaxDyn_noe; }
         else if  (st1 == "coef_arret_pourRelaxDyn_noe_")
          {*(entreePrinc.entree) >> coef_arret_pourRelaxDyn_noe; }
         else if  (st1 == "coef_redemarrage_pourRelaxDyn_noe_")
          {*(entreePrinc.entree) >> coef_redemarrage_pourRelaxDyn_noe; }
         else if  (st1 == "nb_deb_test_amort_cinetique_noe_")
          {*(entreePrinc.entree) >> nb_deb_test_amort_cinetique_noe;};
                  
         if ((st1 != "max_nb_decroit_pourRelaxDyn_noe_")
          && (st1 != "coef_arret_pourRelaxDyn_noe_")
          && (st1 != "coef_redemarrage_pourRelaxDyn_noe_")
          && (st1 != "nb_deb_test_amort_cinetique_noe_") 
          && (st1 != "fi_parametre_au_noeud_amortissement_cinetique_"))
           { cout << "\n erreur en lecture du coef d'arret pour l'amortissement cinetique"
                  << "\n on attendait un mot cle , au lieu de " << st1
                  << "\n Algori::lecture_Parametres( ... ";
             // dans le cas où un comptage du calcul est en cours on l'arrête
             if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
             if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
             Sortie(1);
           };
       };
       // préparation pour une nouvelle donnée
       lecture_effective= true;
     }; 
    //===  contrôle éventuel d'un arrêt sur équilibre statique
    if (strstr(entreePrinc.tablcar,"ARRET_A_EQUILIBRE_STATIQUE_")!=NULL)
     { // on commence par regarder si l'on est en dynamique, si non on met un warning
       if (Evolution_temporelle_du_calcul(typeCalcul) != DYNAMIQUE)
				      cout << "\n *** Attention, on demande un arret sur l'equilibre statique, (cf. para de l'algo) "
							        << " alors que l'algorithme n'est pas dynamique ?????? il y a peut-etre une erreur de parametres " << endl;
				   // --- lecture éventuelle de opt_cal_C_critique
				   int_defaut=0; mot_cle="ARRET_A_EQUILIBRE_STATIQUE_"; min=0;max=2; lecture_effective= false;
				   lecture_effective = entreePrinc.Lecture_un_parametre_int(int_defaut,nom_class_methode,min,max,mot_cle,arret_a_equilibre_statique);
              //  sert aussi en fait pour retirer éventuellement la partie visqueuse numérique
				 };

    //===  contrôle d'un mode débug éventuel
    // cas où le paramètre de mode débug existe
    if (strstr(entreePrinc.tablcar,"mode_debug_=")!=NULL)
				 { int_defaut=0; mot_cle="mode_debug_=";min=1;max=ConstMath::tresgrand;
					  { bool lec_eff = entreePrinc.Lecture_un_parametre_int
                        (int_defaut,nom_class_methode,min,max,mot_cle,mode_debug);
					    lecture_effective = lecture_effective || lec_eff;
       };
     };
    
    //===  contrôle d'un niveau de commentaire spécifique à l'algo
    if (strstr(entreePrinc.tablcar,"permet_affichage_")!=NULL)
     { int_defaut=0; mot_cle="permet_affichage_";min=1;max=ConstMath::tresgrand;
       { bool lec_eff = entreePrinc.Lecture_un_parametre_int
                        (int_defaut,nom_class_methode,min,max,mot_cle,permet_affichage);
         lecture_effective = lecture_effective || lec_eff;
       };
     };

    //===  cas où on veut moduler la précision de la convergence globale via une fct nD
    if (strstr(entreePrinc.tablcar,"modulation_precision_equilibre_=")!=NULL)
     { mot_cle="modulation_precision_equilibre_=";
       { bool lec_eff = entreePrinc.Lecture_mot_cle_et_string
                           (nom_class_methode,mot_cle,nom_modulation_precision);
         lecture_effective = lecture_effective || lec_eff;
       };
     };    

    if ((lecture_effective) && ( !motCle.SimotCle(entreePrinc.tablcar)))
				  {entreePrinc.NouvelleDonnee(); // si on a lue des para sur la ère ligne on passe à la suivante
				   lecture_effective=false; 
				  };
		 };																								
		 // si le paramètre est activé, on définit le vecteur intermédiaire de sauvegarde du résidu statique
   // un pb cependant, la lecture de l'amortissement artificielle n'est éventuellement pas encore fait
   // c'est le cas d'une première lecture du .info, l'amortissement visqueux s'effectue au niveau des paramètres
   // de contrôle donc bien après l'algo, dans ce cas il faut prévoir dans l'initialisation des algo spécifiques
   // le dimensionnement de vglob_stat
		 if (   (arret_a_equilibre_statique ) // le dimensionnement se fera à l'affectation (=)
			  && (type_calcul_visqu_critique || pa.Amort_visco_artificielle()))
			{ if (vglob_stat == NULL) vglob_stat  = new Vecteur();}
		 else // dans le cas on ne veut pas un arret statique spécifique sans viscosité numérique
			{ if (vglob_stat != NULL) {delete vglob_stat; vglob_stat=NULL ;};}; 	
		 // on se positionne éventuellement sur le prochain mot clé si on n'y est pas déjà
		 while ( !motCle.SimotCle(entreePrinc.tablcar)) 
              entreePrinc.NouvelleDonnee();
  }; //-- fin du cas où on a regardé 
 };
  
// création d'un fichier de commande: cas des paramètres spécifiques
// par défaut on indique qu'il n'y a pas de paramètres.
void Algori::Info_com_parametres(UtilLecture& entreePrinc)
 { // écriture dans le fichier de commande
   ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
	  cout << "\n# exemple para facultatifs generaux d'algorithme (rep o ou n defaut) ? ";
	  string rep;
        rep = "_";
        // procédure de lecture avec prise en charge d'un retour chariot
        rep = lect_return_defaut(true,"n");
        // si la taille de rep == 0 cela veut dire que c'est un retour chariot
        if (rep.size()==0)
           {rep = "n";
            cout << "--> valeur par defaut : n "<<endl;
           };
	  if ((rep == "o") || (rep == "O") ||(rep == "0"))
	    { rep=" "; // init
       while ((Minuscules(rep) != "f")&&(Minuscules(rep) != "0"))
        {
         try
          {
           cout 
              << "\n (0 ou f)  (fin)               "
              << "\n (1)  amortissement cinetique "
              << "\n (2)  calcul d'une matrice de viscosite numerique "
              << "\n      pour de l'amortissement critique par exemple "
              << "\n (3)  parametres supplementaires de controle du type d'arret de l'algo "
              << "\n (4)  parametres de controle du mod debug "
              << "\n (5)  def de fichiers externes contenant des contraintes, deformations ..."
              << "\n      preexistantes avant le debut du calcul "
              << "\n (6)  niveau de commentaire specifique aux algos "
              << "\n  ";
           // procédure de lecture avec prise en charge d'un retour chariot
           rep = lect_return_defaut(false,"f");
           if ((Minuscules(rep) == "f") || (Minuscules(rep) == "0"))// sortie directe
             break;
           int num = ChangeEntier(rep);
           bool choix_valide=false;
           if ((num >= 0)&&(num< 7))
                { choix_valide=true; }
           else { cout << "\n Erreur on attendait un entier entre 0 et 6 !!, "
                       << "\n redonnez une bonne valeur"
                       << "\n ou taper f ou 0 pour arreter le programme"; 
                  choix_valide=false;
                }
           switch (num)
            { case 0:  //fin
                { break;} // normalement cela a déjà été filtré avant
              case 1:  // amortissement cinetique
                {
    sort << "\n#---------------------------------------------------------------------"
         << "\n|  parametres communs a differents algorithmes (voir doc specifique) |"
         << "\n#---------------------------------------------------------------------"
         << "\n  #   IMPORTANT: les parametres qui suivent sont à mettre apres ceux specifique de l'algorithme !!! "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  # ---   parametres de controle du calcul de l'amortissement cinetique  --- "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  #  Il est possible de prendre en compte un amortissement cinetique, dont l'objectif  "
         << "\n  #  est d'obtenir une solution statique a l'aide d'un algorithme dynamique."
         << "\n  #  Pour cela on introduit le mot cle: avec_amortissement_cinetique_ suivi "
         << "\n  #  de mots cles eventuels qui precisent des parametres. Ce groupe de parametres "
         << "\n  #  se termine par le mot cle fi_parametre_amortissement_cinetique_ "
         << "\n  #  Chaque parametre est suivi d'un nombre. Chaque parametre + nombre, doit etre sur une ligne a part."
         << "\n  #   deux techniques pour appliquer le test, qui conduit a mettre les vitesses à 0 "
         << "\n  #   1) on regarde si l'|energie cinetique| diminue n1 fois (consécutives) "
         << "\n  #      ou est inferieure n1 fois  fraction * dernier_max_E_cin, si fraction (= test_fraction_energie) est non nulle "
         << "\n  #      NB: le test ne se declanche que lorsque l'on a scrute au moins n1 valeurs de l'energie cinetique "
         << "\n  #   2) on regarde si la moyenne glissante de n2 |energie cinetique| consecutive diminue n1 fois (consecutives) "
         << "\n  #      ou est inferieure n1 fois  fraction * dernier_max_E_cin, si fraction (= test_fraction_energie) est non nulle "
         << "\n  #      NB: le test ne se declanche que lorsque l'on a au moins n2 valeurs consecutives disponibles de l'energie cinetique "
         << "\n  #   NB: le cas 1) est en fait un cas particulier du 2) pour lequel n2 == 1 "                           
         << "\n  #   Les parametres sont: "
         << "\n  #   max_nb_decroit_pourRelaxDyn_ : nombre n2 de diminution de de la moyenne glissante de l'energie cinetique a"
         << "\n  #       (Par exemple si = 3, cela signifie qu'il faut attendre 3 baisses de "
         << "\n  #        l'energie cinetique pour valider la relaxation. Ces 3 baisses peuvent apparaître au "
         << "\n  #        cours de n iterations, n etant plus grand que 3 ! De plus pendant ces en iterations,"
         << "\n  #        on peut egalement avoir des elevations de l'energie cinetique suivi de diminutions."
         << "\n  #        Seules les diminutions sont prises en compte. Apres une relaxation, le compteur de "
         << "\n  #        diminution est remis à 0 "
         << "\n  #       )"
         << "\n  #                                  partir de laquelle on met les vitesses a 0 "
         << "\n  #   inter_nb_entre_relax_        : intervalle autorise entre deux relaxation "
         << "\n  #       (Par exemple si = 10, cela signifie qu'apres une relaxation, une nouvelle relaxation"
         << "\n  #        ne sera accepte qu'apres 10 nouvelles iterations. "
         << "\n  #        Il est egalement possible de remplacer la valeur numerique par le resultat d'une fonction nD "
         << "\n  #        pour cela on indique le mot cle : nom_fct_nD_inter_nb_entre_relax_ "
         << "\n  #        suivi du nom de la fonction nD qui sera appelee a chaque fois que l'on aura besoin       "
         << "\n  #        du parametre inter_nb_entre_relax_  "
         << "\n  #       )"
         << "\n  #   taille_moyenne_glissante_    : valeur du nombre n2 de la moyenne glissante d'energie cinetique "
         << "\n  #   test_fraction_energie_       : si diff de 0, indique la valeur de fraction    " 
         << "\n  #   coef_arret_pourRelaxDyn_     : arret de l'amortissement lorsque l'energie cinetique est"
         << "\n  #                                  inferieur à ce coef * le dernier pic d'energie cinetique "
         << "\n  #   coef_redemarrage_pourRelaxDyn_ : redemarrage de l'amortissement lorsque l'energie cinetique"
         << "\n  #                                    est sup a ce coef * le maxi des pic d'energie enregistre"
         << "\n  #   max_deltaX_pourRelaxDyn_     : limite inferieur de delta x_(t->tdt) a partir de laquelle"
         << "\n  #                                  on arrete le calcul"
         << "\n  #   nb_max_dX_OK_pourRelaxDyn_   : nombre de fois que le critere precedent doit etre satisfait"
         << "\n  #                                  pour que le calcul s'arrete"
         << "\n  #   nb_deb_testfin_pourRelaxDyn_ : le nombre mini de passage dans l'amortissement, a partir duquel on test la fin"
         << "\n  #   nb_deb_test_amort_cinetique_ : le nombre mini d'iteration, a partir duquel on demarre l'algorithme"
         << "\n  #   fi_parametre_amortissement_cinetique_ : mot cle obligatoire indiquant la fin des parametres"
         << "\n  #   exemple de declaration"
         << "\n  #   "
         << "\n  #   avec_amortissement_cinetique_  "
         << "\n  #      max_nb_decroit_pourRelaxDyn_   1 "
         << "\n  #      inter_nb_entre_relax_          1 "
         << "\n  # ou "
         << "\n  #      inter_nb_entre_relax_          nom_fct_nD_inter_nb_entre_relax_   toto  "
         << "\n  #      taille_moyenne_glissante_      1 "
         << "\n  #      coef_arret_pourRelaxDyn_       0. "
         << "\n  #      coef_redemarrage_pourRelaxDyn_ 0.01 "
         << "\n  #      max_deltaX_pourRelaxDyn_       0.4 "
         << "\n  #      nb_max_dX_OK_pourRelaxDyn_     5 "
         << "\n  #      nb_deb_testfin_pourRelaxDyn_   100 "
         << "\n  #      nb_deb_test_amort_cinetique_   100 "
         << "\n  #   fi_parametre_amortissement_cinetique_ "
         << "\n  #   "
         << endl;
                  break;}
              case 2:  // calcul d'une matrice de viscosite numerique
                {
    sort << "\n#---------------------------------------------------------------------"
         << "\n|  parametres communs a differents algorithmes (voir doc specifique) |"
         << "\n#---------------------------------------------------------------------"
         << "\n  #   IMPORTANT: les parametres qui suivent sont à mettre apres ceux specifique de l'algorithme !!! "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  # ---   parametres de controle du calcul de la matrice de viscosite  --- "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  #  Il est possible de prendre en compte de l'amortissement critique numerique "
         << "\n  #   classiquement utilise dans le cas de la relaxation dynamique "
         << "\n  #   "
         << "\n  #   ** PARAMETRES FACULTATIF **   "
         << "\n  #  on indique tout d'abord le mot cle : parametre_calcul_de_la_viscosite_  suivi des parametres associes eventuelles "
         << "\n  #  type_calcul_visqu_critique=  suivi du type de calcul  "
         << "\n  # si type_calcul_visqu_critique= 1 : la matrice C est calculee selon: C_{ii} = c * M_{ii} avec c = 2 * omega_0  "
         << "\n  # si type_calcul_visqu_critique= 2 : la matrice C est calculee selon: C_{ii} = c * M_{ii}   "
         << "\n  #												avec c = sqrt( 4*omega_0^2 - omega_0^4)"
         << "\n  # si type_calcul_visqu_critique= 3 : la matrice C est calculee selon: C_{ii} = c * M_{ii}   "
         << "\n  #												avec c = 2*sqrt(omega_0/(1+omega_0)) "
         << "\n  #                        "
         << "\n  #                                    puis il y a un choix via le parametre  opt_cal_C_critique=  "
         << "\n  #      opt_cal_C_critique= 0 :  (omega_0)^2 \approx  (Delta X^T Delta R_{i(statique)}^n) / (Delta X^T M Delta X) "
         << "\n  #      opt_cal_C_critique= 1 :  (omega_0)^2 \approx  (dot X^T Delta R_{i(statique)}^n) / (dot X^T M dot X)  "
         << "\n  #                        "
         << "\n  #				a) si (omega_0)^2 est négative on met omega_0^2=0 "
         << "\n  #				b) si (omega_0)^2 > f^2 * 4, on pose (omega_0) = f^2 * 4, par defaut: f= 0.9 "
         << "\n  #      on peut changer la valeur de f avec le parametre: f_= suivi de la valeur  "
         << "\n  #                        "
         << "\n  #      il est possible egalement d'amplifier la valeur de la viscosite calculee (avec cependant toujours les limites a) et b) "
         << "\n  #      ampli_visco_=   val_ampli   "
         << "\n  #                        "
         << "\n  #                             Exemple de declaration:  "
         << "\n  ##   parametre_calcul_de_la_viscosite_  type_calcul_visqu_critique= 1 opt_cal_C_critique= 0  f_= 1.3 "
         << "\n  #                       "
         << "\n  #   tous les parametres sont facultatifs: les valeurs par defaut sont:   "
         << "\n  #   opt_cal_C_critique= 1 ; f_= 0.9 ; ampli_visco_= 1.; "
         << "\n  #                        "
         << "\n  #                        "
         << endl;
                break;}
              case 3:  // parametres supplementaires de controle du type d'arret de l'algo
                {
    sort << "\n#---------------------------------------------------------------------"
         << "\n|  parametres communs a differents algorithmes (voir doc specifique) |"
         << "\n#---------------------------------------------------------------------"
         << "\n  #   IMPORTANT: les parametres qui suivent sont à mettre apres ceux specifique de l'algorithme !!! "
         << "\n  #              apres l'amortissement cinetique, apres la def de calcul de  viscosite et les autres qui precedent!!! "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  # ---   parametres de controle du type d'arret                       --- "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  #  ce parametre est utile dans le cas d'un calcul dynamique "
         << "\n  #  on peut egalement indiquer un arret correspondant a un equilibre statique. Dans ce cas "
         << "\n  #  le calcul s'arrete lorsqu'on considere que l'on a atteind un regime statique dont le critere est:"
         << "\n  #   - le residu statique est inferieur au maxi indique dans les parametres generaux de controle"
         << "\n  #   - si l'on utilise de l'amortissement cinetique: le critere d'arret est celui de cet algo"
         << "\n  #   Pour activer cette possibilite d'arret il faut indiquer sur la ligne qui suit le mot cle"
         << "\n  #   ARRET_A_EQUILIBRE_STATIQUE_  suivant d'un nombre "
         << "\n  #    ARRET_A_EQUILIBRE_STATIQUE_ 1  signifie que le residu statique est utilise comme critere"
         << "\n  #    ARRET_A_EQUILIBRE_STATIQUE_ 2  le residu statique et le critere de l'amortissement cinetique "
         << "\n  #                                   sont conjointement utilises, on retient le plus defavorable"
         << "\n  #    ARRET_A_EQUILIBRE_STATIQUE_ 0  pas de prise en compte de ce parametre (valeur par defaut)"
         << "\n  #   "
         << "\n  #    exemple d'utilisation "
         << "\n  #   "
         << "\n  #   ARRET_A_EQUILIBRE_STATIQUE_ 1 "
         << "\n  #   "
         << endl;
                break;}
              case 4:  // parametres de controle du mod debug
                {
     sort << "\n#---------------------------------------------------------------------"
         << "\n|  parametres communs a differents algorithmes (voir doc specifique) |"
         << "\n#---------------------------------------------------------------------"
         << "\n  #   IMPORTANT: les parametres qui suivent sont à mettre apres ceux specifique de l'algorithme !!! "
         << "\n  #              apres les parametres de controle de type d'arret et les autres qui precedent !!! "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  # ---                parametres de controle du mode debug            --- "
         << "\n  # ---------------------------------------------------------------------- "
         << "\n  #   ** PARAMETRES FACULTATIF **   "
         << "\n  #  mode_debug_=  parametre de debug (= 0 par defaut ->  calcul classique) "
         << "\n  #             =  n : a chaque n iteration on sort une visualisation "
         << "\n  ##  mode_debug_= 0   # a mettre apres les para de relaxation           "
         << "\n  #           "
         << "\n  #  modulation_precision_equilibre_=  fonction_nD "
         << "\n  #  la fonction fonction_nD module par multiplication, la precision en cours "
         << "\n  # **important**: ces parametres doivent etre mis chacun sur une ligne differente  "
         << endl;
                break;}
              case 5:  // def de fichiers externes contenant des contraintes, deformations ...
                {
                   sort << "\n  # --------------------------------------------------------------------------------- "
                        << "\n  # ---   lecture sur fichiers externes de grandeurs preexistantes au calcul      --- "
                        << "\n  # --------------------------------------------------------------------------------- "
                        << "\n  #  ce referer a la doc pour la syntaxe dans les fichiers  \n"
                        << "\n   flotExterne  #--- #  mot cle pour introduire la lecture \n";
         
                   string repi= " "; // init
                   rep=" "; // init
                   while ((Minuscules(repi) != "f")&&(Minuscules(repi) != "0"))
                    {
                     try
                       {
                         cout
                            << "\n (0 ou f)  (fin)               "
                            << "\n (1)  nouveau fichier de contraintes aux point d'integration "
                            << "\n (2)  nouveau fichier de deplacement aux noeuds "
                            << "\n  ";
                         // procédure de lecture avec prise en charge d'un retour chariot
                         repi = lect_return_defaut(false,"f");
                         if ((Minuscules(repi) == "f") || (Minuscules(repi) == "0"))// sortie directe
                           break;
                         int numi = ChangeEntier(repi);
                         bool choix_valide=false;
                         if ((numi >= 0)&&(numi<=2))
                              { choix_valide=true; }
                         else { cout << "\n Erreur on attendait un entier entre 0 et 2 !!, "
                                     << "\n redonnez une bonne valeur"
                                     << "\n ou taper f ou 0 pour arreter le programme"; 
                                choix_valide=false;
                              }
                         switch (numi)
                          { case 0:  //fin
                             { break;} // normalement cela a déjà été filtré avant
                            case 1:  // nouveau fichier de contraintes aux point d'integration
                             { string nom;
                               cout << "\n nom du fichier sans extension? (mais qui doit avoir l'extension _cab.isoe ) ";
                               nom = lect_chaine();
                               sort << "\n  1  " << nom << "_cab.isoe   " ;
                               cout << "\n nom du maillage associe au fichier ? ";
                               nom = lect_chaine();
                               sort << nom << endl;
                               break;
                             }
                            case 2:  // nouveau fichier de contraintes aux point d'integration
                             { string nom;
                               cout << "\n nom du fichier sans extension? (mais qui doit avoir l'extension _dpl.points ) ";
                               nom = lect_chaine();
                               sort << "\n  2  " << nom << "_dpl.points   " ;
                               cout << "\n nom du maillage associe au fichier ? ";
                               nom = lect_chaine();
                               sort << nom << endl;
                               break;
                             }
                            default:
                              cout << "\n le cas "<<repi<<" n'est pas encore traite !!, mais cela devrait ce faire sans tarder ";
                           };
                       }
                     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 on attendait un des mots cles proposés !!, "
                               << "\n redonnez une bonne valeur"
                               << "\n ou taper f ou 0 pour sortir "; 
                       };
                    }; //-- fin du while
                break;}
              case 6:  // parametres de controle du mod debug
                {
                 sort << "\n#---------------------------------------------------------------------"
                     << "\n|  parametres communs a differents algorithmes (voir doc specifique) |"
                     << "\n#---------------------------------------------------------------------"
                     << "\n  #   IMPORTANT: les parametres qui suivent sont à mettre apres ceux specifique de l'algorithme !!! "
                     << "\n  #              apres les parametres de controle de type d'arret et les autres qui precedent !!! "
                     << "\n  # ---------------------------------------------------------------------- "
                     << "\n  # ---    parametres de controle du niveau de commentaire             --- "
                     << "\n  # ---------------------------------------------------------------------- "
                     << "\n  #   ** PARAMETRES FACULTATIF **   "
                     << "\n  #  permet_affichage_  3 "
                     << "\n  # **important**: ce parametre doit etre mis  sur une ligne seule  "
                     << endl;
                 break;
                }
              default:
                cout << "\n le cas "<<rep<<" n'est pas encore traite !!, mais cela devrait ce faire sans tarder ";
            };
          }
         catch (ErrSortieFinale)
              // cas d'une direction voulue vers la sortie
              // on relance l'interuption pour le niveau supérieur
            { ErrSortieFinale toto;
              throw (toto);
            }
         catch (...)//(UtilLecture::ErrNouvelleDonnee erreur)
          {  cout << "\n Erreur on attendait un des mots cles proposés !!, "
                  << "\n redonnez une bonne valeur"
                  << "\n ou taper f ou 0 pour sortir "; 
          };
        }; //-- fin du while
	    };
   sort << "\n" << endl;
 };

// initialisation d'un calcul de remonté à des variables secondaires ou d'erreur
// globalement, cette initialisation dépend des paramètres utilisateurs
// ramène true s'il y a initialisation
bool Algori::InitRemont(LesMaillages * lesMail,LesReferences* lesRef, DiversStockage* diversStockage,Charge* charge,
                    LesCondLim* lesCondLim,LesContacts* lesContacts,Resultats* resultats)
   {bool avec_remont=false;

    // on parcourt la liste soustypeDeCalcul si elle n'est pas vide
    if (soustypeDeCalcul->size() > 0)
     {if ((find(soustypeDeCalcul->begin(),soustypeDeCalcul->end(),prevision_visu_sigma) != soustypeDeCalcul->end())
          || (find(soustypeDeCalcul->begin(),soustypeDeCalcul->end(),avec_remonte_erreur) != soustypeDeCalcul->end()))
        {Algori::InitRemontSigma(lesMail,lesRef,diversStockage,charge,lesCondLim,lesContacts,resultats);
         Algori::InitErreur(lesMail,lesRef,diversStockage,charge,lesCondLim,lesContacts,resultats);
         avec_remont=true;avec_remont_sigma=true;}
      if (find(soustypeDeCalcul->begin(),soustypeDeCalcul->end(),prevision_visu_epsilon) != soustypeDeCalcul->end())
        {Algori::InitRemontEps(lesMail,lesRef,diversStockage,charge,lesCondLim,lesContacts,resultats);
         avec_remont=true;avec_remont_eps=true;}
      if (find(soustypeDeCalcul->begin(),soustypeDeCalcul->end(),prevision_visu_erreur) != soustypeDeCalcul->end())
        {if (!avec_remont_sigma)
            {Algori::InitRemontSigma(lesMail,lesRef,diversStockage,charge,lesCondLim,lesContacts,resultats);
             avec_remont_sigma=true;}
          Algori::InitErreur(lesMail,lesRef,diversStockage,charge,lesCondLim,lesContacts,resultats);
          avec_remont=true;avec_erreur=true;}
     };

    return avec_remont;
    };   

// calcul d'une remontée à des variables secondaires ou d'erreur, ceci pour un post traitement de visualisation
// ce calcul dépend des paramètres utilisateurs, il n'est fait qu'à chaque sortie .BI
// ramène true s'il y a un calcul effectif
bool Algori::CalculRemont(LesMaillages * lesMail,OrdreVisu::EnumTypeIncre type_incre,int incre)
   {  bool change=false;
      bool remont_contrainte=false;
      if (pa.SauvegardeAutorisee(incre,temps_derniere_sauvegarde,type_incre))
           { if (avec_remont_sigma) {Algori::RemontSigma(lesMail);change=true;remont_contrainte=true;};
             if (avec_remont_eps) {Algori::RemontEps(lesMail);change=true;};
             if (avec_erreur) {Algori::RemontErreur(lesMail);change=true;};
            }; 
      return change;
    };  

// affichage éventuelle de la matrice de raideur et du second membre
void Algori::Affiche_RaidSM(const Vecteur & vglobin,const Mat_abstraite& matglob) const
 {      // affichage éventuelle de la matrice de raideur et du second membre     
       if (ParaGlob::NiveauImpression() >= 10)
        { string entete = " affichage de la matrice de raideur  ";
          matglob.Affichage_ecran(entete); 
          entete = " affichage du second membre  ";
          vglobin.Affichage_ecran(entete); 
         }                        
   };
   
// algorithme classique de line search
// le line seach agit sur une partie de la puissance et non sur la norme du résidu
// le line search à pour objectif de minimiser une partie de l'acroissement de la puissance
// en jeux, ainsi on peut avoir une minimisation réussi tout en ayant une norme d'équilibre
// qui continue à grandir
// cependant en général les deux sont lié et l'application du line_seach garantit que dans
// la direction de descente choisi on minimise norme (mais ce minimum peut cependant être
// supérieur à l'ancienne norme!)
bool  Algori::Line_search1(Vecteur& sauve_deltadept,double& puis_precedente,
              Vecteur& vres,LesMaillages * lesMail,Vecteur* sol 
              ,int& compteur,Vecteur& sauve_dept_a_tdt,Charge* charge,Vecteur& vglobex
              ,Assemblage& Ass,Vecteur& v_travail,LesCondLim* lesCondLim,Vecteur& vglobaal
              ,LesReferences* lesRef,Vecteur & vglobin,const Nb_assemb& nb_casAssemb
              ,int cas_combi_ddl,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD)

 {  // ------------------------------------------------------------
    // |on regarde s'il faut utiliser l'algorithme de line search |
    // ------------------------------------------------------------
    // nombre d'incrément mini à partir duquel on exécute le line_search
    const int nb_incr_min_line_search = 2;
    // les mini et maxi du facteur de line _search 
    const double alphmin = 0.1; const double alphmax = 1.5; // les mini et maxi permit
    // le nombre de fois que l'algorithme bute sur le mini et maxi
    int nbfois_alphmin = 0; int nbfois_alphmax = 0;
    const int nbfois_alphmin_max = 2; const int nbfois_alphmax_max = 2;
    // nombre maxi de boucle de line_search
    const int iwwmax = 5;
    // limite sur la puissance pour activer le line_search
    const double epsilon_p = 0.01; 
    double ancien_residu_norme = vres.Max_val_abs ();
    double norme_dep = sauve_deltadept.Max_val_abs ();
    double ancien_residuN = (- vres  * (*sol));
    // sauvegarde du déplacement
    sauve_deltadept = *sol;
    // il est nécessaire de calculer le nouveau résidu pour conclure sur l'utilité de mettre
    // en oeuvre l'algorithme 
      // 1 - on sauvegarde les conditions actuelles
      // récupération de l'incrément total de t à tdt avant le déplacement
    lesMail->RecupDepde_tatdt(sauve_dept_a_tdt,nb_casAssemb);
      // 2 - on met en place le nouveau déplacement
    lesMail->PlusDelta_tdt(*sol,nb_casAssemb);
      // 3 - initialisation de la puissance externe puis on la calcul
    vglobex.Zero();
    if (!(charge->ChargeSecondMembre_Ex_mecaSolid
                 (Ass,lesMail,lesRef,vglobex,pa,lesCourbes1D,lesFonctionsnD)))
       { cout << "\n erreur dans l'algo, du a une erreur sur le calcul "
              << " des efforts externes: arret"
              << "\n Algori::Line_search1(...";
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsCalEquilibre.Comptage_en_cours())
                    tempsCalEquilibre.Arret_du_comptage();
         Sortie(1);
       };
     // 4 - idem pour la puissance interne
     vglobin.Zero();
      // 5 - calcul du second membre et des énergies
    // dans le cas d'un calcul inexploitable arrêt de la boucle
    if (!SecondMembreEnerg(lesMail,Ass,vglobin))
    	{ cout << "\n erreur dans l'algo, du a une erreur sur le calcul du residu: arret"
    	       << "\n Algori::Line_search1(...";
       // dans le cas où un comptage du calcul est en cours on l'arrête
       if (tempsCalEquilibre.Comptage_en_cours())
                  tempsCalEquilibre.Arret_du_comptage();
    	  Sortie(1);
    	};
      // 6 - second membre total
    vglobaal = vglobex ;vglobaal += vglobin;
      // 7 - mise en place des conditions limites, vglobal contiend le nouveau résidu
    lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobaal,nb_casAssemb,cas_combi_ddl,vglob_stat);
    double nouvelle_norme = vglobaal.Max_val_abs ();
    
    // calcul d'une partie de la variation de puissance relative au delta ddl
    // cette grandeur sert à mesurer l'écart de normalité entre la direction de delta ddl
    // et la direction du résidu 
    double var_puis  =  - vglobaal * sauve_deltadept; // calcul de  la puissance
    cout << "\n ecart de permendicularité entre deltat ddl et residu = " << var_puis;
    cout << "\n puissance  " << var_puis ;
    cout << "\n  ||dep|| * ||res|| " << ( norme_dep * nouvelle_norme);

    if ((compteur >= nb_incr_min_line_search)  
         && ((Abs(var_puis) >  epsilon_p * norme_dep * nouvelle_norme) 
         || (nouvelle_norme > 0.1 * ancien_residu_norme))) 
 //      ((var_puis >  puis_precedente) || (nouvelle_norme > ancien_residu_norme))) 
     //       || (puis_reelle < 0.))
        {// signifie que le déplacement précédent était trop grand
         // ou pas dans la bonne direction ! 
         // 1) on revient aux conditions du début de l'itération
         v_travail -= sauve_dept_a_tdt;
         lesMail->PlusDelta_tdt(v_travail,nb_casAssemb);
         // 2) on traite la direction
         // si la puissance est négative cela signifie que la matrice n'est pas définie
         // positive, dans ce cas on utilise comme direction de descente 
         // - la direction trouvée ou le résidu
     /*    if (puis_reelle < 0.)
             { // message d'avertissement
               cout << "\n ********** attention puissance négative, changement de descente "; 
             //  ancien_dep = - vres;
               ancien_dep = - ancien_dep;
              } */
    /*        else if (Abs(puis_reelle) > pa.max_puissance)
               // si la puissance est trop grande arrêt des itérations
              {
                non_convergence = true;
                break;
               } */
          // 3) Maintenant application du line search, 
          //calcul d'alpha en utilisant une approximation linéaire de la puissance
          double alph = 1.; // initialisation du coef multiplicateur du delta ddl
          // dans le cas ou la nouvelle et l'ancienne puissance sont très proche
          // on change de alpha
          if (Abs(var_puis- ancien_residuN) < epsilon_p * norme_dep * nouvelle_norme)
             alph = 0.66; // par exemple
          // dans un premier temps on effectue un nombre de  boucle fixé pour alpha
          bool reduc_reussi = false; 
          bool reduc_possible = true;
          for (int iww = 1;(iww<=iwwmax && reduc_possible) ;iww++)
             { // calcul du nouveau résidu
               // mise en place du nouveau déplacement
               // pour éviter de créer un nouveau vecteur qui peut être de grande taille
               // on utilise la place de vglobex pour le déplacement   
               Vecteur& nouveau_dep = vglobex; // place pour un nouveau déplacement
               nouveau_dep = sauve_dept_a_tdt + alph * sauve_deltadept;
               lesMail->TversTdt();
               lesMail->PlusDelta_tdt(nouveau_dep,nb_casAssemb );
               //puissance externe
               vglobex.Zero();
              
               if (!(charge->ChargeSecondMembre_Ex_mecaSolid
                             (Ass,lesMail,lesRef,vglobex,pa,lesCourbes1D,lesFonctionsnD)))
                   { cout << "\n erreur dans l'algo, du a une erreur sur le calcul "
                          << " des efforts externes: arret"
                          << "\n Algori::Line_search1(...";
                     // dans le cas où un comptage du calcul est en cours on l'arrête
                     if (tempsCalEquilibre.Comptage_en_cours())
                                tempsCalEquilibre.Arret_du_comptage();
                     Sortie(1);
                   };
               // puissance interne boucle sur les elements
               vglobin.Zero();
               // calcul du second membre et des énergies
               // dans le cas d'un calcul inexploitable arrêt de la boucle
               if (!SecondMembreEnerg(lesMail,Ass,vglobin))
                 { cout << "\n erreur2 dans l'algo, du a une erreur sur le calcul du residu: arret"
                        << "\n Algori::Line_search1(...";
                   // dans le cas où un comptage du calcul est en cours on l'arrête
                   if (tempsCalEquilibre.Comptage_en_cours())
                     tempsCalEquilibre.Arret_du_comptage();
                   Sortie(1);
                 };
               // second membre total
               vglobaal = vglobex ;vglobaal += vglobin;
               // mise en place des conditions limites
               lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobaal,nb_casAssemb,cas_combi_ddl,vglob_stat);
               // calcul du nouveau coefficient
   //            double nouveau_residuN = vglobal.Max_val_abs ();//alph * (- vglobal * ancien_dep);
   //            double ancien_residuN = ancien_residu.Max_val_abs ();//alph * (- ancien_residu  * ancien_dep);//vres.Norme();
               double nouveau_residuN = alph * (- vglobaal * sauve_deltadept);
               double ancien_alph = alph;
               nouvelle_norme = vglobaal.Max_val_abs ();
               alph = - alph* ancien_residuN/
                       (nouveau_residuN - ancien_residuN);
               // test d'arrêt 
               if ((alph < alphmin) || (alph > alphmax))
                { // troncature d'alpha
                  if (alph < alphmin) 
                   {alph = alphmin;nbfois_alphmin++;
                    if (nbfois_alphmin > nbfois_alphmin_max)
          /*           if (nouvelle_norme < ancien_residu_norme)
                       // soit la norme à diminuée => on prend le nouveau déplacement
                       // et on arrête car dans tous les cas ce n'est pas une bonne convergence
                       // d'être acculé à la borne 
                      { reduc_reussi = true;
                        // on remet au mieux les choses comme au début de la méthode
                        *sol = ancien_alph * sauve_deltadept;
                        vres = vglobal;
                        puis_precedente = Abs(nouveau_residuN);
                        break; 
                       }
                     else */
                     // sinon il n'y a pas en plus de convergence on arrête tout   
                       reduc_possible = false;
                    }
                  if (alph > alphmax) 
                   {alph = alphmax;nbfois_alphmax++;
                    if (nbfois_alphmax > nbfois_alphmax_max)
           /*          if (nouvelle_norme < ancien_residu_norme)
                       // soit la norme à diminuée => on prend le nouveau déplacement
                       // et on arrête car dans tous les cas ce n'est pas une bonne convergence
                       // d'être acculé à la borne 
                      { reduc_reussi = true;
                        // on remet au mieux les choses comme au début de la méthode
                        *sol = ancien_alph * sauve_deltadept;
                        vres = vglobal;
                        puis_precedente = Abs(nouveau_residuN);
                        break; 
                       }
                     else */
                     // sinon il n'y a pas en plus de convergence on arrête tout   
                        reduc_possible = false;
                    }
                 } 
               else // sinon on test la convergence
      //           if (nouvelle_norme <  0.8 * ancien_residu_norme )  
                {if ((Abs(nouveau_residuN) <  0.5 * Abs(puis_precedente) ) 
                    && (nouvelle_norme < ancien_residu_norme)) 
                 { reduc_reussi = true;
                   // on remet au mieux les choses comme au début de la méthode
                   *sol = ancien_alph * sauve_deltadept;
                   vres = vglobaal;
                   puis_precedente = Abs(nouveau_residuN);
                   break; 
                   }
      //           else if (iww==iwwmax) 
                 else if ((iww==iwwmax) && (nouvelle_norme < ancien_residu_norme))
                  {reduc_reussi = true;
                   // on remet au mieux les choses comme au début de la méthode
                   *sol = ancien_alph * sauve_deltadept;
                   vres = vglobaal;
                   puis_precedente =  Abs(nouveau_residuN);
                   break; 
                   }                  
                 }  
               // sinon on incrémente la position à tdt      
               cout << "\n coefficient de line_search = " << alph
                    << " puiss : " << nouveau_residuN << " old_puiss " << ancien_residuN;            
              }
            // si la réduction c'est mal passée alors on laisse la solution
            // identique à avant le line_search
            if (!reduc_reussi)
             {  // mise en place du  déplacement existant avant le line_search
                // ainsi que les autres grandeurs actives
                lesMail->TversTdt();
                v_travail = sauve_dept_a_tdt + sauve_deltadept;
                lesMail->PlusDelta_tdt( v_travail,nb_casAssemb);
                // retour du fait qu'il n'y a pas eu de line search
                return false;
              } 
             else
             {
              // retour du fait que le line search a été effectué
              return true;
              }  
            }
         else // cas ou l'on ne fait pas de line seach
             return false;
   };
                           
    // algorithme de line search qui utilise une approximation cubique du résidu ( (R).<ddl>)
    // en fonction du paramètre lambda de line search ( X+lambda. delta_ddl)
    // la méthode renvoi si oui ou non le line search a été effecué
    // sauve_deltadept : vecteur de travail servant à sauver la valeur calculée de sol
    // vres : correspond en entrée au vecteur résidu
    // sauve_dept_a_tdt : vecteur de travail utilisé pour sauvegarder la totalité de l'incrément de t à t+dt
    // vglobex : puissance externe
    // v_travail : vecteur de travail à la même dimension que sol 
    // vglobal : = & vglobin: vecteur puissance globale
bool  Algori::Line_search2(Vecteur& sauve_deltadept,double& puis_precedente,
              Vecteur& vres,LesMaillages * lesMail,Vecteur* sol 
              ,int& compteur,Vecteur& sauve_dept_a_tdt,Charge* charge,Vecteur& vglobex
              ,Assemblage& Ass,Vecteur& v_travail,LesCondLim* lesCondLim,Vecteur& vglobaal
              ,LesReferences* lesRef,Vecteur & vglobin,const Nb_assemb& nb_casAssemb
              ,int cas_combi_ddl,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD)

 { 
 
 /*  // initialisation des différentes grandeurs nécessaires pour le calcul du polynome cubique complet ou dégénéré
   double f1=      ;// puissance pour lambda = 0 
   double f2=      ;// puissance pour lambda = 1
   double g1=      ;//  */
 
 
    // ------------------------------------------------------------
    // |on regarde s'il faut utiliser l'algorithme de line search |
    // ------------------------------------------------------------
    // nombre d'incrément mini à partir duquel on exécute le line_search
    const int nb_incr_min_line_search = 2;
    // les mini et maxi du facteur de line _search 
    const double alphmin = 0.1; const double alphmax = 1.5; // les mini et maxi permit
    // le nombre de fois que l'algorithme bute sur le mini et maxi
    int nbfois_alphmin = 0; int nbfois_alphmax = 0;
    const int nbfois_alphmin_max = 2; const int nbfois_alphmax_max = 2;
    // nombre maxi de boucle de line_search
    const int iwwmax = 5;
    // limite sur la puissance pour activer le line_search
    const double epsilon_p = 0.01; 
    double ancien_residu_norme = vres.Max_val_abs ();
    double norme_dep = sauve_deltadept.Max_val_abs ();
    double ancien_residuN = (- vres  * (*sol));
    // sauvegarde du déplacement
    sauve_deltadept = *sol;
    // il est nécessaire de calculer le nouveau résidu pour conclure sur l'utilité de mettre
    // en oeuvre l'algorithme 
      // 1 - on sauvegarde les conditions actuelles
      // récupération de l'incrément total de t à tdt avant le déplacement
    lesMail->RecupDepde_tatdt(sauve_dept_a_tdt,nb_casAssemb);
      // 2 - on met en place le nouveau déplacement
    lesMail->PlusDelta_tdt(*sol,nb_casAssemb);
      // 3 - initialisation de la puissance externe puis on la calcul
    vglobex.Zero();
  
    if (!(charge->ChargeSecondMembre_Ex_mecaSolid
                    (Ass,lesMail,lesRef,vglobex,pa,lesCourbes1D,lesFonctionsnD)))
          { cout << "\n erreur dans l'algo, du a une erreur sur le calcul "
                 << " des efforts externes: arret"
                 << "\n Algori::Line_search2(...";
            // dans le cas où un comptage du calcul est en cours on l'arrête
            if (tempsCalEquilibre.Comptage_en_cours())
                       tempsCalEquilibre.Arret_du_comptage();
            Sortie(1);
          };
      // 4 - idem pour la puissance interne
    vglobin.Zero();
      // 5 - calcul du second membre et des énergies
    // dans le cas d'un calcul inexploitable arrêt de la boucle
    if (!SecondMembreEnerg(lesMail,Ass,vglobin))
       { cout << "\n erreur dans l'algo, du a une erreur sur le calcul du residu: arret"
              << "\n Algori::Line_search2(...";
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsCalEquilibre.Comptage_en_cours())
            tempsCalEquilibre.Arret_du_comptage();
         Sortie(1);
       };
      // 6 - second membre total ( vglobal est identique a vglobin)
    vglobaal = vglobex ; vglobaal += vglobin;
      // 7 - mise en place des conditions limites, vglobal contiend le nouveau résidu
    lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobaal,nb_casAssemb,cas_combi_ddl,vglob_stat);
    double nouvelle_norme = vglobaal.Max_val_abs ();
    
    // calcul d'une partie de la variation de puissance relative au delta ddl
    // cette grandeur sert à mesurer l'écart de normalité entre la direction de delta ddl
    // et la direction du résidu 
    double var_puis  =  - vglobaal * sauve_deltadept; // calcul de  la puissance
    cout << "\n ecart de permendicularité entre deltat ddl et residu = " << var_puis;
    cout << "\n puissance  " << var_puis ;
    cout << "\n  ||dep|| * ||res|| " << ( norme_dep * nouvelle_norme);

    if ((compteur >= nb_incr_min_line_search)  
         && ((Abs(var_puis) >  epsilon_p * norme_dep * nouvelle_norme) 
         || (nouvelle_norme > 0.1 * ancien_residu_norme))) 
 //      ((var_puis >  puis_precedente) || (nouvelle_norme > ancien_residu_norme))) 
     //       || (puis_reelle < 0.))
        {// signifie que le déplacement précédent était trop grand
         // ou pas dans la bonne direction ! 
         // 1) on revient aux conditions du début de l'itération
         v_travail -= sauve_dept_a_tdt;
         lesMail->PlusDelta_tdt(v_travail,nb_casAssemb);
         // 2) on traite la direction
         // si la puissance est négative cela signifie que la matrice n'est pas définie
         // positive, dans ce cas on utilise comme direction de descente 
         // - la direction trouvée ou le résidu
     /*    if (puis_reelle < 0.)
             { // message d'avertissement
               cout << "\n ********** attention puissance négative, changement de descente "; 
             //  ancien_dep = - vres;
               ancien_dep = - ancien_dep;
              } */
    /*        else if (Abs(puis_reelle) > pa.max_puissance)
               // si la puissance est trop grande arrêt des itérations
              {
                non_convergence = true;
                break;
               } */
          // 3) Maintenant application du line search, 
          //calcul d'alpha en utilisant une approximation linéaire de la puissance
          double alph = 1.; // initialisation du coef multiplicateur du delta ddl
          // dans le cas ou la nouvelle et l'ancienne puissance sont très proche
          // on change de alpha
          if (Abs(var_puis- ancien_residuN) < epsilon_p * norme_dep * nouvelle_norme)
             alph = 0.66; // par exemple
          // dans un premier temps on effectue un nombre de  boucle fixé pour alpha
          bool reduc_reussi = false; 
          bool reduc_possible = true;
          for (int iww = 1;(iww<=iwwmax && reduc_possible) ;iww++)
             { // calcul du nouveau résidu
               // mise en place du nouveau déplacement
               // pour éviter de créer un nouveau vecteur qui peut être de grande taille
               // on utilise la place de vglobex pour le déplacement   
               Vecteur& nouveau_dep = vglobex; // place pour un nouveau déplacement
               nouveau_dep = sauve_dept_a_tdt + alph * sauve_deltadept;
               lesMail->TversTdt();
               lesMail->PlusDelta_tdt(nouveau_dep,nb_casAssemb );
               //puissance externe
               vglobex.Zero();
              
               if (!(charge->ChargeSecondMembre_Ex_mecaSolid
                             (Ass,lesMail,lesRef,vglobex,pa,lesCourbes1D,lesFonctionsnD)))
                   { cout << "\n erreur dans l'algo, du a une erreur sur le calcul "
                          << " des efforts externes: arret"
                          << "\n Algori::Line_search2(...";
                     // dans le cas où un comptage du calcul est en cours on l'arrête
                     if (tempsCalEquilibre.Comptage_en_cours())
                                tempsCalEquilibre.Arret_du_comptage();
                     Sortie(1);
                   };
               // puissance interne boucle sur les elements
               vglobin.Zero();
               // calcul du second membre et des énergies
               // dans le cas d'un calcul inexploitable arrêt de la boucle
               if (!SecondMembreEnerg(lesMail,Ass,vglobin))
    	            { cout << "\n erreur2 dans l'algo, du a une erreur sur le calcul du residu: arret"
    	                   << "\n Algori::Line_search1(...";
                   // dans le cas où un comptage du calcul est en cours on l'arrête
                   if (tempsCalEquilibre.Comptage_en_cours())
                       tempsCalEquilibre.Arret_du_comptage();
    	              Sortie(1);
    	            };
               // second membre total ( vglobal est identique a vglobin)
               vglobaal = vglobex ;vglobaal += vglobin;
               // mise en place des conditions limites
               lesCondLim->ImposeConLimtdt(lesMail,lesRef,vglobaal,nb_casAssemb,cas_combi_ddl,vglob_stat);
               // calcul du nouveau coefficient
   //            double nouveau_residuN = vglobal.Max_val_abs ();//alph * (- vglobal * ancien_dep);
   //            double ancien_residuN = ancien_residu.Max_val_abs ();//alph * (- ancien_residu  * ancien_dep);//vres.Norme();
               double nouveau_residuN = alph * (- vglobaal * sauve_deltadept);
               double ancien_alph = alph;
               nouvelle_norme = vglobaal.Max_val_abs ();
               alph = - alph* ancien_residuN/
                       (nouveau_residuN - ancien_residuN);
               // test d'arrêt 
               if ((alph < alphmin) || (alph > alphmax))
                { // troncature d'alpha
                  if (alph < alphmin) 
                   {alph = alphmin;nbfois_alphmin++;
                    if (nbfois_alphmin > nbfois_alphmin_max)
          /*           if (nouvelle_norme < ancien_residu_norme)
                       // soit la norme à diminuée => on prend le nouveau déplacement
                       // et on arrête car dans tous les cas ce n'est pas une bonne convergence
                       // d'être acculé à la borne 
                      { reduc_reussi = true;
                        // on remet au mieux les choses comme au début de la méthode
                        *sol = ancien_alph * sauve_deltadept;
                        vres = vglobal;
                        puis_precedente = Abs(nouveau_residuN);
                        break; 
                       }
                     else */
                     // sinon il n'y a pas en plus de convergence on arrête tout   
                       reduc_possible = false;
                    }
                  if (alph > alphmax) 
                   {alph = alphmax;nbfois_alphmax++;
                    if (nbfois_alphmax > nbfois_alphmax_max)
           /*          if (nouvelle_norme < ancien_residu_norme)
                       // soit la norme à diminuée => on prend le nouveau déplacement
                       // et on arrête car dans tous les cas ce n'est pas une bonne convergence
                       // d'être acculé à la borne 
                      { reduc_reussi = true;
                        // on remet au mieux les choses comme au début de la méthode
                        *sol = ancien_alph * sauve_deltadept;
                        vres = vglobal;
                        puis_precedente = Abs(nouveau_residuN);
                        break; 
                       }
                     else */
                     // sinon il n'y a pas en plus de convergence on arrête tout   
                        reduc_possible = false;
                    }
                 } 
               else // sinon on test la convergence
      //           if (nouvelle_norme <  0.8 * ancien_residu_norme )  
                {if ((Abs(nouveau_residuN) <  0.5 * Abs(puis_precedente) ) 
                    && (nouvelle_norme < ancien_residu_norme)) 
                 { reduc_reussi = true;
                   // on remet au mieux les choses comme au début de la méthode
                   *sol = ancien_alph * sauve_deltadept;
                   vres = vglobaal;
                   puis_precedente = Abs(nouveau_residuN);
                   break; 
                   }
      //           else if (iww==iwwmax) 
                 else if ((iww==iwwmax) && (nouvelle_norme < ancien_residu_norme))
                  {reduc_reussi = true;
                   // on remet au mieux les choses comme au début de la méthode
                   *sol = ancien_alph * sauve_deltadept;
                   vres = vglobaal;
                   puis_precedente =  Abs(nouveau_residuN);
                   break; 
                   }                  
                 }  
               // sinon on incrémente la position à tdt      
               cout << "\n coefficient de line_search = " << alph
                    << " puiss : " << nouveau_residuN << " old_puiss " << ancien_residuN;            
              }
            // si la réduction c'est mal passée alors on laisse la solution
            // identique à avant le line_search
            if (!reduc_reussi)
             {  // mise en place du  déplacement existant avant le line_search
                // ainsi que les autres grandeurs actives
                lesMail->TversTdt();
                v_travail = sauve_dept_a_tdt + sauve_deltadept;
                lesMail->PlusDelta_tdt( v_travail,nb_casAssemb);
                // retour du fait qu'il n'y a pas eu de line search
                return false;
              } 
             else
             {
              // retour du fait que le line search a été effectué
              return true;
              }  
            }
         else // cas ou l'on ne fait pas de line seach
             return false;
   };
   
// initialisation des paramètres pour l'amortissement cinétique  (au cas ou on veut plusieurs amortissement)
void Algori::InitialiseAmortissementCinetique()
	{ compteur_decroit_pourRelaxDyn=0;
   dernier_pic = pic_E_cint_t = -1.; 
   max_pic_E_cin = -1.;
   compteur_test_pourRelaxDyn = 0;
   compteur_pic_energie = 0;compteur_pic_energie_noe= 0;
		 // cas de l'utilisation d'une moyenne glissante
		 sauve_E_cin.Change_taille(taille_moyenne_glissante);
		 nb_Ecin_stocker = 0;indice_EC1= 0;
		 moyenne_glissante = 0.;
		 moy_gliss_t = 0.;
	};                         
       
// méthode d'application de l'amortissement cinétique 
// on arrête l'amortissement lorsque E_cin_tdt < E_cin_t * pic_E_cint_t
// on redémarre l'amortissement lorsque  E_cin_tdt > max_pic_E_cin * coef_redemarrage_pourRelaxDyn
    // retour un int :
    //   = 1  : il y a eu amortissement mais le calcul continue
    //   = 0  : pas d'amortissement
    //   = -1 : on peut arrêter le calcul car le critère sur delta_ddl est satisfait
int Algori::AmortissementCinetique(const Vecteur & delta_ddl,const double& coef_mass,const Vecteur& X
                               ,const Mat_abstraite& mat_mass,int icharge,Vecteur& V)
	{int relax_vit = 0;
  nb_depuis_derniere_relax++; // on incrémente à chaque passage
  // dans le cas où inter_nb_entre_relax dépend d'une fonction nD, on calcule sa valeur
  if (fct_nD_inter_nb_entre_relax != NULL)
   {int old_valeur = inter_nb_entre_relax;
    inter_nb_entre_relax = (fct_nD_inter_nb_entre_relax->Valeur_pour_variables_globales())(1);
    if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 2)) || (permet_affichage > 1))
      if (old_valeur != inter_nb_entre_relax)
      cout << "\n >>> changement de inter_nb_entre_relax de " << old_valeur
           << " en ==> "<< inter_nb_entre_relax;
   };
	 
//	 // dans le cas où l'on veut un amortissement individuel aux noeuds, on appelle la méthode spécifique
//  int relax_vit_individuel = 0;
//	 if (amortissement_cinetique_au_noeud && !relax_vit)
//	   { relax_vit_individuel=AmortissementCinetique_individuel_aux_noeuds(delta_ddl,coef_mass,X,mat_mass,icharge,V);};

//*** je ne comprends pas le test que j'ai mis ici: abscon !! 16 oct 2018
  // on ne continue que s'il n'y a pas d'amortissement cinétique en route
  if (!((icharge > nb_deb_test_amort_cinetique_noe) && amortissement_cinetique_au_noeud))
// ---- non la suite n'est pas bonne et converge beaucoup moins bien !!
// je le remplace : on test si on est en amortissement cinétique, et que le nombre d'itérations
// est suffisant
///  if ((icharge > nb_deb_test_amort_cinetique_noe) && (!amortissement_cinetique_au_noeud)
///      && amortissement_cinetique )
// --- fin modif avorté !!!
    {

//  if (!relax_vit_individuel)
////  else
//    {
//	 cout << "\n E_cin_t= " << E_cin_t << " E_cin_tdt= " << E_cin_tdt
//	      << " compteur_decroit_pourRelaxDyn= " << compteur_decroit_pourRelaxDyn;

// remarque: 
// compteur_decroit_pourRelaxDyn: indique combien de fois le test de pic a été effectué avec succés 
//       si == -2 : indique que la procédure de test est gelé
//       si == 0  : indique que la procédure est active, et qu'il n'y a pas eu pour l'instant de pic détecté
//       si > 0   : indique que la procédure est active, et que, de manière consécutive, indique 
//                  le nombre de fois où il y a eu baisse               

      bool relaxation_effectuee = false; // initialisation
      bool relaxation_reactive = false;  // "
      if ((amortissement_cinetique)&&(icharge != 0) && (icharge > nb_deb_test_amort_cinetique))
       { // on regarde si l'énergie cinétique a déjà été calculée, sinon on la calcule
         if (icharge < pa.NbIncrCalEnergie())  // ******* a revoir, car le test n'est peut-être pas suffisant avec tous les algos
             { // initialisation 
               v_int.Change_taille(V.Taille()); v_int.Zero(); 
               // calcul de l'énergie cinétique
               E_cin_tdt = (0.5/coef_mass) * (V * (mat_mass.Prod_vec_mat(V,v_int)));
             };
         // on calcul éventuellement la moyenne glissante
         if (taille_moyenne_glissante > 1)
          { if (nb_Ecin_stocker < taille_moyenne_glissante)
             { // c'est le cas où la moyenne n'est pas finalisée
               nb_Ecin_stocker++; indice_EC1=1;
               sauve_E_cin(nb_Ecin_stocker) = E_cin_tdt;
             }
            else // cas ou la moyenne est finalisée	
             { sauve_E_cin(indice_EC1) = E_cin_tdt;
               indice_EC1++; // mise à jour de l'indice pour le prochain stockage
               // si l'indice dépasse la taille du taille, on remet au début
               if (indice_EC1 > taille_moyenne_glissante) indice_EC1 = 1;
             }; 
            // maintenant on regarde de nouveau si on peut calculer la moyenne
            if (nb_Ecin_stocker == taille_moyenne_glissante)
              {moy_gliss_t = moyenne_glissante; moyenne_glissante = sauve_E_cin.Moyenne();}
            else 	{moy_gliss_t= 0.; moyenne_glissante = 0.;}; // au cas où on met une valeur très grande
          };
         // écriture d'info pour le débug au cas où
         if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6)) || (permet_affichage > 3))
            {if (taille_moyenne_glissante == 1)
               {cout << "\n E_cin_tdt= " << E_cin_tdt << " compteur_decroit_pourRelaxDyn= " << compteur_decroit_pourRelaxDyn
                     << " pic_E_cint_t= " << pic_E_cint_t << endl;
               }
             else
              {cout << "\n E_cin_tdt= " << E_cin_tdt << " compteur_decroit_pourRelaxDyn= " << compteur_decroit_pourRelaxDyn 
                    << "  moyenne_glissante= " << moyenne_glissante	<< " moy_gliss_t= " << moy_gliss_t
                    << " pic_E_cint_t= " << pic_E_cint_t << endl;
              };
            };	 	  
         // --- application de l'algo
         // a) dans le cas où la procédure a été gelée précédemment, on regarde 
         // s'il faut la redémarrer
         if (compteur_decroit_pourRelaxDyn == -2)
           {if (E_cin_tdt > max_pic_E_cin * coef_redemarrage_pourRelaxDyn)
              {compteur_decroit_pourRelaxDyn = 0; // indique que l'amortissement est de nouveau active
               relaxation_reactive = true;
              };
           };
         // b) on regarde pour l'application de l'amortissement
         if (compteur_decroit_pourRelaxDyn != -2) // != -2 veut dire que la procédure est active
           {if (taille_moyenne_glissante == 1)
             { 
           
    //debug ------
    //cout << "\n E_cin_tdt " << E_cin_tdt << " , E_cin_t " << E_cin_t << " pic_E_cint_t " << pic_E_cint_t << endl; 
    // fin debug					  
              if ((E_cin_tdt > E_cin_t) && (E_cin_tdt > pic_E_cint_t)) //pic_E_cint_t)
               { 
           
    //debug ------
    //cout << "\n marquage d'un nouveau pic: E_cin_tdt " << E_cin_tdt << " , E_cin_t " << E_cin_t << " pic_E_cint_t " << pic_E_cint_t << endl; 
    // fin debug					  
           
                 pic_E_cint_t = E_cin_tdt;
               } //X_EcinMax = X; 
              max_pic_E_cin=MaX(pic_E_cint_t,max_pic_E_cin);
             }
            else 
             { if ((moyenne_glissante > moy_gliss_t) && (moyenne_glissante > pic_E_cint_t))//pic_E_cint_t)
                { pic_E_cint_t = moyenne_glissante;
                };// X_EcinMax = X; }
               max_pic_E_cin=MaX(pic_E_cint_t,max_pic_E_cin);
             };

            // on considère qu'il y a décroissance, lorsque soit:
            // 1) on n'est pas en moyenne glissante et on a une décroissance de l'énergie cinétique
            //    ou (si fraction != 0), l'énergie cinétique est inférieure à fraction * dernier_max_E_cin
            // 2) on est en moyenne glissante et on a une décroissance de la moyenne glissante
            //    ou (si fraction != 0), la moyenne glissante est inférieure à fraction * dernier_max_E_cin
            // NB: si le pic_E_cint_t a été updater dans le bloc précédent, le test qui suit ne peut pas être valide
            if (  ((test_fraction_energie == 0.) &&
                     (   ((taille_moyenne_glissante == 1) && (E_cin_tdt < pic_E_cint_t)) // cas sans moyenne glissante
                  || ((taille_moyenne_glissante != 1) && (moyenne_glissante < pic_E_cint_t)) // cas avec moyenne glissante
                          )
                   )
                  ||
                     ((test_fraction_energie != 0.) && // cas où on se réfère à une fraction de l'énergie maxi précédente
                       (   ((taille_moyenne_glissante == 1) && (E_cin_tdt < (test_fraction_energie * pic_E_cint_t))) // cas sans moyenne glissante
                  || ((taille_moyenne_glissante != 1) && (moyenne_glissante < (test_fraction_energie * pic_E_cint_t))) // cas avec moyenne glissante
                       )
                     )
              )
             { 
    //debug ------
    //cout << "\n descente: E_cin_tdt " << E_cin_tdt  << " pic_E_cint_t " << pic_E_cint_t << endl; 
    // fin debug					  
               compteur_decroit_pourRelaxDyn++;
               if ((compteur_decroit_pourRelaxDyn >= max_nb_decroit_pourRelaxDyn)
                   && (nb_depuis_derniere_relax >= inter_nb_entre_relax))
                { // cas ou on met en place l'amortissement
                  V.Zero();relaxation_effectuee=true;
                  nb_depuis_derniere_relax = 0;
                  dernier_pic = pic_E_cint_t; // on sauvegarde
                  if (taille_moyenne_glissante != 1) // on remet tout à 0
                    {nb_Ecin_stocker = 0;indice_EC1= 0;};
                  relax_vit=1;
                  compteur_pic_energie++;
                  compteur_decroit_pourRelaxDyn=0;
           
                  if ( dernier_pic == -1) 
    //					    { if (   ((taille_moyenne_glissante != 1) && ( moyenne_glissante < dernier_pic*coef_arret_pourRelaxDyn))
    //						       || ((taille_moyenne_glissante == 1) && ( E_cin_tdt < dernier_pic*coef_arret_pourRelaxDyn)) )
                    { pic_E_cint_t = 0.; // on remet à 0 le pic
                      E_cin_tdt = 0.;
                    }
                  else 	 
                    { if  ( dernier_pic < max_pic_E_cin*coef_arret_pourRelaxDyn)
    //	     	    if (E_cin_tdt <= pic_E_cint_t*coef_arret_pourRelaxDyn)
                 // cas où on suspend l'amortissement pour la suite
                        {compteur_decroit_pourRelaxDyn=-2;}; // inactive l'amortissement
                      pic_E_cint_t = 0.; // on remet à 0 le pic
                      E_cin_tdt = 0.;
                    };
                  E_cinNoe_tdt.Zero();
                  pic_E_cint_t_noe.Zero();
                };
             }
            else // sinon cela signifie que la diminution n'est pas validée donc 
              // on re-initialise le compteur de décroissance (qui sert à valider des décroissances consécutives)
              {compteur_decroit_pourRelaxDyn=0;};
            // c) on regarde s'il faut arrêter le calcul, relativement au delta_ddl, seulement à partir d'un
            // certain nombre de fois que la méthode à été activé
            // ce qui permet d'éviter des pbs au démarrage, pour lequel l'énergie cinétique démarre forcément
            // de 0 puis petit à petit grandie
            // de plus on ne test le delta_ddl que:
            // 1) si on a passé au moins "un" pic d'énergie cinétique (sinon on peut avoir un deltat ddl = 0 car on est presqu'au
            //     maxi et ce maxi = pas du tout l'équilibre !!)
    // non pour l'instant ...		// 2) si l'énergie cinétique est suffisament faible
            compteur_test_pourRelaxDyn++;
            if (compteur_test_pourRelaxDyn > nb_deb_testfin_pourRelaxDyn)
            { if ((max_pic_E_cin != -1)
                  || // on ajoute le cas très particulier où l'énergie cinétique reste nulle
                  ((E_cin_tdt ==0.) && (E_cin_t==0.)) // c'est le cas où rien ne bouge !!
                 )
              { if (delta_ddl.Max_val_abs() <= max_deltaX_pourRelaxDyn) nb_dX_OK_pourRelaxDyn++;
                if (nb_dX_OK_pourRelaxDyn > nb_max_dX_OK_pourRelaxDyn) 
                  { relax_vit = -1;
                    if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6)) || (permet_affichage > 3))
                      cout << "\n critere en deplacement ok: max_deltaX= " << delta_ddl.Max_val_abs();
                    nb_dX_OK_pourRelaxDyn=0; // on initialise pour une prochaine fois
                      // dans le cas ou l'arrêt réel n'est pas drivé par le déplacement 
                      // mais par exemple sur le résidu 
                  };
              };	
            };
           }; // -- fin du cas où l'on applique l'amortissement
       };
       
       
      // on tente d'afficher que si l'amortissement a effectivement eu lieu
    //	  if (relax_vit && pa.AfficheIncrEnergie())
      if (relaxation_effectuee)
        if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 2))
         {cout << "\n N:" << icharge << "  relaxation cinetique effectuee "
               << " pic_E_cint_t("<<compteur_pic_energie<<") " << dernier_pic 
               << " max_pic_E_cin= " << max_pic_E_cin;
            //--- info pour le  debug
    /*		  if (ParaGlob::NiveauImpression() > 8)
            {if (taille_moyenne_glissante == 1)
        {cout << "\n nb_Ecin_stocker= " << nb_Ecin_stocker << " indice_EC1= " << indice_EC1  << endl;
        }
       else 
        {cout << "\n nb_Ecin_stocker= " << nb_Ecin_stocker << " indice_EC1= " << indice_EC1 
          << "  moyenne_glissante= " << moyenne_glissante	<< " moy_gliss_t= " << moy_gliss_t << endl;
        };
       };*/
          if (compteur_decroit_pourRelaxDyn == -2) 
            cout << " relaxation_gelee " ;
          if (relaxation_reactive)
            cout << " relaxation_re_active ";
          cout << endl;
         };
    };
	 // dans le cas où l'on veut également un amortissement individuel aux noeuds, on appelle la méthode spécifique
  // que si il n'y a pas eu d'amortissement cinétique global
  int relax_vit_individuel = 0;
	 if (amortissement_cinetique_au_noeud && !relax_vit)
	   { relax_vit_individuel=AmortissementCinetique_individuel_aux_noeuds(delta_ddl,coef_mass,X,mat_mass,icharge,V);
      // ici on a relax_vit = 0 donc le retour est uniquement particularisé par relax_vit_individuel
      return relax_vit_individuel;
    }
  else // sinon on est dans le cas où on a seulement un amortissement global et seul relax_vit donne une info à ce sujet
    { return relax_vit;};
  
//  la ligne qui suit est fausse car avec un "ou" le résultat est uniquement 1 ou 0, et en particulier
//  les - sont supprimé !! ce qui pose un pb lorsque relax_vit = -1 !! 
//  return (relax_vit||relax_vit_individuel);
  
	};

// méthode interne d'application de l'amortissement cinétique, exploratoire
// cas d'un amortissement local
// retour un int :
    //   = 1  : il y a eu au moins un amortissement 
    //   = 0  : pas d'amortissement
int Algori::AmortissementCinetique_individuel_aux_noeuds(const Vecteur & delta_ddl,const double& coef_mass,const Vecteur& X
                               ,const Mat_abstraite& mat_mass,int icharge,Vecteur& V)
	{int relax_vit = 0;

	 if ((amortissement_cinetique_au_noeud)&&(icharge != 0) && (icharge > nb_deb_test_amort_cinetique_noe))
	  { // on va faire une boucle sur tous les noeuds
     int nbddl = delta_ddl.Taille(); int dim = ParaGlob::Dimension();
	    int NBN = nbddl / dim;
    
     // mise à jour des tailles
     if (E_cinNoe_tdt.Taille() == 0)
       { E_cinNoe_tdt.Change_taille(NBN,0.);E_cinNoe_t.Change_taille(NBN,0.);
         compteur_decroit_pourRelaxDyn_noe.Change_taille(NBN,0);
         dernier_pic_noe.Change_taille(NBN,-1);
         pic_E_cint_t_noe.Change_taille(NBN,-1);
         max_pic_E_cin_noe.Change_taille(NBN,-1);
       };

     // un compteur de relaxation pour info
     int compteur_relaxation_ok=0; // init
     
		   for (int noe =1; noe <= NBN; noe++)
  		  { // on calcul l'énergie pour le noeud
	       bool relaxation_effectuee = false; // initialisation
	       bool relaxation_reactive = false;  // "
        double& eN_cin_tdt = E_cinNoe_tdt(noe);
        double& eN_cin_t = E_cinNoe_t(noe);
        double& pic_eN_cint_t = pic_E_cint_t_noe(noe);
        double& dernier_pic_eN = dernier_pic_noe(noe);
        eN_cin_tdt = 0.; // init
        int ptint = (noe-1)*dim;
        switch (dim) { case 3: ptint++; eN_cin_tdt += 0.5 * mat_mass(ptint,ptint) * V(ptint)*V(ptint);
                       case 2: ptint++; eN_cin_tdt += 0.5 * mat_mass(ptint,ptint) * V(ptint)*V(ptint);
                       case 1: ptint++; eN_cin_tdt += 0.5 * mat_mass(ptint,ptint) * V(ptint)*V(ptint);
                     };
/*        switch (dim) { case 3: ptint++; eN_cin_tdt += 0.5  * V(ptint)*V(ptint);
                       case 2: ptint++; eN_cin_tdt += 0.5  * V(ptint)*V(ptint);
                       case 1: ptint++; eN_cin_tdt += 0.5  * V(ptint)*V(ptint);
                     };
*/
	       // --- application de l'algo
        // a) dans le cas où la procédure a été gelée précédemment, on regarde 
	       // s'il faut la redémarrer
        double & max_pic_eN_cin = max_pic_E_cin_noe(noe);
        int & compt_decroit_pourRelaxDyn_noe = compteur_decroit_pourRelaxDyn_noe(noe);
	       if (compt_decroit_pourRelaxDyn_noe == -2)
	        {if (eN_cin_tdt > max_pic_eN_cin * coef_redemarrage_pourRelaxDyn_noe)
	         {compt_decroit_pourRelaxDyn_noe = 0; // indique que l'amortissement est de nouveau active
				       relaxation_reactive = true;
				      };
	        };
	       // b) on regarde pour l'application de l'amortissement
  	     if (compt_decroit_pourRelaxDyn_noe != -2) // != -2 veut dire que la procédure est active
					    { 					  
////debug ------
//if (noe==220)
//{ cout << "\n V(i)= "<<V((noe-1)*dim+1)<<" "<<V((noe-1)*dim+2)<<" "<<V((noe-1)*dim+3)<<" "
//       << mat_mass((noe-1)*dim+1,(noe-1)*dim+1)<<" "<<mat_mass((noe-1)*dim+2,(noe-1)*dim+2)
//       <<" "<<mat_mass((noe-1)*dim+3,(noe-1)*dim+3)<<" ";
//  cout << "\n eN_cin_tdt " << eN_cin_tdt << " , eN_cin_t " << eN_cin_t << " pic_eN_cint_t " << pic_eN_cint_t << endl;
// }
// fin debug
					     if ((eN_cin_tdt > eN_cin_t) && (eN_cin_tdt > pic_eN_cint_t)) //pic_eN_cint_t)
					      { 
////debug ------
//if (noe==220)
//cout << "\n marquage d'un nouveau pic: eN_cin_tdt " << eN_cin_tdt << " , eN_cin_t " << eN_cin_t << " pic_eN_cint_t " << pic_eN_cint_t << endl;
//// fin debug					  
							      pic_eN_cint_t = eN_cin_tdt;
           }; //X_EcinMax = X;
						    max_pic_eN_cin=MaX(pic_eN_cint_t,max_pic_eN_cin);
         };

		      // on considère qu'il y a décroissance, lorsque soit:
		      // 1) on n'est pas en moyenne glissante et on a une décroissance de l'énergie cinétique
		      //    ou (si fraction != 0), l'énergie cinétique est inférieure à fraction * dernier_max_E_cin
		      // 2) on est en moyenne glissante et on a une décroissance de la moyenne glissante
		      //    ou (si fraction != 0), la moyenne glissante est inférieure à fraction * dernier_max_E_cin
			     // NB: si le pic_eN_cint_t a été updater dans le bloc précédent, le test qui suit ne peut pas être valide
			     if (  ((test_fraction_energie == 0.) && (eN_cin_tdt < pic_eN_cint_t))
            ||
              ((test_fraction_energie != 0.) && // cas où on se réfère à une fraction de l'énergie maxi précédente
		                   (eN_cin_tdt < (test_fraction_energie * pic_eN_cint_t)))
           )
		       { 
//debug ------
//cout << "\n descente: eN_cin_tdt " << eN_cin_tdt  << " pic_eN_cint_t " << pic_eN_cint_t << endl;
// fin debug					  
 		       compt_decroit_pourRelaxDyn_noe++;
          if (compt_decroit_pourRelaxDyn_noe >= max_nb_decroit_pourRelaxDyn_noe)
				       { // cas ou on met en place l'amortissement
             int depart = (noe-1)*dim+1;
             for (int i=0;i<dim;i++)
               V(depart+i)=0.;
//     V(depart+i)*=0.5; // test: on diminue la vitesse par 2 au lieu de la supprimer
////debug ------
//if (noe==220)
//cout << "\n amortissemen:eN_cin_tdt " << eN_cin_tdt << " pic_eN_cint_t " << pic_eN_cint_t << endl;
//// fin debug					  
             relaxation_effectuee=true;
					        dernier_pic_eN = pic_eN_cint_t; // on sauvegarde
					        relax_vit=1;
		           compteur_pic_energie_noe++;
					        compt_decroit_pourRelaxDyn_noe=0;
             compteur_relaxation_ok++; // pour l'affichage éventuelle
				   
					        if ( dernier_pic_eN == -1)
               { pic_eN_cint_t = 0.; // on remet à 0 le pic
							          eN_cin_tdt = 0.;
//							 eN_cin_tdt *= 0.25;  // test, comme la vitesse est diminuée de 2, ecin est diminuée de 4
               }
					        else 	 
					          { if  ( dernier_pic_eN < max_pic_eN_cin*coef_arret_pourRelaxDyn_noe)
					               {compt_decroit_pourRelaxDyn_noe=-2;}; // inactive l'amortissement
							          pic_eN_cint_t = 0.; // on remet à 0 le pic
							          eN_cin_tdt = 0.;
//							 eN_cin_tdt = 0.25; // test, comme la vitesse est diminuée de 2, ecin est diminuée de 4
               };
           };
         }
		      else // sinon cela signifie que la diminution n'est pas validée donc 
		       { // on re-initialise le compteur de décroissance (qui sert à valider des décroissances consécutives)
		         compt_decroit_pourRelaxDyn_noe=0;
         };
	  
////debug ------
//if (relaxation_effectuee && (noe==220))
//{ cout << "\n relaxation:V(i)= "<<V((noe-1)*dim+1)<<" "<<V((noe-1)*dim+2)<<" "<<V((noe-1)*dim+3)<<" "
//       << mat_mass((noe-1)*dim+1,(noe-1)*dim+1)<<" "<<mat_mass((noe-1)*dim+2,(noe-1)*dim+2)
//       <<" "<<mat_mass((noe-1)*dim+3,(noe-1)*dim+3)<<" ";
//  cout << "\n eN_cin_tdt " << eN_cin_tdt << " , eN_cin_t " << eN_cin_t << " pic_eN_cint_t " << pic_eN_cint_t << endl;
// }
//// fin debug					  
	       // on tente d'afficher que si l'amortissement a effectivement eu lieu
	       if (relaxation_effectuee)
         if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 7)) || (permet_affichage > 5))
	        {cout << "\n relax effectue au noeud " << noe
	              << " pic_eN_cint_t("<<compteur_pic_energie_noe<<") " << dernier_pic_eN
	              << " max_pic_eN_cin= " << max_pic_eN_cin;
	         if (compt_decroit_pourRelaxDyn_noe == -2)
	           cout << " relaxation_gelee " ;
		        if (relaxation_reactive)
		          cout << " relaxation_re_activee ";	 
	         cout << endl;
	        };
      }; // fin de la boucle sur les noeuds
     if (compteur_relaxation_ok>0) // affichage global éventuel
      if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 4)) || (permet_affichage > 2))
         cout << "\n nb de relaxation individuelle effectuee: " << compteur_relaxation_ok;

	  };// -- fin du cas où l'on applique l'amortissement
		
	  return relax_vit;
	};
    
// pilotage de la convergence, ceci en fonction du type de pilotage choisit
// ramène true si le pas de temps a été modifié, false sinon
// arret = true: indique dans le cas d'une non convergence, neanmoins rien n'a changé, il faut donc arrêter
bool Algori::Pilotage_du_temps(Charge* charge,bool& arret)
  { arret = false; // initialisation
    bool modif_pas_de_temps = false;
    // choix en fonction des différents cas
    switch (pa.TypeDePilotage())
     { case  AUCUN_PILOTAGE: // pas de pilotage
       {if (phase_de_convergence >= 0) // cas normal 
        { // on regarde s'il faut augmenter le paramètre de charge
         if (nombre_de_bonnes_convergences > pa.Nb_bonne_convergence())
          { // cas ou on augmente
        //    modif_pas_de_temps = charge->Augmente_inc_charge(pa.Facteur_augmentation());
            // réinitialisation pour la suite
 //           nombre_de_bonnes_convergences = 0;
//				nombre_de_mauvaises_convergences = 0;
           }
         bool etat_charge= charge->Avance();// à ne pas mettre dans l'expression qui suit
         // car si modif_pas_de_temps est true, il ne fera pas charge->avance(), à moins de le mettre en premier
         modif_pas_de_temps = (etat_charge || modif_pas_de_temps ); // increment de charge
         }
       else
        { // cas ou l'incrément précédent n'a pas convergé
     //    modif_pas_de_temps = charge->Diminue_inc_charge(pa.Facteur_diminution());
     //    if (!modif_pas_de_temps)
          { // dans le cas où il n'y a pas de modification du pas de temps on le signale
            cout << "\n pas de pilotage ===> pas de modification possible du pas de temps"; 
            arret = true;
          }
          bool etat_charge= charge->Avance();// à ne pas mettre dans l'expression qui suit
          // car si modif_pas_de_temps est true, il ne fera pas charge->avance(), à moins de le mettre en premier
          modif_pas_de_temps = (etat_charge || modif_pas_de_temps ); // increment de charge
          phase_de_convergence = 0; // initialisation pour repartir d'un bon pied
         }
        break;
        }
       case  PILOTAGE_BASIQUE: case PILOT_GRADIENT : // cas du pilotage basique
		   // et également pour le gradient pour l'instant
        {//cout << "\n phase de convergence " << phase_de_convergence << endl;
         if (phase_de_convergence >= 0) // cas normal
          {// on regarde s'il faut augmenter le paramètre de charge
           if (phase_de_convergence < 2)
            {if (nombre_de_bonnes_convergences > pa.Nb_bonne_convergence())
              { // cas ou on augmente le pas de temps
                bool etat_charge= charge->Augmente_inc_charge(pa.Facteur_augmentation());
                modif_pas_de_temps = (etat_charge || modif_pas_de_temps );
                // réinitialisation pour la suite
                nombre_de_bonnes_convergences = 0;nombre_de_mauvaises_convergences=0;
              };
            }
           else // cas la convergence n'a pas été bonne: phase_de_convergence >= 2
            { // on va diminuer l'incrément
              bool etat_charge= charge->Diminue_inc_charge(pa.Fact_dim_en_mauvaiseConv());
              modif_pas_de_temps = (etat_charge || modif_pas_de_temps );
              // réinitialisation pour la suite
              nombre_de_bonnes_convergences = 0;
      // on laisse le nombre de mauvaises convergence nombre_de_mauvaises_convergences = 0; // mais ça a convergé donc on met à 0
            };
          }
         else
          { // cas ou l'incrément précédent n'a pas convergé, diminution du pas de temps
           // on revient en arrière
           charge->Precedant(false);
           bool etat_charge= charge->Diminue_inc_charge(pa.Facteur_diminution());
           modif_pas_de_temps = (etat_charge || modif_pas_de_temps );
           if (!modif_pas_de_temps)
            { // dans le cas où il n'y a pas de modification du pas de temps on le signale
              // l'arrêt sera traité par le programme appelant
              cout << "\n Pilotage ===> pas de modification possible du pas de temps";
              arret = true;
             };
           phase_de_convergence = 0; // initialisation pour repartir d'un bon pied
          }
         bool etat_charge= charge->Avance();
         modif_pas_de_temps = (etat_charge || modif_pas_de_temps );
         break;
        };
/*       case  PILOT_GRADIENT: // pour l'instant idem cas du pilotage basique
       {if (phase_de_convergence >= 0) // cas normal
        {// on regarde s'il faut augmenter le paramètre de charge
         if (nombre_de_bonnes_convergences > pa.Nb_bonne_convergence())
          { // cas ou on augmente le pas de temps
            modif_pas_de_temps = charge->Augmente_inc_charge(pa.Facteur_augmentation());
            // réinitialisation pour la suite
            nombre_de_bonnes_convergences = 0;
           } 
         } 
        else
        { // cas ou l'incrément précédent n'a pas convergé, diminution du pas de temps
         // on revient en arrière
         charge->Precedant(false);
         modif_pas_de_temps = charge->Diminue_inc_charge(pa.Facteur_diminution());
         if (!modif_pas_de_temps)
          { // dans le cas où il n'y a pas de modification du pas de temps on le signale
            // l'arrêt sera traité par le programme appelant
            cout << "\n Pilotage ===> pas de modification possible du pas de temps"; 
            arret = true; 
           }             
         phase_de_convergence = 0; // initialisation pour repartir d'un bon pied
         }
        charge->Avance(); // increment de charge
        break; 
        }
*/		  
       default:
        { cout << "\n erreur le cas de pilotage " << Nom_TypePilotage(pa.TypeDePilotage())
               << " n'est pas implante !!! "
               << "\n Algori::Pilotage_du_temps(... ";
          // dans le cas où un comptage du calcul est en cours on l'arrête
          if (tempsCalEquilibre.Comptage_en_cours())
                  tempsCalEquilibre.Arret_du_comptage();
          Sortie(1);
        }; 
      }; // fin des différents cas de pilotage  
     // retour de l'information sur la modification ou non du pas de temps
//debug
//cout << "\n Algori::Pilotage_du_temps: phase_de_convergence=" << phase_de_convergence 
//     << " nombre_de_bonnes_convergences= " << nombre_de_bonnes_convergences
//	  << " nombre_de_mauvaises_convergences " << nombre_de_mauvaises_convergences << endl  ;
//fin debug
     cout << flush; // pour que l'on puisse voir quelque chose ! en cours de calcul en particulier s'il y a pb
     return modif_pas_de_temps;
    };

// pilotage de la fin des itération en implicite : utilisation conjointe avec: Pilotage_du_temps()
// ramène vraie si la convergence est correcte, sinon false
bool Algori::Pilotage_fin_iteration_implicite(int compteur)
  {  bool fin_iter = false;
     if (a_converge)
      { if (compteur < pa.Nb_iter_pour_bonne_convergence())
          { // cas où le calcul a convergé en moins de pa.Nb_iter_pour_bonne_convergence
            nombre_de_bonnes_convergences++; // incrémentation de la bonne convergence
            nombre_de_mauvaises_convergences=0;
            phase_de_convergence = 1; //signalement pour Pilotage_du_temps 
           }
         else if (compteur > (pa.Nb_iter_pour_mauvaise_convergence()))
          { // cas où on considère que la convergence n'a pas été bonne et qu'il vaut mieux diminuer l'incrément
            // pour la prochaine fois, on efface également le nombre de bonnes convergences
            nombre_de_bonnes_convergences = 0;
            nombre_de_mauvaises_convergences++;
         	  phase_de_convergence = 2;	
         	}
         else // cas où c'est ni très bon, ni très mauvais, on laisse en l'état
         	{ phase_de_convergence = 1; //signalement pour Pilotage_du_temps
				        nombre_de_bonnes_convergences=0;
				        nombre_de_mauvaises_convergences=0;				 
         	}; 
//		  nombre_de_mauvaises_convergences = 0; // dans tous les cas on a convergé		            
        fin_iter = true; // on signale qu'il y a eu convergence et que l'on a fini les itérations
      }
     else if (compteur > pa.Iterations()) 
      {  nombre_de_mauvaises_convergences++;
			      nombre_de_bonnes_convergences = 0;
         if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
					     { cout << "\n============================================================================"
						           << "\n   ******* NON convergence des iterations d'equilibre ********* "
						           << "\n============================================================================";
          }
				     else 
					     { cout << "\n============================================================================="
						           << "\n *** NON convergence of the balance iteration                           ***** "
						           << "\n=============================================================================";
						    };
         phase_de_convergence = -1; //signalement pour Pilotage_du_temps    
         fin_iter = false; // on signale qu'il n'y a pas convergence
      }
     else if (phase_de_convergence == -2) // cas où c'est Convergence() qui est la cause de l'arret
                                          // pour cause d'évolution divergente des résidus
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
					    { cout << "\n============================================================================"
						          << "\n  *** NON convergence equilibre: evolution divergente des residus ********* "
						          << "\n============================================================================";
         }
				    else 
         { cout << "\n============================================================================="
						          << "\n *** NON convergence of balance: divergente evolution of the residues   ***** "
						          << "\n=============================================================================";
						   };
        fin_iter = false; // on signale qu'il n'y a pas convergence
				    nombre_de_bonnes_convergences=0;
				    nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -3) // cas où c'est Convergence() qui est la cause de l'arret
                                          // pour cause de variation de ddl minimal atteinte
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
          { cout << "\n============================================================================"
                 << "\n  *** NON convergence equilibre: variation de ddl minimal atteinte ********* "
                 << "\n============================================================================";
          }
        else 
          { cout << "\n============================================================================="
                 << "\n *** NON convergence of balance: reaching of the minimum allowed of ddf ***** "
                 << "\n=============================================================================";
          };
        fin_iter = false; // on signale qu'il n'y a pas convergence
				    nombre_de_bonnes_convergences=0;
				    nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -4) // cas où c'est Convergence() qui est la cause de l'arret
                                          // pour cause de variation de ddl maximal atteinte
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
					     { cout << "\n============================================================================"
								         << "\n  *** NON convergence equilibre: variation de ddl maximum atteinte ********* "
								         << "\n============================================================================";
						    }
				    else
          { cout << "\n============================================================================="
								         << "\n *** NON convergence of balance: reaching of the maximum allowed of ddl ***** "
								         << "\n=============================================================================";
						    };
				    nombre_de_bonnes_convergences=0;
				    nombre_de_mauvaises_convergences++;
        fin_iter = false; // on signale qu'il n'y a pas convergence
      }
     else if (phase_de_convergence == -5) // cas d'un arrêt à cause d'un jacobien négatif, qui est géré
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         { cout << "\n============================================================================"
                << "\n  *** NON convergence due a la gestion de jacobien negatif         ********* "
								        << "\n============================================================================";
         }
				    else
					    { cout << "\n============================================================================"
								        << "\n  *** NON convergence due to the management of negative jacobian  ********* "
								        << "\n============================================================================";
						   };
       fin_iter = false; // on signale qu'il n'y a pas convergence
				   nombre_de_bonnes_convergences=0;
				   nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -6) // cas d'un arrêt à cause d'une variation de jacobine trop grande
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         { cout << "\n============================================================================"
								        << "\n  *** NON convergence due a une variation de jacobien trop grande ********* "
								        << "\n============================================================================";
						   }
				    else
					     { cout << "\n============================================================================"
								         << "\n  *** NON convergence due to a too big variation of the jacobian  ********* "
								         << "\n============================================================================";
						    };
        fin_iter = false; // on signale qu'il n'y a pas convergence
				    nombre_de_bonnes_convergences=0;
				    nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -7) // cas d'un arrêt à cause d'une non convergence de la loi de comp
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
			       {if(ParaGlob::Francais())
					       { cout << "\n============================================================================"
								           << "\n  *** NON convergence due a la loi de comportement                 ********* "
								           << "\n============================================================================";
						      }
				       else
					       { cout << "\n============================================================================"
								           << "\n  *** NON convergence due to the material behaviour               ********* "
								           << "\n============================================================================";
						      };
				      };
        fin_iter = false; // on signale qu'il n'y a pas convergence
        nombre_de_bonnes_convergences=0;
        nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -8) // cas d'un arrêt à cause d'une non convergence de la loi de comp
      // mais pour un pb qui couve depuis plusieurs incréments
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         {if(ParaGlob::Francais())
            { cout << "\n============================================================================"
                   << "\n  *** NON convergence du a la loi de comportement                 ********* "
                   << "\n  ***  probablement due a des erreurs sur plusieurs increments    ********* "
                   << "\n============================================================================";
            }
          else 
            { cout << "\n============================================================================"
                   << "\n  *** NON convergence due to the material behaviour               ********* "
                   << "\n  ***  probably due to  errors along several increments           ********* "
                   << "\n============================================================================";
            };
         };
        fin_iter = false; // on signale qu'il n'y a pas convergence
        nombre_de_bonnes_convergences=0;
        nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -9) // cas d'un arrêt à cause d'un pb de la résolution du syteme linéaire
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         {if(ParaGlob::Francais())
            { cout << "\n=============================================================================="
                   << "\n  *** NON convergence due a un pb de resolution du systeme lineaire    ********* "
                   << "\n==============================================================================";
            }
          else 
            { cout << "\n================================================================================="
                   << "\n  *** NON convergence due to the resolution of the global linear system  ********* "
                   << "\n=================================================================================";
            };
         };
         fin_iter = false; // on signale qu'il n'y a pas convergence
         nombre_de_bonnes_convergences=0;
         nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -10) // cas d'un arrêt à cause d'un pb de calcul des efforts extérieurs
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         {if(ParaGlob::Francais())
            { cout << "\n=============================================================================="
                   << "\n  *** NON convergence due a un pb de calcul des efforts exterieurs    ********* "
                   << "\n==============================================================================";
            }
          else
            { cout << "\n================================================================================="
                   << "\n  *** NON convergence due to extern loading problem  ********* "
                   << "\n=================================================================================";
            };
         };
         fin_iter = false; // on signale qu'il n'y a pas convergence
         nombre_de_bonnes_convergences=0;
         nombre_de_mauvaises_convergences++;
      }
     else if (phase_de_convergence == -11) // cas d'un arrêt à cause d'un nan ou infini sur forces ext, int
      { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         {if(ParaGlob::Francais())
            { cout << "\n=============================================================================="
                   << "\n *** NON convergence due a un nan ou infini sur forces int, ext, totales ********* "
                   << "\n==============================================================================";
            }
          else
            { cout << "\n================================================================================="
                   << "\n  *** NON convergence due to nan value   ********* "
                   << "\n=================================================================================";
            };
         };
         fin_iter = false; // on signale qu'il n'y a pas convergence
         nombre_de_bonnes_convergences=0;
         nombre_de_mauvaises_convergences++;
      }
     else // cas d'une erreur de pilotage: on ne sait pas pourquoi il s'est arrêté sans convergence
         { if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
			          {if(ParaGlob::Francais())
					          {cout << "\n cas de pilotage de fin d'iteration inconnues ? ou alors qui suit une erreur signalee precedemment ";
                if (ParaGlob::NiveauImpression() > 5)
                   cout  << "\n   bool Pilotage_fin_iteration_implicite(int compteur) ";
                cout << endl;
               }
				          else {cout << "\n unknown case of management for the end of iteration or, following an error  ";
							             if (ParaGlob::NiveauImpression() > 5)
                      cout  <<  "\n   bool Pilotage_fin_iteration_implicite(int compteur) ";
                    cout << endl;
                   };
			          };
           // si il ne s'agit pas d'une sortie voulue des itérations et incréments
           // on s'arrête
           if (!(pa.EtatSortieEquilibreGlobal()))
             { // dans le cas où un comptage du calcul est en cours on l'arrête
               if (tempsCalEquilibre.Comptage_en_cours())
                  tempsCalEquilibre.Arret_du_comptage();
               Sortie(1);
             };
           // sinon on continue
         }
     return fin_iter;
//debug
if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
 { cout << "\n Algori::Pilotage_fin_iteration_implicite: phase_de_convergence=" << phase_de_convergence
        << " nombre_de_bonnes_convergences= " << nombre_de_bonnes_convergences
	       << " nombre_de_mauvaises_convergences " << nombre_de_mauvaises_convergences << endl ;
 };
//fin debug
		
    };
		 
 // pilotage spécifique de la fin de la relaxation dynamique, en tenant compte des différents critères demandés
 // relax_vit_acce :    
 //   = 1  : il y a eu relaxation mais le calcul continue
 //   = 0  : pas de relaxation
 //   = -1 : on peut arrêter le calcul car le critère sur delta_ddl est satisfait
 // compteur: indique le nombre d'itération ou d'incrément en cours
 // iter_ou_incr : =1: il s'agit d'itération", =0 il s'agit d'incrément 
 // lit arret, puis modifie arret en true ou false = il faut arrêter ou non
 // si true met à jour en interne certains paramètres de gestion de pilotage pour dire que l'on a bien convergé
bool	Algori::Pilotage_fin_relaxation_et_ou_residu(const int& relax_vit_acce,const int& iter_ou_incr,
                                                 const int& compteur, const bool& arretResidu, bool& arret)
  {
    // pour mémoire: arret_a_equilibre_statique indique, si diff de 0, que l'on veut contrôler une convergence vers la solution
		  // statique. Si = 1: le residu statique est utilise comme critere
		  //           si = 2: nécessite que les deux critères précédents soient satisfait                   
	  
    // examen de la convergence éventuelle, utilisant le déplacement et/ou le résidu
		  string nom="l'iteration"; if (!iter_ou_incr) nom = "l'increment";
    switch (relax_vit_acce)
    { case -1: // cas où le critère en déplacement est ok
       {if (!arret_a_equilibre_statique) // cas ou on ne considère pas l'équilibre en résidu
              {arret=true;a_converge = true;
               cout << "\n   critere en deplacement satisfait pour "<< nom << " : "  << compteur
                    << "  max_deltaX= " <<delta_X.Max_val_abs() << endl;
              }
        else if ((arret_a_equilibre_statique==2) && arretResidu) // cas ou  considère l'équilibre en résidu et en amortissement cinétique
              {arret=true;a_converge = true;
               cout << "\n   critere en deplacement et en residu satisfait pour "<< nom << " : "
                    << compteur << "  max_deltaX= " <<delta_X.Max_val_abs()  <<endl;
              }
        
       }
       break;
      case 0: case 1: // cas où le critère en déplacement n'est pas satisfait ou qu'il y a eu une relaxation (donc pas encore de convergence)
       {if ( (arret_a_equilibre_statique==1) && arretResidu) // cas du critère en résidu uniquement
          {arret=true;a_converge = 1;
           cout << "\n   critere equilibre statique satisfait pour "<< nom << " : "  << compteur << " " <<endl;
          }
        else // sinon on est pas avec un contrôle avec le résidu donc on se base uniquement sur le déplacement
         {arret=false;a_converge = false;}
       }
       break;

     default:
       cout << "\n *** erreur: relax_vit_acce= "<<relax_vit_acce
            << "\n ce cas n'est pas prevu... \n Algori::Pilotage_fin_relaxation_et_ou_residu(...";
       Sortie(1);
       break;
   };
   
   
//    if (relax_vit_acce == -1) // en amortissement cinétique, on a convergé
//       {if (!arret_a_equilibre_statique) // cas ou on ne considère pas l'équilibre en résidu
//              {arret=true;a_converge = 1;
//               cout << "\n   critere en deplacement satisfait pour "<< nom << " : "  << compteur << " " <<endl;
//              }
//        else if ((arret_a_equilibre_statique==2) && arretResidu) // cas ou  considère l'équilibre en résidu et en amortissement cinétique
//              {arret=true;a_converge = 1;
//               cout << "\n   critere en deplacement et en residu satisfait pour "<< nom << " : "
//                    << compteur << " " <<endl;
//              }
//
//       }
//    else if ( (arret_a_equilibre_statique==1) && arretResidu) // cas du critère en résidu uniquement
//        {arret=true;a_converge = 1;
//         cout << "\n   critere equilibre statique satisfait pour "<< nom << " : "  << compteur << " " <<endl;
//        };
		  return arret;
		};	    

// pilotage à chaque itération: 
//    - sous ou sur relaxation
//    - limitation de la norme de delta ddl
//    - indication de précision pour les éléments (utilisée pour certain type de gestion d'hourglass)
// en entrée: le maxi de var ddl juste après résolution
// ramène le maxi de var ddl maxDeltaDdl modifié, ainsi que sol modifié éventuellement, 
void Algori::Pilotage_chaque_iteration(Vecteur* sol,double& maxDeltaDdl,const int& itera,LesMaillages * lesMail)
{ // -- sur ou sous relaxation
  double sursousrelax = pa.SurSousRelaxation();
  (*sol) *= sursousrelax; 
  if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 3)) || (permet_affichage > 2))
  if (sursousrelax != 1.)
			{if(ParaGlob::Francais())
			    {cout << "\n intervention d'une sur-sous relaxation:  ";}
			 else {cout << "\n intervention of an over or under relaxation: ";};
			 cout << sursousrelax;
			};	 
  // --limitation de la norme de delta ddl 
  // dans le cas où le maxi de variation de ddl est supérieur à la limite fixée on ramène la valeur
  if (pa.NormeMax_increment() > 0.)
   {// cas d'une limitation globale (tout le vecteur sol change
    if (maxDeltaDdl > pa.NormeMax_increment()) 
     { double facteur=(pa.NormeMax_increment()/maxDeltaDdl);
       if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
				    {if(ParaGlob::Francais())
				         {cout << "\n ===>>>>> intervention de la reduction du vecteur increment de ddl ";}
         else {cout << "\n ===>>>>> intervention of the reduction for the vector increment of dof ";};
			      cout << facteur;		 
				    };
       (*sol) *= facteur;
       // mise à jour  du maxi de variation de ddl  
       maxDeltaDdl = pa.NormeMax_increment();        
     };
   }
  else // ici il s'agit d'un écrétage, seule les composantes supérieures au maxi, sont écrètées
   {// cas d'une limitation globale (tout le vecteur sol change
    if (maxDeltaDdl > Abs(pa.NormeMax_increment())) 
     { double maxa= Abs(pa.NormeMax_increment());
       Vecteur & vect = *sol; // pour simplifier
       int taill=vect.Taille();
       bool modif = false;
       for (int i=1;i<= taill;i++)
         { if (Abs(vect(i)) > maxa) 
            { vect(i) = Signe(vect(i),maxa); modif = true;};
         };
       if (modif) 
         {if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
				       {if(ParaGlob::Francais())
					            {cout << "\n ===>>>>> intervention de la reduction du vecteur increment de ddl, limitation des maxis ";}
					       else {cout << "\n ===>>>>> intervention of the reduction for the vector increment of dof, limitation of the maxis ";};
						     };	 
				     };
          // mise à jour  du maxi de variation de ddl  
          maxDeltaDdl = maxa;        
     };
   }; 
		
	 // -- indication de précision pour les éléments (utilisée pour certain type de gestion d'hourglass)
	 // pour l'instant l'idée est la suivante: si le résidu augmente => on demande un calcul précis sur les éléments
	 // pour le prochain calcul
  int nb_cycle_controle = MaX(2,pa.NbCycleControleResidu());
	 if (itera == 0)
	  { lesMail->Drapeau_preparation_calcul_precis(10);} // par défaut à l'itération 0 on demande un calcul précis
  if (itera > 1)
	  { bool change = false; // indicateur pour la demande d'un calcul précis		
		   if (itera <= nb_cycle_controle)
		     { if (max_var_residu(itera) > max_var_residu(itera-1)) change = true;}
     else { if (max_var_residu(nb_cycle_controle) > max_var_residu(nb_cycle_controle-1)) change = true;};
		   if (change) // dans le cas où le résidu augmente => on demande un calcul précis sur les éléments
		     { lesMail->Drapeau_preparation_calcul_precis(4);}
		   else {lesMail->Drapeau_preparation_calcul_precis(4);}; // sinon pas de modif
   };
};

// pilotage pour l'initialisation de Xtdt à chaque début d'incrément, en implicite
// avec la même direction qu'au pas précédent
// ramène true, si c'est ok pour initialiser Xtdt 
bool Algori::Pilotage_init_Xtdt()
 { return (nombre_de_mauvaises_convergences <= 2); 
  };

 
// pilotage en dynamique implicite, du maxi des déplacements et/ou vitesses
// limitation spécifiquement des maxi en déplacements et/ou vitesses
// il y a calcul éventuel de delta_X, mais c'est considéré comme une variable de travail
void Algori::Pilotage_maxi_X_V(const Vecteur& X_t,Vecteur& X_tdt,const Vecteur& V_t,Vecteur& V_tdt)
{ // --limitation de la norme de delta X
  // dans le cas où le maxi de variation de delta est supérieur à la limite fixée on ramène la valeur
  double normeMax_X_increment = pa.NormeMax_X_increment();
  if ((normeMax_X_increment > 0.)&&( normeMax_X_increment!= ConstMath::tresgrand))
   {// cas d'une limitation globale sur X
    if (vec_trav.Taille() != X_t.Taille()) vec_trav.Change_taille(X_t.Taille());
    vec_trav = X_tdt; vec_trav -= X_t;
    double maxDeltaDdl=vec_trav.Max_val_abs();
    if (maxDeltaDdl > normeMax_X_increment)
     { double facteur=(normeMax_X_increment/maxDeltaDdl);
       if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 2)) || (permet_affichage > 0))
				    {if(ParaGlob::Francais())
				         {cout << "\n ===>>>>> intervention de la reduction du vecteur delta X ";}
         else {cout << "\n ===>>>>> intervention of the reduction for the vector delta X ";};
			      cout << facteur << "( maxDeltaX= " << maxDeltaDdl<<")";
				    };
       vec_trav *= facteur;
       X_tdt = X_t; X_tdt += vec_trav;
     };
   }
  else if ( normeMax_X_increment!= ConstMath::tresgrand)
  // cas négatif ici il s'agit d'un écrétage, seule les composantes supérieures au maxi, sont écrètées
   {// cas d'une limitation globale (tout le vecteur sol change
    int taill=X_t.Taille();
    if (vec_trav.Taille() != X_t.Taille()) vec_trav.Change_taille(X_t.Taille());
    vec_trav = X_tdt; vec_trav -= X_t;
    double maxDeltaDdl=vec_trav.Max_val_abs();
    if (maxDeltaDdl > -normeMax_X_increment)
     { double maxa= -normeMax_X_increment;
       bool modif = false;
       for (int i=1;i<= taill;i++)
         { if (Abs(vec_trav(i)) > maxa)
            { vec_trav(i) = Signe(vec_trav(i),maxa); modif = true;};
         };
       if (modif) 
         {if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 2)) || (permet_affichage > 0))
				       {if(ParaGlob::Francais())
					            {cout << "\n ===>>>>> intervention de la reduction du vecteur delta X, limitation des maxis ";}
					       else {cout << "\n ===>>>>> intervention of the reduction for the vector delta X, limitation of the maxis ";};
						     };	 
          X_tdt = X_t; X_tdt += vec_trav;
				     };
     };
   };
 
// --limitation de la norme de delta V
  // dans le cas où le maxi de variation de delta est supérieur à la limite fixée on ramène la valeur
  double normeMax_V_increment = pa.NormeMax_V_increment();
  if ((normeMax_V_increment > 0.)&&( normeMax_V_increment!= ConstMath::tresgrand))
   {// cas d'une limitation globale sur V
    if (vec_trav.Taille() != V_t.Taille()) vec_trav.Change_taille(V_t.Taille());
    vec_trav = V_tdt; vec_trav -= V_t;
    double maxDeltaDdl=vec_trav.Max_val_abs();
    if (maxDeltaDdl > normeMax_V_increment)
     { double facteur=(normeMax_V_increment/maxDeltaDdl);
       if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 2)) || (permet_affichage > 0))
				    {if(ParaGlob::Francais())
				         {cout << "\n ===>>>>> intervention de la reduction du vecteur delta V ";}
         else {cout << "\n ===>>>>> intervention of the reduction for the vector delta V ";};
			      cout << facteur ;
				    };
       vec_trav *= facteur;
       V_tdt = V_t; V_tdt += vec_trav;
     };
   }
  else if ( normeMax_V_increment!= ConstMath::tresgrand)
    // cas négatif ici il s'agit d'un écrétage, seule les composantes supérieures au maxi, sont écrètées
   {// cas d'une limitation globale (tout le vecteur sol change
    int taill=V_t.Taille();
    if (vec_trav.Taille() != V_t.Taille()) vec_trav.Change_taille(V_t.Taille());
    vec_trav = V_tdt; vec_trav -= V_t;
    double maxDeltaDdl=vec_trav.Max_val_abs();
    if (maxDeltaDdl > -normeMax_V_increment)
     { double maxa= -normeMax_V_increment;
       bool modif = false;
       for (int i=1;i<= taill;i++)
         { if (Abs(vec_trav(i)) > maxa)
            { vec_trav(i) = Signe(vec_trav(i),maxa); modif = true;};
         };
       if (modif) 
         {if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 4)) || (permet_affichage > 0))
				       {if(ParaGlob::Francais())
					            {cout << "\n ===>>>>> intervention de la reduction du vecteur delta V, limitation des maxis ";}
					       else {cout << "\n ===>>>>> intervention of the reduction for the vector delta V, limitation of the maxis ";};
						     };	 
          V_tdt = V_t; V_tdt += vec_trav;
				     };
     };
   }; 
 
};

// passage des grandeurs gérées par l'algorithme de tdt à t
// en particulier les énergies et les puissances
void Algori::TdtversT()
 {  // cas du calcul des énergies, passage des grandeurs de tdt à t
     // les énergies globales
    E_cin_t=E_cin_tdt; E_int_t = E_int_tdt; E_ext_t = E_ext_tdt;
    // la quantité de mouvement globale
    q_mov_t = q_mov_tdt; 
     // les forces internes et externes généralisées
    F_int_t = F_int_tdt;  F_ext_t = F_ext_tdt; 
     // l'énergie visqueuse numérique
    E_visco_numerique_t=E_visco_numerique_tdt;

    // les puissances
    P_acc_t = P_acc_tdt; P_int_t = P_int_tdt; P_ext_t = P_ext_tdt; bilan_P_t = bilan_P_tdt;
  
    // cas particulier de l'amortissement cinétique à chaque noeud
    if (amortissement_cinetique_au_noeud)
       E_cinNoe_t = E_cinNoe_tdt;

 };                                          

// vérification d'une singularité éventuelle de la matrice de raideur de dim: nbddl 
// pour un problème de mécanique, en comparant le nombre de point d'intégration total 
// et le nombre totale de degré de liberté 
void Algori::VerifSingulariteRaideurMeca(int nbddl,const LesMaillages& lesMail) const
	{
   #ifdef UTILISATION_MPI
   // cas d'un  calcul //, seule la matrice du CPU 0 est concernée
   if (ParaGlob::Monde()->rank() != 0)
     return ;
   #endif
   if (ParaGlob::NiveauImpression() != 0)
	   {  // récup du nombre de grandeurs génératrices
//       Enum_ddl enu(SIG11); // init par défaut
     
       // on récupère la liste des problèmes physiques gérés par les éléments de tous les maillages
       const list <Enum_ddl >& li = lesMail.Ddl_representatifs_des_physiques();
       // on boucle sur la liste
       list <Enum_ddl >::const_iterator ili,iliend = li.end();
       int nbpt=0; // init
       for (ili = li.begin();ili != iliend; ili++)
        nbpt += lesMail.NbTotalGrandeursGeneratrices(*ili);
       // affichage
       const double ratio = 1.2; // le pourcentage
       if (nbpt < ratio * nbddl)
         if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 0)) || (permet_affichage > 0))
         	{ cout << "\n *** Attention, le nombre total de grandeurs generatrices aux point d'integrations dans "
         	       << " le probleme: " << nbpt << "  est inferieur "
         	       << " a " << ratio << " * le nombre total de ddl du probleme !! c-a-d: "<< int(ratio * nbddl)
         	       << "\n en consequence il est possible que la matrice de raideur soit naturellement singuliere \n\n ";
               if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 7)) || (permet_affichage > 0))
                 cout << "\n void Algori::VerifSingulariteRaideurMeca(... ";      
         	};
	   };
	};                                         
        
// retouve le numéro de l'incrément sauvegardé trouvé dont le numéro est juste inf ou égale au 
// dernier incrément calculé  - nb_incr_en_arriere.
// ramène false si on n'a pas trouvé d'incrément a-doc
// en retour icharge prend la valeur de l'incrément trouvé (si tout c'est bien passé)
bool Algori::Controle_retour_sur_un_increment_enregistre(int nb_incr_en_arriere,int& icharge)        
{   // tout d'abord on calcul le numéro d'incrément auquel il faut se placer
    List_io <Entier_et_Double>::iterator inc_ok = list_incre_temps_calculer.begin();
    List_io <Entier_et_Double>::iterator inc_fin = list_incre_temps_calculer.end();
    for (int i=1;i<=nb_incr_en_arriere;i++)
     if (inc_ok != inc_fin)
       {inc_ok++;  } // les incréments sont rangés du dernier vers le premier
     else
       {return false;}; // sinon l'opération n'est pas possible 
    // donc inc_ok  représente maintenant le num maxi à partir duquel on va chercher un incrément enregistré
    List_io <Entier_et_Double>::iterator in_sauve =  list_incre_temps_sauvegarder.begin(); 
    List_io <Entier_et_Double>::iterator in_sauve_fin =  list_incre_temps_sauvegarder.end(); 
    // on recherche le premier incrément sauvegardé juste inf à celui que l'on s'est fixé
    while (((*in_sauve).n > (*inc_ok).n)&&(in_sauve != in_sauve_fin))
       in_sauve++;
    // on vérifie que l'incrément est valide sinon on est obligé d'arrêter
    if (in_sauve == in_sauve_fin)
    	{return false;}; // opération non possible
    // sinon c'est ok 
    icharge = (*in_sauve).n;
    return true;
};

// cohérence du type de convergence adopté en fonction de l'algorithme: par exemple on vérifie que l'on n'utilise
// pas une convergence sur l'énergie cinétique avec un algo statique !!
void Algori::Coherence_Algo_typeConvergence()
 { if ((Evolution_temporelle_du_calcul(typeCalcul) != DYNAMIQUE) 
       && (strstr(pa.Norme().nom1.c_str(),"E_cinetique/E_statique"))
		 && (Evolution_temporelle_du_calcul(typeCalcul) != AUCUNE_EVOLUTION))
	    { cout << "\n *** erreur: utilisation de la norme E_cinetique/E_statique avec un calcul non dynamique"
	           << " arret ! ";
       if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5)) || (permet_affichage > 0))
	         cout << "\n Algori::Coherence_Algo_typeConvergence() ";
       // dans le cas où un comptage du calcul est en cours on l'arrête
       if (tempsCalEquilibre.Comptage_en_cours())
           tempsCalEquilibre.Arret_du_comptage();
	      Sortie(1);
	    };
 };