// 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 "AlgoriNewmark.h"
#include "ConstMath.h"

    //------- décomposition en 3 du calcul d'équilibre -------------
    // a priori   : InitAlgorithme  et FinCalcul ne s'appellent qu'une fois, 
    // par contre : CalEquilibre peut s'appeler plusieurs fois, le résultat sera différent si entre deux calcul
    //              certaines variables ont-été changés    


// initialisation
void AlgoriNewmark::InitAlgorithme(ParaGlob * paraGlob,LesMaillages * lesMail,
               LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
               ,VariablesExporter* varExpor
               ,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage,
           Charge* charge,LesCondLim* lesCondLim,LesContacts* lesContacts
           ,Resultats* resultats)
  { // INITIALISATION globale 
    tempsInitialisation.Mise_en_route_du_comptage(); // temps cpu
    Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(DYNA_IMP); // transfert info

  #ifdef UTILISATION_MPI
    // calcul de l'équilibrage initiale par le cpu 0
    if (distribution_CPU_algo.Tableau_element_CPU_en_cours()->Taille() == 0 )
       {distribution_CPU_algo.Calcul_Equilibrage_initiale(lesMail,lesContacts);
        distribution_CPU_algo.Passage_Equilibrage_aux_CPU();
        paraGlob->Init_tableau(distribution_CPU_algo.Tableau_element_CPU_en_cours()
                                ,distribution_CPU_algo.Tab_indique_CPU_en_cours()
                                ,distribution_CPU_algo.Tableau_noeud_CPU_en_cours()
                                ,distribution_CPU_algo.Tab_indique_noeud_CPU_en_cours());
       };
  #endif

    // avant toute chose, au cas où l'algo interviendrait après un autre algo
    // on inactive tous les ddl existants
    lesMail->Inactive_ddl();
    // on regarde s'il s'agit d'un pb purement non méca, dans ce cas il faut cependant initialiser les positions
    // à t et tdt pour le calcul de la métrique associée
    { const list <EnumElemTypeProblem >& type_pb = lesMail->Types_de_problemes();
      bool purement_non_meca = true;
      list <EnumElemTypeProblem >::const_iterator il,ilfin=type_pb.end();
      for (il=type_pb.begin();il != ilfin; il++)
        {switch (*il)
          { case MECA_SOLIDE_DEFORMABLE: case MECA_SOLIDE_INDEFORMABLE: case MECA_FLUIDE:
                 purement_non_meca=false; break;
            default: break; // sinon on ne fait rien
          };
        };
      if (purement_non_meca) // si pas de méca, on initialise les coordonnées à t et tdt avec celles de 0
        {lesMail->Init_Xi_t_et_tdt_de_0();}
      else // sinon a minima on active X1
       { lesMail->Active_un_type_ddl_particulier(X1);
       };
    };

    // cas du chargement, on verifie egalement la bonne adequation des references
    charge->Initialise(lesMail,lesRef,pa,*lesCourbes1D,*lesFonctionsnD);
    // on indique que l'on ne souhaite le temps fin stricte, vrai par défaut en implicite
    charge->Change_temps_fin_non_stricte(0);
    // dans le cas où l'on calcul des contraintes et/ou déformation et/ou un estimateur d'erreur
    // à chaque incrément, initialisation
    tenuXVG.Change_taille(3);tenuXVG(1)=X1;tenuXVG(2)=V1;tenuXVG(3)=GAMMA1;
    bool prepa_avec_remont = Algori::InitRemont(lesMail,lesRef,diversStockage,charge,lesCondLim,lesContacts,resultats);      
    if ( prepa_avec_remont)// remise enservice des ddl du pab 
          {lesMail->Inactive_ddl(); lesMail->Active_un_type_ddl_particulier(X1);};      
    // on crée les ddl d'accélération et de vitesse non actif mais libres
    // car les forces intérieures et extérieures sont les entitées duales
    // des déplacements, qui sont donc les seules grandeurs actives à ce stade
    lesMail->Plus_Les_ddl_Vitesse( HSLIBRE);
    lesMail->Plus_Les_ddl_Acceleration( HSLIBRE);
    // on défini globalement que l'on a une combinaison des ddl X V GAMMA en même temps   
    int cas_combi_ddl=1; 
    //   récup du pas de temps, proposé par l'utilisateur, initialisation 
    // dans le cas où le pas de temps ou le pas de temps maxi dépendent d'un pas critique, on le calcul
    // et on met à jour les pas de temps
    this->Gestion_pas_de_temps(lesMail,1); // 1 signifie qu'il y a initialisation    
    // mise en place éventuelle du bulk viscosity
    lesMail->Init_bulk_viscosity(pa.BulkViscosity(),pa.CoefsBulk());
    // mise a zero de tous les ddl et creation des tableaux a t+dt
    // les ddl de position ne sont pas mis a zero ! ils sont initialise
    // a la position courante
    lesMail->ZeroDdl(true);
    // on vérifie que les noeuds sont bien attachés à un élément sinon on met un warning si niveau > 2
    if (ParaGlob::NiveauImpression() > 2)
      lesMail->AffichageNoeudNonReferencer();
    // init des ddl avec les conditions initials 
    // les conditions limites initiales de vitesse et d'accélération sont prise en compte
    // de manière identiques à des ddl quelconques, ce ne sont pas des ddl fixé !!
    lesCondLim->Initial(lesMail,lesRef,lesCourbes1D,lesFonctionsnD,true,cas_combi_ddl);
    // mise à jour des différents pointeur d'assemblage et activation des ddl
    //  a) pour les déplacements qui sont à ce stade les seuls grandeurs actives 
    //       on définit un nouveau cas d'assemblage pour les Xi
    //       à travers la définition d'une instance de la classe assemblage
    if (Ass1_ == NULL) Ass1_ = new Assemblage(lesMail->InitNouveauCasAssemblage(1));
    Assemblage& Ass1 = *Ass1_;
    lesMail->MiseAJourPointeurAssemblage(Ass1.Nb_cas_assemb());// mise a jour des pointeurs d'assemblage  
    int nbddl_X = lesMail->NbTotalDdlActifs(X1); // nb total de ddl de déplacement
                 // qui est le même pour les accélérations et les vitesses     
    //  b) maintenant le cas des vitesses qui doivent donc être activées
    //       on définit un nouveau cas d'assemblage pour les Vi
    if (Ass2_ == NULL) Ass2_ = new Assemblage(lesMail->InitNouveauCasAssemblage(1)); 
    Assemblage& Ass2 = *Ass2_;
    lesMail->Inactive_un_type_ddl_particulier(X1);  //       on inactive les Xi
    lesMail->Active_un_type_ddl_particulier(V1);    //       on active les Vi
    lesMail->MiseAJourPointeurAssemblage(Ass2.Nb_cas_assemb());        // mise a jour des pointeurs d'assemblage  
    //  c) idem pour les accélérations
    //       on définit le numéro de second membre en cours 
    //       on définit un nouveau cas d'assemblage pour les GAMMAi
    if (Ass3_ == NULL) Ass3_ = new Assemblage(lesMail->InitNouveauCasAssemblage(1)); 
    Assemblage& Ass3 = *Ass3_;
    casAssemb_principal = Ass3_;

    lesMail->Inactive_un_type_ddl_particulier(V1);  //       on inactive les Vi
    lesMail->Active_un_type_ddl_particulier(GAMMA1);    //       on active les GAMMAi
    lesMail->MiseAJourPointeurAssemblage(Ass3.Nb_cas_assemb());        // mise a jour des pointeurs d'assemblage 
    //  d) activation de tous les ddl, maintenant ils peuvent être les 3 actifs
    lesMail->Active_un_type_ddl_particulier(X1);
    lesMail->Active_un_type_ddl_particulier(V1);
    // en fait ces trois pointeurs d'assemblage ne sont utils  que pour la mise en place des conditions
    // limites
    // mise à jour du nombre de cas d'assemblage pour les conditions limites
    // c-a-d le nombre maxi possible (intégrant les autres pb qui sont résolu en // éventuellement)
    lesCondLim->InitNombreCasAssemblage(lesMail->Nb_total_en_cours_de_cas_Assemblage());
    // définition d'un tableau globalisant les numéros d'assemblage de X V gamma
    t_assemb.Change_taille(3);
    t_assemb(1)=Ass1.Nb_cas_assemb();t_assemb(2)=Ass2.Nb_cas_assemb();t_assemb(3)=Ass3.Nb_cas_assemb();
    // récupération des tableaux d'indices généraux des ddl bloqués, y compris les ddls associés
    icas = 1; // pour indiquer au module Tableau_indice que l'on travaille avec l'association X V GAMMA
    li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb,lesRef,charge->Temps_courant(),icas);
    // on définit trois tableau qui serviront à stocker transitoirement les X V GAMMA correspondant au ddl imposés
    int ttsi = li_gene_asso.size();
    X_Bl.Change_taille(ttsi),V_Bl.Change_taille(ttsi),G_Bl.Change_taille(ttsi);

    // choix de la matrice de raideur  du système linéaire en Xi
    Choix_matriciel(nbddl_X,tab_mato,lesMail,lesRef,Ass1.Nb_cas_assemb(),lesCondLim);
    matglob = tab_mato(1); // par défaut c'est la première matrice
    // $$$ essai newton modifié $$$
    //   bool newton_modifie=true;Mat_abstraite matsauv(matglob);
    //  if (newton_modifie) matsauv=new Mat_abstraite(matglob);
   
    // choix de la résolution
    matglob->Change_Choix_resolution(pa.Type_resolution(),pa.Type_preconditionnement());
    // vérification d'une singularité éventuelle de la matrice de raideur 
    VerifSingulariteRaideurMeca(nbddl_X,*lesMail);
    // def vecteurs globaux 
    vglobin.Change_taille(nbddl_X); // puissance interne
    vglobex.Change_taille(nbddl_X); // puissance externe
    vglobaal.Change_taille(nbddl_X,0.); // puissance totale
    // même si le contact n'est pas encore actif, il faut prévoir qu'il le deviendra peut-être !
    if (lesMail->NbEsclave() != 0)
       vcontact.Change_taille(nbddl_X); // puissance de contact

    //  6 vecteur pour une manipulation globale des positions vitesses et accélérations
    // vecteur qui globalise toutes les positions de l'ensemble des noeuds
    X_t.Change_taille(nbddl_X);X_tdt.Change_taille(nbddl_X); delta_X.Change_taille(nbddl_X);
    var_delta_X.Change_taille(nbddl_X);
    // vecteur qui globalise toutes les vitesses de l'ensemble des noeuds
    vitesse_t.Change_taille(nbddl_X);vitesse_tdt.Change_taille(nbddl_X);
    // vecteur qui globalise toutes les accélérations
    acceleration_t.Change_taille(nbddl_X);acceleration_tdt.Change_taille(nbddl_X) ; 
	   inter_tdt.Change_taille(nbddl_X) ;	 
    // vecteur de travail pour la viscosité artificielle
    if (pa.Amort_visco_artificielle()) 
		   { forces_vis_num.Change_taille(nbddl_X);
		     if (Arret_A_Equilibre_Statique())  // si on veut un équilibre statique, on sauvegarde les forces statiques
         { if (vglob_stat != NULL)
             {vglob_stat->Change_taille(vglobaal.Taille());}
           else
             {vglob_stat  = new Vecteur(vglobaal.Taille());}
         };
	    };
    
    // calcul des énergies
    E_cin_tdt = 0.; E_int_t = 0.; E_int_tdt = 0.; // init des différentes énergies
    E_ext_t = 0.; E_ext_tdt = 0.; bilan_E = 0.;   //       "       et du bilan
    F_int_t.Change_taille(nbddl_X); F_ext_t.Change_taille(nbddl_X); // forces généralisées int et ext au pas précédent
    F_int_tdt.Change_taille(nbddl_X); F_ext_tdt.Change_taille(nbddl_X); // forces généralisées int et ext au pas actuel
    residu_final.Change_taille(nbddl_X); // pour la sauvegarde du résidu pour le post-traitement

    // définition d'un pointeur de fonction d'assemblage en fonction du type voulu
    // (symetrique ou non), utilisé dans la prise en compte des efforts extérieurs 
    if (pa.Symetrie_matrice())
      // cas d'un assemblage symétrique
      assembMat = & Assemblage::AssembMatSym;
    else // cas d'un assemblage non symétrique
      assembMat = & Assemblage::AssembMatnonSym;
                                  
    // dans le cas où l'on fait du line search on dimensionne des vecteurs
    // globaux supplémentaires
    if (pa.Line_search())
      { sauve_deltadept = new Vecteur(nbddl_X);
        Vres = new Vecteur(nbddl_X);
        sauve_dept_a_tdt = new Vecteur(nbddl_X);
        v_travail = new  Vecteur(nbddl_X);        
      };

    //  boucle sur les increments de charge
    icharge = 0; // init
   
    // definition des elements de frontiere, ces elements sont utilises pour le contact
    lesMail->CreeElemFront();
    //  calcul éventuel des normales aux noeuds -> init des normales pour t=0
    lesMail->InitNormaleAuxNoeuds(); //utilisé pour la stabilisation des membranes par ex
    // --- init du contact ---
    // doit-être avant la lecture d'un restart, car il y a une initialisation de conteneurs qui est faites
    // qui ensuite est utilisée en restart
    // par exemple il faut initialiser les frontières et la répartition esclave et maître
    // pour préparer la lecture de restart éventuel
    int niveau_substitution = 0; // on intervient sur toutes les matrices
    bool premier_calcul = true;
    if (lesMail->NbEsclave() != 0)
     { // definition des elements de frontiere, ces elements sont utilises pour le contact
       lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_t);
       // initialisation des zones de contacts éventuelles
       lesContacts->Init_contact(*lesMail,*lesRef,lesFonctionsnD);
       // verification qu'il n'y a pas de contact avant le premier increment de charge
       lesContacts->Verification();
       // definition des elements de contact eventuels
       bool nevez_contact = lesContacts->DefElemCont(0.); // au début le déplacement des noeuds est nul
       // mise à jour éventuelle de la matrice de raideur en fonction du contact
       //bool changement_sur_matrice =
       Gestion_stockage_et_renumerotation_avec_contact
                        (premier_calcul,lesMail,nevez_contact,lesCondLim
                        ,lesRef,tab_mato,Ass1.Nb_cas_assemb(),lesContacts,niveau_substitution);
       matglob=tab_mato(1);
     }
    else // sinon il faut éventuellement tenir compte du cht du au CLL
      // NB: c'est pris en compte directement dans le contact, quand celui-ci est actif
     {bool nouvelle_situation_CLL=false;
      Gestion_stockage_et_renumerotation_sans_contact
                          (lesContacts,premier_calcul,lesMail,nouvelle_situation_CLL
                           ,lesCondLim,lesRef,tab_mato,Ass1.Nb_cas_assemb()
                           ,niveau_substitution);
     };
    //--cas de restart et/ou de sauvegarde------------
    // tout d'abord récup du restart si nécessaire
    // dans le cas ou un incrément différent de 0 est demandé -> seconde lecture à l'incrément
    brestart=false;   // booleen qui indique si l'on est en restart ou pas
    if (this->Num_restart() != 0)
        { int cas = 2;
	         // ouverture de base info
	         entreePrinc->Ouverture_base_info("lecture"); 
          this->Lecture_base_info(cas ,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage
                         ,charge,lesCondLim,lesContacts,resultats,(this->Num_restart()));
          icharge = this->Num_restart();//+1;
          //  récup du pas de temps, proposé par l'utilisateur, initialisation 
          this->Gestion_pas_de_temps(lesMail,2); // 2 signifie cas courant
          brestart = true;
          // on oblige les ddls Vi GAMMAi a avoir le même statut que celui des Xi
          // comme les conditions limites cinématiques peuvent être différentes en restart
          // par rapport à celles sauvegardées, on commence par libérer toutes les CL imposées éventuelles
          lesMail->Libere_Ddl_representatifs_des_physiques(LIBRE);
          lesMail->ChangeStatut(cas_combi_ddl,LIBRE);
          // dans le cas d'un calcul axisymétrique on bloque le ddl 3
          if (ParaGlob::AxiSymetrie())
             lesMail->Inactive_un_ddl_particulier(X3);
          // on réinitialise  des ddl avec les conditions initiales because on vient de tout libérer les ddl
          // or dans Initial, il y a des inits à faire au niveau des statuts
//          lesCondLim->Initial(lesMail,lesRef,lesCourbes1D,lesFonctionsnD,true,cas_combi_ddl);

          // on valide l'activité des conditions limites et condition linéaires, pour le temps initial
          // en conformité avec les conditions lues (qui peuvent éventuellement changé / aux calcul qui a donné le .BI)
          lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());

				      li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb,lesRef,charge->Temps_courant(),icas);
				      int ttsi = li_gene_asso.size();
				      X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);

