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

// les différents types de stockage matriciel
#include <complex>
#include "MatLapack.h"
#include "Mat_pleine.h"
#include "MatBand.h"
#include "Mat_creuse_CompCol.h"
#include "MatDiag.h"

#include "ReferenceNE.h"
#include "ReferenceAF.h"
#include "ExceptionsElemMeca.h"
#include "ExceptionsLoiComp.h"
#include "TypeQuelconqueParticulier.h"


//----------  class DeuxString    ---------
 // surcharge de l'operator de lecture /écriture
istream & operator >> (istream & ent, Algori::DeuxString & a)
  { //lecture du type et vérification
    string nom_type; ent >> nom_type;
    if (nom_type != "Algori::DeuxString") Sortie(1);
    // les deux strings
    ent >> a.nomFichier >>  a.nomMaillage;
    return ent;
   }; 
    
ostream & operator << (ostream & sort, const Algori::DeuxString & a)
         { // le type puis les deux strings
           sort << "Algori::DeuxString " << a.nomFichier << " " << a.nomMaillage <<"\n";
           return sort;
          };
           
istream & operator >> (istream & ent , Algori::ListDeuxString & a)
   { ent >> a.infos; return ent;};
ostream & operator << (ostream & sort , const Algori::ListDeuxString &a)
   { sort << a.infos; return sort;};
   
//---------- fin   class DeuxString  ---------

// CONSTRUCTEURS :
Algori::Algori () :  // par defaut
  paraTypeCalcul(),pa(),temps_derniere_sauvegarde(0.),mode_debug(0),permet_affichage(0)
  ,list_incre_temps_sauvegarder(),list_incre_temps_calculer()
  ,tempsDerniereSortieFilCalcul(0.),noms_fichier(),typeFlotExterne()
  ,entreePrinc(NULL),debut_increment(0),visualise(),visualise_maple()
  ,visualise_geomview(),visualise_Gid(),visualise_Gmsh(),toto_masse(NULL)
  ,icharge(0),casAssemb_principal(NULL)
  ,phase_de_convergence(0),nombre_de_bonnes_convergences(0),a_converge(false) // variables de travail
  ,a_converge_iterMoins1(false)
  ,nombre_de_mauvaises_convergences(0)
  ,max_var_residu(),nb_cycle_test_max_var_residu(0),v_int()        // variables de travail 
  ,tab_para_stockage() //,pt_icharge(NULL)
  ,E_cin_0(0.),E_cin_t(0.),E_cin_tdt(0),E_int_t(0),E_int_tdt(0),E_ext_t(0),E_ext_tdt(0),bilan_E(0) // différentes énergies et bilan
  ,q_mov_t(0.),q_mov_tdt(0.)
  ,P_acc_tdt(0.),P_int_tdt(0.),P_ext_tdt(0.),bilan_P_tdt(0.),P_acc_t(0.),P_int_t(0.),P_ext_t(0.),bilan_P_t(0.)
  ,E_visco_numerique_t(0.),E_visco_numerique_tdt(0.),E_bulk(0.),P_bulk(0.)
  ,F_int_t(),F_ext_t(),F_int_tdt(),F_ext_tdt(),F_totale_tdt() // forces généralisées int et ext à t et t+dt
  ,residu_final()

  ,E_cin_ind_t(),E_cin_ind_tdt(),E_int_ind_t(),E_int_ind_tdt(),E_ext_ind_t(),E_ext_ind_tdt(),bilan_ind_E()
  ,q_mov_ind_t(),q_mov_ind_tdt()
  ,P_acc_ind_tdt(),P_int_ind_tdt(),P_ext_ind_tdt(),bilan_P_ind_tdt()
  ,P_acc_ind_t(),P_int_ind_t(),P_ext_ind_t(),bilan_P_ind_t()
  ,E_visco_numerique_ind_t(),E_visco_numerique_ind_tdt()
  ,F_int_ind_t(),F_ext_ind_t(),F_int_ind_tdt(),F_ext_ind_tdt()
  ,X_t(),X_tdt(),delta_X(),var_delta_X(),delta_t_precedent_a_convergence(0.)
  ,vitesse_t(),vitesse_tdt(),acceleration_t(),acceleration_tdt()
  ,E_bulk_ind()
  ,P_bulk_ind()

  ,cas_ass_eps(NULL),Ass_eps(NULL),cas_ass_eps11(NULL),Ass_eps11(NULL)
  ,cas_ass_sig(NULL),Ass_sig(NULL),cas_ass_sig11(NULL),Ass_sig11(NULL)
  ,mat_glob_sigeps(NULL),vglob_sigeps(NULL),vglobal_sigeps(NULL)
  ,nbddl_sigeps11(0),nbddl_sigeps(0) 
  ,cas_ass_err(NULL),Ass_err(NULL),nbddl_err(0)             
  ,avec_remont_sigma(false),avec_remont_eps(false),avec_erreur(false)
  ,energTotal(),energHourglass(0.),energFrottement(),energPenalisation(0.)
  ,energStabiliMembBiel(0.)
  ,deja_lue_entete_parametre(0) 
  //-------- amortissement cinétique
  ,amortissement_cinetique(false),max_nb_decroit_pourRelaxDyn(1),compteur_decroit_pourRelaxDyn(0)
  ,inter_nb_entre_relax(1),fct_nD_inter_nb_entre_relax(NULL),nom_fct_nD_inter_nb_entre_relax("")
  ,nb_depuis_derniere_relax(0)
  ,taille_moyenne_glissante(1),moyenne_glissante(ConstMath::tresgrand), moy_gliss_t(0.)
  ,sauve_E_cin(),indice_EC1(0),nb_Ecin_stocker(0),test_fraction_energie(0.),X_EcinMax()
  ,coef_arret_pourRelaxDyn(0.5),dernier_pic(-1),pic_E_cint_t(-1.)
  ,max_pic_E_cin(-1.),coef_redemarrage_pourRelaxDyn(0.05)
  ,max_deltaX_pourRelaxDyn(0.1)
  ,nb_max_dX_OK_pourRelaxDyn(5),nb_dX_OK_pourRelaxDyn(0)
  ,nb_deb_testfin_pourRelaxDyn(100),nb_deb_test_amort_cinetique(1)
  ,compteur_test_pourRelaxDyn(0),compteur_pic_energie(0)
  // ...... the same pour l'amortissement local
  ,amortissement_cinetique_au_noeud(0),nb_deb_test_amort_cinetique_noe(100)
  ,E_cinNoe_tdt(),E_cinNoe_t(),max_nb_decroit_pourRelaxDyn_noe(1),compteur_decroit_pourRelaxDyn_noe()
  ,coef_arret_pourRelaxDyn_noe(0.5),dernier_pic_noe()
  ,pic_E_cint_t_noe(),max_pic_E_cin_noe(),coef_redemarrage_pourRelaxDyn_noe(0.05)
  ,compteur_pic_energie_noe(0)
  //------- amortissement visqueux critique
  ,opt_cal_C_critique(-1),f_(0.9),ampli_visco(1.)
  ,C_inter(NULL),type_calcul_visqu_critique(0)
  ,F_int_tdt_prec(),F_ext_tdt_prec()
  ,listeVarGlob(),listeVecGlob()

  //----- fin des amortissements ---
  ,arret_a_equilibre_statique(0),vglob_stat(NULL)
  ,volume_total_matiere(0.),vol_total2D_avec_plan_ref(1,Coordonnee(ParaGlob::Dimension()))
  // -- controle précision convergence
  ,nom_modulation_precision(""),modulation_precision(NULL),fct_norme(NULL)
  //---- grandeurs de travail privées ------
  ,vec_trav()
  //------- temps cpu -----
  ,tempsInitialisation(),tempsMiseAjourAlgo(),tempsCalEquilibre()
  ,tempsRaidSmEner(),tempsSecondMembreEnerg(),tempsResolSystemLineaire()
  ,tempsSauvegarde(),tempsSortieFilCalcul(),tempsRaidSmEnerContact()
  ,tempsSecondMembreEnergContact(),temps_CL(),temps_CLL()
  ,temps_lois_comportement(),temps_metrique_K_SM(),temps_chargement()
  ,temps_rech_contact(),tempsCalMasse(),tempsCalViscoNum()
  ,ptalgocombi(NULL)
 #ifndef UTILISATION_MPI
  ,lesTempsCpu(18)
 #else
  ,lesTempsCpu(30),distribution_CPU_algo()
  ,temps_transfert_court_algo(),temps_transfert_long_algo(),temps_attente_algo()
  ,temps_transfert_court_matSm(),temps_transfert_long_matSm(),temps_attente_matSm()
  ,temps_transfert_court_charge(),temps_transfert_long_charge(),temps_attente_charge()
  ,temps_transfert_court_contact(),temps_transfert_long_contact(),temps_attente_contact()
  ,inter_transfer(),inter_transfer2()
  ,buffer_ioBI_MPI(),tailles_buffer_ioBI_MPI(),posi_buffer()
 #endif
  //---------- stockage pour la transmission des grandeurs consultables -----
//  ,stock_compteur(GENERIQUE_UNE_GRANDEUR_GLOBALE),val_stock_compteur(NULL)
//  ,stock_icharge(GENERIQUE_UNE_GRANDEUR_GLOBALE),val_stock_icharge(NULL)

  {  cout << "\n  ce constructeur ne doit pas etre utilise !! \n" ;
     cout << " Algo::Algo () // par defaut " << endl;
     Sortie(1);
  };   
     
// constructeur en fonction du type de calcul
// il y a ici lecture des parametres attaches au type
Algori::Algori (EnumTypeCalcul type,const bool avec_typeDeCal
            ,const list <EnumSousTypeCalcul>& soustype
            ,const list <bool>& avec_soustypeDeCal
            ,UtilLecture& entreePrincipal
            ):
  typeCalcul(type),avec_typeDeCalcul(avec_typeDeCal)
  ,soustypeDeCalcul(&soustype),avec_soustypeDeCalcul(&avec_soustypeDeCal)
  ,paraTypeCalcul()
  ,pa(),temps_derniere_sauvegarde(0.),mode_debug(0),permet_affichage(0)
  ,list_incre_temps_sauvegarder(),list_incre_temps_calculer()
  ,tempsDerniereSortieFilCalcul(0.),noms_fichier(),typeFlotExterne()
  ,entreePrinc(&entreePrincipal),debut_increment(0)
  ,visualise(&entreePrincipal),visualise_maple(&entreePrincipal)
  ,visualise_geomview(&entreePrincipal),visualise_Gid(&entreePrincipal)
  ,visualise_Gmsh(&entreePrincipal),toto_masse(NULL)
  ,icharge(0),casAssemb_principal(NULL)
  ,phase_de_convergence(0),nombre_de_bonnes_convergences(0),a_converge(false) //  variables de travail
  ,a_converge_iterMoins1(false)
  ,nombre_de_mauvaises_convergences(0)
  ,max_var_residu(),nb_cycle_test_max_var_residu(0),v_int()        // variables de travail 
  ,tab_para_stockage()//,pt_icharge(NULL)
  ,E_cin_0(0.),E_cin_t(0),E_cin_tdt(0),E_int_t(0),E_int_tdt(0),E_ext_t(0),E_ext_tdt(0),bilan_E(0) // différentes énergies et bilan
  ,q_mov_t(0.),q_mov_tdt(0.)
  ,P_acc_tdt(0.),P_int_tdt(0.),P_ext_tdt(0.),bilan_P_tdt(0.),P_acc_t(0.),P_int_t(0.),P_ext_t(0.),bilan_P_t(0.)
  ,E_visco_numerique_t(0.),E_visco_numerique_tdt(0.),E_bulk(0.),P_bulk(0.)
  ,F_int_t(),F_ext_t(),F_int_tdt(),F_ext_tdt(),F_totale_tdt() // forces généralisées int et ext à t et t+dt
  ,residu_final()
  ,X_t(),X_tdt(),delta_X(),var_delta_X(),delta_t_precedent_a_convergence(0.)
  ,vitesse_t(),vitesse_tdt(),acceleration_t(),acceleration_tdt()
  ,cas_ass_eps(NULL),Ass_eps(NULL),cas_ass_eps11(NULL),Ass_eps11(NULL)
  ,cas_ass_sig(NULL),Ass_sig(NULL),cas_ass_sig11(NULL),Ass_sig11(NULL)
  ,mat_glob_sigeps(NULL),vglob_sigeps(NULL),vglobal_sigeps(NULL)               
  ,nbddl_sigeps11(0),nbddl_sigeps(0)
  ,cas_ass_err(NULL),Ass_err(NULL),nbddl_err(0)             
  ,avec_remont_sigma(false),avec_remont_eps(false),avec_erreur(false)
  ,energTotal(),energHourglass(0.),energFrottement(),energPenalisation(0.),energStabiliMembBiel(0.)
  ,deja_lue_entete_parametre(0)
  //-------- amortissement cinétique
  ,amortissement_cinetique(false),max_nb_decroit_pourRelaxDyn(1),compteur_decroit_pourRelaxDyn(0)
  ,inter_nb_entre_relax(1),fct_nD_inter_nb_entre_relax(NULL),nom_fct_nD_inter_nb_entre_relax("")
  ,nb_depuis_derniere_relax(0)
  ,taille_moyenne_glissante(1),moyenne_glissante(ConstMath::tresgrand), moy_gliss_t(0.)
  ,sauve_E_cin(),indice_EC1(0),nb_Ecin_stocker(0),test_fraction_energie(0.),X_EcinMax()
  ,coef_arret_pourRelaxDyn(0.5),dernier_pic(-1.),pic_E_cint_t(-1.),max_pic_E_cin(-1.),coef_redemarrage_pourRelaxDyn(0.05)
  ,max_deltaX_pourRelaxDyn(0.1),nb_max_dX_OK_pourRelaxDyn(5),nb_dX_OK_pourRelaxDyn(0)
  ,nb_deb_testfin_pourRelaxDyn(100),nb_deb_test_amort_cinetique(1)
  ,compteur_test_pourRelaxDyn(0),compteur_pic_energie(0)
  // ...... the same pour l'amortissement local
  ,amortissement_cinetique_au_noeud(0),nb_deb_test_amort_cinetique_noe(100)
  ,E_cinNoe_tdt(),E_cinNoe_t(),max_nb_decroit_pourRelaxDyn_noe(1),compteur_decroit_pourRelaxDyn_noe()
  ,coef_arret_pourRelaxDyn_noe(0.5),dernier_pic_noe()
  ,pic_E_cint_t_noe(),max_pic_E_cin_noe(),coef_redemarrage_pourRelaxDyn_noe(0.05)
  ,compteur_pic_energie_noe(0)
  //------- amortissement visqueux critique
  ,opt_cal_C_critique(-1),f_(0.9),ampli_visco(1.)
  ,C_inter(NULL),type_calcul_visqu_critique(0)
  ,F_int_tdt_prec(),F_ext_tdt_prec()
  ,listeVarGlob(),listeVecGlob()

  //----- fin des amortissements ---  
  ,arret_a_equilibre_statique(0),vglob_stat(NULL)
  ,volume_total_matiere(0.),vol_total2D_avec_plan_ref(1,Coordonnee(ParaGlob::Dimension()))
  // -- controle précision convergence
  ,nom_modulation_precision(""),modulation_precision(NULL),fct_norme(NULL)
  //---- grandeurs de travail privées ------
  ,vec_trav()
  //------- temps cpu -----
  ,tempsInitialisation(),tempsMiseAjourAlgo(),tempsCalEquilibre()
  ,tempsRaidSmEner(),tempsSecondMembreEnerg(),tempsResolSystemLineaire()
  ,tempsSauvegarde(),tempsSortieFilCalcul(),tempsRaidSmEnerContact()
  ,tempsSecondMembreEnergContact(),temps_CL(),temps_CLL()
  ,temps_lois_comportement(),temps_metrique_K_SM(),temps_chargement()
  ,temps_rech_contact(),tempsCalMasse(),tempsCalViscoNum()
  ,ptalgocombi(NULL)
  #ifndef UTILISATION_MPI
  ,lesTempsCpu(18)
 #else
  ,lesTempsCpu(30),distribution_CPU_algo()
  ,temps_transfert_court_algo(),temps_transfert_long_algo(),temps_attente_algo()
  ,temps_transfert_court_matSm(),temps_transfert_long_matSm(),temps_attente_matSm()
  ,temps_transfert_court_charge(),temps_transfert_long_charge(),temps_attente_charge()
  ,temps_transfert_court_contact(),temps_transfert_long_contact(),temps_attente_contact()
  ,inter_transfer(),inter_transfer2()
  ,buffer_ioBI_MPI(),tailles_buffer_ioBI_MPI(),posi_buffer()
 #endif

  //---------- stockage pour la transmission des grandeurs consultables -----
//  ,stock_compteur(GENERIQUE_UNE_GRANDEUR_GLOBALE),val_stock_compteur(NULL)
//  ,stock_icharge(GENERIQUE_UNE_GRANDEUR_GLOBALE),val_stock_icharge(NULL)

  { // on met à jour ParaGlob
    ParaGlob::param->ChangeParaAlgoControleActifs(&pa);
    
    // stockage pour la transmission des grandeurs consultables
    // 1) on définie un type quelconque adapté comme conteneur
    {int toto=0;
     Grandeur_scalaire_entier grand_courant(toto);
     TypeQuelconque typQ1(GENERIQUE_UNE_GRANDEUR_GLOBALE,NU_DDL,grand_courant);
 //    // 2) on adapte le stockage
 //    stock_compteur.ChangeGrandeur(typQ1);
 //    val_stock_compteur = (Grandeur_scalaire_entier*) stock_compteur.Grandeur_pointee();
 //    stock_icharge.ChangeGrandeur(typQ1);
 //    val_stock_icharge = (Grandeur_scalaire_entier*) stock_icharge.Grandeur_pointee();
     // 3) on ajoute en global les grandeurs consultables
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,COMPTEUR_ITERATION_ALGO_GLOBAL);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ALGO_GLOBAL_ACTUEL);
    };
    // --- maintenant les énergies
    {double titi=0.; // les énergies sont des doubles
     Grandeur_scalaire_double grand_courant(titi);
     TypeQuelconque typQ1(GENERIQUE_UNE_GRANDEUR_GLOBALE,NU_DDL,grand_courant);
     // et on ajoute
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_CINETIQUE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_EXTERNE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_INTERNE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,PUISSANCE_ACCELERATION);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,PUISSANCE_INTERNE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,PUISSANCE_EXTERNE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,PUISSANCE_BILAN);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_ELASTIQUE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_PLASTIQUE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_VISQUEUSE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_BILAN);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,QUANTITE_MOUVEMENT);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_PENALISATION);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_FROT_ELAST);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_FROT_PLAST);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_FROT_VISQ);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_VISCO_NUMERIQUE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_BULK_VISCOSITY);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_HOURGLASS_);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,ENERGIE_STABILISATION_MEMB_BIEL);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,VOLUME_TOTAL_MATIERE);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXPUISSEXT);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXPUISSINT);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXREACTION);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXRESIDUGLOBAL);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXdeltaX);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXvarDeltaX);
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,MAXvarDdl);
//     //------- debug
//     {void* pt_void =  ParaGlob::param->Mise_a_jour_grandeur_consultable(ENERGIE_PENALISATION);
//       TypeQuelconque* pt_quelc = (TypeQuelconque*) pt_void;
//       Grandeur_scalaire_double& gr
//                      = *((Grandeur_scalaire_double*) pt_quelc->Grandeur_pointee()); // pour simplifier
//     cout << "\n debug*** Algori::Algori "
//          << "\n pt_void= " << pt_void << " gr.ContDouble()= " << gr.ContDouble()
//          << " (*(gr.ConteneurDouble())) = " << (*(gr.ConteneurDouble()))
//          << flush;
//     }
//     // --- fin debug

    switch (ParaGlob::Dimension())
       { case 3: ParaGlob::param->Ajout_grandeur_consultable(&typQ1,VOL_TOTAL2D_AVEC_PLAN_XY);
         case 2: ParaGlob::param->Ajout_grandeur_consultable(&typQ1,VOL_TOTAL2D_AVEC_PLAN_XZ);
         case 1: ParaGlob::param->Ajout_grandeur_consultable(&typQ1,VOL_TOTAL2D_AVEC_PLAN_YZ);
         default: break;
       };
     // au début on met une norme très grande, car a priori on n'a absolument pas convergé !
     Grandeur_scalaire_double& gr
                    = *((Grandeur_scalaire_double*) typQ1.Grandeur_pointee()); // pour simplifier
     *(gr.ConteneurDouble()) = ConstMath::grand;
     ParaGlob::param->Ajout_grandeur_consultable(&typQ1,NORME_CONVERGENCE);
     
    };
  };
        
// constructeur de copie
Algori::Algori (const Algori& algo) :
  typeCalcul(algo.typeCalcul)
  ,avec_typeDeCalcul(algo.avec_typeDeCalcul)
  ,soustypeDeCalcul(algo.soustypeDeCalcul)
  ,avec_soustypeDeCalcul(algo.avec_soustypeDeCalcul)
  ,paraTypeCalcul(algo.paraTypeCalcul),mode_debug(algo.mode_debug)
  ,permet_affichage(algo.permet_affichage)
  ,pa(algo.pa),temps_derniere_sauvegarde(algo.temps_derniere_sauvegarde)
  ,list_incre_temps_sauvegarder(algo.list_incre_temps_sauvegarder)
  ,list_incre_temps_calculer(algo.list_incre_temps_calculer)
  ,tempsDerniereSortieFilCalcul(algo.tempsDerniereSortieFilCalcul)
  ,noms_fichier(algo.noms_fichier)
  ,typeFlotExterne(algo.typeFlotExterne)
  ,entreePrinc(algo.entreePrinc),debut_increment(algo.debut_increment)
  ,visualise(),visualise_maple(),visualise_geomview(),visualise_Gid(),visualise_Gmsh()
  ,toto_masse(NULL),icharge(algo.icharge),casAssemb_principal(algo.casAssemb_principal)
  ,phase_de_convergence(algo.phase_de_convergence)                     //   variables de travail
  ,nombre_de_bonnes_convergences(algo.nombre_de_bonnes_convergences)   //        "
  ,nombre_de_mauvaises_convergences(algo.nombre_de_mauvaises_convergences) //   "
  ,a_converge(algo.a_converge),a_converge_iterMoins1(algo.a_converge_iterMoins1)                                         //        "
  ,max_var_residu(algo.max_var_residu)                                 // variables de travail
  ,nb_cycle_test_max_var_residu(algo.nb_cycle_test_max_var_residu)     // variables de travail 
  ,v_int(algo.v_int)                                                   //      "             
  ,tab_para_stockage(algo.tab_para_stockage) //,pt_icharge(algo.pt_icharge)
  ,E_cin_0(algo.E_cin_0),E_cin_t(algo.E_cin_t),E_cin_tdt(algo.E_cin_tdt)
  ,E_int_t(algo.E_int_t),E_int_tdt(algo.E_int_tdt) // différentes énergies et bilan
  ,E_ext_t(algo.E_ext_t),E_ext_tdt(algo.E_ext_tdt),bilan_E(algo.bilan_E)     //               "
  ,q_mov_t(algo.q_mov_t),q_mov_tdt(algo.q_mov_tdt)
  ,P_acc_tdt(algo.P_acc_tdt),P_int_tdt(algo.P_int_tdt),P_ext_tdt(algo.P_ext_tdt),bilan_P_tdt(algo.bilan_P_tdt)
  ,P_acc_t(algo.P_acc_t),P_int_t(algo.P_int_t),P_ext_t(algo.P_ext_t),bilan_P_t(algo.bilan_P_t)
  ,E_visco_numerique_t(algo.E_visco_numerique_t),E_visco_numerique_tdt(algo.E_visco_numerique_tdt)
  ,E_bulk(algo.E_bulk),P_bulk(algo.P_bulk)
  ,F_int_t(algo.F_int_t),F_ext_t(algo.F_ext_t)          // forces généralisées int et ext à t et t+dt
  ,F_int_tdt(algo.F_int_tdt),F_ext_tdt(algo.F_ext_tdt)  //                  "
  ,F_totale_tdt(algo.F_totale_tdt),residu_final(algo.residu_final)
  ,X_t(algo.X_t),X_tdt(algo.X_tdt),delta_X(algo.delta_X),var_delta_X(algo.var_delta_X)
  ,delta_t_precedent_a_convergence(algo.delta_t_precedent_a_convergence)
  ,vitesse_t(algo.vitesse_t),vitesse_tdt(algo.vitesse_tdt)
  ,acceleration_t(algo.acceleration_t),acceleration_tdt(algo.acceleration_tdt)
  ,cas_ass_eps(NULL),Ass_eps(NULL),cas_ass_eps11(NULL),Ass_eps11(NULL)
  ,cas_ass_sig(NULL),Ass_sig(NULL),cas_ass_sig11(NULL),Ass_sig11(NULL)
  ,mat_glob_sigeps(NULL),vglob_sigeps(NULL),vglobal_sigeps(NULL)               
  ,nbddl_sigeps11(0),nbddl_sigeps(0)              
  ,cas_ass_err(NULL),Ass_err(NULL),nbddl_err(0)             
  ,avec_remont_sigma(false),avec_remont_eps(false),avec_erreur(false)
  ,energTotal(algo.energTotal),energHourglass(algo.energHourglass)
  ,energFrottement(algo.energFrottement),energPenalisation(algo.energPenalisation)
  ,energStabiliMembBiel(algo.energStabiliMembBiel)
  ,deja_lue_entete_parametre(algo.deja_lue_entete_parametre)
  //------- amortissement visqueux critique               
  ,amortissement_cinetique(algo.amortissement_cinetique)
  ,taille_moyenne_glissante(algo.taille_moyenne_glissante)
  ,moyenne_glissante(algo.moyenne_glissante), moy_gliss_t(algo.moy_gliss_t)
  ,sauve_E_cin(algo.sauve_E_cin),indice_EC1(algo.indice_EC1),nb_Ecin_stocker(algo.nb_Ecin_stocker)
  ,test_fraction_energie(algo.test_fraction_energie),X_EcinMax(algo.X_EcinMax)
  ,max_nb_decroit_pourRelaxDyn(algo.max_nb_decroit_pourRelaxDyn)
  ,inter_nb_entre_relax(1),fct_nD_inter_nb_entre_relax(algo.fct_nD_inter_nb_entre_relax)
  ,nom_fct_nD_inter_nb_entre_relax(algo.nom_fct_nD_inter_nb_entre_relax)
  ,nb_depuis_derniere_relax(0)
  ,compteur_decroit_pourRelaxDyn(algo.compteur_decroit_pourRelaxDyn)
  ,coef_arret_pourRelaxDyn(algo.coef_arret_pourRelaxDyn),pic_E_cint_t(algo.pic_E_cint_t)
  ,dernier_pic(algo.dernier_pic)
  ,max_pic_E_cin(algo.max_pic_E_cin),coef_redemarrage_pourRelaxDyn(algo.coef_redemarrage_pourRelaxDyn)
  ,max_deltaX_pourRelaxDyn(algo.max_deltaX_pourRelaxDyn),nb_max_dX_OK_pourRelaxDyn(algo.nb_max_dX_OK_pourRelaxDyn)
  ,nb_dX_OK_pourRelaxDyn(0),nb_deb_testfin_pourRelaxDyn(algo.nb_deb_testfin_pourRelaxDyn)
  ,nb_deb_test_amort_cinetique(algo.nb_deb_test_amort_cinetique)
  ,compteur_test_pourRelaxDyn(algo.compteur_test_pourRelaxDyn),compteur_pic_energie(algo.compteur_pic_energie)
  // ...... the same pour l'amortissement local
  ,amortissement_cinetique_au_noeud(algo.amortissement_cinetique_au_noeud)
  ,nb_deb_test_amort_cinetique_noe(algo.nb_deb_test_amort_cinetique_noe)
  ,E_cinNoe_tdt(algo.E_cinNoe_tdt),E_cinNoe_t(algo.E_cinNoe_t)
  ,max_nb_decroit_pourRelaxDyn_noe(algo.max_nb_decroit_pourRelaxDyn_noe)
  ,compteur_decroit_pourRelaxDyn_noe(algo.compteur_decroit_pourRelaxDyn_noe)
  ,coef_arret_pourRelaxDyn_noe(algo.coef_arret_pourRelaxDyn_noe),dernier_pic_noe(algo.dernier_pic_noe)
  ,pic_E_cint_t_noe(algo.pic_E_cint_t_noe),max_pic_E_cin_noe(algo.max_pic_E_cin_noe)
  ,coef_redemarrage_pourRelaxDyn_noe(algo.coef_redemarrage_pourRelaxDyn_noe)
  ,compteur_pic_energie_noe(algo.compteur_pic_energie_noe)
  //------- amortissement visqueux critique
  ,opt_cal_C_critique(algo.opt_cal_C_critique),f_(algo.f_)
  ,ampli_visco(algo.ampli_visco),C_inter(NULL)
  ,type_calcul_visqu_critique(algo.type_calcul_visqu_critique)
  ,F_int_tdt_prec(algo.F_int_tdt_prec),F_ext_tdt_prec(algo.F_int_tdt_prec)
  ,listeVarGlob(algo.listeVarGlob),listeVecGlob(algo.listeVecGlob)

  //----- fin des amortissements ---  
  ,arret_a_equilibre_statique(algo.arret_a_equilibre_statique),vglob_stat(algo.vglob_stat)
  ,volume_total_matiere(algo.volume_total_matiere),vol_total2D_avec_plan_ref(algo.vol_total2D_avec_plan_ref) 
  // -- controle précision convergence
  ,nom_modulation_precision(algo.nom_modulation_precision),modulation_precision(NULL),fct_norme(NULL)
  //---- grandeurs de travail privées ------
  ,vec_trav()
  //------- temps cpu -----
  ,tempsInitialisation(algo.tempsInitialisation),tempsMiseAjourAlgo(algo.tempsMiseAjourAlgo)
  ,tempsCalEquilibre(algo.tempsCalEquilibre),tempsRaidSmEner(algo.tempsRaidSmEner)
  ,tempsSecondMembreEnerg(algo.tempsSecondMembreEnerg),tempsResolSystemLineaire(algo.tempsResolSystemLineaire)
  ,tempsSauvegarde(algo.tempsSauvegarde),tempsSortieFilCalcul(algo.tempsSortieFilCalcul)
  ,tempsRaidSmEnerContact(algo.tempsRaidSmEnerContact)
  ,tempsSecondMembreEnergContact(algo.tempsSecondMembreEnergContact)
  ,temps_CL(algo.temps_CL),temps_CLL(algo.temps_CLL)
  ,temps_lois_comportement(algo.temps_lois_comportement),temps_metrique_K_SM(algo.temps_metrique_K_SM)
  ,temps_chargement(algo.temps_chargement),temps_rech_contact(algo.temps_rech_contact)
  ,tempsCalMasse(algo.tempsCalMasse),tempsCalViscoNum(algo.tempsCalViscoNum)
  ,ptalgocombi(NULL)
#ifndef UTILISATION_MPI
  ,lesTempsCpu(18)
#else
  ,lesTempsCpu(30),distribution_CPU_algo(algo.distribution_CPU_algo) // distribution de charge sur les CPU
  ,temps_transfert_court_algo(algo.temps_transfert_court_algo),temps_transfert_long_algo(algo.temps_transfert_long_algo)
  ,temps_attente_algo(algo.temps_attente_algo)
  ,temps_transfert_court_matSm(algo.temps_transfert_court_matSm),temps_transfert_long_matSm(algo.temps_transfert_long_matSm)
  ,temps_attente_matSm(algo.temps_attente_matSm)
  ,temps_transfert_court_charge(algo.temps_transfert_court_charge)
  ,temps_transfert_long_charge(algo.temps_transfert_long_charge)
  ,temps_attente_charge(algo.temps_attente_charge)
  ,temps_transfert_court_contact(algo.temps_transfert_court_contact)
  ,temps_transfert_long_contact(algo.temps_transfert_long_contact)
  ,temps_attente_contact(algo.temps_attente_contact)
  ,inter_transfer(),inter_transfer2()
  ,buffer_ioBI_MPI(algo.buffer_ioBI_MPI),tailles_buffer_ioBI_MPI(algo.tailles_buffer_ioBI_MPI)
  ,posi_buffer(algo.posi_buffer)
#endif

  //---------- stockage pour la transmission des grandeurs consultables -----
