// 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 "Charge.h"
#include "ConstMath.h"
#include "MathUtil.h"
#include "ReferenceNE.h"
#include "ReferenceAF.h"


   
// mise en place du chargement impose sur le second membre global
// cas explicite
// la variable relative à la classe assemblage est passé en argument ce qui permet d'utiliser
// éventuellement plusieurs type d'assemblage
// retour false si problème de calcul
bool Charge::ChargeSecondMembre_Ex_mecaSolid
       (Assemblage & assemb,LesMaillages * lesMail,LesReferences* lesRef,Vecteur& vecglob
        ,const ParaAlgoControle & pa,LesCourbes1D* lesCourbes1D,LesFonctions_nD*  lesFonctionsnD)
 {
#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 true;
#endif
  temps_cpu_chargement.Mise_en_route_du_comptage(); // temps cpu
  bool retour = true; // init par défaut
  try
   {
     // pour tous les maillages, pour tous les éléments on effectue une initialisation éventuelle
     int nb_mail = lesMail->NbMaillage();
     for (int imail = 1;imail<=nb_mail;imail++)
      { int nb_ele = lesMail->Nombre_element(imail);
        for (int iele =1;iele<=nb_ele;iele++)
          lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement();
      };
   
     int posi = Id_nom_ddl("X1") -1;
     int dim = ParaGlob::Dimension();
    
     // $$$ --- cas des forces ponctuelles --- $$$
    
     // on regarde si l'espace de travail est axi-symétrique, si oui on utilise une dimension réduite
     int dima=dim; // c'est dima qui est utiliser pour la mise en place des efforts ponctuelles
     if (ParaGlob::AxiSymetrie()) dima--;
     // on parcours le tableau tabPonctuel
     for (int i=1;i<= tabPonctuel.Taille();i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocForces> >& tabPonctuel_i = tabPonctuel(i);
         const ReferenceNE & ref =
                  ((ReferenceNE &) lesRef->Trouve(tabPonctuel_i.NomRef(),tabPonctuel_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabPonctuel_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(tabPonctuel_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabPonctuel_i.Pas_a_prendre_en_compte(Temps_courant())))
          {int reftaille = ref.Taille();
           int cas_assemblage = (assemb.Nb_cas_assemb()).n; // récup du cas d'assemblage
           for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
             { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
               Coordonnee force; // init
               // on traite en fonction du fait que c'est un champ ou non
               switch (tabPonctuel_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      force = (tabPonctuel_i.Vect_de_coordonnee(nn)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      // on récupère les courbes pour les dim composantes
                      { switch (dima)
                         { case 3 : {Courbe1D * pt_courbe_force_nn_3 =
                                       lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,3));
                                     force(3) = pt_courbe_force_nn_3->Valeur(Temps_courant());
                                    }
                           case 2 : {Courbe1D * pt_courbe_force_nn_2 =
                                       lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,2));
                                     force(2) = pt_courbe_force_nn_2->Valeur(Temps_courant());
                                    }
                           case 1 : {Courbe1D * pt_courbe_force_nn_1 =
                                       lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,1));
                                     force(1) = pt_courbe_force_nn_1->Valeur(Temps_courant());
                                    }
                         };
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      force = (tabPonctuel_i.Vect_de_coordonnee(1)); break;
                 };
               // on regarde si jamais on a une dépendance à des grandeurs globales
               // ou au point + grandeurs globales éventuelles
               const string & nomF_charge = tabPonctuel_i.NomF_charge();
               if(nomF_charge.length())
                 {Fonction_nD * pt_fonct = lesFonctionsnD->Trouve(nomF_charge);
                  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 dimension = la dimension de l'espace
                     Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                     if (pt_fonct->NbComposante() == 1)
                      {force *= tava(1);}
                     else if (pt_fonct->NbComposante() == dima)
                      {switch (dima)
                            { case 3 : force(3) *= tava(3);
                              case 2 : force(2) *= tava(2);
                              case 1 : force(1) *= tava(1);
                            };
                      }
                     else
                      {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                            << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef()
                            << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                            << dima << endl;
                       Sortie(1);
                      };
                   }
                  else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M
                   { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima)
                      { // on vérifie également qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace
                       const Coordonnee& M_tdt = noe.Coord2();
                       // 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_tdt);
                       Tableau <double> & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier
                       if (pt_fonct->NbComposante() == 1)
                        {force *= tava(1);}
                       else if (pt_fonct->NbComposante() == dima)
                        {switch (dima)
                              { case 3 : force(3) *= tava(3);
                                case 2 : force(2) *= tava(2);
                                case 1 : force(1) *= tava(1);
                              };
                        }
                       else
                        {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                              << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef()
                              << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : "
                              << dima << endl;
                         Sortie(1);
                        };
                      }
                     else
                      { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima)
                        {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point"
                              << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef()
                              << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales ";
                         cout << "\n fonction nD: "; pt_fonct->Affiche();
                         Sortie(1);
                        };
                      };
                   }
                 };
              
               if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1)) && (noe.UneVariable(X1)))
                 // les forces ponctuelles correspondent aux ddl X1 X2 ...
                 for (int di=1;di<= dima;di++)
                   {// maintenant deux cas: soit avec ou sans courbes de charge
                    int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                    if(pt_courbe == NULL)
                     {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                          // cas d'une grandeur fixe
                       { vecglob (pt_adres) +=  coeff * (force(di));}
                      else // cas d'une fonction nD sans courbe
                       { vecglob (pt_adres) +=  (force(di))* tabPonctuel_i.Echelle_courbe();}
                     }
                    else
                     { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                       double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPonctuel_i.Echelle_courbe();
                       vecglob (pt_adres) +=  coeff_charge * (force(di));
                     };
                   };
             };
          };
       };
    
    
     // $$$ --- cas des torseurs de forces ponctuelles --- $$$
    
     // on parcours le tableau tabTorseurPonct
     for (int i=1;i<= tabTorseurPonct.Taille();i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<PTorseurPonct> >& tabTorseurPonct_i = tabTorseurPonct(i);
         const ReferenceNE & ref =
                  ((ReferenceNE &) lesRef->Trouve(tabTorseurPonct_i.NomRef(),tabTorseurPonct_i.NomMaillage()));
         // dans le cas d'une courbe de charge globale on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabTorseurPonct_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(tabTorseurPonct_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabTorseurPonct_i.Pas_a_prendre_en_compte(Temps_courant())))
          {int reftaille = ref.Taille();
           int cas_assemblage = (assemb.Nb_cas_assemb()).n; // récup du cas d'assemblage
           // on renseigne les positions des points où vont s'appliquer les forces équivalentes au torseur
           Tableau < Coordonnee >& tab_P_i = tab_P(i); // pour simplier
           Tableau < Coordonnee >& t_force_i = t_force(i); // idem
           for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
             { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
               tab_P_i(nn) = noe.Coord2();
             };
           // on calcul les efforts ponctuelles correspondant au torseur
           tabTorseurPonct_i.Force(tab_P_i,t_force_i);
           // maintenant on calcule l'action de chaque force ponctuelle
           for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
             { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
               Coordonnee& force = t_force_i(nn); // pour simplifier
               // on regarde si jamais on a une dépendance globale à des grandeurs globales
               // ou au point + grandeurs globales éventuelles
               const string & nomF_charge = tabTorseurPonct_i.NomF_charge();
               if(nomF_charge.length())
                 {Fonction_nD * pt_fonct = lesFonctionsnD->Trouve(nomF_charge);
                  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 dimension = la dimension de l'espace
                     Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                     if (pt_fonct->NbComposante() == 1)
                      {force *= tava(1);}
                     else if (pt_fonct->NbComposante() == dima)
                      {switch (dima)
                            { case 3 : force(3) *= tava(3);
                              case 2 : force(2) *= tava(2);
                              case 1 : force(1) *= tava(1);
                            };
                      }
                     else
                      {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                            << " cas des torseurs de forces ponctuelles sur la ref: "
                            << tabTorseurPonct_i.NomRef()
                            << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                            << dima << endl;
                       Sortie(1);
                      };
                   }
                  else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M
                   { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima)
                      { // on vérifie également qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace
                       const Coordonnee& M_tdt = noe.Coord2();
                       // 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_tdt);
                       Tableau <double> & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier
                       if (pt_fonct->NbComposante() == 1)
                        {force *= tava(1);}
                       else if (pt_fonct->NbComposante() == dima)
                        {switch (dima)
                              { case 3 : force(3) *= tava(3);
                                case 2 : force(2) *= tava(2);
                                case 1 : force(1) *= tava(1);
                              };
                        }
                       else
                        {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                              << " cas des torseurs de forces ponctuelles sur la ref: "
                              << tabTorseurPonct_i.NomRef()
                              << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : "
                              << dima << endl;
                         Sortie(1);
                        };
                      }
                     else
                      { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima)
                        {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point"
                              << " cas des forces ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef()
                              << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales ";
                         cout << "\n fonction nD: "; pt_fonct->Affiche();
                         Sortie(1);
                        };
                      };
                   }
                 };
              
               if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1)) && (noe.UneVariable(X1)))
                 // les forces ponctuelles correspondent aux ddl X1 X2 ...
                 for (int di=1;di<= dima;di++)
                   {// maintenant deux cas: soit avec ou sans courbes de charge
                    int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                    if(pt_courbe == NULL)
                     {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                          // cas d'une grandeur fixe
                       { vecglob (pt_adres) +=  coeff * (force(di));}
                      else // cas d'une fonction nD sans courbe
                       { vecglob (pt_adres) +=  (force(di))* tabTorseurPonct_i.Echelle_courbe();}
                     }
                    else
                     { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                       double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabTorseurPonct_i.Echelle_courbe();
                       vecglob (pt_adres) +=  coeff_charge * (force(di));
                     };
                   };
             };
          };
       };

     // --- cas des forces surfaciques  ---$$$
     int tabsurfactaille =  tabFsurfac.Taille();
     if ( tabsurfactaille != 0)
      {// on parcours le tableau tabFsurfac dans le repère absolu
       for (int i=1;i<= tabsurfactaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocForces> >& tabFsurfac_i = tabFsurfac(i);
         const ReferenceAF & ref =
               ((ReferenceAF &) lesRef->Trouve(tabFsurfac_i.NomRef(),tabFsurfac_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabFsurfac_i.NomCourbeCharge() != "")
              pt_courbe = lesCourbes1D->Trouve(tabFsurfac_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabFsurfac_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les surfaces associes
           Coordonnee vforce; // init par défaut
           int nbref = ref.Taille();
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               // récupération de la force de référence
               // on traite en fonction du fait que c'est un champ ou non
               switch (tabFsurfac_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      vforce = (tabFsurfac_i.Vect_de_coordonnee(ns)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      // on récupère les courbes pour les dim composantes
                      { switch (dima)
                         { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                       lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,3));
                                     vforce(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                    }
                           case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                       lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,2));
                                     vforce(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                    }
                           case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                       lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,1));
                                     vforce(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                    }
                         };
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      vforce = (tabFsurfac_i.Vect_de_coordonnee(1)); break;
                 };

               // on regarde si jamais on a une dépendance à des grandeurs globales
               const string & nomF_charge = tabFsurfac_i.NomF_charge();
               Fonction_nD * pt_fonct=NULL;
               if(nomF_charge.length())
                 {pt_fonct = lesFonctionsnD->Trouve(tabFsurfac_i.NomF_charge());
                  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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {vforce *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : vforce(3) *= tava(3);
                             case 2 : vforce(2) *= tava(2);
                             case 1 : vforce(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des forces surfaciques sur la ref: "<< tabFsurfac_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                    // on désactive la fonction nD car elle a déjà été utilisée
                    pt_fonct = NULL;
                   };
                 };
              
               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                    { vforce *=  coeff;} // cas d'une grandeur fixe
                  else // sinon il ne reste plus qu'à intégrer l'échelle
                    { vforce *=  tabFsurfac_i.Echelle_courbe();};
                 }
               else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFsurfac_i.Echelle_courbe();
                   vforce *=  coeff_charge;
                 };
               // appel du calcul du second membre correspondant à la charge surfacique
   //            Coordonnee vforce = tabFsurfac_i.Coord() * coeff;
               Vecteur SM = elem.SM_charge_surfacique_E_tdt(vforce,pt_fonct,ref.NumeroFA(ns),pa);
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
               const ElFrontiere* elfront = NULL; // init
               if (ParaGlob::AxiSymetrie())
                 // dans le cas axiSymétrique, les surfaces sont générées par les lignes,
                {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));}
               else // cas normale
                {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
             };
          };
       };
      };
    
     // --- cas des pressions  ---$$$
     int tabPresUniftaille =  tabPresUnif.Taille();
     if ( tabPresUniftaille != 0)
      {// on parcours le tableau tabPresUnif
       for (int i=1;i<= tabPresUniftaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocIntensite> >& tabPresUnif_i = tabPresUnif(i);
         const ReferenceAF & ref =
              ((ReferenceAF &) lesRef->Trouve(tabPresUnif_i.NomRef(),tabPresUnif_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabPresUnif_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(tabPresUnif_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabPresUnif_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les surfaces associes
           int nbref = ref.Taille();
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               double press_ac = 0.; // init
               // récupération de la pression de référence
               // on traite en fonction du fait que c'est un champ ou non
               switch (tabPresUnif_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      press_ac = (tabPresUnif_i.Val(ns)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      { Courbe1D * pt_courbe_force_ns =
                                       lesCourbes1D->Trouve(tabPresUnif_i.Nom_val(ns));
                        press_ac = pt_courbe_force_ns->Valeur(Temps_courant());
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      press_ac = (tabPresUnif_i.Val(1)); break;
                 };

               // on regarde si jamais on a une dépendance à des grandeurs globales
               const string & nomF_charge = tabPresUnif_i.NomF_charge();
               Fonction_nD * pt_fonct = NULL;
               if(nomF_charge.length())
                 {pt_fonct = lesFonctionsnD->Trouve(tabPresUnif_i.NomF_charge());
                  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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {press_ac *= tava(1);}
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas de pression sur la ref: "<< tabPresUnif_i.NomRef()
                           << " en retour, le nombre de composante est different de 1  : "
                           << dima << endl;
                      Sortie(1);
                     };
                    // on désactive la fonction nD car elle a déjà été utilisée
                    pt_fonct = NULL;
                   };
                 };
              
               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                    { press_ac *=  coeff;} // cas d'une grandeur fixe
                  else // sinon il ne reste plus qu'à intégrer l'échelle
                    { press_ac *=  tabPresUnif_i.Echelle_courbe();};
                 }
               else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPresUnif_i.Echelle_courbe();
                   press_ac *=  coeff_charge;
                 };
               // appel du calcul du second membre correspondant à la charge de type pression
         //      double press_ac = (tabPresUnif_i.Val()) * coeff;
               Vecteur SM = elem.SM_charge_pression_E_tdt(press_ac,pt_fonct,ref.NumeroFA(ns),pa);
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
               const ElFrontiere* elfront = NULL; // init
               if (ParaGlob::AxiSymetrie())
                 // dans le cas axiSymétrique, les surfaces sont générées par les lignes,
                {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));}
               else // cas normale
                {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              };
          };
        }
       };
    
     // --- cas des pression unidirectionnelle ---
     // c'est-à-dire des forces qui sont appliquée sur la surface en gardant la même
     // direction localement pendant toute la transformation
     int PresUniDirtaille =  PresUniDir.Taille();
     if ( PresUniDirtaille != 0)
      {// on parcours le tableau PresUniDir
       for (int i=1;i<= PresUniDirtaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocForces> >& PresUniDir_i = PresUniDir(i);
         const ReferenceAF & ref =
             ((ReferenceAF &) lesRef->Trouve(PresUniDir_i.NomRef(),PresUniDir_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(PresUniDir_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(PresUniDir_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(PresUniDir_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les surfaces associes
           int nbref = ref.Taille();
           Coordonnee press_uni; // init par défaut
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               // récupération de la force de référence
               // on traite en fonction du fait que c'est un champ ou non
               switch (PresUniDir_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      press_uni = (PresUniDir_i.Vect_de_coordonnee(ns)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      // on récupère les courbes pour les dim composantes
                      { switch (dima)
                         { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                       lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,3));
                                     press_uni(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                    }
                           case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                       lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,2));
                                     press_uni(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                    }
                           case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                       lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,1));
                                     press_uni(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                    }
                         };
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      press_uni = (PresUniDir_i.Vect_de_coordonnee(1)); break;
                 };

               // on regarde si jamais on a une dépendance à des grandeurs globales
               const string & nomF_charge = PresUniDir_i.NomF_charge();
               Fonction_nD * pt_fonct=NULL;
               if(nomF_charge.length())
                 {pt_fonct = lesFonctionsnD->Trouve(PresUniDir_i.NomF_charge());
                  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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {press_uni *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : press_uni(3) *= tava(3);
                             case 2 : press_uni(2) *= tava(2);
                             case 1 : press_uni(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des des pression unidirectionnelle sur la ref: "<< PresUniDir_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                    // on désactive la fonction nD car elle a déjà été utilisée
                    pt_fonct = NULL;
                   };
                 };
              
               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                    { press_uni *=  coeff;} // cas d'une grandeur fixe
                  else // sinon il ne reste plus qu'à intégrer l'échelle
                    { press_uni *=  PresUniDir_i.Echelle_courbe();};
                 }
               else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresUniDir_i.Echelle_courbe();
                   press_uni *=  coeff_charge;
                 };
               // appel du calcul du second membre correspondant à la charge de type pression
        //       Coordonnee press_uni = (PresUniDir_i.Coord()) * coeff;
               Vecteur SM = elem.SM_charge_presUniDir_E_tdt(press_uni,pt_fonct,ref.NumeroFA(ns),pa);
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
               const ElFrontiere* elfront = NULL; // init
               if (ParaGlob::AxiSymetrie())
                 // dans le cas axiSymétrique, les surfaces sont générées par les lignes,
                {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));}
               else // cas normale
                {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
             };
          };
       };
      };

     // --- cas des forces lineique  ---$$$
     int tablineiquetaille =  tabFlineique.Taille();
     if ( tablineiquetaille != 0)
      {// on parcours le tableau tabFlineique dans le repère absolu
       for (int i=1;i<= tablineiquetaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocForces> >& tabFlineique_i = tabFlineique(i);
         const ReferenceAF & ref =
             ((ReferenceAF &) lesRef->Trouve(tabFlineique_i.NomRef(),tabFlineique_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabFlineique_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(tabFlineique_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabFlineique_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les aretes associes
           int nbref = ref.Taille();
           Coordonnee f_lin; // init par défaut
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               // récupération de la force de référence
               // on traite en fonction du fait que c'est un champ ou non
               switch (tabFlineique_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      f_lin = (tabFlineique_i.Vect_de_coordonnee(ns)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      // on récupère les courbes pour les dim composantes
                      { switch (dima)
                         { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                       lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,3));
                                     f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                    }
                           case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                       lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,2));
                                     f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                    }
                           case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                       lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,1));
                                     f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                    }
                         };
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      f_lin = (tabFlineique_i.Vect_de_coordonnee(1)); break;
                 };

               // on regarde si jamais on a une dépendance à des grandeurs globales
               const string & nomF_charge = tabFlineique_i.NomF_charge();
               Fonction_nD * pt_fonct = NULL;
               if(nomF_charge.length())
                 {pt_fonct = lesFonctionsnD->Trouve(tabFlineique_i.NomF_charge());
                  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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {f_lin *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : f_lin(3) *= tava(3);
                             case 2 : f_lin(2) *= tava(2);
                             case 1 : f_lin(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des des forces lineiques sur la ref: "<< tabFlineique_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                    // on désactive la fonction nD car elle a déjà été utilisée
                    pt_fonct = NULL;
                   };
                 };

               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                    { f_lin *=  coeff;} // cas d'une grandeur fixe
                  else // sinon il ne reste plus qu'à intégrer l'échelle
                    { f_lin *=  tabFlineique_i.Echelle_courbe();};
                 }
               else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineique_i.Echelle_courbe();
                   f_lin *=  coeff_charge;
                 };
               // appel du calcul du second membre correspondant à la charge lineique
         //      Vecteur f_lin = (tabFlineique_i.Coord()) * coeff ;
               Vecteur SM = elem.SM_charge_lineique_E_tdt(f_lin,pt_fonct,ref.NumeroFA(ns),pa);
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée
               const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arete frontière
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              };
          };
        }
       };
    
     // --- cas des forces lineiques suiveuses  ---$$$
     int tablineiqueSuivtaille =  tabFlineiqueSuiv.Taille();
     if ( tablineiqueSuivtaille != 0)
      {// on parcours le tableau tabFlineiqueSuiv des forces exprimées dans le repère global
       for (int i=1;i<= tablineiqueSuivtaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocForces> > & tabFlineiqueSuiv_i = tabFlineiqueSuiv(i);
         const ReferenceAF & ref =
             ((ReferenceAF &) lesRef->Trouve(tabFlineiqueSuiv_i.NomRef(),tabFlineiqueSuiv_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabFlineiqueSuiv_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabFlineiqueSuiv_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les aretes associes
           int nbref = ref.Taille();
           Coordonnee f_lin; // init par défaut
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               // récupération de la force de référence
               // on traite en fonction du fait que c'est un champ ou non
               switch (tabFlineiqueSuiv_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(ns)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      // on récupère les courbes pour les dim composantes
                      { switch (dima)
                         { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                       lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,3));
                                     f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                    }
                           case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                       lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,2));
                                     f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                    }
                           case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                       lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,1));
                                     f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                    }
                         };
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(1)); break;
                 };

               // on regarde si jamais on a une dépendance à des grandeurs globales
               const string & nomF_charge = tabFlineiqueSuiv_i.NomF_charge();
               Fonction_nD * pt_fonct = NULL;
               if(nomF_charge.length())
                 {pt_fonct = lesFonctionsnD->Trouve(tabFlineiqueSuiv_i.NomF_charge());
                  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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {f_lin *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : f_lin(3) *= tava(3);
                             case 2 : f_lin(2) *= tava(2);
                             case 1 : f_lin(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des des forces lineiques suiveuses sur la ref: "<< tabFlineiqueSuiv_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                    // on désactive la fonction nD car elle a déjà été utilisée
                    pt_fonct = NULL;
                   };
                 };
              
               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                    { f_lin *=  coeff;} // cas d'une grandeur fixe
                  else // sinon il ne reste plus qu'à intégrer l'échelle
                    { f_lin *=  tabFlineiqueSuiv_i.Echelle_courbe();};
                 }
               else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineiqueSuiv_i.Echelle_courbe();
                   f_lin *=  coeff_charge;
                 }
              // appel du calcul du second membre correspondant à la charge lineique
       //        Vecteur f_lin = (tabFlineiqueSuiv_i.Coord()) * coeff ;
               Vecteur SM = elem.SM_charge_lineique_Suiv_E_tdt(f_lin,pt_fonct,ref.NumeroFA(ns),pa);
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée
               const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arete frontière
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              };
          };
        }
       };
    
     // --- cas des forces volumiques  ---$$$
     int tabvoltaille =  tabFvol.Taille();
     if ( tabvoltaille != 0)
      {// on parcours le tableau tabFvol dans le repère absolu
       for (int i=1;i<= tabvoltaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<BlocForces> >& tabFvol_i = tabFvol(i);
         const ReferenceNE & ref =
              ((ReferenceNE &) lesRef->Trouve(tabFvol_i.NomRef(),tabFvol_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(tabFvol_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(tabFvol_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(tabFvol_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments
           int nbref = ref.Taille();
           Coordonnee f_vol; // init par défaut
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.Numero(ns));
               // récupération de la force de référence
               // on traite en fonction du fait que c'est un champ ou non
               switch (tabFvol_i.Champ())
                 { case 1: // c'est un champ de valeurs
                      f_vol = (tabFvol_i.Vect_de_coordonnee(ns)); break;
                   case 2: // c'est un champ de fonctions
                      // pour l'instant, le seul cas est celui des fonctions du temps
                      // on récupère les courbes pour les dim composantes
                      { switch (dima)
                         { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                       lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,3));
                                     f_vol(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                    }
                           case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                       lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,2));
                                     f_vol(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                    }
                           case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                       lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,1));
                                     f_vol(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                    }
                         };
                        break;
                      }
                   case 0: default:// ce n'est pas un champ
                      f_vol = (tabFvol_i.Vect_de_coordonnee(1)); break;
                 };

               // on regarde si jamais on a une dépendance à des grandeurs globales
               const string & nomF_charge = tabFvol_i.NomF_charge();
               Fonction_nD * pt_fonct = NULL;
               if(nomF_charge.length())
                 {pt_fonct = lesFonctionsnD->Trouve(tabFvol_i.NomF_charge());
                  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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {f_vol *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : f_vol(3) *= tava(3);
                             case 2 : f_vol(2) *= tava(2);
                             case 1 : f_vol(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des des forces lineiques sur la ref: "<< tabFvol_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                    // on désactive la fonction nD car elle a déjà été utilisée
                    pt_fonct = NULL;
                   };
                 };
              
               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                    { f_vol *=  coeff;} // cas d'une grandeur fixe
                  else // sinon il ne reste plus qu'à intégrer l'échelle
                    { f_vol *=  tabFvol_i.Echelle_courbe();};
                 }
               else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFvol_i.Echelle_courbe();
                   f_vol *=  coeff_charge;
                 };
               // on indique si éventuellement on veut se référer au volume initial
               bool volume_finale = !(tabFvol_i.Attribut() == "sur_volume_initial_");
               // appel du calcul du second membre correspondant à la charge volumique
         //      Vecteur f_vol = (tabFvol_i.Coord()) * coeff;
               Vecteur SM = elem.SM_charge_volumique_E_tdt(f_vol,pt_fonct,pa,volume_finale);
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à l'élément
               assemb.AssemSM (vecglob,SM,elem.TableauDdl(),elem.Tab_noeud()); // assemblage
              };
          };
        }
       }
    
     // --- cas des des pressions hydrostatiques  ---$$$
     int PresHydrotaille =  PresHydro.Taille();
     if ( PresHydrotaille != 0)
      {// on parcours le tableau PresHydro
       for (int i=1;i<= PresHydrotaille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<PHydro> > & PresHydro_i = PresHydro(i);
         const ReferenceAF & ref =
              ((ReferenceAF &) lesRef->Trouve(PresHydro_i.NomRef(),PresHydro_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(PresHydro_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(PresHydro_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(PresHydro_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les surfaces associes
           int nbref = ref.Taille();
           // on récupère le vecteur normal et le point de la surface libre
           const Coordonnee& N = PresHydro_i.Direction_N(); // récup de la direction N normée
           const Coordonnee& A = PresHydro_i.Point_M(); // récup du point de la surface libre
           double  p_hydro_ref = (PresHydro_i.Val(1));
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               double p_hydro = p_hydro_ref;
               // maintenant deux cas: soit avec ou sans courbes de charge
               if(pt_courbe == NULL)
                 { p_hydro *=  coeff;} // cas d'une grandeur fixe
               else
                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                   double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresHydro_i.Echelle_courbe();
                   p_hydro *=  coeff_charge;
                 }
               // appel du calcul du second membre correspondant à la charge de type pression
               bool sans_limitation=false;
               if (PresHydro_i.Attribut() == "sans_limitation_") sans_limitation= true;
               Vecteur SM = elem.SM_charge_hydrostatique_E_tdt
                        (N,p_hydro,ref.NumeroFA(ns),A,pa,sans_limitation) ;
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
               const ElFrontiere* elfront = NULL; // init
               if (ParaGlob::AxiSymetrie())
                 // dans le cas axiSymétrique, les surfaces sont générées par les lignes,
                {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns));}
               else // cas normale
                {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); };// la surface frontière
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              }
          };
        }
       };
    
     // --- cas des des efforts  hydrodynamiques  ---$$$
     int coefHydroDynataille =  coefHydroDyna.Taille();
     if ( coefHydroDynataille != 0)
      {// on parcours le tableau coefHydroDyna
       for (int i=1;i<= coefHydroDynataille;i++)
       { // recup de la reference correspondant au mot cle
         BlocCharge< BlocDdlLim<PHydrodyna> > & coefHydroDyna_i = coefHydroDyna(i);
         const ReferenceAF & ref =
              ((ReferenceAF &) lesRef->Trouve(coefHydroDyna_i.NomRef(),coefHydroDyna_i.NomMaillage()));
         // dans le cas d'une courbe de charge on récupère l'adresse
         Courbe1D * pt_courbe = NULL;
         if(coefHydroDyna_i.NomCourbeCharge() != "")
            pt_courbe = lesCourbes1D->Trouve(coefHydroDyna_i.NomCourbeCharge());
         // pas de calcul, si le temps est inférieur au temps actif
         // ou si d'une part le temps n'est pas actif et qu'au pas précédent
         // il n'était pas également actif on ne fait rien
         if (!(coefHydroDyna_i.Pas_a_prendre_en_compte(Temps_courant())))
          {// Maintenant les éléments et les noeuds ou arretes ou surfaces associes suivant la dimension
           int nbref = ref.Taille();
           // récup du poids volumique poids volumique = ro * g, qui n'est pas soumis à la courbe de charge
           double  poidvol = (coefHydroDyna_i.Val(1));
           for (int ns=1;ns<= nbref;ns++)
             { // récupération de l'élément fini
               Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
               // maintenant deux cas: soit avec ou sans courbes de charge
               double  coeff_charge=1.;
               if(pt_courbe != NULL)
                 { coeff_charge = (pt_courbe->Valeur(Temps_courant())) * coefHydroDyna_i.Echelle_courbe();};
               // appel du calcul du second membre correspondant aux efforts hydrodynamique
               Vecteur SM = elem.SM_charge_hydrodynamique_E_tdt(coefHydroDyna_i.Frot_fluid(),poidvol
                                             ,coefHydroDyna_i.Coef_aero_n(),ref.NumeroFA(ns),coeff_charge
                                             ,coefHydroDyna_i.Coef_aero_t(),pa) ;
               // assemblage du second membre
               // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
               const ElFrontiere* elfront;
               switch (ParaGlob::Dimension())
                { case 3: {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); break;};// la surface frontière
                 case 2: {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); break;};// la ligne frontière
                 case 1: {elfront = elem.Frontiere_points(ref.NumeroFA(ns)); break;};// le point frontière
                 };
               assemb.AssemSM (vecglob,SM,elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              };
          };
        };
       }; //-- fin de hydrodyna
    
     // affichage éventuelle du second membre
     if (ParaGlob::NiveauImpression() >= 10)
       { string entete = " affichage du second membre apres chargement";
         vecglob.Affichage_ecran(entete);
       };
    }
  catch (ErrSortieFinale)
       // cas d'une direction voulue vers la sortie
       // on relance l'interuption pour le niveau supérieur
     { ErrSortieFinale toto;
       throw (toto);
     }
  catch ( ... )
   { if (ParaGlob::NiveauImpression() >= 1)
       {cout << "\n\n  ****!!**  erreur dans le calcul du chargement impose ***";
        if (ParaGlob::NiveauImpression() >= 4)
              cout <<  "\n Charge::ChargeSecondMembre_Ex_mecaSolid(..";
        cout << flush;
       };
      retour = false;
    };
 
    temps_cpu_chargement.Arret_du_comptage();
    return retour;
 };
     
    