//// --- debug
// cout << "\n debug1 AlgoriNewmark::initial( ";
// cout << "\n li_gene_asso.size()= " << li_gene_asso.size() << endl ;
//// --- fin debug


			       // on remet à jour les éléments pour le contact s'il y du contact présumé
			       if (lesMail->NbEsclave() != 0)
           {// on met à jour la boite d'encombrement compte tenue des nouvelles coordonnées
            lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_t);
            bool nevez_contact = lesContacts->Actualisation(0);
            // mise à jour éventuelle de la matrice de raideur en fonction du contact
            //bool changement_sur_matrice =
            Gestion_stockage_et_renumerotation_avec_contact
                             (premier_calcul,lesMail,nevez_contact,lesCondLim
                             ,lesRef,tab_mato,Ass1.Nb_cas_assemb(),lesContacts,niveau_substitution);
            matglob=tab_mato(1);
           }
          else // sinon il faut éventuellement tenir compte du cht du au CLL
            // NB: c'est pris en compte directement dans le contact, quand celui-ci est actif
           {bool nouvelle_situation_CLL=false;
            Gestion_stockage_et_renumerotation_sans_contact
                                (lesContacts,premier_calcul,lesMail,nouvelle_situation_CLL
                                 ,lesCondLim,lesRef,tab_mato,Ass1.Nb_cas_assemb()
                                 ,niveau_substitution);
           };
        };
    // vérif de cohérence pour le contact
    if ((pa.ContactType()) && (lesMail->NbEsclave() == 0)) // là pb
      {cout << "\n *** erreur: il n'y a pas de maillage disponible pour le contact "
            << " la definition d'un type contact possible est donc incoherente "
            << " revoir la mise en donnees !! "<< flush;
       Sortie(1);
      };
    // on regarde s'il y a besoin de sauvegarde
    if (this->Active_sauvegarde() && (ParaGlob::param->TypeCalcul_maitre() == this->typeCalcul) )
     { // si le fichier base_info n'est pas en service on l'ouvre
       entreePrinc->Ouverture_base_info("ecriture");
       // dans le cas ou se n'est pas un restart on sauvegarde l'incrément actuel
       // c'est-à-dire le premier incrément
       // après s'être positionné au début du fichier
       if (this->Num_restart() == 0)
        { entreePrinc->Sort_BI_Positionnement_offset(entreePrinc->position_debut_fichier);
          int cas = 1;
          paraGlob->Ecriture_base_info(*(entreePrinc->Sort_BI()),cas);
          this->Ecriture_base_info
             (cas,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage
              ,charge,lesCondLim,lesContacts,resultats,OrdreVisu::INCRE_0);
         } 
       else
        { // sinon on se place dans le fichier à la position du restart
          // debut_increment a été définit dans algori (classe mère)
          entreePrinc->Sort_BI_Positionnement_offset(debut_increment);
        };
     };
    //--fin cas de restart et/ou de sauvegarde--------

    // ajout d'un conteneur pour les coordonnées à l'itération 0
    {Coordonnee coor(ParaGlob::Dimension()); // un type coordonnee typique
     Grandeur_coordonnee gt(coor); // une grandeur typique de type Grandeur_coordonnee
     // def d'un type quelconque représentatif  à chaque noeud
     TypeQuelconque typQ_gene_int(XI_ITER_0,X1,gt);
     lesMail->AjoutConteneurAuNoeud(typQ_gene_int);
    };

    // choix de la matrice de masse, qui est en fait celle qui correspond au ddl Xi
    mat_masse = Choix_matrice_masse(nbddl_X,mat_masse,lesMail,lesRef
                                     ,Ass1.Nb_cas_assemb(),lesContacts,lesCondLim);

    // ici on utilise un nom différent mais qui pointe sur la même matrice , ceci uniquement pour mettre en évidence
    // que ce n'est pas la matrice masse mais c'est une matrice modifiée ( *= (coef_masse * unSurBetaDeltaTcarre))
    Mat_abstraite& mat_mas_modif = *mat_masse;
    // on calcul la matrice de masse qui est supposée identique dans le temps
    // c'est-à-dire que l'on considère que la masse volumique est constante
    Cal_matrice_masse(lesMail,Ass1,mat_mas_modif,diversStockage,lesRef,X1,lesFonctionsnD);
    // en fait on se sert de la matrice masse divisée par beta delta t carré 
    // d'où le nom de matrice masse modifiée avec éventuellement un coef qui vient de l'amortissement numérique
    mat_mas_modif *= (coef_masse * unSurBetaDeltaTcarre);
    if (ParaGlob::NiveauImpression()>3)
       {cout << "\n init: matrice masse coefficientee initiale : "; mat_mas_modif.MinMaxMoy(true);};
    type_incre = OrdreVisu::PREMIER_INCRE; // pour la visualisation au fil du calcul
	 
    if (amortissement_cinetique)		
        Algori::InitialiseAmortissementCinetique(); // initialisation des compteurs pour l'amortissement au cas ou