//  ,stock_compteur(algo.stock_compteur),val_stock_compteur(NULL)
//  ,stock_icharge(algo.stock_icharge),val_stock_icharge(NULL)

  {  // on met à jour ParaGlob
     ParaGlob::param->ChangeParaAlgoControleActifs(&pa);
       
	    if (vglob_stat != NULL)
	      vglob_stat = new Vecteur(*vglob_stat);
//     // 2) on adapte le stockage
//     val_stock_compteur = (Grandeur_scalaire_entier*) stock_compteur.Grandeur_pointee();
//     val_stock_icharge = (Grandeur_scalaire_entier*) stock_icharge.Grandeur_pointee();

  };
          
// destructeur   
Algori::~Algori ()
  { if (toto_masse != NULL)
      delete toto_masse;
    if (cas_ass_eps != NULL) {delete cas_ass_eps;cas_ass_eps=NULL;} 
    if (Ass_eps != NULL) {delete Ass_eps;Ass_eps=NULL;} 
    if (cas_ass_eps11 != NULL) {delete cas_ass_eps11;cas_ass_eps11=NULL;} 
    if (Ass_eps11 != NULL) {delete Ass_eps11;Ass_eps11=NULL;} 
    if (cas_ass_sig != NULL) {delete cas_ass_sig;cas_ass_sig=NULL;} 
    if (Ass_sig != NULL) {delete Ass_sig;Ass_sig=NULL;} 
    if (cas_ass_sig11 != NULL) {delete cas_ass_sig11;cas_ass_sig11=NULL;} 
    if (Ass_sig11 != NULL) {delete Ass_sig11;Ass_sig11=NULL;} 
    if (cas_ass_err != NULL) {delete cas_ass_err;cas_ass_err=NULL;} 
    if (Ass_err != NULL) {delete Ass_err;Ass_err=NULL;} 
    if (mat_glob_sigeps != NULL) {delete mat_glob_sigeps;mat_glob_sigeps=NULL;} 
    if (vglob_sigeps != NULL) {delete vglob_sigeps;vglob_sigeps=NULL;} 
    if (vglobal_sigeps != NULL) {delete vglobal_sigeps;vglobal_sigeps=NULL;} 
	   if (vglob_stat != NULL) delete vglob_stat;
  };

//lecture des parametres de controle
// et préparation de la sortie des paramètres globaux
void Algori::Lecture(UtilLecture & entreePrinc,ParaGlob & paraglob,LesMaillages&  lesMail)
  { // lecture des paramètres
    pa.Lecture_paraAlgoControle(entreePrinc);
    // mise à jour des paramètres globaux
    paraglob.Change_Nb_diggit_double(1,pa.Nb_diggit_double_calcul());  
    paraglob.Change_Nb_diggit_double(2,pa.Nb_diggit_double_graphique());
    paraglob.Change_Nb_diggit_double(3,pa.Nb_diggit_double_ecran());
    if (paraglob.EtatDeLaLecturePointInfo() <= 0) // cas de la première lecture dans le .info
      Preparation_conteneurs_interne(lesMail);
//    // -- modif concernant le stockage global pour volume déplacé: cas d'élément 2D dans un calcul 3D
//    // on met à jour le tableau de volume déplacé en fonction du nombre de maillage effectif
//    int nbMailMax = lesMail.NbMaillage();
//    // normalement ne fait rien si la taille n'est pas changé
//    if (vol_total2D_avec_plan_ref.Taille() != nbMailMax)
//      {vol_total2D_avec_plan_ref.Change_taille(nbMailMax,Coordonnee(ParaGlob::Dimension()));
//       // on remet à jour le stockage des grandeurs globales
//       if ( pa.CalVolTotalEntreSurfaceEtPlansRef())
//           { listeVarGlob["vol_total2D_avec_plan_yz"]=&(vol_total2D_avec_plan_ref(1)(1));
//             listeVarGlob["vol_total2D_avec_plan_xz"]=&(vol_total2D_avec_plan_ref(1)(2));
//             listeVarGlob["vol_total2D_avec_plan_xy"]=&(vol_total2D_avec_plan_ref(1)(3));
//           };
//       // s'il y a plus de 1 maillage, on définit de nouvelles variables globales
//       if (pa.CalVolTotalEntreSurfaceEtPlansRef())
//         for (int nbMail =2; nbMail<= nbMailMax; nbMail++)
//           { // l'ordre des string est fait pour que le switch ensuite soit ok
//             string nom1("vol_total2D_avec_plan_xy_"+ChangeEntierSTring(nbMail));
//             listeVarGlob[nom1]=&(vol_total2D_avec_plan_ref(nbMail)(3));
//             TypeQuelconque_enum_etendu::Ajout_un_TypeQuelconque_enum_etendu
//                (VOL_ELEM_AVEC_PLAN_REF,nom1,SCALAIRE_DOUBLE);
//             string nom2("vol_total2D_avec_plan_yz_"+ChangeEntierSTring(nbMail));
//             listeVarGlob[nom2]=&(vol_total2D_avec_plan_ref(nbMail)(1));
//             TypeQuelconque_enum_etendu::Ajout_un_TypeQuelconque_enum_etendu
//                (VOL_ELEM_AVEC_PLAN_REF,nom2,SCALAIRE_DOUBLE);
//             string nom3("vol_total2D_avec_plan_xz_"+ChangeEntierSTring(nbMail));
//             listeVarGlob[nom3]=&(vol_total2D_avec_plan_ref(nbMail)(2));
//             TypeQuelconque_enum_etendu::Ajout_un_TypeQuelconque_enum_etendu
//                (VOL_ELEM_AVEC_PLAN_REF,nom3,SCALAIRE_DOUBLE);
//
//             // --- on ajoute les grandeurs au niveau de paraglob
//             // si elles n'existent pas
//             double titi; // une grandeur de travail
//             Grandeur_scalaire_double grand_courant(titi); // idem
//             TypeQuelconque typQ1(GENERIQUE_UNE_GRANDEUR_GLOBALE,NU_DDL,grand_courant); // idem
//             switch (ParaGlob::Dimension())
//                { case 3: if (ParaGlob::param->GrandeurGlobal(nom1) == NULL)
//                             ParaGlob::param->Ajout_grandeur_consultable(&typQ1,nom1);
//                  case 2: if (ParaGlob::param->GrandeurGlobal(nom2) == NULL)
//                             ParaGlob::param->Ajout_grandeur_consultable(&typQ1,nom2);
//                  case 1: if (ParaGlob::param->GrandeurGlobal(nom3) == NULL)
//                             ParaGlob::param->Ajout_grandeur_consultable(&typQ1,nom1);
//                  default: break;
//                };
//           };
//      };
  };
// préparation des conteneurs interne, utilisés après la lecture
// où dans le cas où il n'y a pas de lecture directe
// (pour un sous algo via un algo combiné par exemple)
void Algori::Preparation_conteneurs_interne(LesMaillages&  lesMail)
 { // on teste la présence d'un élément , si oui on considère que tous les autres sont présents
   map < string, const double * , std::less <string> >::iterator imap;
   imap = listeVarGlob.find("energie_cinetique");
   if (imap == listeVarGlob.end()) // la map ne contient pas "energie_cinetique"
     {
        // remplissage de la liste des variables scalaires globales
        // tout d'abord pour les énergies globales
        listeVarGlob["energie_cinetique"]=&E_cin_tdt; listeVarGlob["energie_interne"]=&E_int_tdt;
        listeVarGlob["energie_externe"]=&E_ext_tdt; listeVarGlob["energie_bilan"]=&bilan_E;
        // la quantité de mouvement totale
        listeVarGlob["quantite_mouvement"]=&q_mov_tdt;
        // puis pour les puissances
        listeVarGlob["puissance_acceleration"]=&P_acc_tdt; listeVarGlob["puissance_interne"]=&P_int_tdt;
        listeVarGlob["puissance_externe"]=&P_ext_tdt; listeVarGlob["puissance_bilan"]=&bilan_P_tdt;
 //       // puis les énergies globales pour chaque maillage
 //       listeVarGlob["ener_cin_par_maillage"]=&E_cin_t; listeVarGlob["energie_interne"]=&E_int_t;
 //       listeVarGlob["energie_externe"]=&E_ext_t; listeVarGlob["energie_bilan"]=&bilan_E;
 //       // la quantité de mouvement totale
 //       listeVarGlob["quantite_mouvement"]=&q_mov_t; listeVarGlob["energie_bilan"]=&bilan_E;
 //       // puis pour les puissances
 //       listeVarGlob["puissance_acceleration"]=&P_acc_t; listeVarGlob["puissance_interne"]=&P_int_t;
 //       listeVarGlob["puissance_externe"]=&P_ext_t; listeVarGlob["puissance_bilan"]=&bilan_P_t;
        // pour les énergies internes, elles sont systématiquement calculées
        listeVarGlob["energie_elastique"]=&energTotal.EnergieElastique();
        listeVarGlob["energie_plastique"]=&energTotal.DissipationPlastique();
        listeVarGlob["energie_visqueuse"]=&energTotal.DissipationVisqueuse();
        listeVarGlob["energie_hourglass"]=&energHourglass;
        listeVarGlob["energie_penalisation"]=&energPenalisation;
        listeVarGlob["energie_stabilisation_membrane_biel"]=&energStabiliMembBiel;
        listeVarGlob["energie_frot_elast"]=&energFrottement.EnergieElastique();
        listeVarGlob["energie_frot_plast"]=&energFrottement.DissipationPlastique();
        listeVarGlob["energie_frot_visq"]=&energFrottement.DissipationVisqueuse();
        // cas d'une viscosité numérique on l'a sort systématiquement même si elle est nulle
        // ce qui permet une vérif éventuelle
        listeVarGlob["energie_visco_numerique"]=&E_visco_numerique_tdt;
        // idem pour le bulk viscosity
        listeVarGlob["energie_bulk_viscosity"]=&E_bulk;listeVarGlob["puissance_bulk_viscosity"]=&P_bulk;
        // volume totale de la matière
        listeVarGlob["volume_total_matiere"]=&volume_total_matiere;
        // -- volume déplacé: cas d'élément 2D dans un calcul 3D
        // vol_total2D_avec_plan_ref; // volume total entre la surface et : yz, xz,xy
        // dans une première étape on considère uniquement un seul maillage, ce sera changé éventuellement
        // aprés
        if ( pa.CalVolTotalEntreSurfaceEtPlansRef())
         { listeVarGlob["vol_total2D_avec_plan_yz"]=&(vol_total2D_avec_plan_ref(1)(1));
           listeVarGlob["vol_total2D_avec_plan_xz"]=&(vol_total2D_avec_plan_ref(1)(2));
           listeVarGlob["vol_total2D_avec_plan_xy"]=&(vol_total2D_avec_plan_ref(1)(3));
         };
        // -- cas des temps cpu
        
//   //------- temps cpu ----- on garde des temps spécifiques pour chaque algo
//   //   ainsi en sortie on pourra différencier les temps totaux et les temps partiels
//  Temps_CPU_HZpp tempsInitialisation;          // lesTempsCpu(1)
//  Temps_CPU_HZpp tempsMiseAjourAlgo;           // lesTempsCpu(2)
//  Temps_CPU_HZpp tempsCalEquilibre;            // lesTempsCpu(3)
//  Temps_CPU_HZpp tempsRaidSmEner;              // lesTempsCpu(4)
//  Temps_CPU_HZpp tempsSecondMembreEnerg;       // lesTempsCpu(5)
//  Temps_CPU_HZpp tempsResolSystemLineaire;     // lesTempsCpu(6)
//  Temps_CPU_HZpp tempsSauvegarde;              // lesTempsCpu(7)
//  Temps_CPU_HZpp tempsSortieFilCalcul;         // lesTempsCpu(8)
//  Temps_CPU_HZpp tempsRaidSmEnerContact;       // lesTempsCpu(9)
//  Temps_CPU_HZpp tempsSecondMembreEnergContact; // lesTempsCpu(10)
//  Temps_CPU_HZpp temps_CL;                     // lesTempsCpu(11)
//  Temps_CPU_HZpp temps_CLL;                    // lesTempsCpu(12)
//  Temps_CPU_HZpp temps_lois_comportement;      // lesTempsCpu(13)
//  Temps_CPU_HZpp temps_metrique_K_SM;          // lesTempsCpu(14)
//  Temps_CPU_HZpp temps_chargement;             // lesTempsCpu(15)
//  Temps_CPU_HZpp temps_rech_contact;           // lesTempsCpu(16)
//  Temps_CPU_HZpp tempsCalMasse;                // lesTempsCpu(17)
//  Temps_CPU_HZpp tempsCalViscoNum;             // lesTempsCpu(18)
//#ifdef UTILISATION_MPI
//  // cas d'un calcul parallèle: // passage des infos entre process
//  Temps_CPU_HZpp  temps_transfert_court_algo ;    // lesTempsCpu(19)
//  Temps_CPU_HZpp  temps_transfert_long_algo ;     // lesTempsCpu(20)
//  Temps_CPU_HZpp  temps_attente_algo ;            // lesTempsCpu(21)
//  Temps_CPU_HZpp  temps_transfert_court_matSm ;   // lesTempsCpu(22)
//  Temps_CPU_HZpp  temps_transfert_long_matSm ;    // lesTempsCpu(23)
//  Temps_CPU_HZpp  temps_attente_matSm ;           // lesTempsCpu(24)
//  Temps_CPU_HZpp  temps_transfert_court_charge ;  // lesTempsCpu(25)
//  Temps_CPU_HZpp  temps_transfert_long_charge ;   // lesTempsCpu(26)
//  Temps_CPU_HZpp  temps_attente_charge ;          // lesTempsCpu(27)
//  Temps_CPU_HZpp  temps_transfert_court_contact ; // lesTempsCpu(28)
//  Temps_CPU_HZpp  temps_transfert_long_contact ;  // lesTempsCpu(29)
//  Temps_CPU_HZpp  temps_attente_contact ;         // lesTempsCpu(30)

 //    Tableau <Coordonnee3> lesTempsCpu; // un tableau intermédiaire qui récupère et globalise les temps pour les sorties
 //                                        // via listeVarGlob, mais c'est les variables Temps_CPU_HZpp qui stockent
 //                                        // réellement les temps
        listeVarGlob["tpsU_InitAlgo"]=&lesTempsCpu(1)(1);
 //       listeVarGlob["tps_Syst_InitAlgo"]=&lesTempsCpu(1)(2);
 //       listeVarGlob["tps_Reel_InitAlgo"]=&lesTempsCpu(1)(3);
        listeVarGlob["tpsU_MiseAJourAlgo"]=&lesTempsCpu(2)(1);
 //       listeVarGlob["tps_Syst_MiseAJourAlgo"]=&lesTempsCpu(2)(2);
 //       listeVarGlob["tps_Reel_MiseAJourAlgo"]=&lesTempsCpu(2)(3);
        listeVarGlob["tpsU_CalEquilibre"]=&lesTempsCpu(3)(1);
 //       listeVarGlob["tps_Syst_CalEquilibre"]=&lesTempsCpu(3)(2);
 //       listeVarGlob["tps_Reel_CalEquilibre"]=&lesTempsCpu(3)(3);
        listeVarGlob["tpsU_MatSmLoc"]=&lesTempsCpu(4)(1);
 //       listeVarGlob["tps_Syst_MatSmLoc"]=&lesTempsCpu(4)(2);
 //       listeVarGlob["tps_Reel_MatSmLoc"]=&lesTempsCpu(4)(3);
        listeVarGlob["tpsU_SmLoc"]=&lesTempsCpu(5)(1);
 //       listeVarGlob["tps_Syst_SmLoc"]=&lesTempsCpu(5)(2);
 //       listeVarGlob["tps_Reel_SmLoc"]=&lesTempsCpu(5)(3);
        listeVarGlob["tpsU_ResSystLineaire"]=&lesTempsCpu(6)(1);
 //       listeVarGlob["tps_Syst_ResSystLineaire"]=&lesTempsCpu(6)(2);
 //       listeVarGlob["tps_Reel_ResSystLineaire"]=&lesTempsCpu(6)(3);
        listeVarGlob["tpsU_Sauvegarde"]=&lesTempsCpu(7)(1);
 //       listeVarGlob["tps_Syst_Sauvegarde"]=&lesTempsCpu(7)(2);
 //       listeVarGlob["tps_Reel_Sauvegarde"]=&lesTempsCpu(7)(3);
        listeVarGlob["tpsU_SortieFilCalcul"]=&lesTempsCpu(8)(1);
 //       listeVarGlob["tps_Syst_SortieFilCalcul"]=&lesTempsCpu(8)(2);
 //       listeVarGlob["tps_Reel_SortieFilCalcul"]=&lesTempsCpu(8)(3);
        listeVarGlob["tpsU_contactMatSmLoc"]=&lesTempsCpu(9)(1);
 //       listeVarGlob["tps_Syst_contactMatSmLoc"]=&lesTempsCpu(9)(2);
 //       listeVarGlob["tps_Reel_contactMatSmLoc"]=&lesTempsCpu(9)(3);
        listeVarGlob["tpsU_contactSmLoc"]=&lesTempsCpu(10)(1);
 //       listeVarGlob["tps_Syst_contactSmLoc"]=&lesTempsCpu(10)(2);
 //       listeVarGlob["tps_Reel_contactSmLoc"]=&lesTempsCpu(10)(3);
        // le tableau sert également aux temps pour la mise en place des CL et CLL
        // car la classe LesCondLim n'a pas vocation a gérer les I/O des tps de cpu
        // par contre les temps sont stocké dans LesCondLim, lesTempsCpu sert d'intermédiaire
        listeVarGlob["tpsU_CL"]=&lesTempsCpu(11)(1);
 //       listeVarGlob["tps_Syst_CL"]=&lesTempsCpu(11)(2);
 //       listeVarGlob["tps_Reel_CL"]=&lesTempsCpu(11)(3);
        listeVarGlob["tpsU_CLL"]=&lesTempsCpu(12)(1);
 //       listeVarGlob["tps_Syst_CLL"]=&lesTempsCpu(12)(2);
 //       listeVarGlob["tps_Reel_CLL"]=&lesTempsCpu(12)(3);
        listeVarGlob["tpsU_Lois_comp"]=&lesTempsCpu(13)(1);
        listeVarGlob["tpsU_Metriques"]=&lesTempsCpu(14)(1);
        listeVarGlob["tpsU_chargement"]=&lesTempsCpu(15)(1);
        listeVarGlob["tpsU_rech_contact"]=&lesTempsCpu(16)(1);
        listeVarGlob["tempsCalMasse"]=&lesTempsCpu(17)(1);
        listeVarGlob["tempsCalViscoNum"]=&lesTempsCpu(18)(1);
      #ifdef UTILISATION_MPI
        // cas d'un calcul parallèle
        listeVarGlob["tpsU_transfert_court_inter_cpu_Algo"]=&lesTempsCpu(19)(1);
        listeVarGlob["tpsU_attente_inter_cpu_Algo"]=&lesTempsCpu(20)(1);
        listeVarGlob["tpsU_transfert_long_inter_cpu_Algo"]=&lesTempsCpu(21)(1);
        listeVarGlob["tpsU_transfert_court_inter_cpu_Mat-Smloc"]=&lesTempsCpu(22)(1);
        listeVarGlob["tpsU_attente_inter_cpu_Mat-Smloc"]=&lesTempsCpu(23)(1);
        listeVarGlob["tpsU_transfert_long_inter_cpu-Mat-Smloc"]=&lesTempsCpu(24)(1);

        listeVarGlob["tpsU_transfert_court_inter_cpu_charge"]=&lesTempsCpu(25)(1);
        listeVarGlob["tpsU_transfert_long_inter_cpu_charge"]=&lesTempsCpu(26)(1);
        listeVarGlob["tpsU_attente_inter_cpu_charge"]=&lesTempsCpu(27)(1);
        listeVarGlob["tpsU_transfert_court_inter_cpu_contact"]=&lesTempsCpu(28)(1);
        listeVarGlob["tpsU_transfert_long_inter_cpu_contact"]=&lesTempsCpu(29)(1);
        listeVarGlob["tpsU_attente_inter_cpu_contact"]=&lesTempsCpu(30)(1);

      #endif

        // cas des vecteurs globaux
        int dim = ParaGlob::Dimension();
        if (ParaGlob::AxiSymetrie())
            dim--;
        // def des grandeurs courantes de type coordonnee
        Coordonnee coor(dim); // un type coordonnee typique
        // maintenant on définit une grandeur typique de type Grandeur_coordonnee
        Grandeur_coordonnee gt(coor);
        // def d'un type quelconque représentatif  pour un vecteur force à chaque noeud
        { // pour les forces externes
          TypeQuelconque typQ(FORCE_GENE_EXT,X1,gt);
          listeVecGlob.push_back(typQ);
        };
        { // pour les forces internes
          TypeQuelconque typQ(FORCE_GENE_INT,X1,gt);
          listeVecGlob.push_back(typQ);
        };
        { // pour les forces totales = internes + externes
          TypeQuelconque typQ(FORCE_GENE_TOT,X1,gt);
          listeVecGlob.push_back(typQ);
        };
        { // pour les forces résidu d'équilibre
          TypeQuelconque typQ(RESIDU_GLOBAL,X1,gt);
          listeVecGlob.push_back(typQ);
        };
        { // pour les variations de positions
          TypeQuelconque typQ(DELTA_XI,X1,gt);
          listeVecGlob.push_back(typQ);
        };
        listeVecGlob.sort(); // on ordonne la liste
        listeVecGlob.unique(); // suppression des doublons
       
       // -- modif concernant le stockage global pour volume déplacé: cas d'élément 2D dans un calcul 3D
       // on met à jour le tableau de volume déplacé en fonction du nombre de maillage effectif
       int nbMailMax = lesMail.NbMaillage();
       // normalement ne fait rien si la taille n'est pas changé
       if (vol_total2D_avec_plan_ref.Taille() != nbMailMax)
         {vol_total2D_avec_plan_ref.Change_taille(nbMailMax,Coordonnee(ParaGlob::Dimension()));
          // on remet à jour le stockage des grandeurs globales
          if ( pa.CalVolTotalEntreSurfaceEtPlansRef())
              { listeVarGlob["vol_total2D_avec_plan_yz"]=&(vol_total2D_avec_plan_ref(1)(1));
                listeVarGlob["vol_total2D_avec_plan_xz"]=&(vol_total2D_avec_plan_ref(1)(2));
                listeVarGlob["vol_total2D_avec_plan_xy"]=&(vol_total2D_avec_plan_ref(1)(3));
              };
          // s'il y a plus de 1 maillage, on définit de nouvelles variables globales
          if (pa.CalVolTotalEntreSurfaceEtPlansRef())
            for (int nbMail =2; nbMail<= nbMailMax; nbMail++)
              { // l'ordre des string est fait pour que le switch ensuite soit ok
                string nom1("vol_total2D_avec_plan_xy_"+ChangeEntierSTring(nbMail));
                listeVarGlob[nom1]=&(vol_total2D_avec_plan_ref(nbMail)(3));
                TypeQuelconque_enum_etendu::Ajout_un_TypeQuelconque_enum_etendu
                   (VOL_ELEM_AVEC_PLAN_REF,nom1,SCALAIRE_DOUBLE);
                string nom2("vol_total2D_avec_plan_yz_"+ChangeEntierSTring(nbMail));
                listeVarGlob[nom2]=&(vol_total2D_avec_plan_ref(nbMail)(1));
                TypeQuelconque_enum_etendu::Ajout_un_TypeQuelconque_enum_etendu
                   (VOL_ELEM_AVEC_PLAN_REF,nom2,SCALAIRE_DOUBLE);
                string nom3("vol_total2D_avec_plan_xz_"+ChangeEntierSTring(nbMail));
                listeVarGlob[nom3]=&(vol_total2D_avec_plan_ref(nbMail)(2));
                TypeQuelconque_enum_etendu::Ajout_un_TypeQuelconque_enum_etendu
                   (VOL_ELEM_AVEC_PLAN_REF,nom3,SCALAIRE_DOUBLE);
               
                // --- on ajoute les grandeurs au niveau de paraglob
                // si elles n'existent pas
                double titi=0.; // une grandeur de travail
                Grandeur_scalaire_double grand_courant(titi); // idem
                TypeQuelconque typQ1(GENERIQUE_UNE_GRANDEUR_GLOBALE,NU_DDL,grand_courant); // idem
                switch (ParaGlob::Dimension())
                   { case 3: if (ParaGlob::param->GrandeurGlobal(nom1) == NULL)
                                ParaGlob::param->Ajout_grandeur_consultable(&typQ1,nom1);
                     case 2: if (ParaGlob::param->GrandeurGlobal(nom2) == NULL)
                                ParaGlob::param->Ajout_grandeur_consultable(&typQ1,nom2);
                     case 1: if (ParaGlob::param->GrandeurGlobal(nom3) == NULL)
                                ParaGlob::param->Ajout_grandeur_consultable(&typQ1,nom1);
                     default: break;
                   };
              };
         };
    };
 };

// affichage des parametres
void Algori::Affiche() const// affichage de tous les parametres
  { Affiche1();
    Affiche2();
   };  
void Algori::Affiche1() const// affichage du type de calcul
  { cout << "--> type de calcul  = " << Nom_TypeCalcul(typeCalcul) << '\n';
    cout << paraTypeCalcul.Taille() << " parametre(s) = ";
    int nbcar = ParaGlob::NbdigdoEC();
    for (int i=1; i<= paraTypeCalcul.Taille(); i++)
      cout << setw(nbcar)  << setprecision(nbcar) << paraTypeCalcul(i) << " , ";
    cout << endl;  
  };     
void Algori::Affiche2() const// affichage des parametres de controle
  { pa.Affiche();
     };
      
     