// cas implicite avec également modification de la raideur quand le chargement dépend
// de la géométrie, en fonction des paramètres de pa 
// la variable relative à la classe assemblage est passé en argument ce qui permet d'utiliser
// éventuellement plusieurs type d'assemblage
// de plus il faut prévoir un pointeur d'assemblage pour séparer le cas d'assemblage symétrique ou non
// d'où la variable pointeur de fonction membre d'assemblage
// retour false si problème de calcul
bool Charge::ChargeSMembreRaideur_Im_mecaSolid
                                (Assemblage & assemb,LesMaillages * lesMail,LesReferences* lesRef
                                 ,Vecteur& vecglob, Mat_abstraite& matglob
                                 ,void (Assemblage::* assembMat)  // le pointeur de fonction
                                       (Mat_abstraite & matglob,const Mat_abstraite & matloc,
                                        const DdlElement& tab_ddl,const Tableau<Noeud *>&tab_noeud)
                                 ,const ParaAlgoControle & pa
                                 ,LesCourbes1D* lesCourbes1D,LesFonctions_nD*  lesFonctionsnD)
 {
#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 true;
#endif
  temps_cpu_chargement.Mise_en_route_du_comptage();
  bool retour = true; // init par défaut
  try
   {// pour tous les maillages, pour tous les éléments on effectue une initialisation éventuelle
    int nb_mail = lesMail->NbMaillage();
    for (int imail = 1;imail<=nb_mail;imail++)
     { int nb_ele = lesMail->Nombre_element(imail);
       for (int iele =1;iele<=nb_ele;iele++)
         lesMail->Element_LesMaille(imail,iele).Initialisation_avant_chargement();
     };
  
    // récup
    bool avec_raid = pa.Var_charge_externe();
    int cas_assemblage = (assemb.Nb_cas_assemb()).n; // récup du cas d'assemblage
    int posi = Id_nom_ddl("X1") -1;
    int dim = ParaGlob::Dimension();

    // $$$--- cas des forces ponctuelles ---$$$
    // on regarde si l'espace de travail est axi-symétrique, si oui on utilise une dimension réduite
    int dima=dim; // c'est dima qui est utiliser pour la mise en place des efforts ponctuelles
    if (ParaGlob::AxiSymetrie()) dima--;
    // on parcours le tableau tabPonctuel
    // ici il n'y a pas d'implication sur la raideur
    for (int i=1;i<= tabPonctuel.Taille();i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocForces> >& tabPonctuel_i = tabPonctuel(i);
//tabPonctuel_i.Affiche();cout << endl;
        const ReferenceNE & ref = 
              ((ReferenceNE &) lesRef->Trouve(tabPonctuel_i.NomRef(),tabPonctuel_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabPonctuel_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabPonctuel_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
        // ou si d'une part le temps n'est pas actif et qu'au pas précédent
        // il n'était pas également actif on ne fait rien
        if (!(tabPonctuel_i.Pas_a_prendre_en_compte(Temps_courant())))
         {int reftaille = ref.Taille();
          for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
            { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
              Coordonnee force; // init
              // on traite en fonction du fait que c'est un champ ou non
              switch (tabPonctuel_i.Champ())
                { case 1: // c'est un champ de valeurs
                     force = (tabPonctuel_i.Vect_de_coordonnee(nn)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     // on récupère les courbes pour les dim composantes
                     { switch (dima)
                        { case 3 : {Courbe1D * pt_courbe_force_nn_3 =
                                      lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,3));
                                    force(3) = pt_courbe_force_nn_3->Valeur(Temps_courant());
                                   }
                          case 2 : {Courbe1D * pt_courbe_force_nn_2 =
                                      lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,2));
                                    force(2) = pt_courbe_force_nn_2->Valeur(Temps_courant());
                                   }
                          case 1 : {Courbe1D * pt_courbe_force_nn_1 =
                                      lesCourbes1D->Trouve(tabPonctuel_i.Nom_vect(nn,1));
                                    force(1) = pt_courbe_force_nn_1->Valeur(Temps_courant());
                                   }
                        };
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     force = (tabPonctuel_i.Vect_de_coordonnee(1)); break;
                };
              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabPonctuel_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(nomF_charge);
                 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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {force *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : force(3) *= tava(3);
                             case 2 : force(2) *= tava(2);
                             case 1 : force(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                  }
                 else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M
                  { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima)
                     { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace
                       const Coordonnee& M_tdt = noe.Coord2();
                       Tableau <Coordonnee> tab_M(1,M_tdt);
                       Tableau <double> & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier
                       if (pt_fonct->NbComposante() == 1)
                        {force *= tava(1);}
                       else if (pt_fonct->NbComposante() == dima)
                        {switch (dima)
                              { case 3 : force(3) *= tava(3);
                                case 2 : force(2) *= tava(2);
                                case 1 : force(1) *= tava(1);
////----- debug ---
//cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid ** ";
//cout << "\n tava= " << tava << endl;
//
////--- fin debug ---
                             
                              };
                        }
                       else
                        {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                              << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef()
                              << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : "
                              << dima << endl;
                         Sortie(1);
                        };
                      }
                    else
                     { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima)
                       {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point"
                             << " cas des forces ponctuelles sur la ref: "<< tabPonctuel_i.NomRef()
                             << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales ";
                        cout << "\n fonction nD: "; pt_fonct->Affiche();
                        Sortie(1);
                       };
                     };
                  };
                };
              // les efforts s'appliques sur les ddl de position, il faut donc que X1 soit actif
              // pour le cas d'assemblage considéré et en activité
              if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1))
                    && (noe.UneVariable(X1)))
                // les forces ponctuelles correspondent aux ddl 1 2 ...
                for (int di=1;di<= dima;di++)
                  {// maintenant deux cas: soit avec ou sans courbes de charge
                   int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                   if (pt_courbe == NULL)
                    {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                        // cas d'une grandeur fixe
                      { vecglob (pt_adres) +=  coeff * (force(di));}
                     else // cas d'une fonction nD sans courbe
                      { vecglob (pt_adres) +=  (force(di))* tabPonctuel_i.Echelle_courbe();}
                    }
                   else // cas avec courbe, avec ou sans fonction c'est pareil
                    { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                      double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPonctuel_i.Echelle_courbe();
                      vecglob (pt_adres) +=  coeff_charge * (force(di));
                    };
                  };
            };
         };
      };
      

    // $$$--- cas des torseurs de forces ponctuelles ---$$$
    for (int i=1;i<= tabTorseurPonct.Taille();i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<PTorseurPonct> >& tabTorseurPonct_i = tabTorseurPonct(i);
//tabTorseurPonct_i.Affiche();cout << endl;
        const ReferenceNE & ref =
              ((ReferenceNE &) lesRef->Trouve(tabTorseurPonct_i.NomRef(),tabTorseurPonct_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabTorseurPonct_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabTorseurPonct_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
        // ou si d'une part le temps n'est pas actif et qu'au pas précédent
        // il n'était pas également actif on ne fait rien
        if (!(tabTorseurPonct_i.Pas_a_prendre_en_compte(Temps_courant())))
         {int reftaille = ref.Taille();
          // on renseigne les positions des points où vont s'appliquer les forces équivalentes au torseur
          Tableau < Coordonnee >& tab_P_i = tab_P(i); // pour simplier
          Tableau < Coordonnee >& t_force_i = t_force(i); // idem
          for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
            { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
              tab_P_i(nn) = noe.Coord2();
            };
          // on calcul les efforts ponctuelles correspondant au torseur
          tabTorseurPonct_i.Force(tab_P_i,t_force_i);
          // maintenant on calcule l'action de chaque force ponctuelle
          for (int nn =1; nn<= reftaille;nn++)  // nn = position des  num de noeud
            { Noeud & noe=lesMail->Noeud_LesMaille(ref.Nbmaille(),ref.Numero(nn));
              Coordonnee& force = t_force_i(nn); // pour simplifier
              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabTorseurPonct_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(nomF_charge);
                 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 dimension = la dimension de l'espace
                    Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                    if (pt_fonct->NbComposante() == 1)
                     {force *= tava(1);}
                    else if (pt_fonct->NbComposante() == dima)
                     {switch (dima)
                           { case 3 : force(3) *= tava(3);
                             case 2 : force(2) *= tava(2);
                             case 1 : force(1) *= tava(1);
                           };
                     }
                    else
                     {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                           << " cas des torseurs de forces ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef()
                           << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                           << dima << endl;
                      Sortie(1);
                     };
                  }
                 else if (pt_fonct->Depend_M() != 0) // signifie qu'il y a une dépendance au point M
                  { if ((pt_fonct->Depend_M()<0) && pt_fonct->Nom_variables().Taille() == dima)
                     { // on vérifie qu'en retour on a un scalaire ou un vecteur de dimension = la dimension de l'espace
                       const Coordonnee& M_tdt = noe.Coord2();
                       Tableau <Coordonnee> tab_M(1,M_tdt);
                       Tableau <double> & tava = pt_fonct->Val_FnD_Evoluee(NULL,&tab_M,NULL); // pour simplifier
                       if (pt_fonct->NbComposante() == 1)
                        {force *= tava(1);}
                       else if (pt_fonct->NbComposante() == dima)
                        {switch (dima)
                              { case 3 : force(3) *= tava(3);
                                case 2 : force(2) *= tava(2);
                                case 1 : force(1) *= tava(1);
////----- debug ---
//cout << "\n debug Charge::ChargeSMembreRaideur_Im_mecaSolid ** ";
//cout << "\n tava= " << tava << endl;
//
////--- fin debug ---
                             
                              };
                        }
                       else
                        {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                              << " cas des torseurs de forces ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef()
                              << " en retour, le nombre de composante est different de 1 ou de la dimension de l'espace : "
                              << dima << endl;
                         Sortie(1);
                        };
                      }
                    else
                     { if ((pt_fonct->Depend_M()>0)|| pt_fonct->Nom_variables().Taille() != dima)
                       {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD utilisant les coordonnees du point"
                             << " cas des des torseurs de  ponctuelles sur la ref: "<< tabTorseurPonct_i.NomRef()
                             << " actuellement la fonction doit avoir uniquement comme parametres, les coordonnees du point et eventuellement des variables globales ";
                        cout << "\n fonction nD: "; pt_fonct->Affiche();
                        Sortie(1);
                       };
                     };
                  };
                };
              // les efforts s'appliques sur les ddl de position, il faut donc que X1 soit actif
              // pour le cas d'assemblage considéré et en activité
              if ((noe.Pointeur_assemblage(X1,cas_assemblage) != -1) && (noe.En_service(X1))
                    && (noe.UneVariable(X1)))
                // les forces ponctuelles correspondent aux ddl 1 2 ...
                for (int di=1;di<= dima;di++)
                  {// maintenant deux cas: soit avec ou sans courbes de charge
                   int pt_adres = noe.Pointeur_assemblage(Enum_ddl(posi+di),cas_assemblage);
                   if (pt_courbe == NULL)
                    {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                        // cas d'une grandeur fixe
                      { vecglob (pt_adres) +=  coeff * (force(di));}
                     else // cas d'une fonction nD sans courbe
                      { vecglob (pt_adres) +=  (force(di))* tabTorseurPonct_i.Echelle_courbe();}
                    }
                   else // cas avec courbe, avec ou sans fonction c'est pareil
                    { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                      double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabTorseurPonct_i.Echelle_courbe();
                      vecglob (pt_adres) +=  coeff_charge * (force(di));
                    };
                  };
            };
         };
      };

    // $$$--- cas des forces surfaciques ---$$$
    // il y a une implication éventuelle sur la raideur
    int tabsurfactaille =  tabFsurfac.Taille(); 
    if ( tabsurfactaille != 0)      
     {// on parcours le tableau tabFsurfac dans le repère absolu
      for (int i=1;i<= tabsurfactaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocForces> >& tabFsurfac_i = tabFsurfac(i);
        const ReferenceAF & ref = 
               ((ReferenceAF &) lesRef->Trouve(tabFsurfac_i.NomRef(),tabFsurfac_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabFsurfac_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabFsurfac_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
		      // il n'était pas également actif on ne fait rien
        if (!(tabFsurfac_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments et les surfaces associes
          int nbref = ref.Taille();            
          Coordonnee vforce; // init par défaut
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              // récupération de la force de référence
              // on traite en fonction du fait que c'est un champ ou non
              switch (tabFsurfac_i.Champ())
                { case 1: // c'est un champ de valeurs
                     vforce = (tabFsurfac_i.Vect_de_coordonnee(ns)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     // on récupère les courbes pour les dim composantes
                     { switch (dima)
                        { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                      lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,3));
                                    vforce(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                   }
                          case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                      lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,2));
                                    vforce(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                   }
                          case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                      lesCourbes1D->Trouve(tabFsurfac_i.Nom_vect(ns,1));
                                    vforce(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                   }
                        };
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     vforce = (tabFsurfac_i.Vect_de_coordonnee(1)); break;
                };

              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabFsurfac_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(tabFsurfac_i.NomF_charge());
                 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 dimension = la dimension de l'espace
                   Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                   if (pt_fonct->NbComposante() == 1)
                    {vforce *= tava(1);}
                   else if (pt_fonct->NbComposante() == dima)
                    {switch (dima)
                          { case 3 : vforce(3) *= tava(3);
                            case 2 : vforce(2) *= tava(2);
                            case 1 : vforce(1) *= tava(1);
                          };
                    }
                   else
                    {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                          << " cas des forces surfaciques sur la ref: "<< tabFsurfac_i.NomRef()
                          << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                          << dima << endl;
                     Sortie(1);
                    };
                   // on désactive la fonction nD car elle a déjà été utilisée
                   pt_fonct = NULL;
                  };
                };
             
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                   { vforce *=  coeff;} // cas d'une grandeur fixe
                 else // sinon il ne reste plus qu'à intégrer l'échelle
                   { vforce *=  tabFsurfac_i.Echelle_courbe();};
                }
              else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFsurfac_i.Echelle_courbe();
                  vforce *=  coeff_charge;     
                };
              // -- appel du calcul du second membre et de la raideur éventuelle 
              // correspondant à la charge surfacique
              Element::ResRaid resu  = elem.SMR_charge_surfacique_I(vforce,pt_fonct,ref.NumeroFA(ns),pa);
              //                ((tabFsurfac_i.Coord())*coeff ,ref.NumeroFA(ns),pa);
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
              const ElFrontiere* elfront = NULL;
              if (!ParaGlob::AxiSymetrie())
                 elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal
              else  // cas axisymétrique
                 elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage du second membre
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const());
            };
         };
      };
     };
      
    // $$$--- cas des pressions ---$$$
    int tabPresUniftaille =  tabPresUnif.Taille(); 
    if ( tabPresUniftaille != 0)      
     {// on parcours le tableau tabPresUnif 
      for (int i=1;i<= tabPresUniftaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocIntensite> >& tabPresUnif_i = tabPresUnif(i);
        const ReferenceAF & ref = 
             ((ReferenceAF &) lesRef->Trouve(tabPresUnif_i.NomRef(),tabPresUnif_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabPresUnif_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabPresUnif_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
        // il n'était pas également actif on ne fait rien
        if (!(tabPresUnif_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments et les surfaces associes
          int nbref = ref.Taille();
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              double press_ac = 0.; // init
              // récupération de la pression de référence
              // on traite en fonction du fait que c'est un champ ou non
              switch (tabPresUnif_i.Champ())
                { case 1: // c'est un champ de valeurs
                     press_ac = (tabPresUnif_i.Val(ns)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     { Courbe1D * pt_courbe_force_ns =
                                      lesCourbes1D->Trouve(tabPresUnif_i.Nom_val(ns));
                       press_ac = pt_courbe_force_ns->Valeur(Temps_courant());
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     press_ac = (tabPresUnif_i.Val(1)); break;
                };

              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabPresUnif_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(tabPresUnif_i.NomF_charge());
                 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 dimension = la dimension de l'espace
                   Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                   if (pt_fonct->NbComposante() == 1)
                    {press_ac *= tava(1);}
                   else
                    {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                          << " cas de pression sur la ref: "<< tabPresUnif_i.NomRef()
                          << " en retour, le nombre de composante est different de 1  : "
                          << dima << endl;
                     Sortie(1);
                    };
                   // on désactive la fonction nD car elle a déjà été utilisée
                   pt_fonct = NULL;
                  };
                };
             
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                   { press_ac *=  coeff;} // cas d'une grandeur fixe
                 else // sinon il ne reste plus qu'à intégrer l'échelle
                   { press_ac *=  tabPresUnif_i.Echelle_courbe();};
                }
              else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                { // cas d'une courbe de charge calcul de la valeur du coeff de                 { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabPresUnif_i.Echelle_courbe();
                  press_ac *=  coeff_charge;     
                }
              // appel du calcul du second membre et de la raideur éventuelle 
              // correspondant à la charge de type pression
              Element::ResRaid resu  = elem.SMR_charge_pression_I(press_ac,pt_fonct,ref.NumeroFA(ns),pa);
                      //     (coeff *(tabPresUnif_i.Val()),ref.NumeroFA(ns),pa);
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
              const ElFrontiere* elfront = NULL;
              if (!ParaGlob::AxiSymetrie())
                 elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal
              else  // cas axisymétrique
                 elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); 
             }
         };
       }
      }
      
      
    // $$$---  cas des pression unidirectionnelle ---$$$
    // c'est-à-dire des forces qui sont appliquée sur la surface en gardant la même
    // direction localement pendant toute la transformation
    int PresUniDirtaille =  PresUniDir.Taille(); 
    if ( PresUniDirtaille != 0)      
     {// on parcours le tableau PresUniDir 
      for (int i=1;i<= PresUniDirtaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocForces> >& PresUniDir_i = PresUniDir(i);
        const ReferenceAF & ref = 
            ((ReferenceAF &) lesRef->Trouve(PresUniDir_i.NomRef(),PresUniDir_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(PresUniDir_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(PresUniDir_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
		      // il n'était pas également actif on ne fait rien
        if (!(PresUniDir_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments et les surfaces associes
          int nbref = ref.Taille();
          Coordonnee press_uni; // init par défaut
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              // récupération de la force de référence
              // on traite en fonction du fait que c'est un champ ou non
              switch (PresUniDir_i.Champ())
                { case 1: // c'est un champ de valeurs
                     press_uni = (PresUniDir_i.Vect_de_coordonnee(ns)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     // on récupère les courbes pour les dim composantes
                     { switch (dima)
                        { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                      lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,3));
                                    press_uni(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                   }
                          case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                      lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,2));
                                    press_uni(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                   }
                          case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                      lesCourbes1D->Trouve(PresUniDir_i.Nom_vect(ns,1));
                                    press_uni(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                   }
                        };
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     press_uni = (PresUniDir_i.Vect_de_coordonnee(1)); break;
                };

              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = PresUniDir_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(PresUniDir_i.NomF_charge());
                 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 dimension = la dimension de l'espace
                   Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                   if (pt_fonct->NbComposante() == 1)
                    {press_uni *= tava(1);}
                   else if (pt_fonct->NbComposante() == dima)
                    {switch (dima)
                          { case 3 : press_uni(3) *= tava(3);
                            case 2 : press_uni(2) *= tava(2);
                            case 1 : press_uni(1) *= tava(1);
                          };
                    }
                   else
                    {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                          << " cas des des pression unidirectionnelle sur la ref: "<< PresUniDir_i.NomRef()
                          << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                          << dima << endl;
                     Sortie(1);
                    };
                   // on désactive la fonction nD car elle a déjà été utilisée
                   pt_fonct = NULL;
                  };
                };
             
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                { press_uni *=  coeff;} // cas d'une grandeur fixe
                 else // sinon il ne reste plus qu'à intégrer l'échelle
                   { press_uni *=  PresUniDir_i.Echelle_courbe();};
                }
              else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresUniDir_i.Echelle_courbe();
                  press_uni *=  coeff_charge;     
                }
              // appel du calcul du second membre et de la raideur éventuelle 
              // correspondant à la charge de type pression
        //      Vecteur pres_u = (PresUniDir_i.Coord()) * coeff;
              Element::ResRaid resu  = elem.SMR_charge_presUniDir_I(press_uni,pt_fonct,ref.NumeroFA(ns),pa);
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
              const ElFrontiere* elfront = NULL;
              if (!ParaGlob::AxiSymetrie())
                 elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal
              else  // cas axisymétrique
                 elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); 
             };
         };
       }
      }

    // $$$--- cas des forces lineique ---$$$
    int tablineiquetaille =  tabFlineique.Taille(); 
    if ( tablineiquetaille != 0)      
     {// on parcours le tableau tabFlineique dans le repère absolu
      for (int i=1;i<= tablineiquetaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocForces> >& tabFlineique_i = tabFlineique(i);
        const ReferenceAF & ref = 
            ((ReferenceAF &) lesRef->Trouve(tabFlineique_i.NomRef(),tabFlineique_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabFlineique_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabFlineique_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
		      // il n'était pas également actif on ne fait rien
        if (!(tabFlineique_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments et les aretes associes
          int nbref = ref.Taille();
          Coordonnee f_lin; // init par défaut
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              // récupération de la force de référence
              // on traite en fonction du fait que c'est un champ ou non
              switch (tabFlineique_i.Champ())
                { case 1: // c'est un champ de valeurs
                     f_lin = (tabFlineique_i.Vect_de_coordonnee(ns)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     // on récupère les courbes pour les dim composantes
                     { switch (dima)
                        { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                      lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,3));
                                    f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                   }
                          case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                      lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,2));
                                    f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                   }
                          case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                      lesCourbes1D->Trouve(tabFlineique_i.Nom_vect(ns,1));
                                    f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                   }
                        };
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     f_lin = (tabFlineique_i.Vect_de_coordonnee(1)); break;
                };

              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabFlineique_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(tabFlineique_i.NomF_charge());
                 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 dimension = la dimension de l'espace
                   Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                   if (pt_fonct->NbComposante() == 1)
                    {f_lin *= tava(1);}
                   else if (pt_fonct->NbComposante() == dima)
                    {switch (dima)
                          { case 3 : f_lin(3) *= tava(3);
                            case 2 : f_lin(2) *= tava(2);
                            case 1 : f_lin(1) *= tava(1);
                          };
                    }
                   else
                    {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                          << " cas des des forces lineiques sur la ref: "<< tabFlineique_i.NomRef()
                          << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                          << dima << endl;
                     Sortie(1);
                    };
                   // on désactive la fonction nD car elle a déjà été utilisée
                   pt_fonct = NULL;
                  };
                };
             
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                   { f_lin *=  coeff;} // cas d'une grandeur fixe
                 else // sinon il ne reste plus qu'à intégrer l'échelle
                   { f_lin *=  tabFlineique_i.Echelle_courbe();};
                }
              else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineique_i.Echelle_courbe();
                  f_lin *=  coeff_charge;     
                }
              // appel du calcul du second membre  et de la raideur éventuelle 
              // correspondant à la charge lineique
    //          Vecteur f_lin = (tabFlineique_i.Coord()) * coeff;
              Element::ResRaid resu  = elem.SMR_charge_lineique_I(f_lin,pt_fonct,ref.NumeroFA(ns),pa) ;
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée
              const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arrête frontière
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); 
            };
         };
      };
     };

    // $$$--- cas des forces lineiques suiveuses ---$$$
    int tablineiqueSuivtaille =  tabFlineiqueSuiv.Taille(); 
    if ( tablineiqueSuivtaille != 0)      
     {// on parcours le tableau tabFlineiqueSuiv dans le repère global
      for (int i=1;i<= tablineiqueSuivtaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocForces> > & tabFlineiqueSuiv_i = tabFlineiqueSuiv(i);
        const ReferenceAF & ref = 
           ((ReferenceAF &) lesRef->Trouve(tabFlineiqueSuiv_i.NomRef(),tabFlineiqueSuiv_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabFlineiqueSuiv_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabFlineiqueSuiv_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
		      // il n'était pas également actif on ne fait rien
        if (!(tabFlineiqueSuiv_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments et les aretes associes
          int nbref = ref.Taille();
          Coordonnee f_lin; // init par défaut
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              // récupération de la force de référence
              // on traite en fonction du fait que c'est un champ ou non
              switch (tabFlineiqueSuiv_i.Champ())
                { case 1: // c'est un champ de valeurs
                     f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(ns)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     // on récupère les courbes pour les dim composantes
                     { switch (dima)
                        { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                      lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,3));
                                    f_lin(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                   }
                          case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                      lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,2));
                                    f_lin(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                   }
                          case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                      lesCourbes1D->Trouve(tabFlineiqueSuiv_i.Nom_vect(ns,1));
                                    f_lin(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                   }
                        };
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     f_lin = (tabFlineiqueSuiv_i.Vect_de_coordonnee(1)); break;
                };

              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabFlineiqueSuiv_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(tabFlineiqueSuiv_i.NomF_charge());
                 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 dimension = la dimension de l'espace
                   Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                   if (pt_fonct->NbComposante() == 1)
                    {f_lin *= tava(1);}
                   else if (pt_fonct->NbComposante() == dima)
                    {switch (dima)
                          { case 3 : f_lin(3) *= tava(3);
                            case 2 : f_lin(2) *= tava(2);
                            case 1 : f_lin(1) *= tava(1);
                          };
                    }
                   else
                    {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                          << " cas des des forces lineiques suiveuses sur la ref: "<< tabFlineiqueSuiv_i.NomRef()
                          << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                          << dima << endl;
                     Sortie(1);
                    };
                   // on désactive la fonction nD car elle a déjà été utilisée
                   pt_fonct = NULL;
                  };
                };
             
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                   { f_lin *=  coeff;} // cas d'une grandeur fixe
                 else // sinon il ne reste plus qu'à intégrer l'échelle
                   { f_lin *=  tabFlineiqueSuiv_i.Echelle_courbe();};
                }
              else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFlineiqueSuiv_i.Echelle_courbe();
                  f_lin *=  coeff_charge;     
                }
              // appel du calcul du second membre  et de la raideur éventuelle 
              // correspondant à la charge lineique
    //          Vecteur f_lin = (tabFlineiqueSuiv_i.Coord()) * coeff;
              Element::ResRaid resu  = elem.SMR_charge_lineique_Suiv_I(f_lin,pt_fonct,ref.NumeroFA(ns),pa) ;
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la ligne chargée
              const ElFrontiere* elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // l'arrête frontière
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); 
            };
         };
      };
     };
      
    // $$$---  cas des forces volumiques ---$$$
    int tabvoltaille =  tabFvol.Taille(); 
    if ( tabvoltaille != 0)      
     {// on parcours le tableau tabFvol dans le repère absolu
      for (int i=1;i<= tabvoltaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<BlocForces> >& tabFvol_i = tabFvol(i);
        const ReferenceNE & ref = 
            ((ReferenceNE &) lesRef->Trouve(tabFvol_i.NomRef(),tabFvol_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(tabFvol_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(tabFvol_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
		      // il n'était pas également actif on ne fait rien
        if (!(tabFvol_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments
          int nbref = ref.Taille();
          Coordonnee f_vol; // init par défaut
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.Numero(ns));
              // récupération de la force de référence
              // on traite en fonction du fait que c'est un champ ou non
              switch (tabFvol_i.Champ())
                { case 1: // c'est un champ de valeurs
                     f_vol = (tabFvol_i.Vect_de_coordonnee(ns)); break;
                  case 2: // c'est un champ de fonctions
                     // pour l'instant, le seul cas est celui des fonctions du temps
                     // on récupère les courbes pour les dim composantes
                     { switch (dima)
                        { case 3 : {Courbe1D * pt_courbe_force_ns_3 =
                                      lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,3));
                                    f_vol(3) = pt_courbe_force_ns_3->Valeur(Temps_courant());
                                   }
                          case 2 : {Courbe1D * pt_courbe_force_ns_2 =
                                      lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,2));
                                    f_vol(2) = pt_courbe_force_ns_2->Valeur(Temps_courant());
                                   }
                          case 1 : {Courbe1D * pt_courbe_force_ns_1 =
                                      lesCourbes1D->Trouve(tabFvol_i.Nom_vect(ns,1));
                                    f_vol(1) = pt_courbe_force_ns_1->Valeur(Temps_courant());
                                   }
                        };
                       break;
                     }
                  case 0: default:// ce n'est pas un champ
                     f_vol = (tabFvol_i.Vect_de_coordonnee(1)); break;
                };

              // on regarde si jamais on a une dépendance à des grandeurs globales
              const string & nomF_charge = tabFvol_i.NomF_charge();
              Fonction_nD * pt_fonct = NULL;
              if(nomF_charge.length())
                {pt_fonct = lesFonctionsnD->Trouve(tabFvol_i.NomF_charge());
                 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 dimension = la dimension de l'espace
                   Tableau <double> & tava = pt_fonct->Valeur_pour_variables_globales(); // pour simplifier
                   if (pt_fonct->NbComposante() == 1)
                    {f_vol *= tava(1);}
                   else if (pt_fonct->NbComposante() == dima)
                    {switch (dima)
                          { case 3 : f_vol(3) *= tava(3);
                            case 2 : f_vol(2) *= tava(2);
                            case 1 : f_vol(1) *= tava(1);
                          };
                    }
                   else
                    {cout << "\n *** erreur en retour de l'utilisation d'une fonction nD avec grandeur(s) globale(s)"
                          << " cas des des forces lineiques suiveuses sur la ref: "<< tabFvol_i.NomRef()
                          << " en retour, le nombre de composante est different de 1 et la dimension de l'espace : "
                          << dima << endl;
                     Sortie(1);
                    };
                   // on désactive la fonction nD car elle a déjà été utilisée
                   pt_fonct = NULL;
                  };
                };
             
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                {if (!nomF_charge.length()) // pas de courbe et pas de fonction nD
                   { f_vol *=  coeff;} // cas d'une grandeur fixe
                 else // sinon il ne reste plus qu'à intégrer l'échelle
                   { f_vol *=  tabFvol_i.Echelle_courbe();};
                }
              else //sinon c'est idem avec ou sans fonction nD car c'est déjà intégré
                { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * tabFvol_i.Echelle_courbe();
                  f_vol *=  coeff_charge;     
                }
              // appel du calcul du second membre  et de la raideur éventuelle 
              // correspondant à la charge volumique
              // on indique si éventuellement on veut se référer au volume initial
              bool volume_finale = !(tabFvol_i.Attribut() == "sur_volume_initial_");
          //    Vecteur f_vol = (tabFvol_i.Coord()) * coeff;
              Element::ResRaid resu  = elem.SMR_charge_volumique_I(f_vol,pt_fonct,pa,volume_finale);
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à l'élément
              assemb.AssemSM (vecglob,*(resu.res),elem.TableauDdl(),elem.Tab_noeud()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elem.TableauDdl(),elem.Tab_noeud()); 
             };
         };
       }
      }; 
      
    // $$$---  cas des des pressions hydrostatiques ---$$$
    int PresHydrotaille =  PresHydro.Taille(); 
    if ( PresHydrotaille != 0)      
     {// on parcours le tableau PresHydro 
      for (int i=1;i<= PresHydrotaille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<PHydro> > & PresHydro_i = PresHydro(i);
        const ReferenceAF & ref = 
             ((ReferenceAF &) lesRef->Trouve(PresHydro_i.NomRef(),PresHydro_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(PresHydro_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(PresHydro_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
        // ou si d'une part le temps n'est pas actif et qu'au pas précédent
        // il n'était pas également actif on ne fait rien
        if (!(PresHydro_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// on récupère le vecteur normal et le point de la surface libre
          const Coordonnee& N = PresHydro_i.Direction_N(); // récup de la direction N normée
          const Coordonnee& A = PresHydro_i.Point_M(); // récup du point de la surface libre
          double  p_hydro_ref = (PresHydro_i.Val(1));
          // Maintenant les éléments et les surfaces associes
          int nbref = ref.Taille();
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              double  p_hydro = p_hydro_ref;
              // maintenant deux cas: soit avec ou sans courbes de charge
              if(pt_courbe == NULL)
                { p_hydro *=  coeff;} // cas d'une grandeur fixe
              else
                { // cas d'une courbe de charge calcul de la valeur du coeff de charge
                  double coeff_charge = (pt_courbe->Valeur(Temps_courant())) * PresHydro_i.Echelle_courbe();
                  p_hydro *=  coeff_charge;     
                }
              // appel du calcul du second membre et de la raideur éventuelle 
              // correspondant à la charge de type pression
              bool sans_limitation=false;
              if (PresHydro_i.Attribut() == "sans_limitation_") sans_limitation= true;
              Element::ResRaid resu  = elem.SMR_charge_hydrostatique_I
                       (N,p_hydro,ref.NumeroFA(ns),A,pa,sans_limitation) ;
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
              const ElFrontiere* elfront = NULL;
              if (!ParaGlob::AxiSymetrie())
                 elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); // la surface frontière: cas normal
              else  // cas axisymétrique
                 elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); // la surface arête générant une surface en rotation
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); 
             };
         };
       }
      };
      
   // $$$---  cas des des efforts  hydrodynamiques ---$$$
    int coefHydroDynataille =  coefHydroDyna.Taille(); 
    if ( coefHydroDynataille != 0)      
     {// on parcours le tableau coefHydroDyna 
      for (int i=1;i<= coefHydroDynataille;i++)
      { // recup de la reference correspondant au mot cle
        BlocCharge< BlocDdlLim<PHydrodyna> > & coefHydroDyna_i = coefHydroDyna(i);
        const ReferenceAF & ref = 
             ((ReferenceAF &) lesRef->Trouve(coefHydroDyna_i.NomRef(),coefHydroDyna_i.NomMaillage()));
        // dans le cas d'une courbe de charge on récupère l'adresse
        Courbe1D * pt_courbe = NULL;
        if(coefHydroDyna_i.NomCourbeCharge() != "")
           pt_courbe = lesCourbes1D->Trouve(coefHydroDyna_i.NomCourbeCharge());
        // pas de calcul, si le temps est inférieur au temps actif
		      // ou si d'une part le temps n'est pas actif et qu'au pas précédent
		      // il n'était pas également actif on ne fait rien
        if (!(coefHydroDyna_i.Pas_a_prendre_en_compte(Temps_courant())))
         {// Maintenant les éléments et les surfaces associes
          int nbref = ref.Taille();
          // récup du poids volumique poids volumique = ro * g, qui n'est pas soumis à la courbe de charge     
          double  poidvol = (coefHydroDyna_i.Val(1));  
          for (int ns=1;ns<= nbref;ns++)
            { // récupération de l'élément fini
              Element& elem = lesMail->Element_LesMaille(ref.Nbmaille(), ref.NumeroElem(ns));
              // maintenant deux cas: soit avec ou sans courbes de charge
              double  coeff_charge=1.;
              if(pt_courbe != NULL)
                { coeff_charge = (pt_courbe->Valeur(Temps_courant())) * coefHydroDyna_i.Echelle_courbe();};
              // appel du calcul du second membre correspondant aux efforts hydrodynamique
              Element::ResRaid resu = elem.SMR_charge_hydrodynamique_I(coefHydroDyna_i.Frot_fluid(),poidvol
                                        ,coefHydroDyna_i.Coef_aero_n(),ref.NumeroFA(ns),coeff_charge
                                        ,coefHydroDyna_i.Coef_aero_t(),pa) ;
              // assemblage du second membre
              // il faut utiliser les noeuds et les ddlelement correspondant à la face chargée
              const ElFrontiere* elfront;
              switch (ParaGlob::Dimension())
               { case 3: {elfront = elem.Frontiere_surfacique(ref.NumeroFA(ns)); break;};// la surface frontière
                case 2: {elfront = elem.Frontiere_lineique(ref.NumeroFA(ns)); break;};// la ligne frontière
                case 1: {elfront = elem.Frontiere_points(ref.NumeroFA(ns)); break;};// le point frontière
                };
              assemb.AssemSM (vecglob,*(resu.res),elfront->DdlElem_const(),elfront->TabNoeud_const()); // assemblage
              // assemblage de la raideur si nécessaire
              if (avec_raid)
                 // appel de la fonction adoc, via le pointeur de fonction
                   (assemb.*assembMat) (matglob,*(resu.raid),elfront->DdlElem_const(),elfront->TabNoeud_const()); 
             };
         };
       };
      }; //-- fin de hydrodyna  

      
    // affichage éventuelle de la matrice de raideur et du second membre     
    if (ParaGlob::NiveauImpression() >= 10)
      { string entete = " affichage de la matrice de raideur  apres chargement d'efforts externes ";
        matglob.Affichage_ecran(entete); 
        entete = " affichage du second membre apres chargement d'efforts externes ";
        vecglob.Affichage_ecran(entete); 
      };
  }
  catch (ErrSortieFinale)
       // cas d'une direction voulue vers la sortie
       // on relance l'interuption pour le niveau supérieur
     { ErrSortieFinale toto;
       throw (toto);
     }
  catch ( ... )
   { if (ParaGlob::NiveauImpression() >= 1)
       {cout << "\n *** erreur dans le calcul du chargement impose ";
        if (ParaGlob::NiveauImpression() >= 4)
              cout <<  "\n Charge::ChargeSMembreRaideur_Im_mecaSolid(..";
        cout << flush;
       };
      retour = false;
    };
 
    temps_cpu_chargement.Arret_du_comptage();
    return retour;
 };