//// --- debug
// cout << "\n debug2 AlgoriNewmark::initial( ";
// cout << "\n li_gene_asso.size()= " << li_gene_asso.size() << endl ;
//// --- fin debug
    // mise à jour au cas où
    Algori::MiseAJourAlgoMere(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp,diversStockage
                           ,charge,lesCondLim,lesContacts,resultats);
    tempsInitialisation.Arret_du_comptage(); // temps cpu
  };
  
 
    // mise à jour
void AlgoriNewmark::MiseAJourAlgo
              (ParaGlob * paraGlob,LesMaillages * lesMail
               ,LesReferences* lesRef,LesCourbes1D* lesCourbes1D
               ,LesFonctions_nD* lesFonctionsnD
               ,VariablesExporter* varExpor
               ,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage
               ,Charge* charge,LesCondLim* lesCondLim,LesContacts* lescontacts
               ,Resultats* resultats
              )
  { // INITIALISATION globale
    tempsMiseAjourAlgo.Mise_en_route_du_comptage(); // temps cpu
    Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(DYNA_IMP); // transfert info
    
    // activation des ddl
    lesMail->Inactive_ddl(); // on commence par inactiver tous les ddl
    lesMail->Active_un_type_ddl_particulier(X1); // puis on active les ddl qu'ils faut ici
    lesMail->Active_un_type_ddl_particulier(V1);    //       on active les Vi
    lesMail->Active_un_type_ddl_particulier(GAMMA1);    //       on active les GAMMAi


    // mise à jour au cas où
    Algori::MiseAJourAlgoMere(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,varExpor,lesLoisDeComp
                ,diversStockage,charge,lesCondLim,lescontacts,resultats);
    tempsMiseAjourAlgo.Arret_du_comptage(); // temps cpu
  };

// Calcul de l'équilibre de la pièce en dynamique implicite selon la méthode de newmark
    // calcul