// examine s'il y a convergence ou pas en fonction des parametres
// de controle et du resultat du systeme lineaire
// affiche: indique si l'on veut l'affichage ou non des infos
// itera : correspond au compteur d'itération dans un cas d'équilibre itératif (cal implicite par exemple)
// arrêt: indique qu'il vaut mieux arrêter le calcul
bool Algori::Convergence
   (bool affiche,double last_var_ddl_max,Vecteur& residu,double maxPuissExt,double maxPuissInt,double maxReaction
   ,int itera,bool& arret)
 { string type_norme = pa.Norme().nom1; // pour simplifier
   arret = false; // a priori on ne s'arrête pas
   double maxresidu = residu.Max_val_abs();
          // --- on passe les infos au niveau global
  {ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(MAXPUISSEXT,maxPuissExt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(MAXPUISSINT,maxPuissInt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(MAXREACTION,maxReaction);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(MAXRESIDUGLOBAL,maxresidu);
  };
 
////debug ------
//   int iddl=0;
//   maxresidu = residu.Max_val_abs(iddl);
//   cout << "\n iddl= "<<iddl << " residu_total(iddl)= "<<residu(iddl) ;
////fin debug ------
   
   double pa_precision = pa.Precision(); // initialisation de la précision demandé
   if (modulation_precision != NULL) // cas où en plus on fait une modulation de précision
     {pa_precision *= (modulation_precision->Valeur_pour_variables_globales())(1);}
  
	  double prec_res_mini = 4; // 1 mini
	  // dans le cas où on veut un test de convergence uniquement sur la partie statique  
	  if (   (arret_a_equilibre_statique ) // le dimensionnement se fera à l'affectation (=)
		 && (type_calcul_visqu_critique || pa.Amort_visco_artificielle()))
	   {
    #ifdef MISE_AU_POINT
     if (vglob_stat == NULL)
       { cout << "\n erreur le vecteur vglob_stat n'est pas attribue (=NULL)  "
		            << " on ne peut pas effectuer la convergence avec arret en equilibre statique uniquement "
				          << " revoir les parametres (ou se plaindre officiellement !! ) ";
         cout << "\n Algori::Convergence(..";
		       Sortie(1);
       }; 
    #endif
//	   if (vglob_stat->Taille() != 0) // si la taille est non nulle
//		   // cela veut dire que le vecteur a été affecté, donc on peut l'utiliser
		       maxresidu = vglob_stat->Max_val_abs();
		  };
   int nbcar = ParaGlob::NbdigdoEC();

   // 1) --- affichage des maxis de puissances virtuelles
   if (affiche && (ParaGlob::NiveauImpression() > 1))
    { cout << "\n   max puissance exterieure   =   " << std::fixed << std::setprecision(nbcar)<< maxPuissExt ;
      cout << "\n   max puissance interieurs   =   " << std::fixed << std::setprecision(nbcar)<< maxPuissInt ;
      cout << "\n   max des reactions          =   " << std::fixed << std::setprecision(nbcar)<< maxReaction ;
      cout << "\n   max du residu total        =   " << std::fixed << std::setprecision(nbcar)<< maxresidu ;
    }; 
   // 2) --- affichage de la norme
   double laNorme = 0.; // valeur par défaut

   if (strstr(type_norme.c_str(),"Residu/Reaction_et_Residu"))
      {cout << "\n *** norme exploratoire: ";
		     // on initialise le cas "Residu"
		     laNorme = maxresidu; //init 
		     double maxpuissance =MaX(maxReaction,MaX(maxPuissExt,maxPuissInt));
		     if (maxpuissance > prec_res_mini)
		      // cas de la norme Residu/Reaction
			     {if (maxReaction < 1.e-8)
				       laNorme = maxresidu;
			      else
				       laNorme = maxresidu/maxpuissance;
			      if (affiche &&(ParaGlob::NiveauImpression() > 1))
            cout << "\n**** Norme Residu/Reaction --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
        }
		     else // cas d'un résidu très faible
		      {pa_precision = prec_res_mini/50.;
		       if (affiche&&(ParaGlob::NiveauImpression() > 1))
             cout << "\n**** Norme Residu        -->  " << std::fixed << std::setprecision(nbcar)<< laNorme ;
         cout << "\n: pa_precision est mis au mini de la force admise: "
               << std::fixed << std::setprecision(nbcar)<< pa_precision << endl;
			     };
      } 
   else if (strstr(type_norme.c_str(),"Residu/Reaction"))
      {if (maxReaction < 1.e-8)
         laNorme = maxresidu;
       else
         laNorme = maxresidu/maxReaction;
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** Norme Residu/Reaction --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"Residu/PVExterne"))
     { if (maxPuissExt < 1.e-8)
         laNorme = maxresidu;
       else
         laNorme = maxresidu/maxPuissExt;
       if (affiche&&(ParaGlob::NiveauImpression() > 1))
         cout << "\n*** Norme Residu/PVExterne --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
     } 
   else if (strstr(type_norme.c_str(),"Residu/PVInterne"))
     { if (maxPuissInt < 1.e-8)
         laNorme = maxresidu;
       else
         laNorme = maxresidu/maxPuissInt;
       if (affiche&&(ParaGlob::NiveauImpression() > 1))
         cout << "\n*** Norme Residu/PVInterne --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
     } 
   else if (strstr(type_norme.c_str(),"Residu/maxPVI"))
     { // détermination du maximum des puissances
       double maxpuissance =MaX(maxReaction,MaX(maxPuissExt,maxPuissInt));
       if (maxpuissance < 1.e-8)
         laNorme = maxresidu;
       else
         laNorme = maxresidu/maxpuissance;
       if (affiche&&(ParaGlob::NiveauImpression() > 1))
         cout << "\n*** Norme Residu/maxPVI --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
     } 
   else if (strstr(type_norme.c_str(),"Residu"))
      {laNorme = maxresidu;
       if (affiche&&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** Norme Residu        -->  " << std::fixed << std::setprecision(nbcar)<< laNorme ;
      } 
   else if (strstr(type_norme.c_str(),"min(Res,Res/Reaction)"))
      {if (maxReaction < 1.e-8)
         laNorme = maxresidu;
       else
         laNorme = MiN(maxresidu,maxresidu/maxReaction);
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** Norme min(Res,Res/Reaction) --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"min(Res,Res/MaX(Reaction_et_PVExterne))"))
      {if (MaX(maxReaction,maxPuissExt) < 1.e-8)
         laNorme = maxresidu;
       else
         laNorme = MiN(maxresidu,maxresidu/MaX(maxReaction,maxPuissExt));
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** Norme min(Res,Res/MaX(Reaction_et_PVExterne)) --> "
              << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"E_cinetique/E_statique_ET_ResSurReact"))
      {if (affiche && (ParaGlob::NiveauImpression() > 1))
          { cout << "\n   E_cinetique   =   " << std::fixed << std::setprecision(nbcar)<< E_cin_tdt ;
            cout << "\n   E_ext         =   " << std::fixed << std::setprecision(nbcar)<< -E_ext_tdt ;
            cout << "\n   E_int         =   " << std::fixed << std::setprecision(nbcar)<< E_int_tdt ;
           };
       double max_abs_eners = MaX(Dabs(E_int_tdt),Dabs(E_ext_tdt));
		     if (MaX(max_abs_eners,maxReaction) < 1.e-8)
         laNorme = MaX(E_cin_tdt,maxresidu);
       else if  (maxReaction < 1.e-8)
         laNorme = MaX(E_cin_tdt/max_abs_eners,maxresidu);
       else 
         laNorme = MaX(E_cin_tdt/max_abs_eners,maxresidu/maxReaction);
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** Norme E_cinetique/E_statique_ET_ResSurReact --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"E_cinetique/E_statique_ET_Res/Reac_et_Fext"))
      {if (affiche && (ParaGlob::NiveauImpression() > 1))
          { cout << "\n   E_cinetique   =   " << std::fixed << std::setprecision(nbcar)<< E_cin_tdt ;
            cout << "\n   E_ext         =   " << std::fixed << std::setprecision(nbcar)<< -E_ext_tdt ;
            cout << "\n   E_int         =   " << std::fixed << std::setprecision(nbcar)<< E_int_tdt ;
           };
       double max_abs_eners = MaX(Dabs(E_int_tdt),Dabs(E_ext_tdt));
		     double maxpuiss_reac =MaX(maxReaction,MaX(maxPuissExt,maxPuissInt));
		     if (MaX(max_abs_eners,maxReaction) < 1.e-8)
         laNorme = MaX(E_cin_tdt,maxresidu);
       else if  (maxpuiss_reac < 1.e-8)
         laNorme = MaX(E_cin_tdt/max_abs_eners,maxresidu);
       else 
         laNorme = MaX(E_cin_tdt/max_abs_eners,maxresidu/maxpuiss_reac);
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** E_cinetique/E_statique_ET_Res/Reac_et_Fext --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"E_cin/E_stat_ET_min(Res,Res/Reac_et_Fext)"))
      {if (affiche && (ParaGlob::NiveauImpression() > 1))
          { cout << "\n   E_cinetique   =   " << std::fixed << std::setprecision(nbcar)<< E_cin_tdt ;
            cout << "\n   E_ext         =   " << std::fixed << std::setprecision(nbcar)<< -E_ext_tdt ;
            cout << "\n   E_int         =   " << std::fixed << std::setprecision(nbcar)<< E_int_tdt ;
           };
       double max_abs_eners = MaX(Dabs(E_int_tdt),Dabs(E_ext_tdt));
		     double maxpuiss_reac =MaX(maxReaction,MaX(maxPuissExt,maxPuissInt));
		     if (MaX(max_abs_eners,maxReaction) < 1.e-8)
         laNorme = MaX(E_cin_tdt,maxresidu);
       else if  (maxpuiss_reac < 1.e-8)
         laNorme = MaX(E_cin_tdt/max_abs_eners,maxresidu);
       else 
         laNorme = MaX(E_cin_tdt/max_abs_eners, MiN(maxresidu,maxresidu/maxpuiss_reac));
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** E_cin/E_stat_ET_min(Res,Res/Reac_et_Fext) --> "
              << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"E_cinetique"))
      {if (affiche && (ParaGlob::NiveauImpression() > 1))
          { cout << "\n   E_cinetique   =   " << std::fixed << std::setprecision(nbcar)<< E_cin_tdt ;
            cout << "\n   E_ext         =   " << std::fixed << std::setprecision(nbcar)<< -E_ext_tdt ;
            cout << "\n   E_int         =   " << std::fixed << std::setprecision(nbcar)<< E_int_tdt ;
           };
		     laNorme = E_cin_tdt;
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
          cout << "\n**** Norme E_cinetique --> " << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"Bilan_puissance/Max(puiss)"))
      {// si on est en statique, la puissance d'accélération = 0. normalement donc n'intervient pas
		     if (MaX(MaX(Dabs(P_acc_tdt),Dabs(P_int_tdt)),Dabs(P_ext_tdt)) < 1.e-8)
         laNorme = bilan_P_tdt;
       else
         laNorme = MaX(bilan_P_tdt,bilan_P_tdt/MaX(MaX(Dabs(P_acc_tdt),Dabs(P_int_tdt)),Dabs(P_ext_tdt)));
       if (affiche &&(ParaGlob::NiveauImpression() > 1))
         cout << "\n**** Norme Bilan_puissance/Max(puiss) --> "
              << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   else if (strstr(type_norme.c_str(),"fonction_nD:"))
      { // il s'agit ici de l'utilisation d'une fonction nD
        laNorme = (fct_norme->Valeur_pour_variables_globales())(1);
        if (affiche && (ParaGlob::NiveauImpression() > 1))
          { cout << "\n   E_cinetique   =   " << std::fixed << std::setprecision(nbcar)<< E_cin_tdt ;
            cout << "\n   E_ext         =   " << std::fixed << std::setprecision(nbcar)<< -E_ext_tdt ;
            cout << "\n   E_int         =   " << std::fixed << std::setprecision(nbcar)<< E_int_tdt ;
           };
        if (affiche &&(ParaGlob::NiveauImpression() > 1))
           cout << "\n**** Norme "<< std::fixed << std::setprecision(nbcar)<<pa.Norme().nom2<<" --> "
                << std::fixed << std::setprecision(nbcar)<< laNorme ;
      }
   if (affiche && (modulation_precision != NULL)) // cas où en plus on fait une modulation de précision
     cout << "\n (precision= "<< std::fixed << std::setprecision(nbcar)<<pa_precision<<" ) ";
  
   if (affiche&&(ParaGlob::NiveauImpression() > 1)) cout << endl; 
    
//---- debug
// if (E_cin_tdt == 0)
//   { cout << "\n énergie nulle !! : \n Algori::Convergence(... ";
//	  Sortie(1);};
//---- fin debug
	  

    // 3) --- premier test brut de la convergence      
    if ((laNorme <= pa_precision) && (itera != 0))
    	{a_converge = true;} // cas courant de convergence
    else if ((laNorme <= pa_precision) && (itera == 0))
       // cas particulier de la première itération
       	{ if (maxPuissInt == 0.)
            // si la puissance interne est strictement nulle, là il n'y a rien qui bouge, donc cela veut-dire
            // que tout est bloqué, et  on peut dire que cela a convergé car rien ne bougera après
            {a_converge = true;} 
          else if ((maxPuissInt <= ConstMath::trespetit) && (maxPuissExt > ConstMath::petit))
       	    // on a théoriquement convergé, mais on a une puissance externe non nulle et une puissance interne nulle
       	    // donc on est dans des petits nombres, il vaut mieux refaire un passage avec un calcul de delta ddl 
       		{a_converge = false;}
       	  else  {a_converge = true;};	
       	}						
    else
        {a_converge = false;};
// une partie  exploratoire
    int nbcoup = 1;
    if ((strstr(type_norme.c_str(),"Residu/Reaction_et_Residu")) && (!(laNorme <= pa_precision)) && (itera > nbcoup)
	   && (nombre_de_mauvaises_convergences > nbcoup))
     { cout << "\n norme transitoire = " << std::fixed << std::setprecision(nbcar)
             << pa_precision * PUISSN(10.,nombre_de_mauvaises_convergences-nbcoup )
             << endl;
     if (laNorme < pa_precision * PUISSN(10.,nombre_de_mauvaises_convergences-nbcoup ))
           {a_converge = true;};
     };				
// fin de la  partie  exploratoire
        
    // on met à jour le tableau de l'historique de l'évolution des résidus
	   // au minimum on sauvegarde le résidu actuel et le précédent
    int nb_cycle_controle = MaX(2,pa.NbCycleControleResidu());
    bool cycle_controle = false; // variable de pilotage: a priori le cycle est hors service
	   if (itera > 0)
     {if (itera == 1)
       // dimensionnement le premier incrément, puis mise en place pour les autres incréments
       { max_var_residu.Change_taille(nb_cycle_controle,ConstMath::tresgrand);};
      if (itera <= nb_cycle_controle)
         // cas ou on rempli le tableau max_var_residu sans s'en servir
        { max_var_residu(itera)=maxresidu;}
      else
        { // cas où on rempli les tableaux et on fait une rotation 
          cycle_controle = true; // le cycle est en service
          for (int i=1;i<nb_cycle_controle;i++)
              {max_var_residu(i)=max_var_residu(i+1);}
          max_var_residu(nb_cycle_controle)= maxresidu;
        };
     };	  	  

    // cas où l'on prend en compte l'évolution des var ddl et/ou des max résidu
    double var_mini_ddl = pa.VarMiniDdl();
//    double var_maxi_ddl = pa.VarMaxiDdl();
    int nbmax_cycle_test_var_residu= pa.PlageControleResidu();
    // traitement initial dans le cas iteratif
    if (itera > 0)
     { if ((pa.TypeDePilotage() == PILOT_GRADIENT)
          ||(strstr(type_norme.c_str(),"et_miniVarDdl")) || (strstr(type_norme.c_str(),"et_maxiVarDdl"))
          || (strstr(type_norme.c_str(),"et_VarRes")))
         { // initialisation éventuelle
           if (itera == 1) 
             {nb_cycle_test_max_var_residu=0;}
           if (itera > nb_cycle_controle)
             { //contrôle de la convergence que si ça n'a pas convergé
               if (!a_converge) 
                { // on signale qu'il vaut mieux arrêter car il n'y a pas de décroissance du résidu
                  if (((pa.TypeDePilotage() == PILOT_GRADIENT)||(strstr(type_norme.c_str(),"et_VarRes"))) &&
                      (max_var_residu(nb_cycle_controle) >= max_var_residu(1)))  
                        {if(nb_cycle_test_max_var_residu > nbmax_cycle_test_var_residu)
                               {arret = true; phase_de_convergence =-2;}
                         else  nb_cycle_test_max_var_residu++; }
                  // on signale qu'il vaut mieux arrêter car la variation de ddl est trop faible
                  if (((pa.TypeDePilotage() == PILOT_GRADIENT)||(strstr(type_norme.c_str(),"et_miniVarDdl"))) &&
                      (last_var_ddl_max <= var_mini_ddl))  {
                         arret = true; phase_de_convergence =-3;}
                  // on signale qu'il vaut mieux arrêter car la variation de ddl est trop grande
                  if (((pa.TypeDePilotage() == PILOT_GRADIENT)||(strstr(type_norme.c_str(),"et_maxiVarDdl"))) &&
                      (last_var_ddl_max <= var_mini_ddl))  {
                         arret = true; phase_de_convergence =-4;}
                };
             };
         };
     }; //- fin test sur itera
  
    // cas où on demande de vérifier deux fois le résidu
    if (strstr(type_norme.c_str(),"et_verification_ddouble"))
     {// tout d'abord on ne peut pas s'arrêter si le cycle de convergence
      if (itera > 0) // on n'intervient que pour les itérations > 0
       {if (itera == 1)  // on ne peut rien contrôler, donc par principe on repasse
         {a_converge_iterMoins1 = a_converge; // on sauvegarde pour l'itération qui suie
          a_converge = false;}
        else  // on est sup à 1
           // on applique l'algo
         {if (a_converge_iterMoins1 && a_converge)
            {} // c'est ok on ne fait rien
          else if (a_converge) // mais a_converge_iterMoins1 n'est pas vrai
            {a_converge_iterMoins1 = a_converge; // on le met à vrai
             a_converge = false; // on signal que la convergence n'est pas finie
            };
         };
       }
      else {a_converge = false;} ; // si itera = 0 on refait systématiquement une passe
     };
  
    // --- on passe les infos au niveau global
    Transfert_ParaGlob_NORME_CONVERGENCE(laNorme);
    
    //en dernier lieu on regarde si les puissances sont infinies ou nan
    // si oui on signale une divergence associée
    //  maxPuissExt,double maxPuissInt,double maxReaction
    if ((std::isinf(maxPuissExt))||( std::isnan(maxPuissExt)))
      {if ((ParaGlob::NiveauImpression() > 2) || (permet_affichage> 0))
         cout << "\n ***  maxPuissExt inifni ou nan : divergence imposee *** ";
       a_converge = false; phase_de_convergence=-11;arret = true;
      };
    if ((std::isinf(maxPuissInt))||( std::isnan(maxPuissInt)))
      {if ((ParaGlob::NiveauImpression() > 2) || (permet_affichage> 0))
         cout << "\n ***  maxPuissInt inifni ou nan : divergence imposee *** ";
       a_converge = false; phase_de_convergence=-11;arret = true;
      };
    if ((std::isinf(maxReaction))||( std::isnan(maxReaction)))
      {if ((ParaGlob::NiveauImpression() > 2) || (permet_affichage> 0))
         cout << "\n ***  maxReaction inifni ou nan : divergence imposee *** ";
       a_converge = false; phase_de_convergence=-11;arret = true;
      };

    // retour
    return a_converge;
   }; 
   
// sortie d'info sur un increment de calcul pour les ddl
// nb_casAssemb : donne le numéro du cas d'assemblage
void Algori::InfoIncrementDdl(LesMaillages * lesMail,
          int inSol,double maxDeltaDdl,const Nb_assemb& nb_casAssemb)
   {  // la variation maxi de ddl
      // recherche des Ddl concernes
      int nbNoeud1,nbMaillage1; double val0;
      Ddl ddl1 = lesMail->NoeudIndice(inSol,nbNoeud1,nbMaillage1,val0,nb_casAssemb.n);
      int nbcar = ParaGlob::NbdigdoEC();

      // ecriture des infos
      cout << "\n(assemblage " << nb_casAssemb.n << ") "
           << "delta maxi: " << std::fixed << std::setprecision(nbcar)<< maxDeltaDdl
           <<", ddl " << ddl1.Nom()
           <<", noeud/mail "<< nbNoeud1 << "/" << nbMaillage1 
           <<" , var total:"<< std::fixed << std::setprecision(nbcar)<< ddl1.Valeur() - val0
           << ", sur:" << std::fixed << std::setprecision(nbcar)<< ddl1.Valeur();
      // passage en global
      ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(MAXvarDdl,maxDeltaDdl);
      // on affiche les autres informations du même noeud si elles existent
      Noeud& noe=lesMail->Noeud_LesMaille(nbMaillage1,nbNoeud1); // récup du noeud
      Enum_ddl enumax=ddl1.Id_nom(); 
      int inci= enumax - PremierDdlFamille(enumax) ; // incice de composante du ddl -1
      bool absolue = true; // a priori on sort les infos en absolu
      List_io <Enum_ddl>  lienu = noe.Les_type_de_ddl(absolue); // récup de tous les types de ddl du noeud
      List_io <Enum_ddl>::iterator il,ilfin=lienu.end(); bool autres_valeurs = false;
      for (il=lienu.begin();il!=ilfin;il++) 
       { Enum_ddl enu = *il;
         if ( (enu!=ddl1.Id_nom()) && (noe.En_service(enu)) && ( (enu-PremierDdlFamille(enu))==inci))
          { if(!autres_valeurs) {cout << "\n autres ddl ";};
            cout << Nom_ddl(enu) << " " ;
            if (noe.Tdt())
              {cout << std::fixed << std::setprecision(nbcar)<< noe.Valeur_tdt(enu)
                    << " var " << std::fixed << std::setprecision(nbcar)<< noe.Valeur_tdt(enu)-noe.Valeur_t(enu);
              }
            else
              {cout << std::fixed << std::setprecision(nbcar)<< noe.Valeur_t(enu)
                    << " var " << std::fixed << std::setprecision(nbcar)<< noe.Valeur_t(enu)-noe.Valeur_0(enu);};
            cout << "; ";
            autres_valeurs=true;
           };
         };   
      cout << endl;     
   };
           
// sortie d'info sur un increment de calcul pour les réaction
// dans le cas d'un compteur, par exemple dans le cas implicite
// nb_casAssemb : donne le numéro du cas d'assemblage
void Algori::InfoIncrementReac(LesMaillages * lesMail,int compteur,
          int inreaction,double maxreaction,const Nb_assemb& nb_casAssemb)
   {       int nbNoeud2,nbMaillage2;
    #ifdef UTILISATION_MPI
      // cas d'un  calcul //, pour l'instant seule la (ou les) matrices du CPU 0 sont concernées
      if (ParaGlob::Monde()->rank() != 0)
        return ;
    #endif

      Ddl ddl2 = lesMail->NoeudIndice(inreaction,nbNoeud2,nbMaillage2,nb_casAssemb.n);
      int nbcar = ParaGlob::NbdigdoEC();

      // ecriture des infos
      cout << "\n ITERATION : " << compteur;
      if (ParaGlob::NiveauImpression() > 1)
       {cout << "\n(assemblage " << nb_casAssemb.n << ") "
           << " reaction maxi : " << std::fixed << std::setprecision(nbcar)<< maxreaction <<",  du ddl " << ddl2.Nom()
           <<", du noeud/maillage "<< nbNoeud2 << "/" << nbMaillage2 ;
       };
      cout << endl;     
   };        
      
// sortie d'info sur un increment de calcul pour les réaction
// dans le cas explicite par exemple (pas de compteur par rapport
// au cas précédent)
// nb_casAssemb : donne le numéro du cas d'assemblage
void Algori::InfoIncrementReac(LesMaillages * lesMail,
          int inreaction,double maxreaction,const Nb_assemb& nb_casAssemb)
   {       int nbNoeud2,nbMaillage2;
    #ifdef UTILISATION_MPI
      // cas d'un  calcul //, pour l'instant seule la (ou les) matrices du CPU 0 sont concernées
      if (ParaGlob::Monde()->rank() != 0)
        return ;
    #endif

      Ddl ddl2 = lesMail->NoeudIndice(inreaction,nbNoeud2,nbMaillage2,nb_casAssemb.n);
      int nbcar = ParaGlob::NbdigdoEC();

      // ecriture des infos
      if (ParaGlob::NiveauImpression() > 1)
       {cout << "\n(assemblage " << nb_casAssemb.n << ") "
             << " reaction maxi : " << std::fixed << std::setprecision(nbcar)<< maxreaction <<",  du ddl " << ddl2.Nom()
             <<", du noeud/maillage "<< nbNoeud2 << "/" << nbMaillage2 ;
        cout << endl; 
		  };    
   };        

    
// calcul une approximation des différences énergie en jeux
// affichage éventuelle des ces énergies et du bilan
// V         : ddl de vitesse
// coef_mass : en fait il faut utiliser 1/coef_mass * mat_mass pour la matrice masse, ceci due à la spécificité de l'algo
// les différentes énergie et efforts généralisées sont des variables internes à l'algo
// en sortie: énergie cinétique : E_cin_tdt,  énergie interne : E_int_tdt, énergie externe = E_ext_tdt, bilan : bilan_E
//            idem avec les puissances
//            énergie visqueuse due aux forces visqueuses numériques
// icharge : compteur d'increment
// brestart : booleen qui indique si l'on est en restart ou pas
void Algori::CalEnergieAffichage(const double& coef_mass,const  Vecteur & V,const Mat_abstraite& mat_mass 
                                ,const Vecteur & delta_X,int icharge,bool brestart,const Vecteur & gamma
                                ,const Vecteur & forces_vis_num)
 {
 #ifdef UTILISATION_MPI
  // cas d'un  calcul //, le CPU 0 calcule
  int proc_en_cours = ParaGlob::Monde()->rank();
  // on crée un vecteur intermédiaire pour le passage d'information
  Vecteur  v_val_inter(10);

  if (proc_en_cours == 0)
   {
 #endif
    // -- le calcul est conditionnelle
    if (icharge >= pa.NbIncrCalEnergie())
     {int nbcar = ParaGlob::NbdigdoEC();

       // --- cas des énergies ------
       // initialisation 
       v_int.Change_taille(V.Taille()); v_int.Zero(); 
       // calcul de l'énergie cinétique et de la quantité de mouvement 
       // v_int contient les quantités de mouvement locales à chaque ddl
       E_cin_tdt = (0.5/coef_mass) * (V * (mat_mass.Prod_vec_mat(V,v_int)));
//---- debug
// if (E_cin_tdt == 0)
//   { cout << "\n énergie nulle !! : \n Algori::CalEnergieAffichage(... ";
//	  Sortie(1);};
//---- fin debug
		 
       // calcul de la quantité de mouvement globale
       q_mov_tdt = v_int.Somme();
       // calcul de l'énergie interne dépensée
       v_int.Zero();  // init à 0
       if (!brestart)
        { // cas normal
          v_int += F_int_t; v_int += F_int_tdt; E_int_tdt = E_int_t - 0.5 * (delta_X * v_int); }
       else // cas d'un restart, F_int_t n'est pas défini, on utilise la formule des rectangles plutôt que des trapèzes
        { v_int += F_int_tdt; E_int_tdt = E_int_t - delta_X * v_int; } 
       // calcul de l'énergie externe dépensée
       v_int.Zero();  // init à 0
       if (!brestart)
        { // cas normal
          v_int += F_ext_t; v_int += F_ext_tdt; E_ext_tdt = E_ext_t + 0.5 * (delta_X * v_int);}
       else // cas d'un restart, F_int_t n'est pas défini, on utilise la formule des rectangles plutôt que des trapèzes
        { v_int += F_ext_tdt; E_ext_tdt = E_ext_t + delta_X * v_int;}

       // bilan des énergies (on retire l'énergie cinétique initiale si elle exite)
       bilan_E = E_cin_tdt - E_cin_0 + E_int_tdt - E_ext_tdt;
       // --- cas des puissances ------
       // calcul de la puissance d'accélération
       v_int.Zero(); 
       P_acc_tdt = (1./coef_mass) * (V * (mat_mass.Prod_vec_mat(gamma,v_int)));
       // calcul de la puissance interne
       P_int_tdt = V * F_int_tdt;
       // calcul de la puissance externe
       P_ext_tdt = V * F_ext_tdt;
       // bilan des puissances
       bilan_P_tdt = P_acc_tdt - P_int_tdt - P_ext_tdt;
       // cas des forces visqueuses numériques
       double puissance_visco_num = 0.;
       if (pa.Amort_visco_artificielle())
          {E_visco_numerique_tdt = E_visco_numerique_t + delta_X * forces_vis_num;
           puissance_visco_num = V * forces_vis_num;
           bilan_P_tdt += puissance_visco_num;
           bilan_E += E_visco_numerique_tdt;
          };
       // --- affichage ----
       if (( pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde))
           &&(pa.AfficheIncrEnergie()))
         {// cas des énergies
          cout << "\n energies: cin= "<< std::fixed << std::setprecision(nbcar)<< E_cin_tdt
               <<" int= " << std::fixed << std::setprecision(nbcar)<< E_int_tdt
               << " ext= "<< std::fixed << std::setprecision(nbcar)<< E_ext_tdt
               << " bilan= " << std::fixed << std::setprecision(nbcar)<< bilan_E
               << " qmv= " << std::fixed << std::setprecision(nbcar)<< q_mov_tdt;
               
          // cas des puissances
          cout << "\n puissances: acc= "<< std::fixed << std::setprecision(nbcar)<< P_acc_tdt
               <<" int= " << std::fixed << std::setprecision(nbcar)<< P_int_tdt
               << " ext= " << std::fixed << std::setprecision(nbcar)<< P_ext_tdt
               << " bilan= " << std::fixed << std::setprecision(nbcar)<< bilan_P_tdt;
               
          if (pa.Amort_visco_artificielle())
            cout << "\n energie d'amortissement numerique: "
                 << std::fixed << std::setprecision(nbcar)<< E_visco_numerique_tdt
                 << " puissance d'amortissement numerique: "
                 << std::fixed << std::setprecision(nbcar)<< puissance_visco_num;
         };
     #ifdef UTILISATION_MPI
        v_val_inter(1) = E_cin_tdt;
        v_val_inter(2) = E_ext_tdt;
        v_val_inter(3) = E_int_tdt;
        v_val_inter(4) = E_visco_numerique_t;
        v_val_inter(5) = bilan_E;
        v_val_inter(6) = q_mov_tdt;
        v_val_inter(7) = P_acc_tdt;
        v_val_inter(8) = P_int_tdt;
        v_val_inter(9) = P_ext_tdt;
        v_val_inter(10) = bilan_P_tdt;
     #endif
      };
 #ifdef UTILISATION_MPI
   };
  v_val_inter.Broadcast(0);
  if (proc_en_cours != 0)
   { // récup des infos
    E_cin_tdt = v_val_inter(1);
    E_ext_tdt = v_val_inter(2);
    E_int_tdt = v_val_inter(3);
    E_visco_numerique_t = v_val_inter(4);
    bilan_E = v_val_inter(5);
    q_mov_tdt = v_val_inter(6);
    P_acc_tdt = v_val_inter(7);
    P_int_tdt = v_val_inter(8);
    P_ext_tdt = v_val_inter(9);
    bilan_P_tdt = v_val_inter(10);
   };
 #endif
  // --- on passe les infos au niveau global
  {ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_CINETIQUE,E_cin_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_EXTERNE,E_ext_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_INTERNE,E_int_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_VISCO_NUMERIQUE,E_visco_numerique_t);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_BILAN,bilan_E);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(QUANTITE_MOUVEMENT,q_mov_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(PUISSANCE_ACCELERATION,P_acc_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(PUISSANCE_INTERNE,P_int_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(PUISSANCE_EXTERNE,P_ext_tdt);
   ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(PUISSANCE_BILAN,bilan_P_tdt);
  };

 };

// idem pour un calcul statique, c'est-à-dire sans énergie cinématique
// affichage éventuelle des ces énergies et du bilan
// icharge : compteur d'increment
// brestart : booleen qui indique si l'on est en restart ou pas
void Algori::CalEnergieAffichage(const Vecteur & delta_X,int icharge,bool brestart,const Vecteur & forces_vis_num)
  {
   #ifdef UTILISATION_MPI
    // cas d'un  calcul //, le CPU 0 calcule
    int proc_en_cours = ParaGlob::Monde()->rank();
    // on crée un vecteur intermédiaire pour le passage d'information
    Vecteur  v_val_inter(6);

    if (proc_en_cours == 0)
     {
   #endif

    // le calcul est conditionnelle
    if (icharge >= pa.NbIncrCalEnergie())
      {int nbcar = ParaGlob::NbdigdoEC();
       // initialisation
       v_int.Change_taille(delta_X.Taille()); v_int.Zero(); 
       // calcul de l'énergie interne dépensée       
       if (!brestart)
        { // cas normal
          v_int += F_int_t; v_int += F_int_tdt; E_int_tdt = E_int_t - 0.5 * (delta_X * v_int); }
       else // cas d'un restart, F_int_t n'est pas défini, on utilise la formule des rectangles plutôt que des trapèzes
        { v_int += F_int_tdt; E_int_tdt = E_int_t - delta_X * v_int; } 
       // calcul de l'énergie externe dépensée
       v_int.Zero();  // init à 0
       if (!brestart)
        { // cas normal
          v_int += F_ext_t; v_int += F_ext_tdt; E_ext_tdt = E_ext_t + 0.5 * (delta_X * v_int);}
       else // cas d'un restart, F_ext_t n'est pas défini, on utilise la formule des rectangles plutôt que des trapèzes
        { v_int += F_ext_tdt; E_ext_tdt = E_ext_t + delta_X * v_int;}
       // bilan des énergies (on retire l'énergie cinétique initiale si elle exite)
       bilan_E = E_int_tdt - E_ext_tdt;
       // --- cas des puissances ------
       // recup de l'incrément de temps
	      double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
       double unSurDeltat=0;    
	      if (Abs(deltat) >= ConstMath::trespetit)
         {unSurDeltat = 1./deltat;}
       else
         // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
         { // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
           if (unSurDeltat < 0)
              { cout << "\n le pas de temps est négatif !! "; };
           unSurDeltat = ConstMath::tresgrand;
         };
       // calcul de la puissance interne
       P_int_tdt = unSurDeltat * delta_X * F_int_tdt;
       // calcul de la puissance externe
       P_ext_tdt = unSurDeltat * delta_X * F_ext_tdt;
       // bilan des puissances
       bilan_P_tdt = P_int_tdt + P_ext_tdt;
       
       // cas des forces visqueuses numériques
       double puissance_visco_num = 0.;
       if (pa.Amort_visco_artificielle())
          {E_visco_numerique_tdt = E_visco_numerique_t + delta_X * forces_vis_num;
           puissance_visco_num = unSurDeltat * delta_X * forces_vis_num;
           bilan_P_tdt += puissance_visco_num;
           bilan_E += E_visco_numerique_tdt;
          };
          
       // --- affichage ----
       if (( pa.Vrai_commande_sortie(icharge,temps_derniere_sauvegarde))
           &&(pa.AfficheIncrEnergie()))
         {// cas des énergies
           cout << "\n energies: int= " << std::fixed << std::setprecision(nbcar)<< E_int_tdt
                << " ext= " << std::fixed << std::setprecision(nbcar)<< E_ext_tdt
                << " bilan= " << std::fixed << std::setprecision(nbcar)<< bilan_E;
          if (pa.Amort_visco_artificielle())
              cout << " energie numerique: "<< std::fixed << std::setprecision(nbcar)<< E_visco_numerique_tdt;
          if (pa.BulkViscosity())
             cout << " ener_bulk: " << std::fixed << std::setprecision(nbcar)<< E_bulk
                  << " P_bulk: "<< std::fixed << std::setprecision(nbcar)<< P_bulk;

          // cas des puissances
          cout << "\n puissances(V=delta X): int= " << std::fixed << std::setprecision(nbcar)<< P_int_tdt
               << " ext= "<< std::fixed << std::setprecision(nbcar)<< P_ext_tdt
               << " bilan= " << std::fixed << std::setprecision(nbcar)<< bilan_P_tdt;
         };

    #ifdef UTILISATION_MPI
        v_val_inter(1) = E_ext_tdt;
        v_val_inter(2) = E_int_tdt;
        v_val_inter(3) = bilan_E;
        v_val_inter(4) = P_int_tdt;
        v_val_inter(5) = P_ext_tdt;
        v_val_inter(6) = E_visco_numerique_t;
    #endif
      };
   #ifdef UTILISATION_MPI
     };
    v_val_inter.Broadcast(0);
    if (proc_en_cours != 0)
     { // récup des infos
      E_ext_tdt = v_val_inter(1);
      E_int_tdt = v_val_inter(2);
      bilan_E = v_val_inter(3);
      P_int_tdt = v_val_inter(4);
      P_ext_tdt = v_val_inter(5);
      E_visco_numerique_t = v_val_inter(6);
     };
   #endif
    // --- on passe les infos au niveau global
    {ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_EXTERNE,E_ext_tdt);
     ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_INTERNE,E_int_tdt);
     ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_BILAN,bilan_E);
     ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(PUISSANCE_INTERNE,P_int_tdt);
     ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(PUISSANCE_EXTERNE,P_ext_tdt);
     ParaGlob::param->Mise_a_jour_grandeur_consultable_Scalaire_double(ENERGIE_VISCO_NUMERIQUE,E_visco_numerique_t);
    };
   };

// dans le cas où un grand nombre d'information nécessaire pour l'algorithme
// est stocké de manière externe sur un fichier
// exemple :  calcul d'erreur a partir
// d'un champ  d'information stocké de manière externe sur un fichier
// l'acquisition du (ou des) lieu(s) (flot) où sont stocké ces infos est effectuée
// ce lieu fait parti des paramètres de l'algorithme
void Algori::DefFlotExterne(UtilLecture* entreePrinc)
 {if (ParaGlob::NiveauImpression() >= 4)
     cout << " debut de la lecture eventuelle de fichiers types flotExterne " << endl;
  // il y a lecture si le mot clé : flotExterne est rencontré
  if (strstr(entreePrinc->tablcar,"flotExterne")!=NULL)
   { // on passe en revue tous les cas possibles
     // la lecture s'effectue jusqu'à un nouveau mot clé
     MotCle motCle; // ref aux mots cle
     // 1.  preparation du flot
     entreePrinc->NouvelleDonnee();
     do
     { 
      // 2.  lecture du type de données à lire
      int type;
      *(entreePrinc->entree) >> type ; 

      // ------- cas d'un calcul d'erreur sans calcul d'équilibre -------
      //         fichier des contraintes
      if (type==1) 
        // on balaie la liste des sous_types
        { list <EnumSousTypeCalcul>::const_iterator ili;
          list <EnumSousTypeCalcul>::const_iterator ili_fin = soustypeDeCalcul->end();
          for (ili = soustypeDeCalcul->begin();ili!= ili_fin;ili++)
             if(!(Avec_in(*ili)) && Remonte_in(*ili))
                LectureUnTypeExterne(entreePrinc,type);
         };       
      // ------- cas d'un calcul d'erreur sans calcul d'équilibre -------
      //         fichier des déplacements
      if (type==2)
        // on balaie la liste des sous_types
        { list <EnumSousTypeCalcul>::const_iterator ili;
          list <EnumSousTypeCalcul>::const_iterator ili_fin = soustypeDeCalcul->end();
          for (ili = soustypeDeCalcul->begin();ili!= ili_fin;ili++)
             if(!(Avec_in(*ili)) && Remonte_in(*ili))
                LectureUnTypeExterne(entreePrinc,type);
          }      
       // ------ autre cas à rajouter par la suite -------
       
      // 3.  preparation du flot pour la sortie ou la prochaine lecture
      entreePrinc->NouvelleDonnee();
      }
     while ( !motCle.SimotCle(entreePrinc->tablcar)) ;
    }   
  if (ParaGlob::NiveauImpression() >= 4)
     cout << " fin de la lecture eventuelle de fichiers types flotExterne " << endl;
    };
 