// si tb_combiner est non null -> un tableau de 2 fonctions
//  - la première fct dit si on doit valider ou non le calcul à convergence ok,
//  - la seconde dit si on doit sortir de la boucle ou non à convergence ok
//
// si la validation est effectuée, la sauvegarde pour le post-traitement est également effectuée
//                                 en fonction de la demande de sauvegard,
// sinon pas de sauvegarde pour le post-traitement à moins que l'on a demandé un mode debug
//   qui lui fonctionne indépendamment
void AlgoriNewmark::CalEquilibre(ParaGlob * paraGlob,LesMaillages * lesMail
           ,LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
           ,VariablesExporter* varExpor
           ,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage
           ,Charge* charge,LesCondLim* lesCondLim,LesContacts* lescontacts
           ,Resultats* resultats,Tableau < Fonction_nD* > * tb_combiner)
{
    tempsCalEquilibre.Mise_en_route_du_comptage(); // temps cpu
    Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(DYNA_IMP); // transfert info
    
    // préparation pour les aspects validation du calcul et sortie contrôlée des incréments
    int validation_calcul = 1; // init : par défaut la validation est effective si le calcul converge
    int sortie_boucle_controler_par_fctnD = 0; // init : par défaut la sortie n'est pas contrôler
    
    // on vérifie que le pas de temps est correcte: on propose le deltat actuel qui sera modifié s'il n'est pas bon
    double deltat_actuel = pa.Deltat();
    pa.Modif_Detat_dans_borne(deltat_actuel);
    // on met à jour des variables de travail
//    Assemblage& Ass = *Ass_;
    //Vecteur & vglobal = vglobin; // puissance totale qui ecrase vglobin
    int niveau_substitution = 1; // par défaut on utilise la matrice de raideur matglob = tab_mato(1)
    // ici on utilise un nom différent mais qui pointe sur la même matrice , ceci uniquement pour mettre en évidence
	   // que ce n'est pas la matrice masse mais c'est une matrice modifiée ( *= (coef_masse * unSurBetaDeltaTcarre))
	   Mat_abstraite& mat_mas_modif = *mat_masse;  // il s'agit d'une référence
    Coordonnee3 moy_diag_masse;// init
    if (ParaGlob::NiveauImpression()>3)
       {cout << "\n matrice masse coefficientee initiale : "; moy_diag_masse=mat_mas_modif.MinMaxMoy(true);}
    else if (avec_masse_scaling) {moy_diag_masse=mat_mas_modif.MinMaxMoy(false);};
    Coordonnee3  moy_diag_K; // sert aussi pour le masse scaling
    Assemblage& Ass1 = *Ass1_;    Assemblage& Ass3 = *Ass3_;
    Vecteur vglobal_sauv; // vecteur de sauvegarde éventuel, dans le cas de l'utilisation de matrices de raideur secondaire
    double max_delta_X=0.; //  le maxi du delta X
    double max_var_delta_X=0.; // idem d'une itération à l'autre

    bool HHT = false; if ((*hht) != 0.) { HHT = true;};
    
    // def d'un type générique, utilisé pour le transfert des forces internes, vers les conteneurs noeuds
    int dim = ParaGlob::Dimension();if (ParaGlob::AxiSymetrie()) dim--;
    Coordonnee coor(dim); // un type coordonnee typique
    Grandeur_coordonnee gt(coor); // une grandeur typique de type Grandeur_coordonnee
    // def d'un type quelconque représentatif  pour un vecteur force à chaque noeud
    TypeQuelconque typQ_gene_int(FORCE_GENE_INT,X1,gt);
    
    // on définit un type générique qui sert pour passer aux noeuds les positions à l'itération 0
    TypeQuelconque typQ_XI_ITER_0(XI_ITER_0,X1,gt);

    // tant que la fin du chargement n'est pas atteinte
    // dans le cas du premier chargement on calcul de toute manière, ce qui permet
    // de calculer meme si l'utilisateur indique un increment de charge supérieur
    // au temps final
    bool premier_calcul = true; // utilisé pour la méthode hht entre autre
    int indicCycleContact = 0; // def d'un indicateur donnant la situation dans le  cycle de contact   
    bool arret_pilotage=false; // pour arrêt du calcul au niveau du pilotage
    // un booléen pour uniquement gérer le fait que dans la boucle globale on fait le test après le test du while
    // d'où pour le temps final, on peut sortir de la boucle globale sans avoir convergé, ce qui n'est pas bon
    bool pas_de_convergence_pour_l_instant=1;
    icharge++; // on incrémente le chargement -> donne le num d'incr du prochain incr chargé
    while (   ((!charge->Fin(icharge,!pas_de_convergence_pour_l_instant)) // on n'affiche la sin
                    // que si on a eu convergence
               || pas_de_convergence_pour_l_instant ||(icharge == 1)
              )
           && (charge->Fin(icharge,true)!=2) // si on a dépassé le nombre d'incrément permis on s'arrête dans tous les cas
           && (charge->Fin(icharge,false)!=3) // idem si on a dépassé le nombre d'essai d'incrément permis
                                               // 1er appel avec true: pour affichage et second avec false car c'est déjà affiché
           && (!pa.EtatSortieEquilibreGlobal())
          )
    { double maxPuissExt;   // maxi de la puissance des efforts externes
      double maxPuissInt;   // maxi de la puissance des efforts internes 
      double maxReaction;   // maxi des reactions 
      int inReaction = 0;   // pointeur d'assemblage pour le maxi de reaction
      int inSol =0 ;        // pointeur d'assemblage du maxi de variation de ddl      
      double maxDeltaDdl=0; // // maxi de variation de ddl
      double last_var_ddl_max=0;  // maxi var ddl juste après résolutuion, que l'on sauvegarde d'une it à l'autre
      // initialisation de la variable puissance_précédente d'une itération à l'autre
      double puis_precedente = 0.; 
      bool change_statut = false;  // init des changements de statut,       
      Transfert_ParaGlob_NORME_CONVERGENCE(ConstMath::grand);// on met une norme grande par défaut au début
      // boleen pour savoir si on doit faire un affichage
      bool aff_incr=pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde); // pour simplifier
  
      // mise à jour du calcul éventuel des normales aux noeuds -> mise à jour des normales à t
      // mais ici, on calcule les normales à tdt, et on transfert à t
      // comme on est au début de l'incrément, la géométrie à tdt est identique à celle à t
      // sauf "au premier incrément", si l'algo est un sous algo d'un algo combiné
      // et que l'on suit un précédent algo sur un même pas de temps
      // qui a aboutit à une géométrie à tdt différente de celle de t
      // du coup cela permet d'utiliser la nouvelle géométrie pour ce premier incrément
      lesMail->MiseAjourNormaleAuxNoeuds_de_tdt_vers_T();
      // passage aux noeuds  des vecteurs globaux:  F_INT, F_EXT
      Algori::Passage_aux_noeuds_F_int_t_et_F_ext_t(lesMail);

      // premier contact de l'incrément
		    if (!indicCycleContact) // modification de la charge et du pas de temps qu'au premier passage
                              // mais pas après un changement de statut
		     { // appel du Pilotage
         int sauve_indic_convergence = PhaseDeConvergence(); // on sauvegarde car ce sera
           // changé dans le pilotage qui suit
         bool modif_pas_de_temps = Pilotage_du_temps(charge,arret_pilotage);
         if (arret_pilotage) break; // pilotage -> arret du calcul 
         // gestion du pas de temps et mise à jour des variables associées au pas de temps si besoin
         if (modif_pas_de_temps)
           { this->Gestion_pas_de_temps(lesMail,2); // 2 signifie cas courant
             // mise à jour de la matrice de masse en fonction de la nouvelle
             // valeur de unSurBetaDeltaTcarre et éventuellement de l'amortissement artificielle
             mat_mas_modif *= (coef_masse * unSurBetaDeltaTcarre/unSurBetaDeltaTcarre_o ); 
             if (ParaGlob::NiveauImpression()>3)
                {cout << "\n modif: matrice masse coefficientee : "; moy_diag_masse=mat_mas_modif.MinMaxMoy(true);}
             else if (avec_masse_scaling) {moy_diag_masse=mat_mas_modif.MinMaxMoy(false);};
           }
         else if (avec_masse_scaling) {moy_diag_masse=mat_mas_modif.MinMaxMoy(false);};

		       if (aff_incr)
           {cout << "\n======================================================================"
                 << "\nINCREMENT DE CHARGE : " << icharge 
                 << "  intensite " << charge->IntensiteCharge() 
                 << "  t= " << charge->Temps_courant()
                 << " dt= " << ParaGlob::Variables_de_temps().IncreTempsCourant()
                 << "\n======================================================================";
           };    
         // dans le cas d'un masse scaling avec pilotage, on modifie la masse en conséquence
         // sauf s'il y a eu une modification du pas de temps
         if (avec_masse_scaling)
            Pilotage_masse_scaling(moy_diag_masse(3),moy_diag_K(3),sauve_indic_convergence);

         // -- initialisation des coordonnees et des ddl a tdt en fonctions des
         // ddl imposes et de l'increment du chargement: change_statut sera recalculé ensuite
         lesCondLim->MiseAJour_tdt
            (pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
            ,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
 
 		      // dans le cas où il y a changement de statut des blocages on met à jour
			      // les tableaux d'indices généraux des ddl bloqués, y compris les ddls associés
			      if (change_statut)
				      {li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb
																,lesRef,charge->Temps_courant(),icas);
				       int ttsi = li_gene_asso.size();
				       X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
				      };                                               
         // --- calcul des éléments de contact: (correspond à la définition de la surface  de contact)
			      //   definition ou mise à jour, des elements de contact eventuels
         //  - imposition (en fonction de l'algorithme de contact) de conditions particulières de penetration (nulle par exemple)
         { bool changement_sur_matrice = false; // init par défaut
           int niveau_substitution = 0; // on intervient sur toutes les matrices
           if (pa.ContactType()) // traitement différent lors du premier incément calculé (premier passage) puis des passages courants
             {bool a_changer = false; // init par défaut
              double max_abs_deltaX = delta_X.Max_val_abs();
              if (premier_calcul)
                {a_changer = lescontacts->DefElemCont(max_abs_deltaX);}
              else
                { // on supprime les éléments inactifs testés à l'incr prec dans Actualisation()
                  a_changer = lescontacts->SuppressionDefinitiveElemInactif();
                  a_changer = a_changer || lescontacts->Nouveau(max_abs_deltaX);
                };
              if (premier_calcul || a_changer)
                {changement_sur_matrice = Gestion_stockage_et_renumerotation_avec_contact
                            (premier_calcul,lesMail,a_changer,lesCondLim,lesRef
                             ,tab_mato,Ass1.Nb_cas_assemb(),lescontacts,niveau_substitution);
                 matglob=tab_mato(1);
                };
             }
           else // sinon il faut éventuellement tenir compte du cht du au CLL
              // NB: c'est pris en compte directement dans le contact, quand celui-ci est actif
             {bool nouvelle_situation_CLL=false;
              changement_sur_matrice = Gestion_stockage_et_renumerotation_sans_contact(lescontacts,premier_calcul,lesMail,nouvelle_situation_CLL
                                   ,lesCondLim,lesRef,tab_mato,Ass1.Nb_cas_assemb()
                                   ,niveau_substitution);
             };
           // s'il y a eu un changement de numérotation, les pointeurs d'assemblage dans la matrice
           // masse ne sont plus corrects, on choisit de recalculer la matrice masse (par simplicité)
           if (changement_sur_matrice)
             {Cal_matrice_masse(lesMail,Ass1,mat_mas_modif,diversStockage,lesRef,X1,lesFonctionsnD);
              mat_mas_modif *= (coef_masse * unSurBetaDeltaTcarre);
              if (ParaGlob::NiveauImpression()>3)
               {cout << "\n change pointeur assemblage -> matrice masse coefficientee mise a jour : "; mat_mas_modif.MinMaxMoy(true);};
             };
         };// fin encapsulation du changement éventuel sur les stockages matricielles

       }		
      else
		     {if (arret_pilotage) break; // pilotage -> arret du calcul                
		       cout << "\n============================================================================="
              << "\n ....... re-analyse du contact ........ "
              << "\nINCREMENT DE CHARGE : " << icharge << "  intensite " << charge->IntensiteCharge() 
              << "  t= " << charge->Temps_courant()
              << " dt= " << ParaGlob::Variables_de_temps().IncreTempsCourant()
              << "\n============================================================================";
		     };		 

      lesLoisDeComp->MiseAJour_umat_nbincr(icharge); // init pour les lois Umat éventuelles

      // dans le cas où il y a changement de statut des blocages on met à jour
      // les tableaux d'indices généraux des ddl bloqués, y compris les ddls associés
/// *** dans le cas de contact qui impose la position de ddl, il faudrait en tenir compte en dynamique, au niveau des statuts ***
/// *** donc pour l'instant, ce n'est pas pris en comptes !!! ***

      if (change_statut)
         {li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb
                                             ,lesRef,charge->Temps_courant(),icas);
          int ttsi = li_gene_asso.size();
          X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
          };
      if (pa.ContactType())
		      {Mise_a_jour_Choix_matriciel_contact(tab_mato,Ass1.Nb_cas_assemb(),lescontacts,niveau_substitution);
         matglob=tab_mato(1);
        };

      // récupération (initialisation) des ddl  position, vitesse et accélération 
      lesMail->Vect_loc_vers_glob(TEMPS_t,X1,X_t,X1,Ass1.Nb_cas_assemb());
      lesMail->Vect_loc_vers_glob(TEMPS_t,V1,vitesse_t,V1,Ass1.Nb_cas_assemb());
      lesMail->Vect_loc_vers_glob(TEMPS_t,GAMMA1,acceleration_t,GAMMA1,Ass1.Nb_cas_assemb());
      // récupération au niveau global des ddl locaux à tdt avec conditions limite 
      // pour le vecteur accélération, seules les ddl avec CL sont différents de la précédente
      // récupération
      lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1,Ass1.Nb_cas_assemb());
      lesMail->Vect_loc_vers_glob(TEMPS_tdt,V1,vitesse_tdt,V1,Ass1.Nb_cas_assemb());
                                       
      //  boucle de convergence sur un increment
      Vecteur * sol; // pointeur du vecteur solution
      int compteur; // déclaré à l'extérieur de la boucle car utilisé après la boucle
		    bool arret_convergence = false;
      
      // on démarre avec le compteur à 0 et on sauvegarde la position à t=0
      lesMail->Quelconque_glob_vers_local(X1,X_tdt,typQ_XI_ITER_0,Ass1.Nb_cas_assemb());

      for (compteur = 0; (compteur<= pa.Iterations())&&(!pa.EtatSortieEquilibreGlobal()); compteur++)