// lecture d'un type externe
void Algori::LectureUnTypeExterne(UtilLecture* entreePrinc,const int type)
     { // on incrémente la taille des tableaux
       int taille = typeFlotExterne.Taille(); taille++;
       typeFlotExterne.Change_taille(taille);
       noms_fichier.Change_taille(taille);
       // sauvegarde du type
       typeFlotExterne(taille) = type;
       // on crée une liste intermédiaire
       ListDeuxString la_list; list<DeuxString > & lalist = la_list.infos;
       // deux strings
       DeuxString nom;
       // lecture tant qu'il n'y a pas d'erreur
       while (!entreePrinc->entree->rdstate()) 
        {  // lecture du nom du fichier
           *(entreePrinc->entree) >> nom.nomFichier;
           *(entreePrinc->entree) >> nom.nomMaillage;
           // stockage   
           lalist.push_back(nom);           
         }
       // stokage de la liste
       noms_fichier(taille) = la_list;  
      }; 
   
// ------------METHODES PROTEGEES ----------------:
   
// calcul de la raideur et du second membre pour tous les maillages
// ainsi que des énergies 
// si pb retour false: indique que le calcul c'est mal passé, il y a une erreur
// généré, les résultats: raideur et second membres calculées sont inexploitable
bool Algori::RaidSmEner(LesMaillages * lesMail,Assemblage& Ass,Vecteur & vglobin
                ,Mat_abstraite& matglob )
  {
    tempsRaidSmEner.Mise_en_route_du_comptage(); // comptage cpu
  #ifdef UTILISATION_MPI
    mpi::request reqs1;
    mpi::request reqs2;
    mpi::request reqs3;
    bool premier_passage = true;
    int proc_en_cours = ParaGlob::Monde()->rank();
    int gestion_exception = 0; // pour la gestion
  #endif
    bool retour = true; // init par défaut
    // on doit initialiser les temps globaux de loi de comp et de métrique car
    // on va ajouter ensuite des temps qui eux mêmes globalisent ce qui se passe à chaque pti
//    Temps_CPU_HZpp tps_zero; // un temps d'initialisation
//    temps_lois_comportement = tps_zero;
    Temps_CPU_HZpp temps_debut_lois_comportement; // temps initial
//    temps_metrique_K_SM = tps_zero;
    Temps_CPU_HZpp temps_debut_metrique_K_SM; // temps initial
  #ifdef UTILISATION_MPI
    if (proc_en_cours == 0)
  #endif
    if (permet_affichage >3) cout << "\n -- debut calcul second membre et raideur " << flush;

    energTotal.Inita(0.); // initialisation des énergies à 0
    energHourglass=0.;energStabiliMembBiel=0.;
    volume_total_matiere = 0.; // init
    E_bulk = 0.; // init
    P_bulk=0.; // init
	   if (pa.CalVolTotalEntreSurfaceEtPlansRef())
	     vol_total2D_avec_plan_ref.Inita(Coordonnee(ParaGlob::Dimension()));
      
    // on gère les exceptions éventuelles en mettant le bloc sous surveillance
    int nbMailMax = lesMail->NbMaillage();
    int dernier_nbMail=1; // pour la gestion d'erreur
    list <int >::const_iterator dernier_il; // pour la gestion d'erreur
    try
      {// boucle sur les elements
   #ifndef UTILISATION_MPI
       // cas d'un  calcul mono CPU
       for (int nbMail =1; nbMail<= nbMailMax; nbMail++)
        {int nemax = lesMail->Nombre_element(nbMail);
         for (int ne=1; ne<= nemax;ne++)
            {
   #else
       // cas d'un calcul multi-CPU
       // on va parcourir les éléments associés au cpu
       // on récupère la distribution d'éléments concernant le cpu en cours
       const Tableau < list <int > >&  tab_list_elem_cpu = distribution_CPU_algo.List_element_CPU_en_cours();
       // la taille est identique à nbMailMax, sauf si c'est le cpu 0, là le tableau est vide et il n'y
       // aura pas de boucle, (ou plutôt on sortira directement de la boucle
       // le cas CPU 0 est traité à part
       int nb_mail_distrib = tab_list_elem_cpu.Taille();
       for (int nbMail =1; nbMail<= nb_mail_distrib; nbMail++,dernier_nbMail++)
         {const list <int >& list_elem_cpu= tab_list_elem_cpu(nbMail); // pour simplifier
          // on balaie les éléments nécessaires
          list <int >::const_iterator il,ilfin=list_elem_cpu.end();
          for (il = dernier_il = list_elem_cpu.begin();il != ilfin;il++,dernier_il++)
            { int ne = (*il);
              // on récupère un signal du process 0
              temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
              if (premier_passage) {premier_passage = false;}
              else // on attend que les transferts soient finis
                {if (reqs1.active()) reqs1.wait();
                 if (reqs2.active()) reqs2.wait();
                 if (reqs3.active()) reqs3.wait();
                };
              temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu

   #endif
              //calcul de la raideur local et du residu local
              ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne)); // l'element
              Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
              temps_debut_lois_comportement=el.Temps_lois_comportement(); // init
              temps_debut_metrique_K_SM = el.Temps_metrique_K_SM(); // init
              Element::ResRaid resu = el.Calcul_implicit(pa);
              
              
////---- débug              
//{
//cout << "\n Algori::RaidSmEner debug : ";
//string entete = " affichage de la matrice de raideur avant assemblage";
//matglob.Affichage_ecran(entete); 
//entete = " affichage du second membre (puissance interne) avant assemblage ";
//vglobin.Affichage_ecran(entete); 
//}
//             cout << "\n affichage des matrices locales : element numéro =" << ne << endl;
//             (resu.res)->Affiche();
//             (resu.raid)->Affiche();
// on test la symétrie des matrices locales
//if (!(resu.raid->Symetrie()))
//  resu.raid->AfficheNonSymetries();
//-- fin débug             
            #ifndef UTILISATION_MPI
              // --- assemblage et calcul de grandeur globale, pour mono-CPU
              // assemblage
              Ass.AssemSM (vglobin,*(resu.res),el.TableauDdl(),taN); // du second membre
////---debug---
//double maxPuissInt = vglobin.Max_val_abs();
//if (Dabs(maxPuissInt) > ConstMath::petit)
//  {cout << "\n debug Algori::RaidSmEner maxPuissInt=  "<<maxPuissInt<<endl
//        << "\n residu interne= "<<*(resu.res) << endl ;
//  };
////---fin debug ---
              if (pa.Symetrie_matrice())
                 Ass.AssembMatSym (matglob,*(resu.raid),el.TableauDdl(),taN); // de la raideur
              else
                 Ass.AssembMatnonSym (matglob,*(resu.raid),el.TableauDdl(),taN); // de la raideur
            #else
              // cas d'un calcul parallèle, et CPU != 0
              if (proc_en_cours != 0)
                {tempsRaidSmEner.Arret_du_comptage();
                 temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
                 DeuxEntiers num_el_et_mail(el.Num_elt(),el.Num_maillage());
                 // on transmet les numéros d'élément et de maillage
                 reqs1 = ParaGlob::Monde()->isend(0, 24, num_el_et_mail);
                 // puis on transmets le vecteur résidu
                 reqs2 = resu.res->Ienvoi_MPI(0,25);
                 //puis la matrice
                 reqs3 = resu.raid->Ienvoi_MPI(0, 26);
                 temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
                 tempsRaidSmEner.Mise_en_route_du_comptage();
                };

            #endif
              // globalisation des énergies
              energTotal += el.EnergieTotaleElement();
              energHourglass += el.Energie_Hourglass();
              energStabiliMembBiel += el.Energie_stab_membBiel();
              E_bulk += el.Energie_Bulk();
              P_bulk += el.Puissance_Bulk();
              // globalisation des volumes
              volume_total_matiere += el.Volume();
              if (pa.CalVolTotalEntreSurfaceEtPlansRef())
               { const Coordonnee& volPlan = el.VolumePlan();
                 if (volPlan.Dimension())
                   vol_total2D_avec_plan_ref(nbMail) += volPlan;
              //--- debug
              //  cout << "\n calcul des volumes balayés";
              //  vol_total2D_avec_plan_ref.Affiche();
              //  cout << "\n " << *(listeVarGlob["vol_total2D_avec_plan_yz"]) << " " << *(listeVarGlob["vol_total2D_avec_plan_xz"])
              //       << " " << *(listeVarGlob["vol_total2D_avec_plan_xy"]) << endl;
              // fin débug
               };
              // globalisation des temps cpu loi de comp et métrique
              temps_lois_comportement += el.Temps_lois_comportement();
              temps_lois_comportement -= temps_debut_lois_comportement;
              temps_metrique_K_SM += el.Temps_metrique_K_SM();
              temps_metrique_K_SM -= temps_debut_metrique_K_SM;

////--- debug
//{
//cout << "\n Algori::RaidSmEner debug : ";
//string entete = " affichage de la matrice de raideur (puissance interne) après assemblage";
//matglob.Affichage_ecran(entete); 
//entete = " affichage du second membre (puissance interne) après assemblage ";
//vglobin.Affichage_ecran(entete); 
//}
//// fin debug

            };
        };
      }     
    catch (ErrJacobienNegatif_ElemMeca excep) 
         // cas d'un jacobien négatif "a prendre en compte"
       {
       #ifndef UTILISATION_MPI
         phase_de_convergence = -5;  a_converge = false; retour = false;
       #else
         gestion_exception = 1;
       #endif
       }
    catch (ErrVarJacobienMini_ElemMeca excep) 
         // cas d'une variation de  jacobien trop importante
       {
       #ifndef UTILISATION_MPI
         phase_de_convergence = -6;  a_converge = false; retour = false;
       #else
         gestion_exception = 2;
       #endif
       }
    catch (ErrNonConvergence_loiDeComportement excep)
         // cas d'une erreur survenue à cause d'une non convergence pour la résolution
         // d'une loi de comportement incrémentale
       {
       #ifndef UTILISATION_MPI
         switch (excep.cas)
       	  {case 1: phase_de_convergence = -8; break;
       	   case 0: default : phase_de_convergence = -7; break;
          };
         a_converge = false; retour = false;
       #else
         switch (excep.cas)
          {case 1: gestion_exception = 3; break;
          case 0: default : gestion_exception = 4; break;
         };
       #endif
       }
    catch (ErrCalculFct_nD)
         // cas d'une d'une erreur survenue à cause d'une erreur sur le calcul d'une fct nD
         // au niveau du chargement
       {
       #ifndef UTILISATION_MPI
         phase_de_convergence = -10;  a_converge = false;
         cout << "\n *** erreur exception genere par une fonction nD utilisee en chargement ";
         retour = false;
       #else
         gestion_exception = 5;
       #endif
       }
    catch (ErrSortie)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       {
       #ifndef UTILISATION_MPI
         tempsRaidSmEner.Arret_du_comptage();
         ErrSortie toto;
         throw (toto);
       #else
         gestion_exception = 6;
       #endif
       }
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       {
       #ifndef UTILISATION_MPI
         tempsRaidSmEner.Arret_du_comptage();
         ErrSortieFinale toto;
         throw (toto);
       #else
         gestion_exception = 7;
       #endif
       }
    catch ( ... )
       {
       #ifndef UTILISATION_MPI
         if (ParaGlob::NiveauImpression() >= 1)
           {cout << "\n warning: exception generee par un element mais dont la prise en compte "
                 << " n'est pas prevu !, on ne fait rien et on continue le calcul";
            if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n Algori::RaidSmEner(..";
            };
       #else
         gestion_exception = 8;
       #endif
       }

    tempsRaidSmEner.Arret_du_comptage(); // fin comptage cpu
  #ifdef UTILISATION_MPI
    // tout d'abord on gère le cas des exceptions
    if (proc_en_cours != 0)
      {if (gestion_exception != 0)
        {// on transmets l'info au proc 0 avec le même protocole que pour un résultat
         // on boucle sur tous les derniers éléments qui n'ont pas été calculé
         // ceci pour que proc 0 puisse correctement traiter tous les éléments
         const Tableau < list <int > >&  tab_list_elem_cpu = distribution_CPU_algo.List_element_CPU_en_cours();
         int nb_mail_distrib = tab_list_elem_cpu.Taille();
         for (int nbMail =dernier_nbMail; nbMail<= nb_mail_distrib; nbMail++)
           {const list <int >& list_elem_cpu= tab_list_elem_cpu(nbMail); // pour simplifier
            // on balaie les éléments nécessaires
            list <int >::const_iterator il,ilfin=list_elem_cpu.end();
            for (il = dernier_il;il != ilfin;il++)
              { // on récupère un signal du process 0
                // on considère qu'il y a eu au moins un passage
                tempsRaidSmEner.Arret_du_comptage();
                temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
                if (reqs1.active()) reqs1.wait();
                temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu
                temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
                DeuxEntiers num_el_et_mail(gestion_exception,-1);
                // on transmet les numéros d'élément et de maillage
                reqs1 = ParaGlob::Monde()->isend(0, 24, num_el_et_mail);
                temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
                tempsRaidSmEner.Mise_en_route_du_comptage();
              };
            };
        };
      };
    
    // récupération des grandeurs locales par le proc 0
    if (proc_en_cours == 0)
     { // récup du nombre total d'éléments, cumul sur tous les maillages
       int total_elem = distribution_CPU_algo.NB_total_element();
       // on va boucler sur les éléments récupérés des différents process
       // jusqu'au nombre maxi d'élément
       int nb_elem_deja_calcule = 0; // init
       while (nb_elem_deja_calcule < total_elem)
        { // on récupère un résultat de calcul
          temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
          DeuxEntiers num_el_et_mail;
          mpi::request reqs1 = ParaGlob::Monde()->irecv(mpi::any_source, 24, num_el_et_mail);
          temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
          temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
          mpi::status stat = reqs1.wait(); // on attend que le conteneur soit rempli
          temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu
          temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
          int  ne = num_el_et_mail.un;  // numero d'identification de l'element
          int  nbMail = num_el_et_mail.deux;  // numéro de maillage
          // on regarde s'il s'agit d'une exception
          if (nbMail > 0)
           // cas normal
           {// d'où l'élément
            ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne));
            // récupération des conteneurs ad hoc vecteur et raideur
            int source = stat.source(); // récupération du numéro de la source
            Vecteur * residu = el.Conteneur_Residu();
            mpi::request reqs2 = residu->Irecup_MPI(source, 25);
            Mat_pleine* raideur = el.Conteneur_raideur();
            mpi::request reqs3 = raideur->Irecup_MPI(source, 26);
            temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu

            temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
            reqs2.wait(); // on attend que le conteneur soit rempli
            reqs3.wait(); // on attend que le conteneur soit rempli
            temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu

            tempsRaidSmEner.Mise_en_route_du_comptage(); // comptage cpu
            Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
            // --- assemblage
            Ass.AssemSM (vglobin,*residu,el.TableauDdl(),taN); // du second membre
            if (pa.Symetrie_matrice())
              Ass.AssembMatSym (matglob,*raideur,el.TableauDdl(),taN); // de la raideur
            else
              Ass.AssembMatnonSym (matglob,*raideur,el.TableauDdl(),taN); // de la raideur
            tempsRaidSmEner.Arret_du_comptage(); // fin comptage cpu
           }
          else // la suite concerne le cas d'une exception
           {// s'il y a une exception fatal, cela va arrêter proc 0 et donc tous les autres
            gestion_exception = ne;
            switch (gestion_exception)
             {case 1: phase_de_convergence = -5;  a_converge = false; retour = false; break;
              case 2: phase_de_convergence = -6;  a_converge = false; retour = false; break;
              case 3: phase_de_convergence = -8;  retour = false; break;
              case 4: phase_de_convergence = -7;  retour = false; break;
              case 5: phase_de_convergence = -10;  a_converge = false;
                cout << "\n *** erreur exception genere par une fonction nD utilisee en chargement ";
                retour = false; break;
              case 6: {tempsRaidSmEner.Arret_du_comptage(); ErrSortie toto; throw (toto);}
              case 7: {tempsRaidSmEner.Arret_du_comptage(); ErrSortieFinale toto;; throw (toto);}
              case 8: { if (ParaGlob::NiveauImpression() >= 1)
                         {cout << "\n warning: exception generee par un element mais dont la prise en compte "
                               << " n'est pas prevu !, on ne fait rien et on continue le calcul";
                          if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n Algori::RaidSmEner(..";
                         };break;
                      }
              default:
                cout << "\n erreur en retour d'exception: on a lue gestion_exception = "<< gestion_exception
                     << " la gestion de ce cas n'existe pas !!! on stop l'execution ..."
                     << "\n Algori::RaidSmEner(..";
                Sortie(1);
            }; // fin de la gestion des différents cas d'exception
            // arrivée ici on continue la boucle du while
           tempsRaidSmEner.Arret_du_comptage(); // fin comptage cpu
           };
          // on incrémente le nombre d'élément traité
          nb_elem_deja_calcule++;
        };
     };
    // à la sortie du calcul on transmets la gestion d'excep
    // on transmet à tous les proc l'indicateur calcul_a_effectuer
    broadcast(*ParaGlob::Monde(), gestion_exception, 0);
    // on vérifie que ce n'est pas une sortie due à une exception
    // si c'est ok on continue
    if (gestion_exception != 0)
      retour = false;
  #endif
    if(retour)
      {tempsRaidSmEner.Mise_en_route_du_comptage(); // comptage cpu
       // -- calcul des intégrales éventuelles avec transfert en global
       lesMail->Integration();
  
  #ifdef UTILISATION_MPI
       tempsRaidSmEner.Arret_du_comptage(); // fin comptage cpu
       // récupération de toutes les énergies par le cpu 0
       Algori::Passage_energiesEtVolumes();
  #endif
  
       tempsRaidSmEner.Mise_en_route_du_comptage(); // fin comptage cpu
       // -- on transfert en global les énergies internes
       Transfert_ParaGlob_energies_internesLoisComp();
       Transfert_ParaGlob_energies_hourglass_bulk_stab();
       // idem pour les volumes entre plans
       Transfert_ParaGlob_volume_entre_plans();
    
       // affichage éventuelle de la matrice de raideur et du second membre
 #ifdef UTILISATION_MPI
    if (proc_en_cours == 0)
 #endif
     { if (ParaGlob::NiveauImpression() >= 10)
           { string entete = " affichage de la matrice de raideur (puissance interne) avant conditions limites";
             matglob.Affichage_ecran(entete);
             entete = " affichage du second membre (puissance interne) avant condition limites ";
             vglobin.Affichage_ecran(entete);
           };
       if (permet_affichage >3) cout << "\n -- fin calcul second membre et raideur " << flush;
     };
       tempsRaidSmEner.Arret_du_comptage(); // fin comptage cpu
      }; // fin du cas sans exception
    // retour indiquant si oui ou non il y a eu des exceptions
    return retour ;
  };         

// calcul  du second membre pour tous les maillages ainsi que des énergies totales
bool Algori::SecondMembreEnerg(LesMaillages * lesMail,Assemblage& Ass,Vecteur & vglobin) 
   {tempsSecondMembreEnerg.Mise_en_route_du_comptage(); // comptage cpu
    bool retour = true; // init par défaut
    // on doit initialiser les temps globaux de loi de comp et de métrique car
    // on va ajouter ensuite des temps qui eux mêmes globalisent ce qui se passe à chaque pti
//    Temps_CPU_HZpp tps_zero; // un temps d'initialisation
//    temps_lois_comportement = tps_zero;
    Temps_CPU_HZpp temps_debut_lois_comportement; // temps initial
//    temps_metrique_K_SM = tps_zero;
    Temps_CPU_HZpp temps_debut_metrique_K_SM; // temps initial
    if (permet_affichage >3) cout << "\n -- debut calcul second membre " << flush;

    energTotal.Inita(0.); // initialisation des énergies à 0
    energHourglass=0.;energStabiliMembBiel=0.;
    E_bulk = 0.; // init
    P_bulk=0.; // init
    volume_total_matiere = 0.; // init
    if (pa.CalVolTotalEntreSurfaceEtPlansRef())
      vol_total2D_avec_plan_ref.Inita(Coordonnee(ParaGlob::Dimension()));
        
  #ifdef UTILISATION_MPI
    mpi::request reqs1;
    mpi::request reqs2;
    Vecteur deux_faux_entiers(2);
    bool premier_passage = true;
    int proc_en_cours = ParaGlob::Monde()->rank();
    int gestion_exception = 0; // pour la gestion
  #endif
    // on gère les exceptions éventuelles en mettant le bloc sous surveillance
    int nbMailMax = lesMail->NbMaillage();
    int dernier_nbMail=1; // pour la gestion d'erreur
    list <int >::const_iterator dernier_il; // pour la gestion d'erreur
    try
      {// boucle sur les elements
       int nbMailMax = lesMail->NbMaillage();
      #ifndef UTILISATION_MPI
       // cas d'un  calcul mono CPU
       for (int nbMail =1; nbMail<= nbMailMax; nbMail++,dernier_nbMail++)
         {int nemax = lesMail->Nombre_element(nbMail);
          for (int ne=1; ne<= nemax;ne++)
             {
      #else
        // cas d'un calcul multi-CPU
        // on va parcourir les éléments associés au cpu
        // on récupère la distribution d'éléments concernant le cpu en cours
        const Tableau < list <int > > & tab_list_elem_cpu = distribution_CPU_algo.List_element_CPU_en_cours();
        // la taille est identique à nbMailMax, sauf si c'est le cpu 0, là le tableau est vide et il n'y
        // aura pas de boucle, (ou plutôt on sortira directement de la boucle
        // le cas CPU 0 est traité à part
        int nb_mail_distrib = tab_list_elem_cpu.Taille();
        if (proc_en_cours != 0)
         for (int nbMail =1; nbMail<= nb_mail_distrib; nbMail++)
          {const list <int >& list_elem_cpu= tab_list_elem_cpu(nbMail); // pour simplifier
           // on balaie les éléments nécessaires
           list <int >::const_iterator il,ilfin=list_elem_cpu.end();
           for (il = dernier_il = list_elem_cpu.begin();il != ilfin;il++,dernier_il++)
             { int ne = (*il);
               // on récupère un signal du process 0
               temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
               if (!premier_passage) // si ce n'est pas le premier passage
                // on vérifie que le premier conteneur a bien été lue
                 {if (reqs1.active()) reqs1.wait();
                  if (reqs2.active()) reqs2.wait();
                 };
               temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu

       #endif
              //calcul  du residu local
              ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne)); // l'element
              Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
              Temps_CPU_HZpp temps_debut_lois_comportement=el.Temps_lois_comportement(); // init
              Temps_CPU_HZpp temps_debut_metrique_K_SM = el.Temps_metrique_K_SM(); // init
              Vecteur* res = el.CalculResidu_tdt(pa);
////--- debug
//cout << "\n debug .... Algori::SecondMembreEnerg( "
//         << "\n CPU= "<< ParaGlob::Monde()->rank();
//    // vecteur résidu
//cout <<  "\n res  = " <<    *res ;
//
////--- fin debug

////--- debug
//if (res->Max_val_abs() > 10)
//  { cout << "\n debug .... Algori::SecondMembreEnerg( "
//         << "\n **** res->Max_val_abs() > 10 ";
//    // on recalcule pour le debug le vecteur résidu
//    res = el.CalculResidu_tdt(pa);
//  };
//
////--- fin debug
            #ifndef UTILISATION_MPI
              // --- assemblage et calcul de grandeur globale, pour mono-CPU
              // assemblage
              Ass.AssemSM (vglobin,*res,el.TableauDdl(),taN); // du second membre
            #else
              // cas d'un calcul parallèle, et CPU != 0
              tempsSecondMembreEnerg.Arret_du_comptage();
              temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
//              DeuxEntiers num_el_et_mail(el.Num_elt(),el.Num_maillage());
              deux_faux_entiers(1) = (double ) el.Num_elt();
              deux_faux_entiers(2) = (double ) el.Num_maillage();
              // on transmet les numéros d'élément et de maillage
              reqs1 = deux_faux_entiers.Ienvoi_MPI(0, 24111);
//              reqs1 = ParaGlob::Monde()->isend(0, 24, num_el_et_mail);
              temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
              // puis on transmets le vecteur résidu
              temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
              temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu
              temps_transfert_long_matSm.Mise_en_route_du_comptage(); // comptage cpu
              reqs2 = res->Ienvoi_MPI(0,25111);
              if (premier_passage) {premier_passage = false;};
              temps_transfert_long_matSm.Arret_du_comptage(); // fin comptage cpu
              tempsSecondMembreEnerg.Mise_en_route_du_comptage();
            #endif
////--- debug
//  { cout << "\n debug .... Algori::SecondMembreEnerg( ";
//    // on recalcule pour le debug le vecteur résidu
//    cout << "\n second membre local = " << (*res);
////    res = el.CalculResidu_tdt(pa);
//  };

//--- fin debug
              // globalisation des énergies
              energTotal += el.EnergieTotaleElement();
              energHourglass += el.Energie_Hourglass();
              energStabiliMembBiel += el.Energie_stab_membBiel();
              E_bulk += el.Energie_Bulk();
              P_bulk += el.Puissance_Bulk();
              // globalisation des volumes
              volume_total_matiere += el.Volume();
              if (pa.CalVolTotalEntreSurfaceEtPlansRef())
                { const Coordonnee& volPlan = el.VolumePlan();
                  if (volPlan.Dimension())
                   vol_total2D_avec_plan_ref(nbMail) += volPlan;
                };
              // globalisation des temps cpu loi de comp et métrique
              temps_lois_comportement += el.Temps_lois_comportement();
              temps_lois_comportement -= temps_debut_lois_comportement;
              temps_metrique_K_SM += el.Temps_metrique_K_SM();
              temps_metrique_K_SM -= temps_debut_metrique_K_SM;
              };
          };
      }
    catch (ErrJacobienNegatif_ElemMeca excep) 
         // cas d'un jacobien négatif "a prendre en compte"
       {
       #ifndef UTILISATION_MPI
         phase_de_convergence = -5;  a_converge = false; retour = false;
       #else
         gestion_exception = 1;
       #endif
       }
    catch (ErrVarJacobienMini_ElemMeca excep)
         // cas d'une variation de  jacobien trop importante
       {
       #ifndef UTILISATION_MPI
         phase_de_convergence = -6;  a_converge = false; retour = false;
       #else
         gestion_exception = 2;
       #endif
       }
    catch (ErrNonConvergence_loiDeComportement excep)
         // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
         // d'une loi de comportement incrémentale
       {
       #ifndef UTILISATION_MPI
         cout << "\n non convergence sur une loi de comportement ";
         switch (excep.cas)
          {case 1: phase_de_convergence = -8; break;
           case 0: default : phase_de_convergence = -7; break;
          };
         a_converge = false; retour = false;
       #else
         switch (excep.cas)
          {case 1: gestion_exception = 3; break;
          case 0: default : gestion_exception = 4; break;
         };
       #endif
       }
    catch (ErrCalculFct_nD)
         // cas d'une d'une erreur survenue à cause d'une erreur sur le calcul d'une fct nD
         // au niveau du chargement
       {
       #ifndef UTILISATION_MPI
         phase_de_convergence = -10;  a_converge = false;
         cout << "\n *** erreur exception genere par une fonction nD utilisee en chargement ";
         retour = false;
       #else
         gestion_exception = 5;
       #endif
       }
    catch (ErrSortie)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
      {
       #ifndef UTILISATION_MPI
         tempsSecondMembreEnerg.Arret_du_comptage();
         ErrSortie toto;
         throw (toto);
       #else
         gestion_exception = 6;
       #endif
       }
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
          {
          #ifndef UTILISATION_MPI
            tempsSecondMembreEnerg.Arret_du_comptage();
            ErrSortieFinale toto;
            throw (toto);
          #else
            gestion_exception = 7;
          #endif
          }
    catch ( ... )
          {
          #ifndef UTILISATION_MPI
            if (ParaGlob::NiveauImpression() >= 1)
              {cout << "\n warning: exception generee par un element mais dont la prise en compte "
                    << " n'est pas prevu !, on ne fait rien et on continue le calcul";
               if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n Algori::SecondMembreEnerg(..";
               };
          #else
            gestion_exception = 8;
          #endif
          }


  tempsSecondMembreEnerg.Arret_du_comptage(); // fin comptage cpu
#ifdef UTILISATION_MPI
  // tout d'abord on gère le cas des exceptions
  if (proc_en_cours != 0)
    {if (gestion_exception != 0)
      {// on transmets l'info au proc 0 avec le même protocole que pour un résultat
       // on boucle sur tous les derniers éléments qui n'ont pas été calculé
       // ceci pour que proc 0 puisse correctement traiter tous les éléments
       const Tableau < list <int > >&  tab_list_elem_cpu = distribution_CPU_algo.List_element_CPU_en_cours();
       int nb_mail_distrib = tab_list_elem_cpu.Taille();
       for (int nbMail =dernier_nbMail; nbMail<= nb_mail_distrib; nbMail++)
         {const list <int >& list_elem_cpu= tab_list_elem_cpu(nbMail); // pour simplifier
          // on balaie les éléments nécessaires
          list <int >::const_iterator il,ilfin=list_elem_cpu.end();
          for (il = dernier_il;il != ilfin;il++)
            { // on récupère un signal du process 0
              // on considère qu'il y a eu au moins un passage
              tempsSecondMembreEnerg.Arret_du_comptage();
              temps_attente_matSm.Mise_en_route_du_comptage(); // comptage cpu
              if (reqs1.active()) reqs1.wait();
              temps_attente_matSm.Arret_du_comptage(); // fin comptage cpu
              temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
              deux_faux_entiers(1) = gestion_exception;
              deux_faux_entiers(2) = -1;
              // on transmet au proc 0
              reqs1 = deux_faux_entiers.Ienvoi_MPI(0, 24111);
              temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
              tempsSecondMembreEnerg.Mise_en_route_du_comptage();
            };
          };
      };
    };

    if (proc_en_cours == 0)
      { tempsSecondMembreEnerg.Mise_en_route_du_comptage(); // comptage cpu
        // récup du nombre total d'éléments, cumul sur tous les maillages
        int total_elem = distribution_CPU_algo.NB_total_element();
        // on va boucler sur les éléments récupérés des différents process
        // jusqu'au nombre maxi d'élément
        int nb_elem_deja_calcule = 0; // init
        while (nb_elem_deja_calcule < total_elem)
         { tempsSecondMembreEnerg.Arret_du_comptage(); // fin comptage cpu
            // on récupère un résultat de calcul
           temps_transfert_court_matSm.Mise_en_route_du_comptage(); // comptage cpu
//           DeuxEntiers num_el_et_mail;
//           mpi::request reqs1 = ParaGlob::Monde()->irecv(mpi::any_source, 24, num_el_et_mail);
           reqs1 = deux_faux_entiers.Irecup_MPI(mpi::any_source, 24111);
           mpi::status stat = reqs1.wait(); // on attend que le conteneur soit rempli
           temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
           tempsSecondMembreEnerg.Mise_en_route_du_comptage(); // comptage cpu

//           temps_attente.Mise_en_route_du_comptage(); // comptage cpu
//           mpi::status stat = reqs1.wait(); // on attend que le conteneur soit rempli
//           temps_attente.Arret_du_comptage(); // fin comptage cpu
           int  ne = (int) deux_faux_entiers(1); //num_el_et_mail.un;  // numero d'identification de l'element
           int  nbMail = (int) deux_faux_entiers(2);//num_el_et_mail.deux;  // numéro de maillage
           if (nbMail > 0)
             // cas normal
             {
              // d'où l'élément
              ElemMeca & el = *((ElemMeca*) &lesMail->Element_LesMaille(nbMail,ne));
              // récupération des conteneurs ad hoc vecteur et raideur
              int source = stat.source(); // récupération du numéro de la source
              Vecteur * residu = el.Conteneur_Residu();
              tempsSecondMembreEnerg.Arret_du_comptage();
              temps_transfert_long_matSm.Mise_en_route_du_comptage(); // comptage cpu
              mpi::request reqs2 = residu->Irecup_MPI(source, 25111);
              reqs2.wait(); // on attend que le conteneur soit rempli
              temps_transfert_long_matSm.Arret_du_comptage(); // fin comptage cpu
              tempsSecondMembreEnerg.Mise_en_route_du_comptage();
              Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
            #ifdef MISE_AU_POINT
              if (ParaGlob::NiveauImpression() > 9)
               {cout << "\n proc0: , retour du proc "<<source ;
                cout << "\n Second membre local: elem= " << ne
                     << ", maillage= " << nbMail;
                cout << "\n second membre: " << (*residu);
              };
           #endif

              // --- assemblage
              Ass.AssemSM (vglobin,*residu,el.TableauDdl(),taN); // du second membre
             }
           else // la suite concerne le cas d'une exception
            {// s'il y a une exception fatal, cela va arrêter proc 0 et donc tous les autres
             gestion_exception = ne;
             switch (gestion_exception)
              {case 1: phase_de_convergence = -5;  a_converge = false; retour = false; break;
               case 2: phase_de_convergence = -6;  a_converge = false; retour = false; break;
               case 3: phase_de_convergence = -8;  retour = false; break;
               case 4: phase_de_convergence = -7;  retour = false; break;
               case 5: phase_de_convergence = -10;  a_converge = false;
                 cout << "\n *** erreur exception genere par une fonction nD utilisee en chargement ";
                 retour = false; break;
               case 6: {tempsSecondMembreEnerg.Arret_du_comptage(); ErrSortie toto; throw (toto);}
               case 7: {tempsSecondMembreEnerg.Arret_du_comptage(); ErrSortieFinale toto;; throw (toto);}
               case 8: { if (ParaGlob::NiveauImpression() >= 1)
                          {cout << "\n warning: exception generee par un element mais dont la prise en compte "
                                << " n'est pas prevu !, on ne fait rien et on continue le calcul";
                           if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n Algori::SecondMembreEnerg(..";
                          };break;
                       }
               default:
                 cout << "\n erreur en retour d'exception: on a lue gestion_exception = "<< gestion_exception
                      << " la gestion de ce cas n'existe pas !!! on stop l'execution ..."
                      << "\n Algori::SecondMembreEnerg(..";
                 Sortie(1);
             }; // fin de la gestion des différents cas d'exception
             // arrivée ici on continue la boucle du while
            tempsSecondMembreEnerg.Arret_du_comptage(); // fin comptage cpu
            };
           // on incrémente le nombre d'élément traité
           nb_elem_deja_calcule++;
         };
        tempsSecondMembreEnerg.Arret_du_comptage(); // fin comptage cpu
      };
    // à la sortie du calcul on transmets la gestion d'excep
    // on transmet à tous les proc l'indicateur calcul_a_effectuer
    broadcast(*ParaGlob::Monde(), gestion_exception, 0);
    // on vérifie que ce n'est pas une sortie due à une exception
    // si c'est ok on continue
    if (gestion_exception != 0)
      retour = false;
  #endif
    if(retour)
      {tempsSecondMembreEnerg.Mise_en_route_du_comptage(); // comptage cpu
       // -- calcul des intégrales éventuelles avec transfert en global
       lesMail->Integration();
  
  #ifdef UTILISATION_MPI
       tempsSecondMembreEnerg.Arret_du_comptage(); // fin comptage cpu
       // récupération de toutes les énergies par le cpu 0
       Algori::Passage_energiesEtVolumes();
  #endif
  
    tempsSecondMembreEnerg.Mise_en_route_du_comptage();
    // -- on transfert en global les énergies internes
    Transfert_ParaGlob_energies_internesLoisComp();
    Transfert_ParaGlob_energies_hourglass_bulk_stab();
    // idem pour les volumes entre plans
    Transfert_ParaGlob_volume_entre_plans();

    // affichage éventuelle  du second membre
 #ifdef UTILISATION_MPI
    if (proc_en_cours == 0)
 #endif
     {if (ParaGlob::NiveauImpression() >= 10)
        { string entete = " affichage du second membre (puissance interne) avant conditions limites ";
          vglobin.Affichage_ecran(entete);
        };
      if (permet_affichage >3) cout << "\n -- fin calcul second membre " << flush;
     };
       tempsSecondMembreEnerg.Arret_du_comptage(); // fin comptage cpu
       }; // fin du cas sans exception

    // retour indiquant si oui ou non il y a eu des exceptions
    return retour ;
   };
   
// Mise à jour  de la raideur et du second membre pour le contact
// calcul des énergies spécifiques au contact
// si pb retour false: indique que le calcul c'est mal passé, il y a une erreur
// généré, les résultats: raideur et second membres calculées sont inexploitables
bool Algori::RaidSmEnerContact(LesContacts * lesCont,Assemblage& Ass,Vecteur & vglobin
                ,Mat_abstraite& matglob )
  {
    tempsRaidSmEnerContact.Mise_en_route_du_comptage(); // temps cpu
    energFrottement.Inita(0.); // initialisation des énergies de frottement à 0
    energPenalisation = 0.;
    if (permet_affichage >3) cout << "\n -- debut calcul second membre et raideur pour le  contact " << flush;

  #ifdef UTILISATION_MPI
    mpi::request reqs1;
    mpi::request reqs2;
    mpi::request reqs3;
    bool premier_passage = true;

    // mise à jour de la liste des contacts pour proc 0: il s'agit des infos minimales
    // qui permettent ensuite de faire un assemblage
    lesCont->Mise_a_jour_liste_contacts_actif_interProc();
    int proc_en_cours = ParaGlob::Monde()->rank();
    // en // ce sont les proc i>0 qui gèrent les éléments de contact
    if (proc_en_cours != 0)
     {
  #endif
    // on récupère la liste des éléments de contact
    LaLIST <ElContact>& listElContact = lesCont->LesElementsDeContact(); 
    LaLIST <ElContact>::iterator il,ilfin=listElContact.end();
    // on gère les exceptions éventuelles en mettant le bloc sous surveillance
    try 
      {// boucle sur les elements de contact
       for (il=listElContact.begin();il != ilfin; il++)
         { 
////------ debug
//        cout << "\n debug Algori::RaidSmEnerContact ";
//// ------ fin debug
          if ((*il).Actif()) // on n'intervient que si le contact est actif
          {	//calcul de la raideur locale et du residu local
      //------ debug
      //        (*il).Affiche();
      // ------ fin debug		  			 
         #ifdef UTILISATION_MPI
            // on récupère un signal du process 0
            temps_attente_contact.Mise_en_route_du_comptage(); // comptage cpu
            if (premier_passage) {premier_passage = false;}
            else // on regarde l'activité , car au début avant le balayage de tous les éléments des tableaux
               {if (reqs1.active()) reqs1.wait(); // les requests ne sont pas alimentés
                if (reqs2.active()) reqs2.wait(); // car aucun transfert n'a encore été effectué
                if (reqs3.active()) reqs3.wait(); //
               };
           temps_attente_contact.Arret_du_comptage();
         #endif
         
           ElContact& elcontact = (*il); // pour simplifier
           Tableau<Noeud *>& taN = elcontact.TabNoeud_pour_assemblage(); // tableau de noeuds 
           // calcul effectif de la raideur et résidu locals
           Element::ResRaid resu = elcontact.SM_K_charge_contact();
//   cout << "numéro =" << ne << endl;
//------ debug
//              (resu.res)->Affiche();
//              (resu.raid)->Affiche(); //1,8,1,1,8,1);
//              (resu.raid)->Affiche1(1,6,1,1,6,1);
        //      (int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd)
// ------ fin debug
         #ifndef UTILISATION_MPI
           // assemblage
           // dans le cas où le retour est un pointeur nul, cela signifie qu'il n'y a pas de second membre et raideur calculés
           if (resu.res != NULL)
            {  Ass.AssemSM (vglobin,*(resu.res),elcontact.TableauDdlCont(),taN); // du second membre
             if (pa.Symetrie_matrice())
              Ass.AssembMatSym (matglob,*(resu.raid),elcontact.TableauDdlCont(),taN); // de la raideur
             else
              Ass.AssembMatnonSym (matglob,*(resu.raid),elcontact.TableauDdlCont(),taN); // de la raideur
            };
         #else
           // cas d'un calcul parallèle, et CPU != 0
           {tempsRaidSmEnerContact.Arret_du_comptage();
            // on va transmettre les infos permettant au proc 0 de récupérer
            // l'élément de contact ad hoc et de préparer la place
            temps_transfert_court_contact.Mise_en_route_du_comptage(); // comptage cpu
            // on commence par transmettre les infos de l'élément de contact
            QuatreEntiers nums;
            nums.un = elcontact.Esclave()->Num_Mail();
            nums.deux = elcontact.Esclave()->Num_noeud();
            nums.trois = elcontact.Elfront()->NumUnique();
            if (resu.res != NULL)
              {nums.quatre = resu.res->Taille();}
            else // sinon resu.res == NULL, on envoie un signal via une taille négative
              {nums.quatre = -1;}

            // on transmet les infos permettants de récupérer l'élément de contact
            reqs1 = ParaGlob::Monde()->isend(0, 240, nums);
            if (resu.res != NULL)
              {// puis on transmets le vecteur résidu
               reqs2 = resu.res->Ienvoi_MPI(0,250);
               //puis la matrice
               reqs3 = resu.raid->Ienvoi_MPI(0, 260);
               };
            temps_transfert_court_contact.Arret_du_comptage(); // fin comptage cpu
            tempsRaidSmEnerContact.Mise_en_route_du_comptage();
           };

         #endif

           // globalisation des énergies dues aux frottements
           energFrottement += elcontact.EnergieFrottement(); 
           energPenalisation += elcontact.EnergiePenalisation();
          }; //-- fin du test sur l'activité
         }; //-- fin de la boucle sur les éléments de contact
      } //-- fin du try
    catch (ErrSortie)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortie toto;
         throw (toto);
       }
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortieFinale toto;
         tempsRaidSmEnerContact.Arret_du_comptage(); // temps cpu
         throw (toto);
       }
    catch ( ... )
       { if (ParaGlob::NiveauImpression() >= 1) 
           {cout << "\n warning: exception generee par un element de contact mais dont la prise en compte "
                 << " n'est pas prevu !, on ne fait rien et on continue le calcul";
            if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n Algori::RaidSmEnerContact(..";
            };
        }
  #ifdef UTILISATION_MPI
     } // fin du cas des proc i>0
  #endif

        
        
        
    // affichage éventuelle de la matrice de raideur et du second membre     
//------ debug
//        cout << "\n debug Algori::RaidSmEnerContact ";
//         string entete = " affichage du second membre  uniquement du au  contact ";
//         vglobin.Affichage_ecran(entete);
// ------ fin debug
    if (ParaGlob::NiveauImpression() >= 10)
        { string entete = " affichage de la matrice de raideur (puissance interne) apres le contact ";
          matglob.Affichage_ecran(entete); 
          entete = " affichage du second membre  uniquement du au  contact ";
          vglobin.Affichage_ecran(entete); 
         }
    // -- on transfert en global les énergies
    Transfert_ParaGlob_energies_contact();
    // affichage éventuelle de la force maxi de contact
    // pas nécessaire ici car c'est fait dans le calcul des réactions de contact
//    lesCont->Forces_contact_maxi(true);
    // idem pour les gap N et tangentiel
//    lesCont->Gap_contact_maxi(true);
    
    if (permet_affichage >3) cout << "\n -- fin calcul second mmembre pour le  contact " << flush;

    if (permet_affichage >3) cout << "\n -- fin calcul second membre et raideur pour le  contact " << flush;

    tempsRaidSmEnerContact.Arret_du_comptage(); // temps cpu
    
#ifdef UTILISATION_MPI
    // récupération des grandeurs locales par le proc 0
    if (ParaGlob::Monde()->rank() == 0)
      { tempsRaidSmEnerContact.Mise_en_route_du_comptage(); // comptage cpu
        // récup du nombre total d'éléments de contact, cumulé sur tous les maillages
        int total_elem_contact  = lesCont->Nb_actifs();
        // on va boucler sur les éléments récupérés des différents process
////---debug
//cout << "\n debug Algori::RaidSmEnerContact "
//     << "\n total_elem_contact actif = "<<total_elem_contact;
////---debug
        // jusqu'au nombre maxi d'élément
        int nb_elem_deja_calcule = 0; // init
        while (nb_elem_deja_calcule < total_elem_contact)
         { // on récupère un résultat de calcul
           temps_transfert_court_contact.Mise_en_route_du_comptage(); // comptage cpu
           QuatreEntiers nums;
           mpi::request reqs1 = ParaGlob::Monde()->irecv(mpi::any_source, 240, nums);
           temps_transfert_court_contact.Arret_du_comptage(); // fin comptage cpu
           temps_attente_contact.Mise_en_route_du_comptage(); // comptage cpu
           mpi::status stat = reqs1.wait(); // on attend que le conteneur soit rempli
           temps_attente_contact.Arret_du_comptage(); // fin comptage cpu
           int numMailEsclave = nums.un;
           int numNoeudEsclave = nums.deux;
           int numUnique = nums.trois;
           int taille_res = nums.quatre;
           // on récupére l'élément de contact
           ElContact*  elcontact = lesCont->RecupElContactActif
                  (numMailEsclave,numNoeudEsclave,numUnique);
           if (elcontact == NULL )
             {cout << "\n erreur *** en calcul parallele, le proc 0 ne parvient pas a recuperer"
                    << " l'element de contact pour les infos suivantes: "
                    << "\n numMailEsclave= "<< numMailEsclave
                    << " numNoeudEsclave= "<< numNoeudEsclave
                    << " numUnique= "<< numUnique
                    << "\n la suite n'est pas possible ";
              Sortie(1);
             };
           // récupération des conteneurs ad hoc vecteur et raideur
           int source = stat.source(); // récupération du numéro de la source
           
           if (nums.quatre > 0) // s'il est négatif => il n'y a pas eu de résultat venant de l'élément de contact
             {// récup uniquement des conteneurs raideurs et résidu (pas forcément remplis, mais de la bonne taille)
              Element::ResRaid  resRaid = elcontact->Conteneur_ResRaid();
              Vecteur * residu = resRaid.res;
              if (residu->Taille() != taille_res )
                {cout << "\n erreur *** en calcul parallele, le proc 0 ne parvient pas a recuperer"
                       << " la bonne taille pour le résidu de contact avec les infos suivantes: "
                       << "\n numMailEsclave= "<< numMailEsclave
                       << " numNoeudEsclave= "<< numNoeudEsclave
                       << " numUnique= "<< numUnique
                       << "\n taille_residu transmis : "<<taille_res
                       << " taille actuelle pour le proc 0 "<< residu->Taille()
                       << "\n ce n'est pas normal ==> erreur "
                       << "\n on arrete ";
                 Sortie(1);
                };

              temps_transfert_court_contact.Mise_en_route_du_comptage(); // comptage cpu
              mpi::request reqs2 = residu->Irecup_MPI(source, 250);
              Mat_pleine* raideur = resRaid.raid;
              mpi::request reqs3 = raideur->Irecup_MPI(source, 260);
              temps_transfert_court_contact.Arret_du_comptage(); // fin comptage cpu

              temps_attente_contact.Mise_en_route_du_comptage(); // comptage cpu
              reqs2.wait(); // on attend que le conteneur soit rempli
              reqs3.wait(); // on attend que le conteneur soit rempli
              temps_attente_contact.Arret_du_comptage(); // fin comptage cpu
              tempsRaidSmEnerContact.Mise_en_route_du_comptage(); // temps cpu
              // assemblage
              Tableau<Noeud *>& taN = elcontact->TabNoeud_pour_assemblage(); // tableau de noeuds
              Ass.AssemSM (vglobin,*(residu),elcontact->TableauDdlCont(),taN); // du second membre
              if (pa.Symetrie_matrice())
                Ass.AssembMatSym (matglob,*(raideur),elcontact->TableauDdlCont(),taN); // de la raideur
              else
                Ass.AssembMatnonSym (matglob,*(raideur),elcontact->TableauDdlCont(),taN); // de la raideur
              tempsRaidSmEnerContact.Arret_du_comptage(); // fin comptage cpu
             };
           // on incrémente le nombre d'élément traité
           nb_elem_deja_calcule++;
         };
        tempsRaidSmEnerContact.Arret_du_comptage(); // fin comptage cpu
        //
      };
      
    // récupération de toutes les énergies par le cpu 0
    tempsRaidSmEnerContact.Arret_du_comptage(); // fin comptage cpu
    Algori::Passage_energiesContact();
    tempsRaidSmEnerContact.Mise_en_route_du_comptage();
#endif

    tempsRaidSmEnerContact.Mise_en_route_du_comptage();
    // -- on transfert en global les énergies internes
    Transfert_ParaGlob_energies_contact();
    tempsRaidSmEnerContact.Arret_du_comptage(); // fin comptage cpu
    
    // retour indiquant que tout c'est bien passé
    return true ;                           
  };

// cas du contact: calcul  et mise à jour du second membre ainsi que des énergies spécifiques  
bool Algori::SecondMembreEnergContact(LesContacts * lesCont,Assemblage& Ass,Vecteur & vglobin
                                      ,bool aff_iteration) 
  { tempsSecondMembreEnergContact.Mise_en_route_du_comptage(); // temps cpu
    energFrottement.Inita(0.); // initialisation des énergies de frottement à 0
    energPenalisation = 0.;
    volume_total_matiere = 0.; // init
    if (permet_affichage >3) cout << "\n -- debut calcul second membre pour le  contact " << flush;

  #ifdef UTILISATION_MPI
    mpi::request reqs1;
    mpi::request reqs2;
    bool premier_passage = true;

    int proc_en_cours = ParaGlob::Monde()->rank();
    // mise à jour de la liste des contacts pour proc 0: il s'agit des infos minimales
    // qui permettent ensuite de faire un assemblage
    lesCont->Mise_a_jour_liste_contacts_actif_interProc();

    // en // ce sont les proc i>0 qui gèrent les éléments de contact
    if (proc_en_cours != 0)
     {
  #endif
    // on récupère la liste des éléments de contact
    LaLIST <ElContact>& listElContact = lesCont->LesElementsDeContact(); 
    LaLIST <ElContact>::iterator il,ilfin=listElContact.end();
    // on gère les exceptions éventuelles en mettant le bloc sous surveillance
    try 
      {// boucle sur les elements de contact
       for (il=listElContact.begin();il != ilfin; il++)
         { if ((*il).Actif()) // on n'intervient que si le contact est actif
            {
          #ifdef UTILISATION_MPI
             // on récupère un signal du process 0
             temps_attente_contact.Mise_en_route_du_comptage(); // comptage cpu
             if (premier_passage) {premier_passage = false;}
             else // on regarde l'activité , car au début avant le balayage de tous les éléments des tableaux
                {if (reqs1.active()) reqs1.wait(); // les requests ne sont pas alimentés
                 if (reqs2.active()) reqs2.wait(); // car aucun transfert n'a encore été effectué
                };
            temps_attente_contact.Arret_du_comptage();
          #endif
           
              //calcul du residu local
              ElContact& elcontact = (*il); // pour simplifier
              Tableau<Noeud *>& taN = elcontact.TabNoeud_pour_assemblage(); // tableau de noeuds 
              //calcul  du residu local
              Vecteur* res = elcontact.SM_charge_contact();
           #ifndef UTILISATION_MPI
              // assemblage
              Ass.AssemSM (vglobin,*res,elcontact.TableauDdlCont(),taN); // du second membre
           #else
             // cas d'un calcul parallèle, et CPU != 0
             {tempsSecondMembreEnergContact.Arret_du_comptage();
              // on va transmettre les infos permettant au proc 0 de récupérer
              // l'élément de contact ad hoc et de préparer la place
              temps_transfert_long_contact.Mise_en_route_du_comptage(); // comptage cpu
              // on commence par transmettre les infos de l'élément de contact
              QuatreEntiers nums;
              nums.un = elcontact.Esclave()->Num_Mail();
              nums.deux = elcontact.Esclave()->Num_noeud();
              nums.trois = elcontact.Elfront()->NumUnique();
              nums.quatre = res->Taille();
              // on transmet les infos permettants de récupérer l'élément de contact
              reqs1 = ParaGlob::Monde()->isend(0, 240, nums);
              // puis on transmets le vecteur résidu
              reqs2 = res->Ienvoi_MPI(0,250);
              temps_transfert_court_matSm.Arret_du_comptage(); // fin comptage cpu
              tempsRaidSmEnerContact.Mise_en_route_du_comptage();
             };

           #endif

              // globalisation des énergies dues aux frottements
              energFrottement += elcontact.EnergieFrottement(); 
              energPenalisation += elcontact.EnergiePenalisation(); 
            }; //-- fin du test sur l'activité
         }; //-- fin de la boucle sur les éléments de contact
      } //-- fin du try
    catch (ErrSortie)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortie toto;
         throw (toto);
       }
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortieFinale toto;
         tempsSecondMembreEnergContact.Arret_du_comptage(); // temps cpu
         throw (toto);
       }
    catch ( ... )
       { if (ParaGlob::NiveauImpression() >= 1) 
           {cout << "\n warning: exception generee par un element de contact mais dont la prise en compte "
                 << " n'est pas prevu !, on ne fait rien et on continue le calcul";
            if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n Algori::SecondMembreEnergContact(..";
           };
        }                

   #ifdef UTILISATION_MPI
      } // fin du cas des proc i>0
   #endif

    // affichage éventuelle  du second membre
    if (ParaGlob::NiveauImpression() >= 10)
     { string entete = " affichage du second membre   uniquement du au  contact  ";
       vglobin.Affichage_ecran(entete); 
     };

    tempsSecondMembreEnergContact.Arret_du_comptage(); // temps cpu
          
#ifdef UTILISATION_MPI
    // récupération des grandeurs locales par le proc 0
    if (ParaGlob::Monde()->rank() == 0)
      { tempsSecondMembreEnergContact.Mise_en_route_du_comptage(); // comptage cpu
        // récup du nombre total d'éléments de contact, cumulé sur tous les maillages
        int total_elem_contact  = lesCont->Nb_actifs();
        // on va boucler sur les éléments récupérés des différents process
        // jusqu'au nombre maxi d'élément
        int nb_elem_deja_calcule = 0; // init
        while (nb_elem_deja_calcule < total_elem_contact)
         { // on récupère un résultat de calcul
           temps_transfert_court_contact.Mise_en_route_du_comptage(); // comptage cpu
           QuatreEntiers nums;
           mpi::request reqs1 = ParaGlob::Monde()->irecv(mpi::any_source, 240, nums);
           temps_transfert_court_contact.Arret_du_comptage(); // fin comptage cpu
           temps_attente_contact.Mise_en_route_du_comptage(); // comptage cpu
           mpi::status stat = reqs1.wait(); // on attend que le conteneur soit rempli
           temps_attente_contact.Arret_du_comptage(); // fin comptage cpu
           int numMailEsclave = nums.un;
           int numNoeudEsclave = nums.deux;
           int numUnique = nums.trois;
           int taille_res = nums.quatre;
           // on récupére l'élément de contact
           ElContact*  elcontact = lesCont->RecupElContactActif
                  (numMailEsclave,numNoeudEsclave,numUnique);
           if (elcontact == NULL )
             {cout << "\n erreur *** en calcul parallele, le proc 0 ne parvient pas a recuperer"
                    << " l'element de contact pour les infos suivantes: "
                    << "\n numMailEsclave= "<< numMailEsclave
                    << " numNoeudEsclave= "<< numNoeudEsclave
                    << " numUnique= "<< numUnique
                    << "\n la suite n'est pas possible ";
              Sortie(1);
             };
           // récupération des conteneurs ad hoc vecteur et raideur
           int source = stat.source(); // récupération du numéro de la source
           
           // récup uniquement du conteneur résidu (pas forcément rempli, mais de la bonne taille)
           Vecteur* residu = elcontact->Conteneur_Residu();

           if (residu->Taille() != taille_res )
             {cout << "\n erreur *** en calcul parallele, le proc 0 ne parvient pas a recuperer"
                    << " la bonne taille pour le résidu de contact avec les infos suivantes: "
                    << "\n numMailEsclave= "<< numMailEsclave
                    << " numNoeudEsclave= "<< numNoeudEsclave
                    << " numUnique= "<< numUnique
                    << "\n taille_residu transmis : "<<taille_res
                    << " taille actuelle pour le proc 0 "<< residu->Taille()
                    << "\n ce n'est pas normal ==> erreur "
                    << "\n on arrete ";
              Sortie(1);
             };

           temps_transfert_court_contact.Mise_en_route_du_comptage(); // comptage cpu
           mpi::request reqs2 = residu->Irecup_MPI(source, 250);
           temps_transfert_court_contact.Arret_du_comptage(); // fin comptage cpu

           temps_attente_contact.Mise_en_route_du_comptage(); // comptage cpu
           reqs2.wait(); // on attend que le conteneur soit rempli
           temps_attente_contact.Arret_du_comptage(); // fin comptage cpu
           tempsSecondMembreEnergContact.Mise_en_route_du_comptage(); // temps cpu
           // assemblage
           Tableau<Noeud *>& taN = elcontact->TabNoeud_pour_assemblage(); // tableau de noeuds
           Ass.AssemSM (vglobin,*(residu),elcontact->TableauDdlCont(),taN); // du second membre
           tempsSecondMembreEnergContact.Arret_du_comptage(); // fin comptage cpu
           // on incrémente le nombre d'élément traité
           nb_elem_deja_calcule++;
         };
        tempsSecondMembreEnergContact.Arret_du_comptage(); // fin comptage cpu
        //
      };
      
    // récupération de toutes les énergies par le cpu 0
    tempsSecondMembreEnergContact.Arret_du_comptage(); // fin comptage cpu
    Algori::Passage_energiesContact();
    tempsSecondMembreEnergContact.Mise_en_route_du_comptage();
#endif

    tempsSecondMembreEnergContact.Mise_en_route_du_comptage();
    // -- on transfert en global les énergies internes
    Transfert_ParaGlob_energies_contact();
    tempsSecondMembreEnergContact.Arret_du_comptage(); // fin comptage cpu

    // affichage éventuelle de la force maxi de contact
    lesCont->Forces_contact_maxi(aff_iteration);
    // idem pour les gap N et tangentiel
    lesCont->Gap_contact_maxi(aff_iteration);
    if (permet_affichage >3) cout << "\n -- fin calcul second mmembre pour le  contact " << flush;


    tempsSecondMembreEnergContact.Arret_du_comptage(); // temps cpu
    // retour indiquant que tout c'est bien passé
    return true ;                           
   };
          