//---//\\//\\// début de la boucle sur les itérations d'équilibres //\\//\\//
        {// initialisation de la matrice et du second membre
         matglob->Initialise (0.);
         vglobin.Zero();
         vglobex.Zero();
         if (pa.ContactType())
            vcontact.Zero();
         vglobaal.Zero(); // puissance totale
         Transfert_ParaGlob_COMPTEUR_ITERATION_ALGO_GLOBAL(compteur);
         
         // prise en compte du cas particulier ou l'utilisateur demande
         // une mise à jour des conditions limites à chaque itération
         if (cL_a_chaque_iteration)
           {// -- initialisation des coordonnees et des ddl a tdt en fonctions des
            // ddl imposes et de l'increment du chargement: change_statut sera recalculé ensuite
            lesCondLim->MiseAJour_tdt
                   (pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
                   ,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
            // dans le cas où il y a changement de statut des blocages on met à jour
            // les tableaux d'indices généraux des ddl bloqués, y compris les ddls associés
            if (change_statut)
             {li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb
                   ,lesRef,charge->Temps_courant(),icas);
              int ttsi = li_gene_asso.size();
              X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
             };
            // on récupère les nouvelles positions globalement de manière à pouvoir calculer le delta_X pour le contact
            lesMail->Vect_loc_vers_glob(TEMPS_tdt,X1,X_tdt,X1,Ass1.Nb_cas_assemb());
         //        delta_X = X_tdt; delta_X -= X_t; // X_tdt - X_t
            Algori::Cal_Transfert_delta_et_var_X(max_delta_X,max_var_delta_X);
           };
         
         // renseigne les variables définies par l'utilisateur via les valeurs déjà calculées par Herezh
         Algori::Passage_de_grandeurs_globales_vers_noeuds_pour_variables_globales(lesMail,varExpor,Ass1.Nb_cas_assemb(),*lesRef);
         varExpor->RenseigneVarUtilisateur(*lesMail,*lesRef);
         lesMail->CalStatistique(); // calcul éventuel de statistiques
         // affichage ou non de l'itération
         bool aff_iteration = (pa.freq_affich_iter() > 0) ?
                (aff_incr && (compteur % pa.freq_affich_iter()==0) &&(compteur!=0)) : false ;
		       lesMail->Force_Ddl_aux_noeuds_a_une_valeur(R_X1,0.0,TEMPS_tdt,true); // mise à 0 des ddl de réactions, qui sont uniquement des sorties
		       lesMail->Force_Ddl_etendu_aux_noeuds_a_zero(Ddl_enum_etendu::Tab_FN_FT()); // idem pour les composantes normales et tangentielles
         //  mise en place des conditions linéaires (ne comprend pas les conditions linéaires de contact éventuelles)
         lesCondLim->MiseAJour_condilineaire_tdt      
              (pa.Multiplicateur(),lesMail,charge->Increment_de_Temps(),lesRef,charge->Temps_courant()
               ,lesCourbes1D,lesFonctionsnD,charge->MultCharge(),change_statut,cas_combi_ddl);
         // il faut éventuellement tenir compte du cht du au CLL
         // NB: c'est pris en compte directement dans le contact, quand celui-ci est actif
         {bool premier_calcul=false;
          bool changement_sur_matrice = Gestion_stockage_et_renumerotation_sans_contact
                                   (lescontacts,premier_calcul,lesMail,change_statut
                                 ,lesCondLim,lesRef,tab_mato,Ass1.Nb_cas_assemb()
                                 ,niveau_substitution);
          if (changement_sur_matrice)
            {Cal_matrice_masse(lesMail,Ass1,mat_mas_modif,diversStockage,lesRef,X1,lesFonctionsnD);
             mat_mas_modif *= (coef_masse * unSurBetaDeltaTcarre);
             if ((aff_iteration)&&(ParaGlob::NiveauImpression()>3))
               {cout << "\n change pointeur assemblage -> matrice masse coefficientee mise a jour : "; mat_mas_modif.MinMaxMoy(true);};
            };
         };
         
         // utilisation d'un comportement tangent simplifié si nécessaire
         if (compteur <= pa.Init_comp_tangent_simple() )
           lesLoisDeComp->Loi_simplifie(true);
         else 
           lesLoisDeComp->Loi_simplifie(false);
         // on récupère les accélérations calculées au pas précédent
         lesMail->Vect_loc_vers_glob(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1,Ass1.Nb_cas_assemb());
         // les autres grandeurs n'ont pas changées depuis la fin de la boucle précédente
         
         // maintenant on met les conditions limites sur les ddls bloqués secondaires c-a-d associés
         // aux ddl bloqués par l'utilisateur, leur calcul dépend de l'algorithme d'où un calcul global
         list <LesCondLim::Gene_asso>::iterator ie,iefin=li_gene_asso.end();; // def d'un iterator adoc
         int ih=1; // indice
         for(ie=li_gene_asso.begin(),ih=1;ie!=iefin;ie++,ih++)
           // comme les valeurs des X V Gamma vont être écrasé par le calcul global, on utilise
           // des conteneurs intermédiaires
           {//trois cas
            LesCondLim::Gene_asso & s = (*ie); // pour simplifier
            int ix=s.pointe(1); // début des Xi 
            int iv=s.pointe(2); // début des Vi 
            int ig=s.pointe(3); // début des gammai
            if (PremierDdlFamille(s.ty_prin) == X1) 
            // cas ou les Xi sont imposés, on calcul Vi et Gammai
               { X_Bl(ih) = X_tdt(ix);
                 G_Bl(ih)= unSurBetaDeltaTcarre*(X_tdt(ix)-X_t(ix)) 
                                       - unsurbetadeltat * vitesse_t(iv)
                                       - unmoinsdeuxbetasurdeuxbeta * acceleration_t(ig);
                 V_Bl(ih) = vitesse_t(iv) + deltatUnMoinGamma * acceleration_t(iv)
                                            + deltatGamma * G_Bl(ih);                 
                }
            else if (PremierDdlFamille(s.ty_prin) == V1) 
            // cas ou les Vi sont imposés, calcul des Xi et Gammai
               { V_Bl(ih) = vitesse_tdt(iv);
                 G_Bl(ih) = unSurGammaDeltat * (vitesse_tdt(iv)-vitesse_t(iv))
                                                           - unMoinGammasurgamma*acceleration_t(ig);
                 X_Bl(ih) = X_t(ix) + delta_t * vitesse_t(iv) 
                                + zero5deltatcarreUnMoinDeuxBeta * acceleration_t(ig)
                                + zero5deltatcarredeuxBeta * G_Bl(ih);
                }
            else if (PremierDdlFamille(s.ty_prin) == GAMMA1) 
            // cas ou les gammai sont imposés, calcul des Vi et Xi
               { G_Bl(ih) = acceleration_tdt(ig);
                 V_Bl(ih) = vitesse_t(iv) + deltatUnMoinGamma * acceleration_t(iv)
                                            + deltatGamma * acceleration_tdt(ig);
                 X_Bl(ih) = X_t(ix) + delta_t * vitesse_t(iv) 
                                + zero5deltatcarreUnMoinDeuxBeta * acceleration_t(ig)
                                + zero5deltatcarredeuxBeta * acceleration_tdt(ig);
                }
           };
            
         // -- calcul global des nouvelles vitesses et accélérations, valable en fait que pour
         // les ddl libres   
         // calcul du champ de vitesse pour les ddl, valable en fait pour les ddl sans CL
         vitesse_tdt = vitesse_t + deltatUnMoinGamma * acceleration_t
                                 + deltatGamma * acceleration_tdt;                                 
         // calcul du champ de position à t+dt des ddl, valable en fait pour les ddl sans CL
         X_tdt = X_t + delta_t * vitesse_t 
                  + zero5deltatcarreUnMoinDeuxBeta  * acceleration_t
                  + zero5deltatcarredeuxBeta * acceleration_tdt;

         // 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 travaille
         Pilotage_maxi_X_V(X_t,X_tdt,vitesse_t,vitesse_tdt);
         // -- maintenant on met réellement en place les CL a partir de la sauvegarde
         for(ie=li_gene_asso.begin(),ih=1;ie!=iefin;ie++,ih++)
           {LesCondLim::Gene_asso & s = (*ie); // pour simplifier
            int ix=s.pointe(1); // début des Xi 
            int iv=s.pointe(2); // début des Vi 
            int ig=s.pointe(3); // début des gammai
            X_tdt(ix) = X_Bl(ih);
            vitesse_tdt(iv) = V_Bl(ih);
            acceleration_tdt(ig) = G_Bl(ih); 
           };
         
         // passage des valeurs calculées aux niveaux des maillages
         // on met ainsi en place des conditions limites sur les ddls secondaires
         // pour les gammai, comme seules les gammis avec CL ont été modifié 
         // c'est idem au niveau local bien que tout le vecteur soit transporté
         lesMail->Vect_glob_vers_local(TEMPS_tdt,X1,X_tdt,X1,Ass1.Nb_cas_assemb());
         lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1,Ass1.Nb_cas_assemb());
         lesMail->Vect_glob_vers_local(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1,Ass1.Nb_cas_assemb());
         lesLoisDeComp->MiseAJour_umat_nbiter(compteur); // init pour les lois Umat éventuelles
         // actualisation des conditions de contact qui peuvent changer la largeur de bande, 
			      // quand un noeud glisse d'une facette sur une voisine, peut changer la position du noeud
			      // qui est projeté sur la facette dans le cas de l'algorithme cinématique
         { bool changement_sur_matrice = false; // init par défaut
           int niveau_substitution = 0; // on intervient sur toutes les matrices
           if (pa.ContactType())
             { int niveau_substitution = 0; // on intervient sur toutes les matrices
               bool a_changer = false; // init
               if (compteur != 0)
                  a_changer = lescontacts->Actualisation(1);
               // mise à jour éventuelle des matrices de raideur en fonction du contact
               if (a_changer)
                {//bool changement_sur_matrice =
                 Gestion_stockage_et_renumerotation_avec_contact
                            (premier_calcul,lesMail,a_changer,lesCondLim,lesRef
                             ,tab_mato,Ass1.Nb_cas_assemb(),lescontacts,niveau_substitution);
                 matglob=tab_mato(1);
                };
             }
           else // sinon il faut éventuellement tenir compte du cht du au CLL
              // NB: c'est pris en compte directement dans le contact, quand celui-ci est actif
             {bool nouvelle_situation_CLL=false;
              changement_sur_matrice = Gestion_stockage_et_renumerotation_sans_contact
                                   (lescontacts,premier_calcul,lesMail,nouvelle_situation_CLL
                                   ,lesCondLim,lesRef,tab_mato,Ass1.Nb_cas_assemb()
                                   ,niveau_substitution);
             };
           // s'il y a eu un changement de numérotation, les pointeurs d'assemblage dans la matrice
           // masse ne sont plus corrects, on choisit de recalculer la matrice masse (par simplicité)
           if (changement_sur_matrice)
             {Cal_matrice_masse(lesMail,Ass1,mat_mas_modif,diversStockage,lesRef,X1,lesFonctionsnD);
              mat_mas_modif *= (coef_masse * unSurBetaDeltaTcarre);
              if ((aff_iteration)&&(ParaGlob::NiveauImpression()>3))
               {cout << "\n change pointeur assemblage -> matrice masse coefficientee mise a jour : "; mat_mas_modif.MinMaxMoy(true);};
             };
         };// fin encapsulation du changement éventuel sur les stockages matricielles
         
         // calcul des puissances internes et externe a partir des nouvelles positions
         // des nouvelles vitesses et pour les nouvelles accélérations imposées par les CL
         // -- mise en place du chargement impose sur le second membre
         // -- et éventuellement sur la raideur en fonction de sur_raideur
         // mise en place du chargement impose, c-a-d calcul de la puissance externe
         // si pb on sort de la boucle
         if (!(charge->ChargeSMembreRaideur_Im_mecaSolid
           (Ass1,lesMail,lesRef,vglobex,(*matglob),assembMat,pa,lesCourbes1D,lesFonctionsnD)))
           { Change_PhaseDeConvergence(-10);break;};
         // -- appel du calcul de la raideur et du second membre, énergies
         // dans le cas d'un calcul inexploitable arrêt de la boucle
         if (!RaidSmEner(lesMail,Ass1,vglobin,(*matglob))) break;
         if (ParaGlob::NiveauImpression()>3)
           {cout << "\n raideur : ";moy_diag_K= matglob->MinMaxMoy(true);}
         else if (avec_masse_scaling) {moy_diag_K= matglob->MinMaxMoy(false);};
			      // cas des efforts de contact, contribution second membre et matrice,
			      // suivant le type de contact on utilise ou on recalcule les reactions de contact éventuelles (contact et frottement)
         if (pa.ContactType()) // et des énergies développées pendant le contact
            {// dans le cas où le calcul est inexploitable (pb dans le calcul) arrêt de la boucle
             if (!RaidSmEnerContact(lescontacts,Ass1,vcontact,(*matglob))) break;
             if (ParaGlob::NiveauImpression()>3)
               {cout << "\n raideur apres contact : "; matglob->MinMaxMoy(true);}
            };
         // calcul des maxi des puissances internes
         maxPuissInt = vglobin.Max_val_abs();
         F_int_tdt = vglobin; // sauvegarde des forces généralisées intérieures
         
         if (pa.ContactType())
              vglobex += vcontact;
         maxPuissExt = vglobex.Max_val_abs();
         F_ext_tdt = vglobex; // sauvegarde des forces généralisées extérieures
         // second membre total
         vglobaal += vglobex ;vglobaal += vglobin ;

         