// choix de la matrice de raideur du système linéaire:
//  1) Si le pointeur est non null, on vérifie les tailles, si c'est différent
//     la matrice est reconstruite
//  2) Si le pointeur est null, il y a création d'une matrice
// lesCondLim est utilisé pour modifier la largeur de bande par exemple
// en fonction des conditions limites linéaire
// lescontacts: si différent de NULL indique qu'il y a du contact
// NB: la dimension du tableau dépend du paramétrage lu dans ParaAlgoGlobal
Tableau <Mat_abstraite* >  Algori::Choix_matriciel
            (int nbddl,Tableau <Mat_abstraite* >& tab_matglob,LesMaillages * lesMail
             ,LesReferences* lesRef,const Nb_assemb& nb_casAssemb,LesCondLim* lesCondLim
             , LesContacts* lescontacts )
 {
    // on s'occupe tout d'abord du contact (extérieur à la boucle éventuelle sur
    // les matrices )
    int demicontact=0,totalcontact=0,cumulecontact=0;
    bool largeur_bande_contact_sup_1 = false; // init
    if (lescontacts != NULL)
        lescontacts->Largeur_Bande(demicontact,totalcontact,nb_casAssemb,cumulecontact);

    // tout d'abord le choix de la matrice de stockage
  
    Tableau<Enum_matrice> tab_enu_matrice_sec; // init
    tab_enu_matrice_sec.Init_from_list(pa.Type_matrice_secondaire());
    Tableau<Enum_type_resolution_matri> tab_enu_type_resolution_matri;
    tab_enu_type_resolution_matri.Init_from_list(pa.Type_resolution_secondaire());
    Tableau<Enum_preconditionnement> tab_enu_preconditionnement;
    tab_enu_preconditionnement.Init_from_list(pa.Type_preconditionnement_secondaire());
    Tableau<int> tab_nb_iter_nondirecte_secondaire;
    tab_nb_iter_nondirecte_secondaire.Init_from_list(pa.Nb_iter_nondirecte_secondaire());
    Tableau<double> tab_tolerance_secondaire;
    tab_tolerance_secondaire.Init_from_list(pa.Tolerance_secondaire());
  
    int nb_mat = 1+tab_enu_matrice_sec.Taille(); // nombre de matrices
  
    // dans le cas où le tableau serait non nulle on regarde s'il a la bonne dimension
    if (tab_matglob.Taille() > nb_mat)
      // là il faut supprimer les matrices pointées en trop
      {int taille = tab_matglob.Taille();
       for (int i=nb_mat+1;i<=taille;i++)
         if (tab_matglob(i)!= NULL)
            delete tab_matglob(i);
       // maintenant on peut changer la taille, en moins donc normalement
       // on ne modifie pas les pointeurs déjà existants
       tab_matglob.Change_taille(nb_mat);
      }
    else if (tab_matglob.Taille() < nb_mat)
      // s'il est trop petit il faut augmenter la taille: on met les nouveaux pointeurs à NULL
      tab_matglob.Change_taille(nb_mat,NULL);
    // sinon il a la bonne taille et le traitement sera fait en fonction
    // de l'existence ou non d'une matrice associée au pointeur

    // on va boucler sur le nombre de matrice
    for (int inum_mat = 1; inum_mat <= nb_mat; inum_mat++)
     { // init des valeurs par défaut de la première matrice
       Enum_matrice enumatrice = pa.Type_Matrice(); // init
       Enum_type_resolution_matri enumResolutionMatri = pa.Type_resolution();
       Enum_preconditionnement enumPreconditionnement = pa.Type_preconditionnement();
       // modif au cas où il y a des matrices secondaire
       if (inum_mat > 1)
         {enumatrice = tab_enu_matrice_sec(inum_mat-1);
          if (tab_enu_type_resolution_matri.Taille() >= (inum_mat-1))
             enumResolutionMatri = tab_enu_type_resolution_matri(inum_mat-1);
          if (tab_enu_preconditionnement.Taille() >= (inum_mat-1))
             enumPreconditionnement = tab_enu_preconditionnement(inum_mat-1);
         };
       int  nbddl_MPI = nbddl; // init
       bool cal_parallele_et_rang_non_nul = false;
       #ifdef UTILISATION_MPI
         // cas d'un  calcul //, seule la matrice du CPU 0 est initialisée
         // à la taille nécessaire pour une résolution globale
         // pour les autres CPU on définit une matrice de taille unitaire
         // qui ne sera pas utilisée mais qui sera valide
         if (ParaGlob::Monde()->rank() != 0)
           {nbddl_MPI = 1;
            cal_parallele_et_rang_non_nul = true;
           };
       #endif

       switch (enumatrice)
       { case CARREE : case RECTANGLE : // ici comme les tailles sont identiques rectangle
                                        // et CARRE -> même traitement
          {// cas s'une matrice carrée
           if (tab_matglob(inum_mat) == NULL)
             {tab_matglob(inum_mat) = new Mat_pleine (nbddl_MPI,nbddl_MPI);}
           // sinon on regarde si la taille est bonne
           else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
             // dans ce cas on change de matrice
             { delete tab_matglob(inum_mat);
               tab_matglob(inum_mat) = new Mat_pleine (nbddl_MPI,nbddl_MPI);
             };
           tab_matglob(inum_mat)->Change_Choix_resolution
                      (enumResolutionMatri,enumPreconditionnement);
           // sinon il n'y a rien n'a faire au niveau des tailles
           break;
          }
         case CARREE_LAPACK :
          {// cas s'une matrice carrée a priori quelconque
           if (tab_matglob(inum_mat) == NULL)
             {tab_matglob(inum_mat) = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.
                            ,enumResolutionMatri,enumPreconditionnement);
             }
           // sinon on regarde si la taille est bonne
           else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
             // dans ce cas on change de matrice
             { delete tab_matglob(inum_mat);
               tab_matglob(inum_mat) = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.
                            ,enumResolutionMatri,enumPreconditionnement);
             };
           // sinon il n'y a rien n'a faire au niveau des tailles
           break;
          }
         case CARREE_SYMETRIQUE_LAPACK :
          {// cas s'une matrice carrée symétrique
           bool symetrique = true;
           if (tab_matglob(inum_mat) == NULL)
             {tab_matglob(inum_mat) = new MatLapack (nbddl_MPI,nbddl_MPI,symetrique,0.
                            ,enumResolutionMatri,enumPreconditionnement);
             }
           // sinon on regarde si la taille est bonne
           else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
             // dans ce cas on change de matrice
             { delete tab_matglob(inum_mat);
               tab_matglob(inum_mat) = new MatLapack (nbddl_MPI,nbddl_MPI,symetrique,0.
                            ,enumResolutionMatri,enumPreconditionnement);
             };
           // sinon il n'y a rien n'a faire au niveau des tailles
           break;
          }
         case RECTANGLE_LAPACK :
          {// cas s'une matrice rectangulaire
           if (tab_matglob(inum_mat) == NULL)
             {tab_matglob(inum_mat) = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.
                            ,enumResolutionMatri,enumPreconditionnement);
             }
           // sinon on regarde si la taille est bonne
           else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
             // dans ce cas on change de matrice
             { delete tab_matglob(inum_mat);
               tab_matglob(inum_mat) = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.
                            ,enumResolutionMatri,enumPreconditionnement);
             };
           // sinon il n'y a rien n'a faire au niveau des tailles
           break;
          }
         case BANDE_SYMETRIQUE :
          { //  calcul de la largeur de bande effective
            int demi=0,total=0;
            lesMail->Largeur_Bande(demi,total,nb_casAssemb);
            total = std::min(total,nbddl_MPI);
            demi = std::min(demi,nbddl_MPI);
            // --- on stocke les largeurs initiales de ddl
            { // on redimensionne si nécessaire
              int nombre_cas_assemb = tab_para_stockage.Taille();
              if (nb_casAssemb.n  > nombre_cas_assemb)
                tab_para_stockage.Change_taille(nb_casAssemb.n);
              // maintenant on stocke
              ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
              parstock.nbass = nb_casAssemb;
              parstock.demi = demi;
              parstock.total = total;
            };
            
            // dans le cas non // ou // pour le CPU 0, on précise la dim
            // ici on garde nbddl au lien de nbddl_MPI: ils sont identiques
            if (!cal_parallele_et_rang_non_nul)
              { if (ParaGlob::NiveauImpression() >= 5)
                  cout << "\n raideur: largeur de bande initiale pour les elements (en ddl) " << total << " ";
                if (lescontacts!= NULL)
                 {// calcul de la largeur de bande effective en tenant compte du contact
                  int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
                  // on n'intervient que si la largeur en noeud augmente
//                  if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
                  if (largeur_bande_contact_sup_1)
                    {totalc = std::min(totalc,nbddl);
                     demic = std::min(demic,nbddl);
                     cumule = std::min(cumule,nbddl);
                     if (ParaGlob::NiveauImpression() >= 5)
                         cout << "\n augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << " ";
                     demi += cumule; // c'est donc toute la largeur hors diag que l'on ajoute
                     total += 2 * cumule;
                     // on limite au nombres total de ddl
                     total = std::min(total,nbddl);
                     demi = std::min(demi,nbddl);
                    }
                 };
                // cas de l'influence des conditions linéaires d'entrée
                int demic=0,totalc=0,cumule=0;
                // on n'intervient que si la largeur en noeud augmente
                if (lesCondLim->Largeur_Bande(demic,totalc,nb_casAssemb,lesMail,lesRef,cumule))
                 {totalc = std::min(totalc,nbddl);
                  demic = std::min(demic,nbddl);
                  cumule = std::min(cumule,nbddl);
                  if (ParaGlob::NiveauImpression() >= 5)
                    cout << "\n augmentation de la largeur de bande pour les conditions lineaires (en ddl) " << cumule << " ";
       //  ---- les deux lignes qui suivent sont fausses, car dans le cas du produit de deux matrices bandes la largeur de bande du produit
       // = non pas le max, mais la "somme" des deux largeurs de bandes (parties hors diagonales),
       // d'une manière plus précise cf. J Chevalier 1993, soit une matrice bande avec deux bandes de part et autre de la diagonale, différentes
       // A_ij (a,b) telle que  A_ij = 0 si j<i-a ou j>i+b, donc a est la bande inf (hors diag) et b est la bande sup (hors diag)
       // soient A(a,b) et B(c,d) alors A*B est de type (a+B,b+d)
       // or pour mettre en place les CL (actuellement)
       // on utilise une matrice de rotation et on fait S^T * K * S , tel S est de type (a,b) et donc S^T est de type (b,a) donc
       // si on suppose que K est de type (l,l) et bien S^T * K * S sera de type (l+a+b,l+a+b) donc c'est a+b qui compte (hors diag)
            
       //           demi = MaX(demi,demic);
       //           total = MaX(total,totalc);
                  demi += cumule; // c'est donc toute la largeur hors diag que l'on ajoute
                  total += 2 * cumule;
                  // on limite au nombres total de ddl
                  total = std::min(total,nbddl);
                  demi = std::min(demi,nbddl);
                 };

                if (ParaGlob::NiveauImpression() >= 5)
                  cout << "\n largeur de bande finale (en ddl) " << total << " "
                       << "\n stockage symetrique d'ou utilisation de la demi largeur= " << demi << " ";
                  };
              
            // création éventuelle et initialisation a zero
            if (tab_matglob(inum_mat) == NULL)
              {tab_matglob(inum_mat) = new MatBand (BANDE_SYMETRIQUE,demi,nbddl_MPI,0);
              }
            // sinon on regarde si la taille est bonne
            else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
              // dans ce cas on change de matrice
              { delete tab_matglob(inum_mat);
                tab_matglob(inum_mat) = new MatBand (BANDE_SYMETRIQUE,demi,nbddl_MPI,0);
              }
            else // sinon on se contente d'initialiser à 0.
              tab_matglob(inum_mat)->Initialise(0.);
            // on spécifie le type de résolution
            tab_matglob(inum_mat)->Change_Choix_resolution
                      (enumResolutionMatri,enumPreconditionnement);
            // on vérifie que l'assemblage choisit est bien de type symétrique
            if (pa.Symetrie_matrice() != true)
             { cout << "Erreur de type d'assemblage, le type de matrice choisit est symetrique"
                    << " alors que l'assemblage n'est pas symetrique !! \n";
               cout << " Algori::Choix_matriciel(..)";
               // dans le cas où un comptage du calcul est en cours on l'arrête
               if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
               if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
               if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
               Sortie(1);
             };
          }
          break;
         case BANDE_NON_SYMETRIQUE_LAPACK : case BANDE_SYMETRIQUE_LAPACK: // le choix est fait au moment de l'allocation
          { //  calcul de la largeur de bande effective
            int demi=0,total=0;
            lesMail->Largeur_Bande(demi,total,nb_casAssemb);
            total = std::min(total,2*nbddl_MPI-1);
            demi = std::min(demi,nbddl_MPI);
            // --- on stocke les largeurs initiales de ddl
            { // on redimensionne si nécessaire
              int nombre_cas_assemb = tab_para_stockage.Taille();
              if (nb_casAssemb.n  > nombre_cas_assemb)
                tab_para_stockage.Change_taille(nb_casAssemb.n);
              // maintenant on stocke
              ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
              parstock.nbass = nb_casAssemb;
              parstock.demi = demi;
              parstock.total = total;
            };
            
            // dans le cas non // ou // pour le CPU 0, on précise la dim
            // ici on garde nbddl au lien de nbddl_MPI: ils sont identiques
            if (!cal_parallele_et_rang_non_nul)
               {if (ParaGlob::NiveauImpression() >= 5)
                  cout << "\n raideur: largeur de bande initiale pour les elements (en ddl) " << total << " ";
                if (lescontacts!= NULL)
                  {// calcul de la largeur de bande effective en tenant compte du contact
                   int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
                   // on n'intervient que si la largeur en noeud augmente
//                   if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
                   if (largeur_bande_contact_sup_1)
                    {totalc = std::min(totalc,2*nbddl-1);
                     demic = std::min(demic,nbddl);
                     cumule = std::min(cumule,nbddl);
                     if (ParaGlob::NiveauImpression() >= 5)
                       cout << "\n augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << " ";
                     demi += cumule; // c'est donc toute la largeur hors diag que l'on ajoute
                     total += 2 * cumule;
                     // on limite au nombres total de ddl
                     total = std::min(total,2*nbddl-1);
                     demi = std::min(demi,nbddl);
                    };
                  };
                // cas de l'influence des conditions linéaires d'entrée
                int demic=0,totalc=0,cumule=0;
                // on n'intervient que si la largeur en noeud augmente
                if (lesCondLim->Largeur_Bande(demic,totalc,nb_casAssemb,lesMail,lesRef,cumule))
                 {totalc = std::min(totalc,2*nbddl-1);
                  demic = std::min(demic,nbddl);
                  cumule = std::min(cumule,nbddl);
                  if (ParaGlob::NiveauImpression() >= 5)
                    cout << "\n augmentation de largeur de bande pour les conditions lineaires (en ddl) " << totalc << " ";
       //  ---- les deux lignes qui suivent sont faussent, car dans le cas du produit de deux matrices bandes la largeur de bande du produit
       // = non pas le max, mais la "somme" des deux largeurs de bandes (parties hors diagonales),
       // d'une manière plus précise cf. J Chevalier 1993, soit une matrice bande avec deux bandes de part et autre de la diagonale, différentes
       // A_ij (a,b) telle que  A_ij = 0 si j<i-a ou j>i+b, donc a est la bande inf (hors diag) et b est la bande sup (hors diag)
       // soient A(a,b) et B(c,d) alors A*B est de type (a+c,b+d)
       // or pour mettre en place les CL (actuellement)
       // on utilise une matrice de rotation et on fait S^T * K * S , tel S est de type (a,b) et donc S^T est de type (b,a) donc
       // si on suppose que K est de type (l,l) et bien S^T * K * S sera de type (l+a+b,l+a+b) donc c'est a+b qui compte (hors diag)
            
       //           demi = MaX(demi,demic);
       //           total = MaX(total,totalc);
                  demi += cumule; // c'est donc toute la largeur hors diag que l'on ajoute
                  total += 2 * cumule;
                  // on limite au nombres total de ddl
                  total = std::min(total,2*nbddl-1);
                  demi = std::min(demi,nbddl);
                 };
                if (ParaGlob::NiveauImpression() >= 5)
                  cout << "\n largeur de bande finale (en ddl) " << total << " ";
               };
             
            // création éventuelle et initialisation a zero
            if (tab_matglob(inum_mat) == NULL)
              { if (enumatrice == BANDE_SYMETRIQUE_LAPACK)
                  { tab_matglob(inum_mat) = new MatLapack (enumatrice,demi,nbddl_MPI,true,0.0
                                       ,enumResolutionMatri,enumPreconditionnement);}
                else
                  { tab_matglob(inum_mat) = new MatLapack (enumatrice,total,nbddl_MPI,false,0.0
                                       ,enumResolutionMatri,enumPreconditionnement);};
              }
            // sinon on regarde si la taille est bonne
            else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
              // dans ce cas on change de matrice
              { delete tab_matglob(inum_mat);
                if (enumatrice == BANDE_SYMETRIQUE_LAPACK)
                  { tab_matglob(inum_mat) = new MatLapack (enumatrice,demi,nbddl_MPI,true,0.0
                                       ,enumResolutionMatri,enumPreconditionnement);}
                else
                  { tab_matglob(inum_mat) = new MatLapack (enumatrice,total,nbddl_MPI,false,0.0
                                       ,enumResolutionMatri,enumPreconditionnement);};
              }
            else // sinon on se contente d'initialiser à 0.
              tab_matglob(inum_mat)->Initialise(0.);
            if (enumatrice == BANDE_SYMETRIQUE_LAPACK) // si on a choisit un stockage symétrique
              // on vérifie que l'assemblage choisit est bien de type symétrique
              if (pa.Symetrie_matrice() != true)
               { cout << "Erreur de type d'assemblage, le type de matrice choisit est symetrique"
                      << " alors que l'assemblage n'est pas symetrique !! \n";
                 cout << " Algori::Choix_matriciel(..)";
                 // dans le cas où un comptage du calcul est en cours on l'arrête
                 if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                 if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                 if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                 Sortie(1);
               };
          }
          break;
       case CREUSE_COMPRESSEE_COLONNE :
         { //  cas d'une matrice creuse compressée par colonne
           // tout d'abord on demande l'ensemble des tables de connexion
           // pour les éléments
           Tableau < Tableau <int> > petites_matrices;
           lesMail->Table_connexion(petites_matrices,nb_casAssemb);
           // pour l'instant le cas avec contact n'est pas traité
           if (lescontacts!= NULL)
             { cout << "\n *** cas avec contact non actuellement traite "
                    << "\n Algori::Choix_matriciel(... ) \n";
               // dans le cas où un comptage du calcul est en cours on l'arrête
               if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
               if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
               if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
               Sortie(1);
             };
       // !!! pour l'instant on ne prend pas en compte les conditions linéaires initiales
       // on vérifie 
           if (lesCondLim->ExisteCondiLimite())
             { cout << "\n*** Erreur de choix de matrice de masse, il y a des conditions lineaires intiales, ce qui conduit "
              << "  a intervenir sur le stockage. Or le cas CREUSE_COMPRESSEE_COLONNE ne prend pas en compte pour l'instant "
                      << " la presence de conditions limites  !! \n";
               if (ParaGlob::NiveauImpression() >= 5)
                 cout << "\n Algori::Choix_matriciel(... ) \n";
               // dans le cas où un comptage du calcul est en cours on l'arrête
               if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
               if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
               if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
               Sortie(1);
             };
           // puis on définit la matrice en conséquence
           if (tab_matglob(inum_mat) == NULL)
             {tab_matglob(inum_mat) = new Mat_creuse_CompCol(nbddl_MPI,nbddl_MPI, petites_matrices);
             }
           // sinon on regarde si la taille est bonne
           else if ((tab_matglob(inum_mat)->Nb_ligne()!= nbddl_MPI) || (tab_matglob(inum_mat)->Nb_colonne()!= nbddl_MPI))
             // dans ce cas on change de matrice
             { delete tab_matglob(inum_mat);
               tab_matglob(inum_mat) = new Mat_creuse_CompCol(nbddl_MPI,nbddl_MPI, petites_matrices);
             };
           // on spécifie le type de résolution
           tab_matglob(inum_mat)->Change_Choix_resolution
                      (enumResolutionMatri,enumPreconditionnement);
           // sinon il n'y a rien n'a faire au niveau des tailles
         }
         break;
       default :
         cout << "\nErreur : valeur incorrecte ou non implanté du type  de matrice retenue, type  !\n"
              << " = " << enumatrice ;
         cout << "\n Algori::Choix_matriciel(... ) \n";
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
         if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
         if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
         Sortie(1);
       };
       if (ParaGlob::NiveauImpression() >= 5)
         cout << "\n stockage raideur" << enumatrice << " taille (entiers) :" << tab_matglob(inum_mat)->Place() << endl;

     };
    return tab_matglob;
  };        
    
// Mise à jour éventuelle de la taille de la matrice de raideur ou masse du système linéaire
// en fonction du contact ou en fonction de "nouvelle_largeur_imposee" dans ce cas, lescontacts
// n'est pas utilisé ! donc la méthode peut servir dans ce cas sans le contact
//
// la matrice *matglob: doit déjà existée
// en retour le pointeur désigne toujours la même matrice qu'en entrée
// avec les  même caractéristique de fonctionnement que l'ancienne (mêm résolution, stockage ...)
// NB: la dimension du tableau dépend du paramétrage lu dans ParaAlgoGlobal
// 1) si niveau_substitution =0 : -> mise à jour de toutes les matrices
// 2) si niveau_substitution = i : -> uniquement mise à jour de la matrices i
// par contre la dimension de tab_matglob est toujours mis à jour si c'est nécessaire
// si nouvelle_largeur_imposee : n'est pas nulle, cela veut-dire qu'il y a eu une nouvelle numérotation
//        et une nouvelle mise en place des pointeurs d'assemblage
// alors c'est directement cette valeur de ddl (et non de noeud) qui est imposée sur la matrice
//  nouvelle_largeur_imposee.un = la largeur totale
//  nouvelle_largeur_imposee.deux = la demie largeur
//  nouvelles_largeur_en_ddl.trois = la demie largeur maximale pour la partie éléments finis
//                                   uniquement (sans les CLL)
void Algori::Mise_a_jour_Choix_matriciel_contact
            (Tableau <Mat_abstraite* >& tab_matglob,const Nb_assemb& nb_casAssemb
            , LesContacts* lescontacts
            ,int niveau_substitution ,TroisEntiers* nouvelle_largeur_imposee)
 {
    // on s'occupe tout d'abord du contact (extérieur à la boucle éventuelle sur
    // les matrices )
    int demicontact=0,totalcontact=0,cumulecontact=0;
    bool largeur_bande_contact_sup_1 = false; // init
    if (lescontacts != NULL)
        largeur_bande_contact_sup_1 = lescontacts->Largeur_Bande(demicontact,totalcontact,nb_casAssemb,cumulecontact);
    
  #ifdef UTILISATION_MPI
    int proc_en_cours = ParaGlob::Monde()->rank();
    // la suite ne concerne que le proc 0 qui gère la résolution globale
    if (proc_en_cours == 0)
     {
  #endif
 
    // tout d'abord le choix de la matrice de stockage
    Tableau<Enum_matrice> tab_enu_matrice_sec; // init
    tab_enu_matrice_sec.Init_from_list(pa.Type_matrice_secondaire());
    Tableau<Enum_type_resolution_matri> tab_enu_type_resolution_matri;
    tab_enu_type_resolution_matri.Init_from_list(pa.Type_resolution_secondaire());
    Tableau<Enum_preconditionnement> tab_enu_preconditionnement;
    tab_enu_preconditionnement.Init_from_list(pa.Type_preconditionnement_secondaire());
    Tableau<int> tab_nb_iter_nondirecte_secondaire;
    tab_nb_iter_nondirecte_secondaire.Init_from_list(pa.Nb_iter_nondirecte_secondaire());
    Tableau<double> tab_tolerance_secondaire;
    tab_tolerance_secondaire.Init_from_list(pa.Tolerance_secondaire());
  
    int nb_mat = 1+tab_enu_matrice_sec.Taille(); // nombre de matrices
  
    // dans le cas où le tableau serait non nulle on regarde s'il a la bonne dimension
    if (tab_matglob.Taille() > nb_mat)
      // là il faut supprimer les matrices pointées en trop
      {int taille = tab_matglob.Taille();
       for (int i=nb_mat+1;i<=taille;i++)
         if (tab_matglob(i)!= NULL)
            delete tab_matglob(i);
       // maintenant on peut changer la taille, en moins donc normalement
       // on ne modifie pas les pointeurs déjà existants
       tab_matglob.Change_taille(nb_mat);
      }
    else if (tab_matglob.Taille() < nb_mat)
      // s'il est trop petit il faut augmenter la taille: on met les nouveaux pointeurs à NULL
      tab_matglob.Change_taille(nb_mat,NULL);
    // sinon il a la bonne taille on suppose que les matrices pointées sont de bonnes dimensions
  
    int i_deb = 1; int i_fin = nb_mat+1; // initi par défaut des bornes
    if (niveau_substitution != 0)
      { i_deb = niveau_substitution; i_fin = i_deb+1;};
  
    // on va boucler sur le nombre de matrice
    for (int inum_mat = i_deb; inum_mat < i_fin; inum_mat++)
     { // init des valeurs par défaut de la première matrice
       Mat_abstraite*  matglob = tab_matglob(inum_mat);
       Enum_matrice enumatrice = pa.Type_Matrice(); // init
       Enum_type_resolution_matri enumResolutionMatri = pa.Type_resolution();
       Enum_preconditionnement enumPreconditionnement = pa.Type_preconditionnement();
       // modif au cas où il y a des matrices secondaire
       if (inum_mat > 1)
         {enumatrice = tab_enu_matrice_sec(inum_mat-1);
          if (tab_enu_type_resolution_matri.Taille() > (inum_mat-1))
             enumResolutionMatri = tab_enu_type_resolution_matri(inum_mat-1);
          if (tab_enu_preconditionnement.Taille() > (inum_mat-1))
             enumPreconditionnement = tab_enu_preconditionnement(inum_mat-1);
         };

       switch (enumatrice)
        { case CARREE : break; // pas de modification
          case RECTANGLE : break; // pas de modification
          case CARREE_LAPACK : break; // pas de modification
          case CARREE_SYMETRIQUE_LAPACK : break; // pas de modification
          case RECTANGLE_LAPACK : break; // pas de modification
          case BANDE_SYMETRIQUE : case BANDE_SYMETRIQUE_LAPACK:
            { // -- récup des anciennes tailles
              int nbddl = matglob->Nb_ligne();
              ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
              int ancien_demi = matglob->Nb_colonne();
              // calcul de la largeur de bande effective en tenant compte du contact
              int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
              int demi = parstock.demi; // init
              if (nouvelle_largeur_imposee == NULL) // cas où il faut calculer
               {// on n'intervient que si la largeur en noeud augmente
//                if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
                if (largeur_bande_contact_sup_1)
                  {// la mise à jour dépend du type de contact
                   if (pa.ContactType() == 1) // contact cinématique, sans multiplicateur ni pénalisation
                     // on doit augmenter la lb en fonction de cumule, si la largeur en noeud est différente
                     // de 1
                     {cumule = std::min(cumule,nbddl);
                       demi = parstock.demi + cumule; // c'est donc toute la largeur hors diag que l'on ajoute
                       if (demi > ancien_demi)
                         if (ParaGlob::NiveauImpression() >= 5)
                              cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << "old_1/2lb: "<<ancien_demi << " new_1/2lb: "<< demi <<" ";
                     }
                   else if ((pa.ContactType() == 2)|| (pa.ContactType() == 3)
                            ||(pa.ContactType() == 4) )// cas d'un contact par pénalisation ou type 3 *** exploratoire
                     // dans ce cas on considère uniquement le LB maxi
                     { demic = std::min(demic,nbddl);
                       int demi_inter = std::max(parstock.demi, demic);
                       if (demi_inter > ancien_demi)
                         if ((ParaGlob::NiveauImpression() >= 5) && (demic > demi))
                         cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << (demic-demi)  << "old_1/2lb: "<<ancien_demi << " new_1/2lb: "<< demi_inter <<" ";
                       demi = demi_inter;
                     };
                   };
                }
              else // sinon on applique la demande
                { demi = nouvelle_largeur_imposee->deux;
                  // et  on met à jour : uniquement la partie relative aux maillages
                  ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
                  parstock.nbass = nb_casAssemb;
                  parstock.demi = nouvelle_largeur_imposee->trois;
                  parstock.total = 2 * nouvelle_largeur_imposee->trois -1;
                };
              // on limite au nombres total de ddl
              demi = std::min(demi,nbddl);

              // on change la taille, si elle doit augmenter, ou si elle est inférieur d'au moins 5% par rapport à la taille actuelle
              if ( (demi > matglob->Nb_colonne()) ||
                  (demi < 0.95 * matglob->Nb_colonne()))
                   // on change la taille
                { if (ParaGlob::NiveauImpression() >= 5)
                    cout << "\n augmentation effective: LB/2+1 finale(en ddl) " << demi  << " ";
                  switch (enumatrice)
                   { case BANDE_SYMETRIQUE:
                       ((MatBand*)matglob)->Change_taille(demi,nbddl );
                       break;
                     case BANDE_SYMETRIQUE_LAPACK:
                       ((MatLapack*)matglob)->Change_taille(nbddl,demi);
                       break;
                     default:
                       cout << "\n erreur non prevue Algori::Mise_a_jour_Choix_matriciel_contact "<<endl;
                       // dans le cas où un comptage du calcul est en cours on l'arrête
                       if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                       if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                       if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                       Sortie(1);
                   };
                };
              break;
            }
          case BANDE_NON_SYMETRIQUE : case BANDE_NON_SYMETRIQUE_LAPACK : 
            { // -- récup des anciennes tailles
              int nbddl = matglob->Nb_ligne();
              int ancien_lb = matglob->Nb_colonne();
              ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
                  //int demi = matglob->Nb_colonne();
              // calcul de la largeur de bande effective en tenant compte du contact
              int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
              int total = parstock.total; // init
              if (nouvelle_largeur_imposee == NULL) // cas où il faut calculer
               {// on n'intervient que si la largeur en noeud augmente
//                if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
                if (largeur_bande_contact_sup_1)
                  {// la mise à jour dépend du type de contact
                   if (pa.ContactType() == 1) // contact cinématique, sans multiplicateur ni pénalisation
                    // on doit augmenter la lb en fonction de cumule
                    { cumule = std::min(cumule,nbddl);
                      int total_inter = parstock.total + 2 * cumule;
                      if (total_inter > ancien_lb)
                        if (ParaGlob::NiveauImpression() >= 5)
                             cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << "old_lb: "<<ancien_lb << " new_lb: "<< total_inter <<" ";
                      total = total_inter;
                    }
                   
                   else if ((pa.ContactType() == 2)|| (pa.ContactType() == 3)
                            || (pa.ContactType() == 4))// cas d'un contact par pénalisation ou type 3 *** exploratoire
                    // dans ce cas on considère uniquement le LB maxi
                    { totalc = std::min(totalc,2*nbddl-1);
                      int total_inter = std::max(parstock.total, totalc);
                      if (total_inter > ancien_lb)
                         if ((ParaGlob::NiveauImpression() >= 5) && (totalc > total))
                             cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << (totalc-total)  << "old_lb: "<<ancien_lb << " new_lb: "<< total_inter << " ";
                      total = total_inter;
                   };
                  };
                 }
              else // sinon on applique la demande
                { total = nouvelle_largeur_imposee->un;
                  // et  on met à jour : uniquement la partie relative aux maillages
                  ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
                  parstock.nbass = nb_casAssemb;
                  parstock.demi = nouvelle_largeur_imposee->trois;
                  parstock.total = 2 * nouvelle_largeur_imposee->trois -1;
                };
             // on limite au nombres total de ddl
              total = std::min(total,2*nbddl-1);

              // on change la taille, si elle doit augmenter, ou si elle est inférieur d'au moins 5% par rapport à la taille actuelle
              if ( (total > matglob->Nb_colonne()) ||
                  (total < 0.95 * matglob->Nb_colonne()))
                       // on change la taille
                { if (ParaGlob::NiveauImpression() >= 5)
                            cout << "\n augmentation effective: LB finale(en ddl) " << total  << " ";
                  switch (enumatrice)
                   { case BANDE_NON_SYMETRIQUE:
                       ((MatBand*)matglob)->Change_taille(total,nbddl );
                       break;
                     case BANDE_NON_SYMETRIQUE_LAPACK:
                       ((MatLapack*)matglob)->Change_taille(nbddl,total);
                       break;
                    
                     default:
                       cout << "\n erreur2 non prevue Algori::Mise_a_jour_Choix_matriciel_contact "<<endl;
                       // dans le cas où un comptage du calcul est en cours on l'arrête
                       if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                       if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                       if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                       Sortie(1);
                   };
                };	
              break;
            }
          case CREUSE_COMPRESSEE_COLONNE :
            { cout << "\n *** cas avec contact non actuellement traite "
                   << "\n Algori::Mise_a_jour_Choix_matriciel_contact(... ) \n";
              // dans le cas où un comptage du calcul est en cours on l'arrête
              if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
              if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
              if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
              Sortie(1);
              break;
            }
          default :
             cout << "\nErreur : valeur incorrecte ou non implante du type  de matrice retenue, type  !\n"
                  << " = " << enumatrice ;
             cout << "\n Algori::Mise_a_jour_Choix_matriciel_contact(... ) \n";
             // dans le cas où un comptage du calcul est en cours on l'arrête
             if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
             if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
             if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
             Sortie(1);
        };
       if (ParaGlob::NiveauImpression() >= 5)
          cout  << "\n stockage raideur " << enumatrice << " taille (entiers) :"
                << matglob->Place() << endl;
     };
  #ifdef UTILISATION_MPI
     }; // fin du cas proc 0
  #endif

     
  };
     