//---- je ne sais pas s'il faut ou non mettre cette partie après l'introduction de l'amortissement numérique et/OU HHT ?????
         // calcul des reactions de contact pour les noeuds esclaves
         // dans le repere absolu ( pour la sortie des infos sur le contact)
         // et test s'il y a decollement de noeud en contact (pour les contacts actifs)
			// mais ici, il n'y a pas de modification des éléments de contact (elles ont été faites dans  lescontacts->Actualisation())
         bool decol=false;
         if (pa.ContactType()) 
			        lescontacts->CalculReaction(vglobin,decol,Ass1.Nb_cas_assemb(),aff_iteration);
         //  - definition éventuelle de conditions limites linéaires de contact, en fonction des éléments du contact existant à ce niveau 
			//  (certains contacts (par pénalisation par exemple) ne produise pas de conditions linéaires)
         list <Condilineaire>* listCondLine=NULL;
         if (pa.ContactType())
			          listCondLine = &(lescontacts->ConditionLin(Ass1.Nb_cas_assemb()));
//------ fin de la partie pour laquelle je m'interroge --------------------

         // dans le cas où l'on utilise de l'amortissement numérique, calcul de la part visqueuse dans SM et K
         if (pa.Amort_visco_artificielle())
           { if (Arret_A_Equilibre_Statique())  // si on veut un équilibre statique, on sauvegarde les forces statiques
				           (*vglob_stat) = (vglobaal);
			          if (pa.CoefRayleighMasse()!=0.)
              { forces_vis_num = (nuAphaPrimeunsurcoef_masse * betaDelta_t_2) * vitesse_tdt;
                vglobaal -= mat_mas_modif.Prod_vec_mat(forces_vis_num);
               }
             if (pa.CoefRayleighRaideur()!=0.)
              { forces_vis_num = nuBetaPrime * vitesse_tdt;
                vglobaal -= matglob->Prod_vec_mat(forces_vis_num);
                (*matglob) *= coef_raideur;
                if (ParaGlob::NiveauImpression()>3)
                  {cout << "\n raideur apres viscosite Rayleigh : "; matglob->MinMaxMoy(true);
                   cout << " fact augmentation: "<<coef_raideur;
                  }
               }
           };

         // atténuation HHT    
         if (HHT) // cas de la méthode hht
          { vglobal_n_inter =  vglobaal; // sauvegarde de vglobal
            if (!premier_calcul)
             { // pour le second membre : on a vglobal_n = (*hht)* vglobal_n_inter 
               vglobaal *= (1.+(*hht)); vglobaal -= vglobal_n;
               // pour la raideur
               (*matglob) *= (1.+(*hht));
               if (ParaGlob::NiveauImpression()>3)
                 {cout << "\n raideur apres hht : "; matglob->MinMaxMoy(true);}
              }
          };

         // prise en compte des forces inertielles
         // -- pour le second membre
         inter_tdt = acceleration_tdt * (unsurcoef_masse * betaDelta_t_2  ); // on tient compte que la masse est * 1/(beta Delta T2)
         vglobaal -=  mat_mas_modif.Prod_mat_vec(inter_tdt) ;
         // -- pour la raideur
         (*matglob) += mat_mas_modif;
         if (ParaGlob::NiveauImpression()>3)
           {cout << "\n raideur avec masse : "; matglob->MinMaxMoy(true);}
   
         // initialisation des sauvegardes sur second membre (uniquement pour les Xi)
         lesCondLim->InitSauve(Ass1.Nb_cas_assemb());

         //--- prise en compte des conditions linéaire 
         // on récupère les réactions avant changement de repère et calcul des torseurs de réaction
         lesCondLim->ReacAvantCHrepere(vglobaal,lesMail,lesRef,Ass1.Nb_cas_assemb(),cas_combi_ddl);
         // -->> expression de la raideur et du second membre dans un nouveau repere
         lesCondLim->CoLinCHrepere_int((*matglob),vglobaal,Ass1.Nb_cas_assemb(),vglob_stat);
         if (pa.ContactType()==1)
			        lesCondLim->CoLinCHrepere_ext((*matglob),vglobaal,*listCondLine,Ass1.Nb_cas_assemb(),vglob_stat);
         // sauvegarde des reactions pour les  ddl bloques (simplement)