// choix de la matrice de la matrice de masse globale
// ou mise à jour de la matrice masse si elle existe
Mat_abstraite*  Algori::Choix_matrice_masse
            (int nbddl,Mat_abstraite* mat_mass,LesMaillages * lesMail,LesReferences* lesRef
             ,const Nb_assemb& nb_casAssemb,LesContacts* ,LesCondLim* lesCondLim)
{
  int  nbddl_MPI = nbddl; // init
  bool cal_parallele_et_rang_non_nul = false;
  #ifdef UTILISATION_MPI
    // cas d'un  calcul //, seule la matrice du CPU 0 est initialisée
    // à la taille nécessaire pour une résolution globale
    // pour les autres CPU on définit une matrice de taille unitaire
    // qui ne sera pas utilisée mais qui sera valide
    if (ParaGlob::Monde()->rank() != 0)
      {nbddl_MPI = 1;
       cal_parallele_et_rang_non_nul = true;
      };
  #endif

  // tout d'abord le choix de la matrice de stockage
  // s'il n'y a pas de condition linéaires classiques, la largeur de bande ne doit pas changer
  if ( ((pa.Type_calcul_masse() == MASSE_DIAG_COEF_EGAUX) || (pa.Type_calcul_masse() == MASSE_DIAG_COEF_VAR))
        && (!(lesCondLim->ExisteCondiLimite()))
//        && (pa.ContactType() != 1) // non ce sera géré dans la mise à jour, pour éviter
//                                   // un stockage en bande systématique
     )
    { // on utilise une matrice diagonale
      if (mat_mass == NULL)
        { mat_mass = new MatDiag (nbddl_MPI);
          if ((ParaGlob::NiveauImpression() > 2)&&(!cal_parallele_et_rang_non_nul))
            cout << "\n stockage initial matrice masse --> diagonal, nbddl= " << nbddl_MPI ;
        }
      else // sinon on redimentionne ou on change si nécessaire
       { if (mat_mass->Type_matrice() !=  DIAGONALE)
              { delete mat_mass; mat_mass = new MatDiag (nbddl_MPI);} // cas où on change de matrice
         else { MatDiag * mat_mass_diag = (MatDiag*) mat_mass; // cas où on redimentionne la matrice
                mat_mass_diag->Change_taille(nbddl_MPI);
              };
         if ((ParaGlob::NiveauImpression() > 2)&&(!cal_parallele_et_rang_non_nul))
            cout << "\n changement stockage matrice masse --> diagonal, nbddl= " << nbddl_MPI ;
       };
      // --- on stocke les largeurs initiales de ddl
      // on redimensionne si nécessaire
      int nombre_cas_assemb = tab_para_stockage.Taille();
      if (nb_casAssemb.n  > nombre_cas_assemb) 
           tab_para_stockage.Change_taille(nb_casAssemb.n);
      // maintenant on stocke
      ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
      parstock.nbass = nb_casAssemb;
      parstock.demi = parstock.total = 1; // tout est diagonal
    }
  else if ( ((pa.Type_calcul_masse() == MASSE_DIAG_COEF_EGAUX) ||
             (pa.Type_calcul_masse() == MASSE_DIAG_COEF_VAR))
            &&
             ((lesCondLim->ExisteCondiLimite())
              || (pa.ContactType() == 1) // on part avec l'idée d'une largeur de 1 en noeud
               // ce qui est le cas d'un contact solide - déformable
               // sera ensuite éventuellement augmenté
              // || contact_type_1_et_maitre_deformable
              )
//        && (pa.Type_Matrice() == DIAGONALE)
          )
    { // dans le cas ou on devrait utiliser une matrice diagonale, on va arbitrairement utiliser une matrice bande
      // mais avec une largeur uniquement déterminé par les conditions linéaires
      // en particulier, si les CLL => un stockage diagonale en noeuds, on utilise une matrice de bande = dim espace
      //  calcul de la largeur de bande effective due à  l'influence des conditions linéaires d'entrée
      // --- on stocke les largeurs initiales de ddl
      // on redimensionne si nécessaire
      int nombre_cas_assemb = tab_para_stockage.Taille();
      if (nb_casAssemb.n  > nombre_cas_assemb) 
           tab_para_stockage.Change_taille(nb_casAssemb.n);
      // maintenant on stocke
      ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
      parstock.nbass = nb_casAssemb;
      parstock.demi = parstock.total = 1; // tout est diagonal
     
      int demic=0,totalc=0,cumule=0;
      int demi=0,total=0;
      // on n'intervient que si la largeur en noeud augmente
      // et que l'on est sur le CPU 0
      if (   (lesCondLim->Largeur_Bande(demic,totalc,nb_casAssemb,lesMail,lesRef,cumule))
          && (!cal_parallele_et_rang_non_nul)
         )
        { demi = 1 + cumule; // c'est donc toute la largeur hors diag que l'on ajoute
          total = 1 + 2 * cumule;
        }
      else // sinon on est dans le cas où on peut utiliser une matrice bande = dim espace (ex: si dim 3: tridiagonale)
        { int dim_espace = ParaGlob::Dimension();
          demi = dim_espace;
          total = 1 + 2 * (dim_espace-1);
        };
      // on limite au nombres total de ddl
      total = std::min(total,nbddl_MPI);
      demi = std::min(demi,nbddl_MPI);
      // si la matrice actuelle est de type diagonale il faut la supprimer et définir une
      // matrice bande
      bool change_stockage = false;
      if (mat_mass != NULL)
       if (mat_mass->Type_matrice() ==  DIAGONALE)
        { delete mat_mass; mat_mass = NULL;change_stockage = true;}; // cas où on change de stockage
     
      // maintenant on crée une matrice en fonction de la demande
      switch (pa.Type_Matrice())
        { case CARREE :
          {// cas d'une matrice carrée, là a priori on n'a pas besoin d'augmenter la largeur, on est au maxi
           if (mat_mass == NULL)
             { mat_mass = new Mat_pleine (nbddl_MPI,nbddl_MPI);}
           else // sinon on redimentionne ou on change si nécessaire
            { Mat_pleine * mat_mass_d = (Mat_pleine*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Initialise (nbddl_MPI,nbddl_MPI,0.);
            };
           break;
          }
         case RECTANGLE :
          {// cas s'une matrice rectangulaire, pour l'instant idem que le cas précédent
           // car pour l'instant on ne voit pas très bien à quoi cela servirait d'être différent
           if (mat_mass == NULL)
             { mat_mass = new Mat_pleine (nbddl_MPI,nbddl_MPI);}
           else // sinon on redimentionne ou on change si nécessaire
            { Mat_pleine * mat_mass_d = (Mat_pleine*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Initialise (nbddl_MPI,nbddl_MPI,0.);
            };
           break;
          }
         case CARREE_LAPACK :
          {// cas s'une matrice carrée
           if (mat_mass == NULL)
             { mat_mass = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
           else // sinon on redimentionne ou on change si nécessaire
            { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI,0.);
            };
           break;
          }
         case CARREE_SYMETRIQUE_LAPACK :
          {// cas s'une matrice carrée symétrique
           if (mat_mass == NULL)
             { mat_mass = new MatLapack (nbddl_MPI,nbddl_MPI,true,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
           else // sinon on redimentionne ou on change si nécessaire
            { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI,0.);
            };
           break;
          }
         case RECTANGLE_LAPACK :
          {// cas s'une matrice rectangulaire, pour l'instant idem que le cas précédent
           if (mat_mass == NULL)
             { mat_mass = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
           else // sinon on redimentionne ou on change si nécessaire
            { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI,0.);
            };
           break;
          }
         case BANDE_SYMETRIQUE :
          {// initialisation a zero
           if (mat_mass == NULL)
             { mat_mass = new MatBand (BANDE_SYMETRIQUE,demi,nbddl_MPI,0);}
           else // sinon on redimentionne ou on change si nécessaire
            { MatBand * mat_mass_d = (MatBand*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Change_taille(demi,nbddl_MPI);
              mat_mass_d->Initialise (0.);
            };
           break;
          }
         case BANDE_NON_SYMETRIQUE_LAPACK : case BANDE_SYMETRIQUE_LAPACK: // le choix est fait au moment de l'allocation
          {// initialisation a zero
           if (pa.Type_Matrice() == BANDE_SYMETRIQUE_LAPACK)
            { if (mat_mass == NULL)
               { mat_mass = new MatLapack(pa.Type_Matrice(),demi,nbddl_MPI,true,0.0
                                    ,pa.Type_resolution(),pa.Type_preconditionnement());}
              else // sinon on redimentionne  si nécessaire
               { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                 mat_mass_d->Change_taille(demi,nbddl_MPI,0.0);
               };
            }
           else 
            { if (mat_mass == NULL)
               { mat_mass = new MatLapack(pa.Type_Matrice(),total,nbddl_MPI,false,0.0
                                   ,pa.Type_resolution(),pa.Type_preconditionnement());}
              else // sinon on redimentionne  si nécessaire
               { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                       mat_mass_d->Change_taille(total,nbddl_MPI,0.0);
               };
            };
           break;
          }
         case CREUSE_COMPRESSEE_COLONNE :
          { //  cas d'une matrice creuse compressée par colonne
            // tout d'abord on demande l'ensemble des tables de connexion
            // pour les éléments
            Tableau < Tableau <int> > petites_matrices;
            lesMail->Table_connexion(petites_matrices,nb_casAssemb);
            // puis on définit la matrice en conséquence
           
            if (mat_mass == NULL)
             { mat_mass = new Mat_creuse_CompCol(nbddl_MPI,nbddl_MPI, petites_matrices);}
            else // sinon on redimentionne ou on change si nécessaire
             { if (mat_mass->Type_matrice() !=  CREUSE_COMPRESSEE_COLONNE)
                    { delete mat_mass; mat_mass = new Mat_creuse_CompCol(nbddl_MPI,nbddl_MPI, petites_matrices);} // cas où on change de matrice
               else {Mat_creuse_CompCol * mat_mass_d = (Mat_creuse_CompCol*) mat_mass; // cas où on redimentionne la matrice
                     mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI, petites_matrices);
                    };
             };
            // pour l'instant on ne prend pas en compte les conditions linéaires initiales qui augmente la largeur
            // diagonale par bloc de dim
            // par contre on considère que par défaut la diagonale par bloc existe !!
            if (change_stockage)
             { cout << "\n*** Erreur de choix de matrice de masse, il y a des conditions lineaires intiales,"
                    << " qui vont changer le stockage a l'execution. Or le cas CREUSE_COMPRESSEE_COLONNE ne "
                    << " prend pas en compte pour l'instant ce cas de changement   !! \n";
               if (ParaGlob::NiveauImpression() >= 5)
                 cout << " Algori::Choix_matrice_masse(..)";
               // dans le cas où un comptage du calcul est en cours on l'arrête
               if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
               if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
               if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
               Sortie(1);
             };
            break;
          }
         default :
            cout << "\nErreur : valeur incorrecte ou non implante du type  de matrice retenue, type  !\n"
                 << " = " << pa.Type_Matrice() ;
            cout << "\n Algori::Algori::Choix_matrice_masse(... ) \n";
            // dans le cas où un comptage du calcul est en cours on l'arrête
            if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
            if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
            if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
            Sortie(1);
        };
      if ((ParaGlob::NiveauImpression() > 2)&&(!cal_parallele_et_rang_non_nul))
         { cout << "\n"; if (change_stockage) cout << " changement de ";
           cout << " stockage matrice masse -->" << Nom_matrice(mat_mass->Type_matrice())
                << " nb_ligne= "<< mat_mass->Nb_ligne() <<" nb_col= "<< mat_mass->Nb_colonne() << endl;
         };
    }
  else
    // sinon on se trouve dans le cas de matrice voulue non diagonale car dont le calcul des coefs ne conduit pas
    // à une matrice a priori diagonale (exemple: matrice consistante)
    {switch (pa.Type_Matrice())
     { case CARREE :
       {// cas d'une matrice carrée, là a priori on n'a pas besoin d'augmenter la largeur, on est au maxi
        if (mat_mass == NULL)
          { mat_mass = new Mat_pleine (nbddl_MPI,nbddl_MPI);}
        else // sinon on redimentionne ou on change si nécessaire
         { Mat_pleine * mat_mass_d = (Mat_pleine*) mat_mass; // cas où on redimentionne la matrice
           mat_mass_d->Initialise (nbddl_MPI,nbddl_MPI,0.);
         };
        break;
       }
      case RECTANGLE :
       {// cas s'une matrice rectangulaire, pour l'instant idem que le cas précédent
        // car pour l'instant on ne voit pas très bien à quoi cela servirait d'être différent
        if (mat_mass == NULL)
          { mat_mass = new Mat_pleine (nbddl_MPI,nbddl_MPI);}
        else // sinon on redimentionne ou on change si nécessaire
         { Mat_pleine * mat_mass_d = (Mat_pleine*) mat_mass; // cas où on redimentionne la matrice
           mat_mass_d->Initialise (nbddl_MPI,nbddl_MPI,0.);
         };
        break;
       }
      case CARREE_LAPACK :
       {// cas s'une matrice carrée
        if (mat_mass == NULL)
          { mat_mass = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
        else // sinon on redimentionne ou on change si nécessaire
         { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
           mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI,0.);
         };
        break;
       }
      case CARREE_SYMETRIQUE_LAPACK :
       {// cas s'une matrice carrée symétrique
        if (mat_mass == NULL)
          { mat_mass = new MatLapack (nbddl_MPI,nbddl_MPI,true,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
        else // sinon on redimentionne ou on change si nécessaire
         { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
           mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI,0.);
         };
        break;
       }
      case RECTANGLE_LAPACK :
       {// cas s'une matrice rectangulaire, pour l'instant idem que le cas précédent
        if (mat_mass == NULL)
          { mat_mass = new MatLapack (nbddl_MPI,nbddl_MPI,false,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
        else // sinon on redimentionne ou on change si nécessaire
         { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
           mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI,0.);
         };
        break;
       }
      case BANDE_SYMETRIQUE :
       {//  calcul de la largeur de bande effective
        int demi=0,total=0;
        lesMail->Largeur_Bande(demi,total,nb_casAssemb);
        total = std::min(total,nbddl_MPI);
        demi = std::min(demi,nbddl_MPI);
        // --- on stocke les largeurs initiales de ddl
        // on redimensionne si nécessaire
        int nombre_cas_assemb = tab_para_stockage.Taille();
        if (nb_casAssemb.n  > nombre_cas_assemb) 
           tab_para_stockage.Change_taille(nb_casAssemb.n);
        // maintenant on stocke
        ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
        parstock.nbass = nb_casAssemb;
        parstock.demi = demi;
        parstock.total = total;
        // l'influence des CLL n'intervient que sur le CPU 0
        if (!cal_parallele_et_rang_non_nul)
          {if (ParaGlob::NiveauImpression() >= 5)
               cout << "\n matrice masse: largeur de bande pour les elements (en ddl) " << total << " ";
           // cas de l'influence des conditions linéaires d'entrée
           int demic=0,totalc=0,cumule=0;
           // on n'intervient que si la largeur en noeud augmente
           if (lesCondLim->Largeur_Bande(demic,totalc,nb_casAssemb,lesMail,lesRef,cumule))
             {totalc = std::min(totalc,nbddl_MPI);
              demic = std::min(demic,nbddl_MPI);
              cumule = std::min(cumule,nbddl_MPI);
              if (ParaGlob::NiveauImpression() >= 5)
                   cout << "\n matrice masse: largeur de bande pour les conditions lineaires (en ddl) " << totalc << " ";
              demi += cumule; // c'est donc toute la largeur hors diag que l'on ajoute
              total += 2 * cumule;
              // on limite au nombres total de ddl
              total = std::min(total,nbddl_MPI);
              demi = std::min(demi,nbddl_MPI);
             };

           if (ParaGlob::NiveauImpression() >= 5)
                cout << "\n matrice masse: largeur de bande finale (en ddl) " << total << " "
                     << "\n stockage symetrique d'ou utilisation de la demi largeur= " << demi << " ";
          };
        // initialisation a zero
        if (mat_mass == NULL)
          { mat_mass = new MatBand (BANDE_SYMETRIQUE,demi,nbddl_MPI,0);}
        else // sinon on redimentionne ou on change si nécessaire
         { MatBand * mat_mass_d = (MatBand*) mat_mass; // cas où on redimentionne la matrice
           mat_mass_d->Change_taille(demi,nbddl_MPI);
           mat_mass_d->Initialise (0.);
         };
//        // on vérifie que l'assemblage choisit est bien de type symétrique
//        if (pa.Symetrie_matrice() != true)
//         { cout << "Erreur de type d'assemblage, le type de matrice choisit est non symetrique"
//                << " alors que l'assemblage n'est pas symetrique !! \n";
//           cout << " Algori::Choix_matrice_masse(..)";
//           Sortie(1);
//         };
        break;
       }
      case BANDE_NON_SYMETRIQUE_LAPACK : case BANDE_SYMETRIQUE_LAPACK: // le choix est fait au moment de l'allocation
       {//  calcul de la largeur de bande effective
        int demi=0,total=0;
        lesMail->Largeur_Bande(demi,total,nb_casAssemb);
        total = std::min(total,2*nbddl_MPI-1);
        demi = std::min(demi,nbddl_MPI);
        // --- on stocke les largeurs initiales de ddl
        // on redimensionne si nécessaire
        int nombre_cas_assemb = tab_para_stockage.Taille();
        if (nb_casAssemb.n  > nombre_cas_assemb) 
           tab_para_stockage.Change_taille(nb_casAssemb.n);
        // maintenant on stocke
        ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
        parstock.nbass = nb_casAssemb;
        parstock.demi = demi;
        parstock.total = total;
        // cas de l'influence des conditions linéaires d'entrée  
        int demic=0,totalc=0,cumule=0;
        // on n'intervient que si la largeur en noeud augmente
        // et que l'on est sur le CPU 0
        if (   (lesCondLim->Largeur_Bande(demic,totalc,nb_casAssemb,lesMail,lesRef,cumule))
            && (!cal_parallele_et_rang_non_nul)
           )
          {totalc = std::min(totalc,2*nbddl_MPI-1);
           demic = std::min(demic,nbddl_MPI);
           cumule = std::min(cumule,nbddl_MPI);
           demi += cumule; // c'est donc toute la largeur hors diag que l'on ajoute
           total += 2 * cumule;
           if (ParaGlob::NiveauImpression() >= 5)
             cout << "\n matrice masse: augmentation de la largeur de bande pour les conditions lineaires (en ddl) " << cumule   << " ";
           // on limite au nombres total de ddl
           total = std::min(total,2*nbddl_MPI-1);
           demi = std::min(demi,nbddl_MPI);
          };
        if ((ParaGlob::NiveauImpression() >= 5)&& (!cal_parallele_et_rang_non_nul))
            cout << "\n matrice masse: largeur de bande finale (en ddl) " << total << " ";
        // initialisation a zero
        if (pa.Type_Matrice() == BANDE_SYMETRIQUE_LAPACK)
         { if (mat_mass == NULL)
            { mat_mass = new MatLapack(pa.Type_Matrice(),demi,nbddl_MPI,true,0.0
                                 ,pa.Type_resolution(),pa.Type_preconditionnement());}
           else // sinon on redimentionne  si nécessaire
            { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
              mat_mass_d->Change_taille(demi,nbddl_MPI,0.0);
            };
         }
        else 
         { if (mat_mass == NULL)
            { mat_mass = new MatLapack(pa.Type_Matrice(),total,nbddl_MPI,false,0.0
                                ,pa.Type_resolution(),pa.Type_preconditionnement());}
           else // sinon on redimentionne  si nécessaire
            { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                    mat_mass_d->Change_taille(total,nbddl_MPI,0.0);
            };
         };
        break;
       }
      case CREUSE_COMPRESSEE_COLONNE :
       { //  cas d'une matrice creuse compressée par colonne
         // tout d'abord on demande l'ensemble des tables de connexion
         // pour les éléments
         Tableau < Tableau <int> > petites_matrices;
         lesMail->Table_connexion(petites_matrices,nb_casAssemb);
         // puis on définit la matrice en conséquence
        
         if (mat_mass == NULL)
          { mat_mass = new Mat_creuse_CompCol(nbddl_MPI,nbddl_MPI, petites_matrices);}
         else // sinon on redimentionne ou on change si nécessaire
         { if (mat_mass->Type_matrice() !=  CREUSE_COMPRESSEE_COLONNE)
                { delete mat_mass; mat_mass = new Mat_creuse_CompCol(nbddl_MPI,nbddl_MPI, petites_matrices);} // cas où on change de matrice
           else {Mat_creuse_CompCol * mat_mass_d = (Mat_creuse_CompCol*) mat_mass; // cas où on redimentionne la matrice
                 mat_mass_d->Change_taille(nbddl_MPI,nbddl_MPI, petites_matrices);
                };
         };
      // !!! pour l'instant on ne prend pas en compte les conditions linéaires initiales
      // on vérifie 
         if (lesCondLim->ExisteCondiLimite())
          { cout << "\n*** Erreur de choix de matrice de masse, il y a des conditions lineaires intiales, ce qui conduit "
                 << "  a intervenir sur le stockage. Or le cas CREUSE_COMPRESSEE_COLONNE ne prend pas en compte pour l'instant "
                 << " la presence de conditions limites  !! \n";
            if (ParaGlob::NiveauImpression() >= 5)
              cout << " Algori::Choix_matrice_masse(..)";
            // dans le cas où un comptage du calcul est en cours on l'arrête
            if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
            if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
            if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
            Sortie(1);
          };
         break;
       }
      default :
         cout << "\nErreur : valeur incorrecte ou non implante du type  de matrice retenue, type  !\n"
              << " = " << pa.Type_Matrice() ;
         cout << "\n Algori::Algori::Choix_matrice_masse(... ) \n";
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
         if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
         if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
         Sortie(1);
     };
    };
	 if ((ParaGlob::NiveauImpression() >= 5)&& (!cal_parallele_et_rang_non_nul))
	    cout << "\n Taille du stockage masse " << mat_mass->Type_matrice() << " (taille en entiers) :" 
	         << mat_mass->Place() << endl;
  return mat_mass;
};        

// mise à jour éventuelle du type et/ou de la taille de la matrice de masse
// retourne un pointeur sur la nouvelle matrice, s'il y a eu une modification
//             l'ancienne est supprimée
// sinon retour d'un pointeur NULL
Mat_abstraite*  Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite
        (int nbddl,Mat_abstraite* mat_mass,LesMaillages * lesMail,LesReferences* lesRef
         ,const Nb_assemb& nb_casAssemb, LesContacts* lescontacts)
{

//// préparation pour le cas particulier : pa.ContactType() == 1
//  int contact_type_1_et_maitre_deformable = false;
//  if (pa.ContactType() == 1)
//   {// on veut savoir s'il y a des contacts autres que solide-déformables
//    if (lescontacts->Maitres_avec_deplacement_imposer())
//      // cas où tous les surfaces maîtres sont fixées
//      {contact_type_1_et_maitre_deformable = true;};
//      // dans ce cas on peut utiliser une tridiagonale !
//   };
  bool change_stockage = false; // init pour le retour

  if (( pa.ContactType() == 2) || (pa.ContactType() == 0)
      || (pa.ContactType() == 3)|| (pa.ContactType() == 4))
      
     // cas sans contact ou d'une pénalisation en explicite
     // il n'y a rien n'a faire
     {return NULL;}
  else // sinon on statut en fonction du type existant de matrice
   { // .. on va utiliser des infos de contact, on va en préliminaire
     // .. récupérer ces infos
     int demicontact=0,totalcontact=0,cumulecontact=0;
     bool largeur_bande_contact_sup_1 = lescontacts->Largeur_Bande(demicontact,totalcontact,nb_casAssemb,cumulecontact);
  #ifdef UTILISATION_MPI
    int proc_en_cours = ParaGlob::Monde()->rank();
    // la suite ne concerne que le proc 0 qui gère la résolution globale
    if (proc_en_cours != 0)
      {return NULL;} // retour par défaut
    // pour le proc 0 on continue
  #endif
   
     // en fait il n'y a que deux cas très différents:
     // soit on est en diagonal, et on doit peut-être passer dans un autre type
     // de stockage
     // soit ce n'est pas du diagonal, et on regarde s'il faut agrandir
     if (mat_mass->Type_matrice() == DIAGONALE)
      {// si on a un contact de
       // type 1 c-a-d qui induit des conditions linéaires, on va analyser
       if (pa.ContactType() == 1)
        { // dans le cas ou on devrait utiliser une matrice diagonale, on va utiliser la matrice demandée
          // mais avec une largeur uniquement déterminé par les conditions linéaires éventuelles de contact
          // en particulier, si les CLL => un stockage diagonale en noeuds, on utilise une matrice de bande = dim espace
          //  calcul de la largeur de bande effective due à  l'influence des conditions linéaires de contact
          // --- on stocke les largeurs initiales de ddl
          // on redimensionne si nécessaire
          int nombre_cas_assemb = tab_para_stockage.Taille();
          if (nb_casAssemb.n  > nombre_cas_assemb)
              tab_para_stockage.Change_taille(nb_casAssemb.n);
          // maintenant on stocke
          ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
          parstock.nbass = nb_casAssemb;
          parstock.demi = parstock.total = 1; // tout est diagonal
        
          int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
          int demi=0,total=0;
//          int ancien_lb = mat_mass->Nb_colonne();
       
          //-- dans tous les cas, à minima il faut passer en tridiagonale
//          if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
          if (largeur_bande_contact_sup_1)
             { demi = 1 + cumule; // c'est donc toute la largeur hors diag que l'on ajoute
               total = 1 + 2 * cumule;
             }
          else // sinon on est dans le cas où on peut utiliser une matrice bande = dim espace (ex: si dim 3: tridiagonale)
             { int dim_espace = ParaGlob::Dimension();
               demi = dim_espace;
               total = 1 + 2 * (dim_espace-1);
             };
          // on limite au nombres total de ddl
          total = std::min(total,nbddl);
          demi = std::min(demi,nbddl);
          if (ParaGlob::NiveauImpression() >= 2)
                 cout << "\n augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << " new_lb: "<< total <<" ";
          // comme la matrice actuelle est de type diagonale il faut la supprimer et définir une autre
          if (mat_mass != NULL)
           if (mat_mass->Type_matrice() ==  DIAGONALE)
            { delete mat_mass; mat_mass = NULL;change_stockage = true;}; // cas où on change de stockage
     
          // maintenant on crée une matrice en fonction de la demande
          switch (pa.Type_Matrice())
            { case CARREE :
              {// cas d'une matrice carrée, là a priori on n'a pas besoin d'augmenter la largeur, on est au maxi
               if (mat_mass == NULL)
                 { mat_mass = new Mat_pleine (nbddl,nbddl);}
               else // sinon on redimentionne ou on change si nécessaire
                { Mat_pleine * mat_mass_d = (Mat_pleine*) mat_mass; // cas où on redimentionne la matrice
                  mat_mass_d->Initialise (nbddl,nbddl,0.);
                };
               break;
              }
             case RECTANGLE :
              {// cas s'une matrice rectangulaire, pour l'instant idem que le cas précédent
               // car pour l'instant on ne voit pas très bien à quoi cela servirait d'être différent
               if (mat_mass == NULL)
                 { mat_mass = new Mat_pleine (nbddl,nbddl);}
               else // sinon on redimentionne ou on change si nécessaire
                { Mat_pleine * mat_mass_d = (Mat_pleine*) mat_mass; // cas où on redimentionne la matrice
                  mat_mass_d->Initialise (nbddl,nbddl,0.);
                };
               break;
              }
             case CARREE_LAPACK :
              {// cas s'une matrice carrée
               if (mat_mass == NULL)
                 { mat_mass = new MatLapack (nbddl,nbddl,false,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
               else // sinon on redimentionne ou on change si nécessaire
                { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                  mat_mass_d->Change_taille(nbddl,nbddl,0.);
                };
               break;
              }
             case CARREE_SYMETRIQUE_LAPACK :
              {// cas s'une matrice carrée symétrique
               if (mat_mass == NULL)
                 { mat_mass = new MatLapack (nbddl,nbddl,true,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
               else // sinon on redimentionne ou on change si nécessaire
                { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                  mat_mass_d->Change_taille(nbddl,nbddl,0.);
                };
               break;
              }
             case RECTANGLE_LAPACK :
              {// cas s'une matrice rectangulaire, pour l'instant idem que le cas précédent
               if (mat_mass == NULL)
                 { mat_mass = new MatLapack (nbddl,nbddl,false,0.,pa.Type_resolution(),pa.Type_preconditionnement());}
               else // sinon on redimentionne ou on change si nécessaire
                { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                  mat_mass_d->Change_taille(nbddl,nbddl,0.);
                };
               break;
              }
             case BANDE_SYMETRIQUE :
              {// initialisation a zero
               if (mat_mass == NULL)
                 { mat_mass = new MatBand (BANDE_SYMETRIQUE,demi,nbddl,0);}
               else // sinon on redimentionne ou on change si nécessaire
                { MatBand * mat_mass_d = (MatBand*) mat_mass; // cas où on redimentionne la matrice
                  mat_mass_d->Change_taille(demi,nbddl);
                  mat_mass_d->Initialise (0.);
                };
               break;
              }
             case BANDE_NON_SYMETRIQUE_LAPACK : case BANDE_SYMETRIQUE_LAPACK: // le choix est fait au moment de l'allocation
              {// initialisation a zero
               if (pa.Type_Matrice() == BANDE_SYMETRIQUE_LAPACK)
                { if (mat_mass == NULL)
                   { mat_mass = new MatLapack(pa.Type_Matrice(),demi,nbddl,true,0.0
                                        ,pa.Type_resolution(),pa.Type_preconditionnement());}
                  else // sinon on redimentionne  si nécessaire
                   { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                     mat_mass_d->Change_taille(demi,nbddl,0.0);
                   };
                }
               else
                { if (mat_mass == NULL)
                   { mat_mass = new MatLapack(pa.Type_Matrice(),total,nbddl,false,0.0
                                       ,pa.Type_resolution(),pa.Type_preconditionnement());}
                  else // sinon on redimentionne  si nécessaire
                   { MatLapack * mat_mass_d = (MatLapack*) mat_mass; // cas où on redimentionne la matrice
                           mat_mass_d->Change_taille(total,nbddl,0.0);
                   };
                };
               break;
              }
             case CREUSE_COMPRESSEE_COLONNE :
              { //  cas d'une matrice creuse compressée par colonne
                // tout d'abord on demande l'ensemble des tables de connexion
                // pour les éléments
                Tableau < Tableau <int> > petites_matrices;
                lesMail->Table_connexion(petites_matrices,nb_casAssemb);
                // puis on définit la matrice en conséquence
               
                if (mat_mass == NULL)
                 { mat_mass = new Mat_creuse_CompCol(nbddl,nbddl, petites_matrices);}
                else // sinon on redimentionne ou on change si nécessaire
                 { if (mat_mass->Type_matrice() !=  CREUSE_COMPRESSEE_COLONNE)
                        { delete mat_mass; mat_mass = new Mat_creuse_CompCol(nbddl,nbddl, petites_matrices);} // cas où on change de matrice
                   else {Mat_creuse_CompCol * mat_mass_d = (Mat_creuse_CompCol*) mat_mass; // cas où on redimentionne la matrice
                         mat_mass_d->Change_taille(nbddl,nbddl, petites_matrices);
                        };
                 };
                // pour l'instant on ne prend pas en compte les conditions linéaires initiales qui augmente la largeur
                // diagonale par bloc de dim
                // par contre on considère que par défaut la diagonale par bloc existe !!
                if (change_stockage)
                 { cout << "\n*** Erreur de choix de matrice de masse, il y a des conditions lineaires intiales,"
                        << " qui vont changer le stockage a l'execution. Or le cas CREUSE_COMPRESSEE_COLONNE ne "
                        << " prend pas en compte pour l'instant ce cas de changement   !! \n";
                   if (ParaGlob::NiveauImpression() >= 2)
                     cout << " Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite(..)";
                   // dans le cas où un comptage du calcul est en cours on l'arrête
                   if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                   if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                   if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                   Sortie(1);
                 };
                break;
              }
             default :
                cout << "\nErreur : valeur incorrecte ou non implante du type  de matrice retenue, type  !\n"
                     << " = " << pa.Type_Matrice() ;
                cout << "\n Algori::Algori::Choix_matrice_masse(... ) \n";
                // dans le cas où un comptage du calcul est en cours on l'arrête
                if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                Sortie(1);
            };
          if (ParaGlob::NiveauImpression() > 2)
             { cout << "\n"; if (change_stockage) cout << " changement de ";
               cout << " stockage matrice masse -->" << Nom_matrice(mat_mass->Type_matrice())
                    << " nb_ligne= "<< mat_mass->Nb_ligne() <<" nb_col= "<< mat_mass->Nb_colonne() << endl;
             };
        }
       else // sinon pour l'instant les autres cas de contact ne sont pas traités
         {cout << "\n*** Erreur pour la mise a jour de la matrice de masse "
               << " avec du contact, le type de contact = " <<pa.ContactType()
               << " n'est pas pris en compte  \n";
          if (ParaGlob::NiveauImpression() >= 5)
              cout << " Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite(..)";
         };
      }
     else // sinon on peut augmenter la taille directement
      {Enum_matrice enumatrice = mat_mass->Type_matrice();// init

       switch (enumatrice)
        { case CARREE : break; // pas de modification
          case RECTANGLE : break; // pas de modification
          case CARREE_LAPACK : break; // pas de modification
          case CARREE_SYMETRIQUE_LAPACK : break; // pas de modification
          case RECTANGLE_LAPACK : break; // pas de modification
          case BANDE_SYMETRIQUE : case BANDE_SYMETRIQUE_LAPACK:
            { // -- récup des anciennes tailles
              int nbddl = mat_mass->Nb_ligne();
              ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
              int ancien_demi = mat_mass->Nb_colonne();
              // calcul de la largeur de bande effective en tenant compte du contact
              int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
              int demi = parstock.demi; // init
              // on n'intervient que si la largeur en noeud augmente
//              if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
              if (largeur_bande_contact_sup_1)
                {// la mise à jour dépend du type de contact
                 if (pa.ContactType() == 1) // contact cinématique, sans multiplicateur ni pénalisation
                   // on doit augmenter la lb en fonction de cumule, si la largeur en noeud est différente
                   // de 1
                   {cumule = std::min(cumule,nbddl);
                     demi = parstock.demi + cumule; // c'est donc toute la largeur hors diag que l'on ajoute
                     if (demi > ancien_demi)
                       if (ParaGlob::NiveauImpression() >= 5)
                            cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << "old_1/2lb: "<<ancien_demi << " new_1/2lb: "<< demi <<" ";
                   }
                 else if ((pa.ContactType() == 2) // cas d'un contact par pénalisation
                          || (pa.ContactType() == 3)|| (pa.ContactType() == 4)
                         )
                   // dans ce cas on considère uniquement le LB maxi
                   { demic = std::min(demic,nbddl);
                     int demi_inter = std::max(parstock.demi, demic);
                     if (demi_inter > ancien_demi)
                       if ((ParaGlob::NiveauImpression() >= 5) && (demic > demi))
                       cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << (demic-demi)  << "old_1/2lb: "<<ancien_demi << " new_1/2lb: "<< demi_inter <<" ";
                     demi = demi_inter;
                   };
                 };
              // on limite au nombres total de ddl et pas inférieur à la dimension de l'espace 
              // car on a au minimum une largeur de 1 noeud donc 3 ddl
              // ( à moins qu'il n'y ait plus de contact) mais pour l'instant on n'en tient pas compte ....
              demi = std::max(std::min(demi,nbddl),ParaGlob::Dimension());

              // on change la taille, si elle doit augmenter, ou si elle est inférieur d'au moins 5% par rapport à la taille actuelle
              if ( (demi > mat_mass->Nb_colonne()) ||
                  (demi < 0.95 * mat_mass->Nb_colonne()))
                   // on change la taille
                { if (ParaGlob::NiveauImpression() >= 5)
                    cout << "\n augmentation effective: LB/2+1 finale(en ddl) " << demi  << " ";
                  switch (enumatrice)
                   { case BANDE_SYMETRIQUE:
                       ((MatBand*)mat_mass)->Change_taille(demi,nbddl );
                       break;
                     case BANDE_SYMETRIQUE_LAPACK:
                       ((MatLapack*)mat_mass)->Change_taille(nbddl,demi);
                       break;
                     default:
                       cout << "\n erreur non prevue Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite "<<endl;
                       // dans le cas où un comptage du calcul est en cours on l'arrête
                       if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                       if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                       if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                       Sortie(1);
                   };
                };
              break;
            }
          case BANDE_NON_SYMETRIQUE : case BANDE_NON_SYMETRIQUE_LAPACK :
            { // -- récup des anciennes tailles
              int nbddl = mat_mass->Nb_ligne();
              int ancien_lb = mat_mass->Nb_colonne();
              ParaStockage & parstock = (tab_para_stockage(nb_casAssemb.n));
                  //int demi = matglob->Nb_colonne();
              // calcul de la largeur de bande effective en tenant compte du contact
              int demic=demicontact,totalc=totalcontact,cumule=cumulecontact;
              int total = parstock.total; // init
              // on n'intervient que si la largeur en noeud augmente
//              if (lescontacts->Largeur_Bande(demic,totalc,nb_casAssemb,cumule))
              if (largeur_bande_contact_sup_1)
                {// la mise à jour dépend du type de contact
                 if (pa.ContactType() == 1) // contact cinématique, sans multiplicateur ni pénalisation
                  // on doit augmenter la lb en fonction de cumule
                  { cumule = std::min(cumule,nbddl);
                    int total_inter = parstock.total + 2 * cumule;
                    if (total_inter > ancien_lb)
                      if (ParaGlob::NiveauImpression() >= 5)
                           cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << cumule  << "old_lb: "<<ancien_lb << " new_lb: "<< total_inter <<" ";
                    total = total_inter;
                  }
                 
                 else if ((pa.ContactType() == 2) // cas d'un contact par pénalisation
                         || (pa.ContactType() == 3)|| (pa.ContactType() == 4)
                         )
                  // dans ce cas on considère uniquement le LB maxi
                  { totalc = std::min(totalc,2*nbddl-1);
                    int total_inter = std::max(parstock.total, totalc);
                    if (total_inter > ancien_lb)
                       if ((ParaGlob::NiveauImpression() >= 5) && (totalc > total))
                           cout << "\n propositon d'augmentation de la largeur de bande pour les conditions de contact (en ddl) " << (totalc-total)  << "old_lb: "<<ancien_lb << " new_lb: "<< total_inter << " ";
                    total = total_inter;
                 };
                };
              // on limite au nombres total de ddl
              total = std::min(total,2*nbddl-1);

              // on change la taille, si elle doit augmenter, ou si elle est inférieur d'au moins 5% par rapport à la taille actuelle
              if ( (total > mat_mass->Nb_colonne()) ||
                  (total < 0.95 * mat_mass->Nb_colonne()))
                       // on change la taille
                { if (ParaGlob::NiveauImpression() >= 5)
                            cout << "\n augmentation effective: LB finale(en ddl) " << total  << " ";
                  switch (enumatrice)
                   { case BANDE_NON_SYMETRIQUE:
                       ((MatBand*)mat_mass)->Change_taille(total,nbddl );
                       break;
                     case BANDE_NON_SYMETRIQUE_LAPACK:
                       ((MatLapack*)mat_mass)->Change_taille(nbddl,total);
                       break;
                    
                     default:
                       cout << "\n erreur2 non prevue Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite "<<endl;
                       // dans le cas où un comptage du calcul est en cours on l'arrête
                       if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
                       if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
                       if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
                       Sortie(1);
                   };
                };
              break;
            }
          case CREUSE_COMPRESSEE_COLONNE :
            { cout << "\n *** cas avec contact non actuellement traite "
                   << "\n Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite(... ) \n";
              // dans le cas où un comptage du calcul est en cours on l'arrête
              if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
              if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
              if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
              Sortie(1);
              break;
            }
          default :
             cout << "\nErreur : valeur incorrecte ou non implante du type  de matrice retenue, type  !\n"
                  << " = " << enumatrice ;
             cout << "\n Algori::Mise_a_jour_type_et_taille_matrice_masse_en_explicite(... ) \n";
             // dans le cas où un comptage du calcul est en cours on l'arrête
             if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
             if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
             if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
             Sortie(1);
        };
       if (ParaGlob::NiveauImpression() >= 5)
          cout  << "\n stockage matrice masse " << enumatrice << " taille (entiers) :"
                << mat_mass->Place() << endl;

      }; // fin changement matrice

      
    };
 
  if (ParaGlob::NiveauImpression() >= 7)
     cout << "\n Taille du stockage masse " << mat_mass->Type_matrice() << " (taille en entiers) :"
          << mat_mass->Place() << endl;
  // retour
  if (change_stockage)
       {return mat_mass;}
  else {return NULL;};

  
};

    // calcul la matrice de masse 
    // N_ddl : donne le ddl où doit être mis la masse
void Algori::Cal_matrice_masse(LesMaillages * lesMail,Assemblage& Ass
                               ,Mat_abstraite& mat_mass
                               ,const DiversStockage* divStoc,LesReferences* lesRef
                               ,const Enum_ddl & N_ddl ,LesFonctions_nD*  lesFonctionsnD)
{   tempsCalMasse.Mise_en_route_du_comptage();
    bool use_diagonale = false; // par défaut on n'utilise pas de stockage diagonal
    if (    (mat_mass.Type_matrice() == DIAGONALE)
	       || (pa.Type_calcul_masse() == MASSE_DIAG_COEF_EGAUX)
        || (pa.Type_calcul_masse() == MASSE_DIAG_COEF_VAR)
      )
			    {use_diagonale = true;}
		  else if (   (mat_mass.Type_matrice() == DIAGONALE)
	            && (pa.Type_calcul_masse() == MASSE_CONSISTANTE)
            )
       { cout << "\n *** erreur, le stockage globale est diagonale et on veut "
              << " egalement une matrice masse consistante !! c'est incoherent. ";
         if (ParaGlob::NiveauImpression() > 5)
            cout << "\n Algori::Cal_matrice_masse( ";
         Sortie(1);
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
         if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
         if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
       };
 
		  if (use_diagonale)	  
       // on définit un vecteur identique à la diagonale si c'est la première fois
       { if (toto_masse == NULL)
            toto_masse = new Vecteur(mat_mass.Nb_ligne());
          else
            // on initialise le vecteur
            toto_masse->Zero();
       };
    // on initialise la matrice masse à 0
    mat_mass.Initialise(0.);
    
    // boucle sur les elements
    for (int nbMail =1; nbMail<= lesMail->NbMaillage(); nbMail++)
     for (int ne=1; ne<= lesMail->Nombre_element(nbMail);ne++)
        { 
          //calcul de la matrice masse  local
          Element & el = lesMail->Element_LesMaille(nbMail,ne); // l'element
          Tableau<Noeud *>& taN = el.Tab_noeud(); // tableau de noeuds de l'el
          Mat_pleine * mat_mass_loc = el.CalculMatriceMasse (pa.Type_calcul_masse());
     //    cout << "numéro =" << ne << endl;
     //     mat_mass_loc->Affiche();
    //      mat_mass_loc->Affiche1(1,8,1,1,8,1);
          // assemblage
          // il faut distinguer si l'on a une matrice globale diagonale ou si l'on a une matrice pleine
				      // et ... il faut tenir compte du type de calcul de matrice masse locale
				      if (use_diagonale)
//              if (mat_mass.Type_matrice() == DIAGONALE)
					       {  // récupération de la première ligne
						         Vecteur* vectloc = &(mat_mass_loc->Ligne_set(1)) ;
						         // ensuite on fait l'assemblage suivant un vecteur comme pour le second membre
						         Ass.AssemSM (*toto_masse,*vectloc,el.TableauDdl(),taN);
            }
          else // cas d'une matrice pleine, donc consistante
            { if (pa.Symetrie_matrice())
                   Ass.AssembMatSym (mat_mass,*mat_mass_loc,el.TableauDdl(),taN); // de la masse
              else
                   Ass.AssembMatnonSym (mat_mass,*mat_mass_loc,el.TableauDdl(),taN); // de la masse
					       };
        };
    // ---- cas des masses concentrées -----
    Ajout_masses_ponctuelles(lesMail,Ass,mat_mass,divStoc,lesRef,N_ddl,lesFonctionsnD);
    // ---- fin cas des masses concentrées -----
    
            
    // dans le cas où la masse est diagonale on transfert du vecteur assemblé à la matrice
    // diagonale
    if ((use_diagonale) && (mat_mass.Type_matrice() == DIAGONALE))
      { MatDiag & mat_d = *((MatDiag*) & mat_mass);
        mat_d += (*toto_masse);
      }
		  else if (use_diagonale)  
		    { // cas d'une matrice masse concentrée sur la diagonale, mais d'un stockage non diagonale
			     MatDiag mat_inter(mat_mass.Nb_ligne());
				    mat_inter = (*toto_masse);
			     mat_mass += mat_inter;
      };
		  // sinon l'assemblage est déjà ok	  
	  
    // affichage éventuelle de la matrice de masse       
    if (ParaGlob::NiveauImpression() >= 10)
     { string entete = " affichage de la matrice masse ";
       mat_mass.Affichage_ecran(entete); 
     };
    tempsCalMasse.Arret_du_comptage();
};         


// uniquement ajout dans la matrice de masse des masses ponctuelles  si elles existent
// N_ddl : donne le ddl où doit être mis la masse
// par exemple, sert pour la relaxation dynamique
void Algori::Ajout_masses_ponctuelles(LesMaillages * lesMail,Assemblage& Ass,Mat_abstraite& mat_masse
                               ,const DiversStockage* divStoc,LesReferences* lesRef
                               ,const Enum_ddl & N_ddl,LesFonctions_nD*  lesFonctionsnD)
 {
  #ifdef UTILISATION_MPI
    // cas d'un  calcul //, seule la matrice du CPU 0 est concernée
    if (ParaGlob::Monde()->rank() != 0)
      return ;
  #endif
  tempsCalMasse.Mise_en_route_du_comptage();
  // prise en compte de masses ponctuelles
  int dima = ParaGlob::Dimension();
  // dans le cas où on est en axisymétrie il faut diminuer de 1
  if (ParaGlob::AxiSymetrie())
    dima--;
//  int nbddl_X = mat_masse.Nb_ligne();
  // ---- cas des masses concentrées -----
  // on récupère les masses concentrées
  const Tableau < BlocDdlLim<BlocScal_ou_fctnD> > & Tab_M_A = divStoc->TabMasseAddi();
  int cas_assemblage = (Ass.Nb_cas_assemb()).n; // récup du cas d'assemblage
  int posi = N_ddl -1;
  // on balaie le tableau
  int  Tab_M_ATaille = Tab_M_A.Taille();
  for (int i=1;i<= Tab_M_ATaille;i++)
   { // on récup la référence
     const ReferenceNE & ref =
     ((ReferenceNE &) lesRef->Trouve(Tab_M_A(i).NomRef(),Tab_M_A(i).NomMaillage()));

     // recup de la masse: soit c'est une valeur fixe soit une fonction nD
     double val_mass = 0; // init par défaut
     Fonction_nD * pt_fonct = NULL; // init
     Coordonnee mass_noeu_differencie(dima);
     int cas_masse = 0; // = 1 -> masse fixe sur tous les axes
                        // = 2 -> masse différenciée sur chaque axe
                        // = 0 -> n'a pas encore de valeur
     string nom_fct_nD("_"); // ""
     if ((Tab_M_A(i)).Val() != NULL) // cas d'une valeur fixe
       { val_mass = *((Tab_M_A(i)).Val()); // recup de la valeur
         cas_masse = 1;
       }
     else if ((Tab_M_A(i)).Fct_nD() != NULL) // cas d'une fonction nD
       { nom_fct_nD = *((Tab_M_A(i)).Fct_nD());
         // on récupère la fonction
         if (lesFonctionsnD->Existe(nom_fct_nD))
           {pt_fonct  = lesFonctionsnD->Trouve(nom_fct_nD);
            if (pt_fonct->Nom_variables().Taille() == 0) // cas où il n'y a que des variables globales
             { // on vérifie qu'en retour on a un scalaire ou un vecteur de dima = la dima de l'espace
               Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
               if (pt_fonct->NbComposante() == 1)
                {val_mass = tava(1);cas_masse=1;}
               else if (pt_fonct->NbComposante() == dima)
                {switch (dima)
                      { case 3 : mass_noeu_differencie(3) = tava(3);
                        case 2 : mass_noeu_differencie(2) = tava(2);
                        case 1 : mass_noeu_differencie(1) = tava(1);
                        cas_masse=2;
                      };
                }
               else
                {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                      << " cas des masses ponctuelles sur la ref: "<< Tab_M_A(i).NomRef()
                      << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                      << dima << endl;
                 Sortie(1);
                };
             }; // fin du cas où la fonction ne dépend que de variables globales
           }
         else // sinon c'est un problème
           { cout << "\n *** erreur dans la definition des masses ponctuelles via la fonction nD "
                  << nom_fct_nD << ", a priori cette fonction n'est pas disponible!! ";
             Sortie(1);
           };
      }
     else // sinon il y a un pb
       { cout << "\n *** erreur de definition des masses ponctuelles aux noeuds"
              << " les infos lues ne sont pas correctes: ";
         (Tab_M_A(i)).Affiche();
         Sortie(1);
       };
  
     // maintenant on va balayer les noeuds de la ref
     int reftaille = ref.Taille();
     // on boucle sur les noeuds de la référence
     for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
      { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
        if ((noe.Pointeur_assemblage(N_ddl,cas_assemblage) != -1)
           && (noe.En_service(N_ddl)) && (noe.UneVariable(N_ddl)))
       {// les masses ponctuelles correspondent aux ddl N_ddl N_ddl+1  ...
        switch (cas_masse)
          {case 0: // il s'agit de l'utilisation d'une fonction nD qui ne doit dépendre que de M
             {if ((pt_fonct->Depend_M() < 0) || (pt_fonct->Depend_Mt() < 0)
                 || (pt_fonct->Depend_M0() < 0)
                ) // signifie que la fonction dépend exclusivement d'un des 3 temps: tdt, t ou 0
                { Coordonnee M(dima); // init
                
                  // on vérifie également qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace
                  if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima)
                       M = noe.Coord2();
                  if ((pt_fonct->Depend_Mt()<0) && pt_fonct->Nom_variables().Taille() == dima)
                       M = noe.Coord1();
                  if ((pt_fonct->Depend_M0()<0) && pt_fonct->Nom_variables().Taille() == dima)
                       M = noe.Coord0();

                    // comme pt_fonct->Depend_M()<0, cela signifie que la fonction dépend que de M
                    // donc dans les paramètres d'appel on ne transmet que M
                    Tableau <Coordonnee> tab_M(1,M);
                    Tableau <double> & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     { double& val_masse_fct = tava(1);
                       for (int di=1;di<= dima;di++)
                        { int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                          mat_masse (pt_adres,pt_adres) +=  val_masse_fct;
                        };
                     }
                    else if (pt_fonct->NbComposante() == dima)
                     { int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi),cas_assemblage);
                     
                     
                       switch (dima)
                           { case 3 : mat_masse(pt_adres+3,pt_adres+3) +=  tava(3);
                             case 2 : mat_masse(pt_adres+2,pt_adres+2) +=  tava(2);
                             case 1 : mat_masse(pt_adres+1,pt_adres+1) +=  tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD "
                           << " cas des masses ponctuelles sur la ref: "<< Tab_M_A(i).NomRef()
                           << " en retour, le nombre de composante: " << pt_fonct->NbComposante()
                           << " est different de 1 ou de la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                 }
              else
                {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD "
                      << " cas des masses ponctuelles sur la ref: "<< Tab_M_A(i).NomRef()
                      << " la fonction ne doit dependre que de M : valeur courante, ou a t "
                      << " ou au temps initial "
                      << dima ;
                 cout << "\n fonction nD: "; pt_fonct->Affiche();
                 Sortie(1);
                };
             }
             break;
           case 1: // la masse a déjà été déterminé, et elle est identique sur tous les axes
             for (int di=1;di<= dima;di++)
               { int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                 mat_masse(pt_adres,pt_adres) +=  val_mass;
               };
             break;
           case 2: // la masse a déjà été déterminé, et elle est différentes sur tous les axes
             for (int di=1;di<= dima;di++)
               { int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                 mat_masse(pt_adres,pt_adres) +=  mass_noeu_differencie(di);
               };
             break;
           default:
             cout << "\n *** erreur dans l'ajout de masses ponctuelles "
                  << " cas non prevu \n Algori::Ajout_masses_ponctuelles(...";
             Sortie(1);
             break;
          };
       };
      };
   };
       // ---- fin cas des masses concentrées -----
  tempsCalMasse.Arret_du_comptage();
 };

// définition de la matrice de viscosité numérique, dans le cas d'une approximation
// diagonale de la raideur, dans un calcul explicit 
// la matrice masse doit être diagonale 
// inita indique si l'on est dans une phase d'initialisation (creation de la matrice) ou non 
//  dans ce dernier cas il y a modification éventuelle des valeurs de C
// ramène un pointeur sur la matrice visqueuse instancié ou nul si aucune instantiation 
Mat_abstraite* Algori::Cal_mat_visqueux_num_expli(const Mat_abstraite& mat_mass,Mat_abstraite* mat_C_pt
                                                  ,const Vecteur & delta_X,bool inita, const Vecteur & vitesse)
 {
   #ifdef UTILISATION_MPI
     // cas d'un  calcul //, seule la matrice du CPU 0 est concernée
     if (ParaGlob::Monde()->rank() != 0)
       return mat_C_pt;
   #endif
   tempsCalViscoNum.Mise_en_route_du_comptage();
   // on regarde s'il ne faut pas redimentionner la matrice: 
	  // si inadéquation: on supprime la matrice C existante
	  if (inita)
	   {if (mat_C_pt != NULL)
			   { if ((mat_mass.Type_matrice() != mat_C_pt->Type_matrice()) ||
					       (mat_mass.Nb_ligne() != mat_C_pt->Nb_ligne()))
				    delete mat_C_pt;
			   }; 
		   // maintenant création éventuelle 
	    if (mat_C_pt == NULL)
	         {mat_C_pt = mat_mass.NouvelElement();}; // amortissement du même type que la masse
	    //idem avec la matrice de travail C_inter
	    if (C_inter != NULL)
			   { if ((mat_mass.Type_matrice() != C_inter->Type_matrice()) ||
					       (mat_mass.Nb_ligne() != C_inter->Nb_ligne()))
				    delete C_inter;
			   }; 
		   // maintenant création éventuelle
		   if (C_inter == NULL)
	      C_inter = mat_mass.NouvelElement(); // amortissement du même type et de même valeur que la masse
		   // dimensionnement des forces int et ext mémorisée à l'itération précédente
		   F_int_tdt_prec.Change_taille(mat_C_pt->Nb_ligne());		
		   F_ext_tdt_prec.Change_taille(mat_C_pt->Nb_ligne());
		  };
			 
	  // on s'occupe maintenant des valeurs		 
   switch (pa.Amort_visco_artificielle())
   	{case 1 : // cas d'un amortissement de Rayleigh classique
   	  {if (inita)
   	    {if ((ParaGlob::NiveauImpression() >= 0)&&((pa.CoefRayleighMasse()!= 1.)
   	  	     ||(pa.CoefRayleighRaideur()!=0.)))
           { cout << "\n ** ATTENTION : dans le cas d'un algorithme explicite avec un amortissement "
                  << " de Rayleigh classique la matrice d'amortissement numerique [C] n'est construite"
                  << "  qu'avec la matrice masse selon : [C] = nu [M], ainsi le coefficient de"
                  << " Rayleigh sur la raideur est impose a 0 et sur la masse a 1";
           };      
         // on définit la matrice visqueuse en multipliant par la viscosité
			      (*mat_C_pt) = mat_mass;  // assignation à la valeur de la masse
         (*mat_C_pt) *= pa.Visco_artificielle();
   	  	 };
   	   break;	
   	  }
   	 case 2 : // cas d'un amortissement en fonction de l'amortissement critique approché
		         // première méthode
   	  { if (inita)
   		    {if (mat_mass.Type_matrice() != DIAGONALE)
 	   	 	  { cout << "\n erreur, la matrice masse doit etre diagonale ";
 	   	      if (ParaGlob::NiveauImpression()>5)
 	   	        cout << "\n Algori::Cal_mat_visqueux_num(...";
            // dans le cas où un comptage du calcul est en cours on l'arrête
            if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
            if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
            if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
 	   	 	    Sortie(1);
          };
   		    }; // fin initialisation
           
        // calcul des termes de la matrice visqueuse
        MatDiag & m_masse = *((MatDiag*) & mat_mass);
        MatDiag & m_C = *((MatDiag*)  mat_C_pt);
        int mtaille = m_masse.Nb_ligne();
        double eta = pa.Visco_artificielle();
        double coef_limite_C = pa.BorneMaxiSurCfonctionDeLaMasse();
        for (int i=1;i<=mtaille;i++)
        	{ double deltaxi=delta_X(i);
        	  if (Dabs(deltaxi) <= ConstMath::petit) deltaxi = ConstMath::petit;
        	  double m_masseii=m_masse(i,i);
        	  double m_Cii =2.*sqrt(Dabs(m_masseii *(F_int_tdt(i)-F_int_t(i)+F_ext_tdt(i)-F_ext_t(i))/deltaxi));
        	  // on limite la valeur 
        	  if (m_Cii > coef_limite_C * m_masseii)
        	      m_Cii = coef_limite_C * m_masseii;
        	  m_C(i,i)=eta * m_Cii; // on définit directement chaque terme
        	};
   	    break;	
   	  }
   	 case 3 : // cas d'un amortissement en fonction de l'amortissement critique approché
		         // deuxième méthode à l'aide du quotient de Rayleight sur les forces internes
					      // C = eta * 2. sqrt(lambda_0) * M,  avec lambda_0 = DX^T * (delta_F_int) / DX^T M DX
					      // le delta_F_int est supposé = K_tangent DX
		    {
         // --- calcul des termes de la matrice visqueuse
		       // on ne modifie la valeur de la matrice visqueuse, que s'il y a eu un déplacement, sinon on garde la valeur précédament stockée
		       if (delta_X.Norme() > ConstMath::trespetit)
		        { double lambda_0 = ((delta_X * F_int_tdt) - (delta_X * F_int_t)) / mat_mass.vectT_mat_vec(delta_X , delta_X );
			         (*mat_C_pt) = mat_mass;  // assignation à la valeur de la masse
				        // prise en compte des limitations
				        double coef = 2. * sqrt(lambda_0);
				        double coef_limite_C = pa.BorneMaxiSurCfonctionDeLaMasse();
				        if (coef > coef_limite_C)
				            coef = coef_limite_C;
            (*mat_C_pt) *= (pa.Visco_artificielle() * coef);
		        };  
   	     break;	
   	  }
       //   ----------------- cas de l'amortissement critique utilisé par exemple avec les algo de relaxation dynamique ----------------
		   case 4 : 
		     CalculEnContinuMatriceViscositeCritique(mat_mass,*mat_C_pt,delta_X,vitesse);
		   break;																		
       // dans les autres cas on ne fait rien   	
   	};
    
   tempsCalViscoNum.Arret_du_comptage();
   return mat_C_pt; // retour de la matrice 	
 };
                               
// mise à jour du calcul de la viscosité critique en continu
// on ajoute la viscosité critique ce qui permet de cumuler avec de la viscosité numérique
void Algori::CalculEnContinuMatriceViscositeCritique(const Mat_abstraite& mat_mass,Mat_abstraite& mat_C_pt
                                                             ,const Vecteur & delta_X, const Vecteur & vitesse)
 {
    #ifdef UTILISATION_MPI
      // cas d'un  calcul //, seule la matrice du CPU 0 est concernée
      if (ParaGlob::Monde()->rank() != 0)
        return ;
    #endif
    tempsCalViscoNum.Mise_en_route_du_comptage();
    // ---- on commence par calculer (omega_0)^2 qui représente la plus basse fréquence
	   double omega02=0.; // init
    switch (opt_cal_C_critique)
   	{case -1 : // version test : c= 0
		    { omega02 = 0.;  // pas de prise en compte de viscosité
			     break;
		    };	 
   	 case 0 : // utilisation d'une approximation d'un quotient de Rayleigh
        // (omega_0)^2 \approx  (dot X^T Delta R_{i(statique)}^n) / (dot X^T M dot X)
		    { // on ne calcule que s'il y a eu une vitesse non nulle, sinon on ne fait rien
		      if (vitesse.Norme() > ConstMath::trespetit)
		         omega02 = -((vitesse * F_int_tdt) - (vitesse * F_int_t) 
				                 + (vitesse * F_ext_tdt) - (vitesse * F_ext_t)) / mat_mass.vectT_mat_vec(vitesse , vitesse );
			     break;
		    };	 
     case 1 : // utilisation d'une approximation d'un quotient de Rayleigh
        // (omega_0)^2 \approx  (Delta X^T Delta R_{i(statique)}^n) / (Delta X^T M Delta X)
		    { // on ne calcule que s'il y a eu un déplacement, sinon on ne fait rien
		      if (delta_X.Norme() > ConstMath::trespetit)
           omega02 = -((delta_X * F_int_tdt) - (delta_X * F_int_t)
								             + (delta_X * F_ext_tdt) - (delta_X * F_ext_t)) / mat_mass.vectT_mat_vec(delta_X , delta_X );
			     break;
      };	 
   	 case 2 : // utilisation d'une approximation d'un quotient de Rayleigh
		         // (omega_0)^2 \approx  (dot X^T Delta R_{i(statique)}^n) / (dot X^T M dot X)
		    { // on ne calcule que s'il y a eu une vitesse non nulle, sinon on ne fait rien
		      if (vitesse.Norme() > ConstMath::trespetit)
		         omega02 = -((vitesse * F_int_tdt) - (vitesse * F_int_tdt_prec) 
				                 + (vitesse * F_ext_tdt) - (vitesse * F_ext_tdt_prec)) / mat_mass.vectT_mat_vec(vitesse , vitesse );
			     break;
		    };	 
     default :
		     cout << "*** erreur le cas du coefficient opt_cal_C_critique= " << opt_cal_C_critique
					       << " n'est pas pris actuellement en compte car cela signifie qu'il n'est pas defini !!, "
            << "\n *** peut-etre revoir les pararametres de l'algorithme ***  ,  "
				       	<< "\n Algori::CalculEnContinuMatriceViscositeCritique(..." << endl;
       // dans le cas où un comptage du calcul est en cours on l'arrête
       if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
       if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
       if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
			    Sortie(1);
		  };
    // on regarde pour les limitations de la valeur de omega_0^2 et on en déduit la valeur d'omega0
	   double omega0 = 0.;
    int nbcar = ParaGlob::NbdigdoEC();

    //.... debug 	
    if (ParaGlob::NiveauImpression() >= 8)
      cout << "\n omega02= " << std::fixed << std::setprecision(nbcar)<< omega02
           << " llF_int_tdtll= " << std::fixed << std::setprecision(nbcar)<< F_int_tdt.Norme()
           << " llF_int_tll= " << std::fixed << std::setprecision(nbcar)<< F_int_t.Norme() << endl;
    //.... fin debug 	 
	 
	   if (omega02 < 0.)  
	        {  omega0 = 0.; }
	   else if (omega02 > (f_ * f_ * 4.))
	        {  omega0 = f_ * 2.; }
	   else {  omega0 = sqrt(omega02);};
	 
	   // ---- maintenant on calcul une viscosité critique globale
	   double viscosite = 0.;
	   switch (type_calcul_visqu_critique)
	   { case 1: // cas de la méthode historique d'Underwood C_{ii} += c * M_{ii} et c = 2 * omega_0
		     { viscosite = (omega0 * 2.);
			      break;
			    }
	     case 2: // Pajand: la matrice C est calculee selon: C_{ii} = c * M_{ii}, avec c = sqrt( 4*omega_0^2 - omega_0^4)
		     { double coef = 4.* PUISSN(omega0, 2) - PUISSN(omega0, 4);
			      if (coef < 0.)
			        { if (ParaGlob::NiveauImpression() > 7)
				           cout << "\n *** attention la grandeur ( 4*omega_0^2 - omega_0^4)= "
                    << coef << "  est negative, on la met a 0 ";
					        coef = 0.;
				       };
			      viscosite = sqrt(coef);	  	  
			      break;
			    }
	     case 3: // la matrice C est calculee selon: C_{ii} = c * M_{ii}, avec c = 2*sqrt(omega_0/(1+omega_0))
		     { double coef = omega0/(1.+ omega0);
			      if (coef < 0.)
			        { if (ParaGlob::NiveauImpression() > 7)
				           cout << "\n *** attention la grandeur omega_0/(1+omega_0)= "
                    << coef << "  est negative, on la met a 0 ";
					        coef = 0.;
				       }
			      viscosite = sqrt(coef);	  	  
			      break;
       }
		    default:
				   { cout << "\n **** erreur le cas type_calcul_visqu_critique= " << type_calcul_visqu_critique 
					         << " n'est pas pris actuellement en compte car cela signifie qu'il n'est pas defini !!, "
              << "\n *** peut-etre revoir les pararametres de l'algorithme ***  ,  "
              << "\n Algori::CalculEnContinuMatriceViscositeCritique( ... ";
         // dans le cas où un comptage du calcul est en cours on l'arrête
         if (tempsCalEquilibre.Comptage_en_cours()) tempsCalEquilibre.Arret_du_comptage();
         if (tempsInitialisation.Comptage_en_cours()) tempsInitialisation.Arret_du_comptage();
         if (tempsMiseAjourAlgo.Comptage_en_cours()) tempsMiseAjourAlgo.Arret_du_comptage();
         Sortie(1);
       };			
	   };

	   // affichage éventuelle de la viscosité      
	   if (ParaGlob::NiveauImpression() >= 7)
         { cout << "\n viscosite c= " << std::fixed << std::setprecision(nbcar)<< viscosite << endl ;
           cout << endl;
         };
		
	 // ---- maintenant on va calculer la matrice visqueuse:  C_{ii} += c * M_{ii} 

//-- cas ou on cumulait, abandonné actuellement
//	 (*C_inter) = mat_mass;
//	 (*C_inter) *= viscosite * ampli_visco;
//	 mat_C_pt += (*C_inter);
//-- cas ou on cumulait, abandonné actuellement
	 
	 
	   mat_C_pt = mat_mass;
	   mat_C_pt *= viscosite * ampli_visco;	 
    tempsCalViscoNum.Arret_du_comptage();
 };
                               

// calcul de l'impact d'une viscosité artificielle sur la raideur et sur le second membre
// l'algo travaille avec les vecteurs: delta_X, var_delta_X, et vitesse stockés dans Algori
// l'algo modifie mat_glob et calcul les forces visqueuses
void Algori::Cal_mat_visqueux_num_stat(Mat_abstraite& mat_glob,Vecteur& forces_vis_num)
 {
   #ifdef UTILISATION_MPI
     // cas d'un  calcul //, seule la matrice du CPU 0 est concernée
     if (ParaGlob::Monde()->rank() != 0)
       return ;
   #endif
   tempsCalViscoNum.Mise_en_route_du_comptage();
   // le fonctionnement n'est pas identique au cas explicite
   switch (pa.Amort_visco_artificielle())
    {case 1 : // cas d'un amortissement de Rayleigh classique
      {// ici on n'utilise que la matrice de raideur qui est la seule connue en statique
       // K -> K*(1.-nu/delta_t)
       // recup de l'incrément de temps
       double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
       double unSurDeltat=0;
       if (Abs(deltat) >= ConstMath::trespetit)
        {unSurDeltat = 1./deltat;}
       else
        // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
        { // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
          if (unSurDeltat < 0)
           { cout << "\n le pas de temps est négatif !! "; };
          unSurDeltat = ConstMath::tresgrand;
        };
       
//CoefRayleighRaideur()
       
       // on calcul les forces visqueuses
       // ici il s'agit de force qui dépendent de la variation de delta_X d'une itération à l'autre
       // cela ne doit pas être très stable mais ... à voir
       // mais à convergence, les forces ne sont plus actives puisque var_delta_X = 0
       mat_glob.Prod_mat_vec(
//           ((-pa.Visco_artificielle()*unSurDeltat)*var_delta_X + (-pa.CoefRayleighRaideur()*unSurDeltat)*delta_X)
           ((pa.Visco_artificielle()*var_delta_X )+ (pa.CoefRayleighRaideur()*unSurDeltat)*delta_X)
           ,forces_vis_num);
       // calcul de la nouvelle matrice de raideur qui tient compte des forces visqueuse
//       mat_glob *= (1.+pa.Visco_artificielle()*unSurDeltat+pa.CoefRayleighRaideur()*unSurDeltat);
       mat_glob *= (1.-pa.Visco_artificielle()-pa.CoefRayleighRaideur()*unSurDeltat);
       break;
      }
       // dans les autres cas on ne fait rien
     default:
       if (ParaGlob::NiveauImpression()>0)
         cout << "\n attention on cherche a utiliser de l'amortissement visqueux de type "
              << pa.Amort_visco_artificielle()
              << " celui-ci n'est pas disponible en statique, on n'en tient pas compte ! " << flush;
    };
  tempsCalViscoNum.Arret_du_comptage();
 };

                               
                               
                               