// en fait ne sert à rien, car les réactions maxi sont calculées dans condlin, mais comme c'est peut-être un peu spécial ici on laisse
// mais dans les autres algos c'est supprimé !!!!!
         // sauvegarde des reactions aux ddl bloque (uniquement pour les Xi)
         lesCondLim->ReacApresCHrepere(vglobaal,lesMail,lesRef,Ass1.Nb_cas_assemb(),cas_combi_ddl);

         // prise en compte des conditions imposée sur la matrice et second membres (donc dans le nouveau repère)
         lesCondLim->ImposeConLimtdt(lesMail,lesRef,(*matglob),vglobaal,Ass1.Nb_cas_assemb()
			                            ,cas_combi_ddl,vglob_stat);
         // blocage de toutes les  conditions lineaires, quelque soit leur origines
         lesCondLim->CoLinBlocage((*matglob),vglobaal,Ass1.Nb_cas_assemb(),vglob_stat);
         // calcul du maxi des reactions
         maxReaction = lesCondLim->MaxEffort(inReaction,Ass1.Nb_cas_assemb());

         // sortie d'info sur l'increment concernant les réactions
         if (compteur != 0) 
           if (aff_iteration) InfoIncrementReac(lesMail,compteur,inReaction,maxReaction,Ass1.Nb_cas_assemb());
         // ---dans le cas du mode debug on sort éventuellement les infos au fil du calcul (un peu bricolé)
         if ((mode_debug > 0)||(pa.EtatSortieEtatActuelDansCVisu()))
           {bool a_sortir_debug=false; // pour éviter une division par 0 du test (compteur % mode_debug == 0)
            if (pa.EtatSortieEtatActuelDansCVisu()) {a_sortir_debug=true;}
            else if (compteur % mode_debug == 0) {a_sortir_debug=true;};
            if (a_sortir_debug)
         	  {// on passe les déplacements de tdt à t
         	   lesMail->Vect_glob_vers_local(TEMPS_t,X1,X_tdt,X1,Ass1.Nb_cas_assemb());
             tempsCalEquilibre.Arret_du_comptage(); // on arrête le compteur pour la sortie
             // visualisation éventuelle au fil du calcul: essentiellement la déformée
             VisuAuFilDuCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
                          ,lesCondLim,lescontacts,resultats,type_incre,(icharge*1000000+compteur));
             tempsCalEquilibre.Mise_en_route_du_comptage(); // on remet en route le compteur
         	   // on remet les choses dans l'ordre initial
         	   lesMail->Vect_glob_vers_local(TEMPS_t,X1,X_t,X1,Ass1.Nb_cas_assemb());
         	  };
           };
			      //--- fin du mode debug -------
         // calcul des énergies et affichage des balances
         Algori::Cal_Transfert_delta_et_var_X(max_delta_X,max_var_delta_X);
         CalEnergieAffichage(coef_masse * unSurBetaDeltaTcarre,vitesse_tdt,mat_mas_modif
                            ,delta_X,icharge,brestart,acceleration_tdt,forces_vis_num);
         if (permet_affichage > 3)
          { cout << "\n --- |max_var_DeltaDdl|= "<< max_var_delta_X
                 << " , |max_deltaDdl|= " << max_delta_X << flush;};

         if (icharge==1)// dans le cas du premier incrément on considère que la balance vaut l'énergie
            // cinétique initiale, car vu que l'on ne met pas de CL à t=0, E_cin_0 est difficile à calculer
            {E_cin_0 = E_cin_tdt - bilan_E + E_int_tdt - E_ext_tdt; 
            };

         // test de convergence sur un increment
         if (Convergence(aff_iteration,last_var_ddl_max,vglobaal,maxPuissExt,maxPuissInt,maxReaction,compteur,arret_convergence))
          { // sortie des itérations sauf si l'on est en loi simplifiée
            if (lesLoisDeComp->Test_loi_simplife() )
              // cas d'une loi simplifié on remet normal
              lesLoisDeComp->Loi_simplifie(false);
            else 
              // cas normal, arrêt
              break;
          }
         else if (arret_convergence) {break;} // cas ou la méthode Convergence() demande l'arret
         // sinon on continue    
         // pour le pilotage, sauvegarde du résidu
         if (pa.Line_search()) (*Vres) = vglobaal;
         // ici sol en fait = vecglob qui est ecrase par la resolution
         // ici l'inconnue est l'accélération 
         // la résolution nous donne la variation de ddl, qui est ensuite multipliée par 1/(beta delta t carré) 
         // pour avoir la variation d'accélération
         // la raideur du système est M/(beta deltaT2) + K, 
         
         // ---- resolution ----
         tempsResolSystemLineaire.Mise_en_route_du_comptage(); // temps cpu
         bool erreur_resolution_syst_lineaire = false; // init
         int nb_matrice_secondaire = tab_mato.Taille(); // = 1 par défaut, mais on peut en avoir d'autre
         int niveau_substitution = 1; // par défaut on utilise la matrice de raideur matglob = tab_mato(1)
         while (niveau_substitution <= nb_matrice_secondaire)
          { // on sauvegarde éventuellement le second membre
            if (nb_matrice_secondaire > 1) // cela veut dire que l'on est suceptible de faire plusieurs boucles
             { if (niveau_substitution == 1)
                 // init au premier passage
                 {vglobal_sauv = vglobaal;
                  matglob->Transfert_vers_mat(*tab_mato(niveau_substitution+1));
                 }
               else
                 {vglobaal = vglobal_sauv; // récup du second membre
                  // s'il y a une matrice de substitution supplémentaire, on sauvegarde 
                  if (nb_matrice_secondaire > niveau_substitution)
                     tab_mato(niveau_substitution)->Transfert_vers_mat
                                         (*tab_mato(niveau_substitution+1));
                  if (ParaGlob::NiveauImpression() >= 3)
                   { cout <<  "\n  nouvelle resolution du systeme avec la matrice "<<niveau_substitution
                          << " de substitution " << endl;
                   };
                 };
             };
            residu_final = vglobaal; // sauvegarde pour le post-traitement
            try
             {// ici sol en fait = vecglob qui est ecrase par la resolution
//              sol = &(matglob->Resol_systID
//                 (vglobal,pa.Tolerance(),pa.Nb_iter_nondirecte(),pa.Nb_vect_restart()));
              sol = &(tab_mato(niveau_substitution)->Resol_systID
                 (vglobaal,pa.Tolerance(),pa.Nb_iter_nondirecte(),pa.Nb_vect_restart()));
              (*sol) *= unSurBetaDeltaTcarre ;  // variation de l'accélération
              // retour des ddl dans les reperes generaux, dans le cas où ils ont ete modifie
               // par des conditions linéaires
              lesCondLim->RepInitiaux( *sol,Ass1.Nb_cas_assemb());
              // affichage éventuelle du vecteur solution : deltat_ddl     
              if (ParaGlob::NiveauImpression() >= 10) 
                { string entete = " affichage du vecteur solution (delta_ddl) ";
                  sol->Affichage_ecran(entete);
                };
   //---debug
   //cout << "\ndelta ddl "; for(int i=1;i<=3;i++) { cout << "\nnoe:"<<i<< " ";for (int j=1;j<=3;j++)  cout << (*sol)((i-1)*3+j)<<" ";};
   ////--fin-debug
   //      cout << "\n résolution "; sol->Affiche();  // --- pour le debug
              erreur_resolution_syst_lineaire = false;
              break; // on sort du while
             }
            catch (ErrResolve_system_lineaire excep)
                 // cas d'une d'une erreur survenue pendant la résolution
               { erreur_resolution_syst_lineaire = true; // on prépare l'arrêt
               } // car il faut néanmoins effacer les marques avant l'arrêt
            catch (ErrSortieFinale)
                 // cas d'une direction voulue vers la sortie
                 // on relance l'interuption pour le niveau supérieur
               { ErrSortieFinale toto;
                 tempsResolSystemLineaire.Arret_du_comptage(); // temps cpu
                 throw (toto);
               }
            catch ( ... )  
               { erreur_resolution_syst_lineaire = true; // on prépare l'arrêt
               }
           // il y a eu un pb de résolution
           niveau_substitution++;
          };
         tempsResolSystemLineaire.Arret_du_comptage(); // temps cpu

         // effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée
         lesCondLim->EffMarque();
         if (pa.ContactType()) lescontacts->EffMarque();
         
         if (erreur_resolution_syst_lineaire)
           {Change_PhaseDeConvergence(-9); // on signale une divergence due à la résolution
            break; // arrêt si on avait détecté une erreur à la résolution
           };

         // calcul du maxi de variation de ddl
         maxDeltaDdl = sol->Max_val_abs(inSol); last_var_ddl_max=maxDeltaDdl;       
         double maxDeltatDdl_signe=(*sol)(inSol); // récupération de la grandeur signée      
         // pilotage à chaque itération:  ramène: maxDeltaDdl,sol modifiés éventuellement et insol
         Pilotage_chaque_iteration(sol,maxDeltaDdl,compteur,lesMail);
         // sortie d'info sur l'increment concernant les variations de ddl 
         if ((aff_iteration)&&(ParaGlob::NiveauImpression() > 1))
            InfoIncrementDdl(lesMail,inSol,maxDeltatDdl_signe,Ass3.Nb_cas_assemb());
         
         // suite du pilotage
         // ------------------------------------------------------------
         // |on regarde s'il faut utiliser l'algorithme de line search |
         // ------------------------------------------------------------
         bool avec_line_search = false;
         if (pa.Line_search())
           { // la méthode de line_search incrémente la solution 
             // même s'il n'y a pas de mise en oeuvre de la méthode
             avec_line_search = Line_search1
                        (*sauve_deltadept,puis_precedente,*Vres,lesMail,sol 
                        ,compteur,*sauve_dept_a_tdt,charge,vglobex,Ass1,*v_travail,
                        lesCondLim,vglobaal,lesRef,vglobaal,Ass3.Nb_cas_assemb()
                        ,cas_combi_ddl,lesCourbes1D,lesFonctionsnD);
            }
         else      
           //cas ou l'on ne fait pas de line seach
           { // actualisation des ddl actifs a t+dt
             // à chaque incrément, ici il s'agit de l'accélération
             lesMail->PlusDelta_tdt(*sol,Ass3.Nb_cas_assemb()); 
           };
         
        }
//---//\\//\\// fin de la boucle sur les itérations d'équilibres //\\//\\//
      // effacement du marquage de ddl bloque du au conditions lineaire imposée par l'entrée
      // car ici on n'est pas passé par cette endroit si l'on a convergé
      lesCondLim->EffMarque();
		    if (pa.ContactType()) lescontacts->EffMarque();
      
      // mise à jour des indicateurs contrôlés par le tableau: *tb_combiner
      if (tb_combiner != NULL) // cas d'un contrôle via des fonctions nD
       {if ((*tb_combiner)(1) != NULL)
           validation_calcul = (*tb_combiner)(1)->Valeur_pour_variables_globales()(1);
        if ((*tb_combiner)(2) != NULL)
           sortie_boucle_controler_par_fctnD = (*tb_combiner)(2)->Valeur_pour_variables_globales()(1);
       };

      // gestion de la fin des itérations
      if ((!Pilotage_fin_iteration_implicite(compteur))||(pa.EtatSortieEquilibreGlobal()))
          { // cas d'une non convergence
            pas_de_convergence_pour_l_instant = 1;
            // comme on incrémente pas les ddl on remet cependant les ddl
            //  et les grandeurs actives de tdt aux valeurs de ceux de t
            lesMail->TversTdt();   
            lescontacts->TversTdt(); // idem pour les contacts éventuels
				        // on remet indicCycleContact à 0 pour que le pilotage sur le temps se refasse
				        indicCycleContact = 0;
            // ------ cas particulier où on a une divergence qui demande de remonter sur plus d'un incrément
            if (Algori::PhaseDeConvergence() == -8)
              { int nb_incr_en_arriere = 3; // !!!!! nombre actuellement arbitraire -> par la suite mettre dans les para
                if (Controle_retour_sur_un_increment_enregistre(nb_incr_en_arriere,icharge))
                  { // cas ou on a réussi à trouver un incrément sauvegardé adoc = maintenant icharge
                    int cas=2;
                    this->Lecture_base_info(cas ,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage
                                           ,charge,lesCondLim,lescontacts,resultats,icharge);
                    // comme les conditions limites cinématiques peuvent être différentes en restart
                    // on libére toutes les CL imposées éventuelles
                    lesMail->Libere_Ddl_representatifs_des_physiques(LIBRE);
                    lesMail->ChangeStatut(cas_combi_ddl,LIBRE);
                    // dans le cas d'un calcul axisymétrique on bloque le ddl 3
                    if (ParaGlob::AxiSymetrie())
                       lesMail->Inactive_un_ddl_particulier(X3);
                    // on valide l'activité des conditions limites et condition linéaires, pour le temps initial
                    // en conformité avec les conditions lues (qui peuvent éventuellement changé / aux calcul qui a donné le .BI)
                    lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());
                    li_gene_asso = lesCondLim->Tableau_indice (lesMail,t_assemb,lesRef,charge->Temps_courant(),icas);
                    int ttsi = li_gene_asso.size();
                    X_Bl.Change_taille(ttsi);V_Bl.Change_taille(ttsi);G_Bl.Change_taille(ttsi);
                    // on remet à jour les éléments pour le contact s'il y du contact présumé
                    if (pa.ContactType())
                       lesMail->Mise_a_jour_boite_encombrement_elem_front(TEMPS_t);
                              brestart = true; // on signale que l'on repars avec un restart
                    brestart = true; // on signale que l'on repars avec un restart
                  };
                // sinon on ne fait rien, on se contente du pilotage de divergence existant  
              };        
            // ------ fin cas particulier où on a une divergence qui demande de remonter sur plus d'un incrément
         }
      else   
        {// --- sinon calcul correcte
         pas_de_convergence_pour_l_instant = 0;
         if (permet_affichage > 0)
          { cout << "\n --- |max_var_DeltaDdl|= "<< max_var_delta_X
                 << " , |max_deltaDdl|= " << max_delta_X << flush;};
         if (pa.ContactType()) 
           { bool nevez_contact = lescontacts->Actualisation(0); // actualisation du contact en fonction du dernier incrément
             // réexamen du contact pour voir s'il n'y a pas de nouveau element de contact
			          // en fait on  fera au plus deux passages supplémentaire, sinon la boucle peut être infini,
				         // à la fin du second passage, on regarde s'il y a décollement, si oui on relâche et on refait un passage
				         // sinon on valide
				         //I)--------------
				         if (indicCycleContact == 0 )
				          { nevez_contact = nevez_contact || lescontacts->Nouveau(lesMail->Max_var_dep_t_a_tdt());
                if (nevez_contact)
                  {indicCycleContact=1;} // on a de nouveau contact on refait le deuxième cycle
					           else 
					             // sinon, on n'a pas de nouveau contact, on regarde s'il y a du relachement	 
					             { indicCycleContact=2;}; 
              }
				         else if (indicCycleContact == 1)  
				          {indicCycleContact=2;} // pour regarder le relachement
				         else                              
					         {indicCycleContact=0;}; // pas de newcontact, ni de relachement donc c'est ok
				         //II)---------------
             bool relachement_contact = false;
             if (indicCycleContact == 2)
              {relachement_contact = lescontacts->RelachementNoeudcolle();
               if (relachement_contact)
                {indicCycleContact=3;} // pour une dernière boucle d'équilibre
               else
                {indicCycleContact=0;};// pas de relachement donc ok
              };
            
             if (   (nevez_contact || relachement_contact)
                )
                 {int niveau_substitution = 0; // on intervient sur toutes les matrices
                  bool nouveau = (nevez_contact || relachement_contact);
                  bool changement_sur_matrice = Gestion_stockage_et_renumerotation_avec_contact
                          (premier_calcul,lesMail,nouveau,lesCondLim
                          ,lesRef,tab_mato,Ass1.Nb_cas_assemb(),lescontacts,niveau_substitution);
                  matglob=tab_mato(1);
                  // s'il y a eu un changement de numérotation, les pointeurs d'assemblage dans la matrice
                  // masse ne sont plus corrects, on choisit de recalculer la matrice masse (par simplicité)
                  if (changement_sur_matrice)
                    {Cal_matrice_masse(lesMail,Ass1,mat_mas_modif,diversStockage,lesRef,X1,lesFonctionsnD);
                     mat_mas_modif *= (coef_masse * unSurBetaDeltaTcarre);
                     if (ParaGlob::NiveauImpression()>3)
                      {cout << "\n change pointeur assemblage -> matrice masse coefficientee mise a jour : "; mat_mas_modif.MinMaxMoy(true);};
                    };
                 };
			        }
			      else
			        // concerne ici: soit le cas où c'est un problème sans contact (=contact non activé)
			        { indicCycleContact = 0; }; //pour le prochain increment 
			  
			      // cas où on n'a pas de contact, ou cas où le contact est ok, (indicCycleContact == 0), on valide
			      if (!(pa.ContactType()) || (indicCycleContact == 0))
			        { // calcul éventuelle de l'amortissement cinétique
			          int relax_vit_acce = AmortissementCinetique(delta_X,coef_masse * unSurBetaDeltaTcarre
                                                   ,X_tdt,mat_mas_modif,icharge,vitesse_tdt);
			          // s'il y a amortissement cinétique il faut re-updater les vitesses
			          if (Abs(relax_vit_acce) == 1) {lesMail->Vect_glob_vers_local(TEMPS_tdt,V1,vitesse_tdt,V1,Ass1.Nb_cas_assemb());}
             delta_prec_X = delta_X; // on sauvegarde le delta_X qui a convergé
             
             // si on est sans validation, on ne fait rien, sinon on valide l'incrément
             // avec sauvegarde éventuelle
             if (validation_calcul)
               {//  passage des accélérations calculées aux niveaux des maillages
                lesMail->Vect_glob_vers_local(TEMPS_tdt,GAMMA1,acceleration_tdt,GAMMA1,Ass1.Nb_cas_assemb());
                // actualisation des ddl et des grandeurs actives de t+dt vers t
                lesMail->TdtversT();
                lescontacts->TdtversT();
                Algori::TdtversT();
                // on valide l'activité des conditions limites  et condition linéaires
                lesCondLim->Validation_blocage (lesRef,charge->Temps_courant());
                // sauvegarde de l'incrément si nécessaire
                tempsCalEquilibre.Arret_du_comptage(); // on arrête le compteur pour la sortie
                Ecriture_base_info(2,lesMail,lesRef,lesCourbes1D,lesFonctionsnD
                ,lesLoisDeComp,diversStockage,charge,lesCondLim,lescontacts
                ,resultats,type_incre,icharge);
                // enregistrement du num d'incrément et du temps correspondant
                list_incre_temps_calculer.push_front(Entier_et_Double(icharge,pa.Variables_de_temps().TempsCourant()));
                // visualisation éventuelle au fil du calcul
                VisuAuFilDuCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
                             ,lesCondLim,lescontacts,resultats,type_incre,icharge);
                tempsCalEquilibre.Mise_en_route_du_comptage(); // on remet en route le compteur
               }; // -- fin du test: if (validation_calcul)
                
			          // cas de la méthode hht
			          premier_calcul = false;// -- on enregistre le fait que l'on n'est plus au premier passage
			          // -- sauvegarde de la puissance interne
			          if (HHT) vglobal_n = (*hht)*vglobal_n_inter;
			          //s'il y a remonté des sigma et/ou def aux noeuds et/ou calcul d'erreur
			          bool change =false;  // calcul que s'il y a eu initialisation
			          if(prepa_avec_remont) {change = Algori::CalculRemont(lesMail,type_incre,icharge);};
			          if (change) // dans le cas d'une remonté il faut réactiver les bon ddls
			          	{lesMail->Inactive_ddl(); lesMail->Active_un_type_ddl_particulier(tenuXVG);};
			          brestart = false; // dans le cas où l'on était en restart, on passe l'indicateur en cas courant
			          // test de fin de calcul effectue dans charge via : charge->Fin()
			          if (relax_vit_acce == -1)  // cas d'une relaxation cinétique active
//*** là je pense qu'il faudrait revoir le cas relaxation cinétique avec newmark .... a tester plus en détail !!				
			          	{arret_convergence=true;	};  // arrêt du calcul car deltaddl satisfait le critère de relaxation cinétique
             if (validation_calcul)
               {icharge++;
                Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
               };
			        }; //- fin du cas où l'incrément est validé pour l'équilibre avec le contact
        }; // -- fin du cas ou le calcul converge
      // cas particulier où la sortie de la boucle est pilotée
      if (sortie_boucle_controler_par_fctnD && (!pas_de_convergence_pour_l_instant))
        break;
    }; // -- fin du while sur les incréments de charge
    // on remet à jour le nombre d'incréments qui ont convergés:
    if (validation_calcul)
      {icharge--;
       Transfert_ParaGlob_COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL(icharge);
      }
    else // si on ne valide pas le calcul, on reinitialise la charge
     // c-a-d l'avancement en temps, incrément et facteur multiplicatif
     // de manière à avoir les mêmes conditions de départ pour le prochain calcul
      { charge->Precedant(true);}  ;
    tempsCalEquilibre.Arret_du_comptage(); // temps cpu    
};

// Calcul de l'équilibre de la pièce en dynamique implicite selon la méthode de newmark
void AlgoriNewmark::FinCalcul(ParaGlob * paraGlob,LesMaillages * lesMail,
               LesReferences* lesRef,LesCourbes1D* lesCourbes1D,LesFonctions_nD* lesFonctionsnD
               ,VariablesExporter* varExpor
               ,LesLoisDeComp* lesLoisDeComp,DiversStockage* diversStockage,
           Charge* charge,LesCondLim* lesCondLim,LesContacts* lesContacts
           ,Resultats* resultats)
  { 
    // passage finale dans le cas d'une visualisation  au fil du calcul  
    type_incre = OrdreVisu::DERNIER_INCRE;  
    Transfert_ParaGlob_ALGO_GLOBAL_ACTUEL(DYNA_IMP); // transfert info
    VisuAuFilDuCalcul(paraGlob,lesMail,lesRef,lesCourbes1D,lesFonctionsnD,lesLoisDeComp,diversStockage,charge
                          ,lesCondLim,lesContacts,resultats,type_incre,icharge);
    // enregistrement du num d'incrément et du temps correspondant                 
    list_incre_temps_calculer.push_front(Entier_et_Double(icharge,pa.Variables_de_temps().TempsCourant()));                 
    // sauvegarde de l'incrément si nécessaire
    Ecriture_base_info(2,lesMail,lesRef,lesCourbes1D,lesFonctionsnD
             ,lesLoisDeComp,diversStockage,charge,lesCondLim,lesContacts
             ,resultats,type_incre,icharge);
  };

