// 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 "ElContact.h"
#include "Droite.h"
#include "ConstMath.h"
#include "MathUtil.h"
#include "ElemMeca.h"
#include "Util.h"
#include <cmath>
#include "TypeConsTens.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"

        using namespace std;

//---------- variables statiques --------------
list <DdlElement> ElContact::list_Ddl_global; // liste de tous les DdlElements des éléments de contact
list <Vecteur> ElContact::list_SM; // list  de second membre local: sert pour tous les éléments de contact
list <Mat_pleine> ElContact::list_raideur; // list des raideurs locales:     "     "        "         "
// stockage du maximum de distance tolérée entre noeud à tdt et le projeté, sert pour éliminer les contacts aberrants
double ElContact::dep_max=1.e15; // au début un nombre très grand par défaut
int ElContact::niveau_commentaire=0; // init par défaut -> est mis à jour par les_contacts
Fonction_nD * ElContact::fct_pilotage_contact4=NULL ; //  pour le pilotage du type de contact 4


// stockage transitoire pour les quelconques vraiment quelconque
Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_penalisationPenetration;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_penalisationPenetration;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penalisationPenetration;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_contact_maxi;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_penetration_contact_maxi;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_contact_maxi;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_force_contact_noeud_maxi;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_force_contact_noeud_maxi;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_force_contact_noeud_maxi;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_penalisationTangentielle;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_penalisationTangentielle;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penalisationTangentielle;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_tangentielle_contact_maxi;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_tangentielle_contact_maxi;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_tangentielle_contact_maxi;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_tangentielle_borne_regularisation;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_tangentielle_borne_regularisation;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_tangentielle_borne_regularisation;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_force_tangentielle_noeud_maxi;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_force_tangentielle_noeud_maxi;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_force_tangentielle_noeud_maxi;

Tableau <const TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_const_fct_nD_niveau_commentaire;
Tableau < TypeQuelconque *  > ElContact::Fct_nD_contact::tqi_fct_nD_niveau_commentaire;
Tableau <int>  ElContact::Fct_nD_contact::t_num_ordre_fct_nD_niveau_commentaire;


//----------- fin variables statiques ----------

///------- particularité class Fct_nD_contact ------------------
ElContact::Fct_nD_contact::Fct_nD_contact():
    fct_nD_penalisationPenetration(NULL)
    ,fct_nD_penetration_contact_maxi(NULL)
    ,fct_nD_penetration_borne_regularisation(NULL)
    ,fct_nD_force_contact_noeud_maxi(NULL)
    ,fct_nD_penalisationTangentielle(NULL)
    ,fct_nD_tangentielle_contact_maxi(NULL)
    ,fct_nD_tangentielle_borne_regularisation(NULL)
    ,fct_nD_force_tangentielle_noeud_maxi(NULL)
    ,fct_niveau_commentaire(NULL)
      {};
   
ElContact::Fct_nD_contact::Fct_nD_contact(const Fct_nD_contact& a) :
    fct_nD_penalisationPenetration(a.fct_nD_penalisationPenetration)
    ,fct_nD_penetration_contact_maxi(a.fct_nD_penetration_contact_maxi)
    ,fct_nD_penetration_borne_regularisation(a.fct_nD_penetration_borne_regularisation)
    ,fct_nD_force_contact_noeud_maxi(a.fct_nD_force_contact_noeud_maxi)
    ,fct_nD_penalisationTangentielle(a.fct_nD_penalisationTangentielle)
    ,fct_nD_tangentielle_contact_maxi(a.fct_nD_tangentielle_contact_maxi)
    ,fct_nD_tangentielle_borne_regularisation(a.fct_nD_tangentielle_borne_regularisation)
    ,fct_nD_force_tangentielle_noeud_maxi(a.fct_nD_force_tangentielle_noeud_maxi)
    ,fct_niveau_commentaire(a.fct_niveau_commentaire)
      {};

ElContact::Fct_nD_contact::~Fct_nD_contact()
  {};
  
ElContact::Fct_nD_contact& ElContact::Fct_nD_contact::operator= (const ElContact::Fct_nD_contact& a)
  {fct_nD_penalisationPenetration = a.fct_nD_penalisationPenetration;
   fct_nD_penetration_contact_maxi = a.fct_nD_penetration_contact_maxi;
   fct_nD_penetration_borne_regularisation = a.fct_nD_penetration_borne_regularisation;
   fct_nD_force_contact_noeud_maxi = a.fct_nD_force_contact_noeud_maxi;
   fct_nD_penalisationTangentielle = a.fct_nD_penalisationTangentielle;
   fct_nD_tangentielle_contact_maxi = a.fct_nD_tangentielle_contact_maxi;
   fct_nD_tangentielle_borne_regularisation = a.fct_nD_tangentielle_borne_regularisation;
   fct_nD_force_tangentielle_noeud_maxi = a.fct_nD_force_tangentielle_noeud_maxi;
   fct_niveau_commentaire = a.fct_niveau_commentaire;
   return *this;
  
  };
  
  // initialisation des conteneurs statique des fonction nD
  void ElContact::Fct_nD_contact::Init_conteneur_statique()
    {// on initialise chaque fct nD
     if (fct_nD_penalisationPenetration != NULL)
        Definition_conteneurs_static_TypeQuelconque
         (fct_nD_penalisationPenetration
          ,tqi_fct_nD_penalisationPenetration,tqi_const_fct_nD_penalisationPenetration
          ,t_num_ordre_fct_nD_penalisationPenetration);
     if (fct_nD_penetration_contact_maxi != NULL)
        Definition_conteneurs_static_TypeQuelconque
          (fct_nD_penetration_contact_maxi
           ,tqi_fct_nD_penetration_contact_maxi,tqi_const_fct_nD_penetration_contact_maxi
           ,t_num_ordre_fct_nD_penetration_contact_maxi);
     if (fct_nD_penetration_borne_regularisation != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_nD_penetration_borne_regularisation
          ,tqi_fct_nD_penetration_borne_regularisation,tqi_const_fct_nD_penetration_borne_regularisation
          ,t_num_ordre_fct_nD_penetration_borne_regularisation);
     if (fct_nD_force_contact_noeud_maxi != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_nD_force_contact_noeud_maxi
          ,tqi_fct_nD_force_contact_noeud_maxi,tqi_const_fct_nD_force_contact_noeud_maxi
          ,t_num_ordre_fct_nD_force_contact_noeud_maxi);
     if (fct_nD_penalisationTangentielle != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_nD_penalisationTangentielle
          ,tqi_fct_nD_penalisationTangentielle,tqi_const_fct_nD_penalisationTangentielle
          ,t_num_ordre_fct_nD_penalisationTangentielle);
     if (fct_nD_tangentielle_contact_maxi != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_nD_tangentielle_contact_maxi
          ,tqi_fct_nD_tangentielle_contact_maxi,tqi_const_fct_nD_tangentielle_contact_maxi
          ,t_num_ordre_fct_nD_tangentielle_contact_maxi);
     if (fct_nD_tangentielle_borne_regularisation != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_nD_tangentielle_borne_regularisation
          ,tqi_fct_nD_tangentielle_borne_regularisation,tqi_const_fct_nD_tangentielle_borne_regularisation
          ,t_num_ordre_fct_nD_tangentielle_borne_regularisation);
     if (fct_nD_force_tangentielle_noeud_maxi != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_nD_force_tangentielle_noeud_maxi
          ,tqi_fct_nD_force_tangentielle_noeud_maxi,tqi_const_fct_nD_force_tangentielle_noeud_maxi
          ,t_num_ordre_fct_nD_force_tangentielle_noeud_maxi);
     if (fct_niveau_commentaire != NULL)
           Definition_conteneurs_static_TypeQuelconque
         (fct_niveau_commentaire
          ,tqi_fct_nD_niveau_commentaire,tqi_const_fct_nD_niveau_commentaire
          ,t_num_ordre_fct_nD_niveau_commentaire);

    };

  // définition des conteneurs static des TypeQuelconque
  void ElContact::Fct_nD_contact::Definition_conteneurs_static_TypeQuelconque
       (Fonction_nD * pt_fonct,Tableau < TypeQuelconque *  >& tqi,Tableau < const TypeQuelconque *  >& tqii
        ,Tableau <int>&  t_num_ordre )
   {// on commence par récupérer le conteneurs des grandeurs à fournir
    const Tableau <EnumTypeQuelconque>& tab_queconque = pt_fonct->Tab_enu_quelconque();
    //  on dimensionne
    int tail_tab_queconque = tab_queconque.Taille();
    tqi.Change_taille(tail_tab_queconque);
    t_num_ordre.Change_taille(tail_tab_queconque);
    Grandeur_scalaire_entier grand_courant_entier(0); // par défaut pour la création des conteneurs quelconques
    Grandeur_scalaire_double grand_courant_double(0.); // par défaut pour la création des conteneurs quelconques
    // grandeurs de travail
    int dim = ParaGlob::Dimension();
    Coordonnee  inter(dim);
    Grandeur_coordonnee   grand_coor_courant(inter);

    // on va balayer les grandeurs quelconques
    for (int i=1;i<= tail_tab_queconque;i++)
    { bool trouver = false;
      EnumTypeQuelconque enu = tab_queconque(i);

      tqi(i) = NULL; // init
      t_num_ordre(i) = 1; //***** pour l'instant uniquement la première coordonnée s'il s'agit d'un Coordonnee
       // *** à abonder par la suite !!!!!!!!!!
      
      // on regarde tout d'abord les grandeurs spécifiques à l'élément contact
      switch (enu)
        { // il y a des grandeurs vectorielles qui pour l'instant ne sont pas prises en compte
          // cf. fct nD
          // NORMALE_CONTACT, GLISSEMENT_CONTACT ,PENETRATION_CONTACT,FORCE_CONTACT,
          case CONTACT_NB_DECOL:
            { tqi(i) = new TypeQuelconque(CONTACT_NB_DECOL,X1,grand_courant_entier);
              trouver = true;
              break;
            }
          case CONTACT_PENALISATION_N:
            { tqi(i) = new TypeQuelconque(CONTACT_PENALISATION_N,X1,grand_courant_double);
              trouver = true;
              break;
            }
          case CONTACT_PENALISATION_T:
            { tqi(i) = new TypeQuelconque(CONTACT_PENALISATION_T,X1,grand_courant_double);
              trouver = true;
              break;
            }
          case CONTACT_NB_PENET:
            { tqi(i) = new TypeQuelconque(CONTACT_NB_PENET,X1,grand_courant_entier);
              trouver = true;
              break;
            }
          case CONTACT_CAS_SOLIDE:
            { tqi(i) = new TypeQuelconque(CONTACT_CAS_SOLIDE,X1,grand_courant_entier);
              trouver = true;
              break;
            }
          case CONTACT_ENERG_GLISSE_ELAS:
            { tqi(i) = new TypeQuelconque(CONTACT_ENERG_GLISSE_ELAS,X1,grand_courant_double);
              trouver = true;
              break;
            }
          case CONTACT_ENERG_GLISSE_PLAS:
            { tqi(i) = new TypeQuelconque(CONTACT_ENERG_GLISSE_PLAS,X1,grand_courant_double);
              trouver = true;
              break;
            }
          case CONTACT_ENERG_GLISSE_VISQ:
            { tqi(i) = new TypeQuelconque(CONTACT_ENERG_GLISSE_VISQ,X1,grand_courant_double);
              trouver = true;
              break;
            }
          case CONTACT_ENERG_PENAL:
            { tqi(i) = new TypeQuelconque(CONTACT_ENERG_PENAL,X1,grand_courant_double);
              trouver = true;
              break;
            }
          case NOEUD_PROJECTILE_EN_CONTACT:
            { tqi(i) = new TypeQuelconque(NOEUD_PROJECTILE_EN_CONTACT,X1,grand_courant_double);
              trouver = true;
              break;
             }
          case NOEUD_FACETTE_EN_CONTACT:
            { tqi(i) = new TypeQuelconque(NOEUD_FACETTE_EN_CONTACT,X1,grand_courant_double);
//              Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*tqi(i)).Grandeur_pointee()));
              trouver = true;
              break;
            }
          case NUM_NOEUD:
            { tqi(i) = new TypeQuelconque(NUM_NOEUD,X1,grand_courant_entier);
              trouver = true;
              break;
            }
          case NUM_MAIL_NOEUD:
            { tqi(i) = new TypeQuelconque(NUM_MAIL_NOEUD,X1,grand_courant_entier);
              trouver = true;
              break;
            }
          case POSITION_GEOMETRIQUE:
            {tqi(i) = new TypeQuelconque(POSITION_GEOMETRIQUE,X1,grand_coor_courant);
             trouver = true;
              break;
             }
          case POSITION_GEOMETRIQUE_t:
            {tqi(i) = new TypeQuelconque(POSITION_GEOMETRIQUE_t,X1,grand_coor_courant);
             trouver = true;
              break;
             }
          case POSITION_GEOMETRIQUE_t0:
            {tqi(i) = new TypeQuelconque(POSITION_GEOMETRIQUE_t0,X1,grand_coor_courant);
             trouver = true;
              break;
             }
          // *** pour l'instant la suite n'est pas opérationelle car il s'agit de vecteur
          // dont les composantes n'ont pas d'équivalent en ddl_enum_etendu
          //  il faudrait les définir si on veut pouvoir s'en servir
          case PENETRATION_CONTACT:
             {tqi(i) = new TypeQuelconque(PENETRATION_CONTACT,X1,grand_coor_courant);
              trouver = true;
               break;
              }
          case PENETRATION_CONTACT_T:
             {tqi(i) = new TypeQuelconque(PENETRATION_CONTACT_T,X1,grand_coor_courant);
              trouver = true;
               break;
              }
          case GLISSEMENT_CONTACT:
             {tqi(i) = new TypeQuelconque(GLISSEMENT_CONTACT,X1,grand_coor_courant);
              trouver = true;
               break;
              }
          case GLISSEMENT_CONTACT_T:
             {tqi(i) = new TypeQuelconque(GLISSEMENT_CONTACT_T,X1,grand_coor_courant);
              trouver = true;
               break;
              }
          case NORMALE_CONTACT:
             {tqi(i) = new TypeQuelconque(NORMALE_CONTACT,X1,grand_coor_courant);
              trouver = true;
               break;
              }
          case FORCE_CONTACT:
             {tqi(i) = new TypeQuelconque(FORCE_CONTACT,X1,grand_coor_courant);
              trouver = true;
               break;
              }
          case FORCE_CONTACT_T:
           {tqi(i) = new TypeQuelconque(POSITION_GEOMETRIQUE,X1,grand_coor_courant);
            trouver = true;
             break;
            }
          case CONTACT_COLLANT:
           { tqi(i) = new TypeQuelconque(CONTACT_COLLANT,X1,grand_courant_entier);
             trouver = true;
             break;
           }
          case NUM_ZONE_CONTACT:
           { tqi(i) = new TypeQuelconque(NUM_ZONE_CONTACT,X1,grand_courant_entier);
             trouver = true;
             break;
           }
        case NUM_ELEMENT:
         { tqi(i) = new TypeQuelconque(NUM_ELEMENT,X1,grand_courant_entier);
           trouver = true;
           break;
         }
        case NUM_MAIL_ELEM:
         { tqi(i) = new TypeQuelconque(NUM_MAIL_ELEM,X1,grand_courant_entier);
           trouver = true;
           break;
         }
         //*** fin

          default:
            trouver = false;
            break;
       };

      // si on n'a rien trouvé
      if (!trouver)
        {cout << "\n erreur***: la grandeur  " << NomTypeQuelconque(enu)
              << " n'existe pas les elements de contact "
              << " la fonction nD " << pt_fonct->NomFonction()
              << " ne peut pas etre renseignee "
              << "\n ElContact::Fct_nD_contact::Definition_conteneurs_static_TypeQuelconque(.." << flush;
         pt_fonct->Affiche();
         Sortie(1);
        };

    };
    // tableau de const intermédiaire pour l'appel de la fonction
    tqii.Change_taille(tail_tab_queconque);
    for (int i=1;i<= tail_tab_queconque;i++)
      tqii(i) = tqi(i);
   };

///------- fin particularité class Fct_nD_contact ------------------


    // CONSTRUCTEURS :
// par defaut
ElContact::ElContact () :
  tabNoeud(),tabNoeud_t()
  ,Mtdt(ParaGlob::Dimension()),Mt(0) // Mt de dim 0, au début pour dire qu'il n'est pas activé
  ,M_noeud_tdt_avant_projection(ParaGlob::Dimension())
  ,phi_theta_0()
  ,force_contact(),force_contact_t(),tabForce_cont(),tabForce_cont_t()
  ,F_N_max(0.),F_T_max(0.),F_N_max_t(0.),F_T_max_t(0.)
  ,ddlElement_assemblage(NULL),raideur(NULL),residu(NULL)
  ,energie_penalisation(0.),energie_frottement()
  ,elfront(NULL),elfront_t(NULL),noeud(NULL),noeud_t(NULL)
  ,num_zone_contact(0),normale_lisser(0)
  ,nb_decol_t(0),nb_decol_tdt(0),gap_t(0.),gap_tdt(0.),nb_pene_t(0.),nb_pene_tdt(0.)
  ,dep_T_t(0.),dep_T_tdt(0.)
  ,mult_pene_t(1.),mult_pene_tdt(1.)
  ,mult_tang_t(1.),mult_tang_tdt(1.)
  ,actif(1),type_trajectoire_t(0),type_trajectoire_tdt(0)
  ,cas_solide(0),cas_collant(0),nb_change_frontiere(0)
  ,tabNoeud_pour_assemblage(),tab_posi_esclave()
  ,nb_posi_esclave_stocker_t(0),nb_posi_esclave_stocker_tdt(0)
  ,indice_stockage_glissant_t(1),indice_stockage_glissant_tdt(1)
  ,penalisation(0.),penalisation_tangentielle(0.)
  ,penalisation_t(0.),penalisation_tangentielle_t(0.)
  // utilisation éventuelle de fonction nD
  ,fct_nD_contact()
  
  ,N(),dep_tangentiel()
  // pour le contact 4, pour le calcul de la pénalisation avec une moyenne glissante
  ,val_penal(10),pt_dans_val_penal(1)
  {};
    // fonction d'un pointeur d'element frontiere et d'un pointeur de noeud
  // du fait éventuel qu'il peut-être collant ou pas
ElContact::ElContact ( const Front * elf, const Noeud * noe, Fct_nD_contact & fct_contact_,int collant):
  tabNoeud(),tabNoeud_t()
  ,Mtdt(ParaGlob::Dimension()),Mt(0) // Mt de dim 0, au début pour dire qu'il n'est pas activé
  ,M_noeud_tdt_avant_projection(ParaGlob::Dimension())
  ,phi_theta_0()
  ,force_contact(ParaGlob::Dimension()),tabForce_cont()
  ,force_contact_t(ParaGlob::Dimension()),tabForce_cont_t()
  ,F_N_max(0.),F_T_max(0.),F_N_max_t(0.),F_T_max_t(0.)
  ,ddlElement_assemblage(NULL),raideur(NULL),residu(NULL)
  ,energie_penalisation(0.),energie_frottement()
  ,elfront(NULL),elfront_t(NULL),noeud(NULL),noeud_t(NULL)
  ,num_zone_contact(0),normale_lisser(0)
  ,nb_decol_t(0),nb_decol_tdt(0),gap_t(0.),gap_tdt(0.),nb_pene_t(0.),nb_pene_tdt(0.)
  ,dep_T_t(0.),dep_T_tdt(0.)
  ,mult_pene_t(1.),mult_pene_tdt(1.),mult_tang_t(1.),mult_tang_tdt(1.)
  ,actif(1),type_trajectoire_t(0),type_trajectoire_tdt(0)
  ,cas_solide(0),cas_collant(collant),nb_change_frontiere(0)
  ,tabNoeud_pour_assemblage(),tab_posi_esclave()
  ,nb_posi_esclave_stocker_t(0),nb_posi_esclave_stocker_tdt(0)
  ,indice_stockage_glissant_t(1),indice_stockage_glissant_tdt(1)
  ,penalisation(0.),penalisation_tangentielle(0.)
  ,penalisation_t(0.),penalisation_tangentielle_t(0.)
  // utilisation éventuelle de fonction nD
  ,fct_nD_contact(fct_contact_)

  ,N(),dep_tangentiel()
  // pour le contact 4, pour le calcul de la pénalisation avec une moyenne glissante
  ,val_penal(10),pt_dans_val_penal(1)
  {  noeud = ( Noeud *) noe; // le noeud esclave
     noeud_t = noeud;
     elfront = new Front(*elf);
     elfront_t=elfront;
     // on cree un nouvelle element generique frontiere identique
     // pour stocker les info specifique aux intersections et tangence
     // par contre les info generique (noeuds etc) sont identiques
     Construction_TabNoeud(); // on construit le tableau de noeud global
     tabNoeud_t = tabNoeud; // init
     tabForce_cont.Change_taille(tabNoeud.Taille()-1);
     tabForce_cont_t.Change_taille(tabNoeud.Taille()-1);
   
  };
// de copie
ElContact::ElContact ( const ElContact & a): 
  tabNoeud(a.tabNoeud),tabNoeud_t(a.tabNoeud_t)
  ,Mtdt(a.Mtdt),Mt(a.Mt),M_noeud_tdt_avant_projection(a.M_noeud_tdt_avant_projection)
  ,phi_theta_0(a.phi_theta_0)
  ,force_contact(a.force_contact),tabForce_cont(a.tabForce_cont)
  ,force_contact_t(a.force_contact_t),tabForce_cont_t(a.tabForce_cont_t)
  ,F_N_max(a.F_N_max),F_T_max(a.F_T_max)
  ,F_N_max_t(a.F_N_max_t),F_T_max_t(a.F_T_max_t)
  ,ddlElement_assemblage(a.ddlElement_assemblage)
  ,raideur(a.raideur),residu(a.residu)
  ,energie_penalisation(a.energie_penalisation)
  ,energie_frottement(a.energie_frottement)
  ,elfront(NULL),elfront_t(NULL),noeud(NULL),noeud_t(NULL)
  ,num_zone_contact(a.num_zone_contact),normale_lisser(a.normale_lisser)
  ,nb_decol_t(a.nb_decol_t),nb_decol_tdt(a.nb_decol_tdt)
  ,gap_t(a.gap_t),gap_tdt(a.gap_tdt),nb_pene_t(a.nb_pene_t),nb_pene_tdt(a.nb_pene_tdt)
  ,dep_T_t(a.dep_T_t),dep_T_tdt(a.dep_T_tdt)
  ,mult_pene_t(a.mult_pene_t),mult_pene_tdt(a.mult_pene_tdt)
  ,mult_tang_t(a.mult_tang_t),mult_tang_tdt(a.mult_tang_tdt)
  ,actif(a.actif),type_trajectoire_t(a.type_trajectoire_t)
  ,type_trajectoire_tdt(a.type_trajectoire_tdt)
  ,cas_solide(a.cas_solide),cas_collant(a.cas_collant),nb_change_frontiere(a.nb_change_frontiere)
  ,tabNoeud_pour_assemblage(a.tabNoeud_pour_assemblage),tab_posi_esclave(a.tab_posi_esclave)
  ,nb_posi_esclave_stocker_t(a.nb_posi_esclave_stocker_t),nb_posi_esclave_stocker_tdt(a.nb_posi_esclave_stocker_tdt)
  ,indice_stockage_glissant_t(a.indice_stockage_glissant_t),indice_stockage_glissant_tdt(a.indice_stockage_glissant_t)
  ,penalisation(a.penalisation),penalisation_tangentielle(a.penalisation_tangentielle)
  ,penalisation_t(a.penalisation_t),penalisation_tangentielle_t(a.penalisation_tangentielle_t)
  // utilisation éventuelle de fonction nD
  ,fct_nD_contact(a.fct_nD_contact)

  ,N(a.N),dep_tangentiel(a.dep_tangentiel)
  // pour le contact 4, pour le calcul de la pénalisation avec une moyenne glissante
  ,val_penal(a.val_penal),pt_dans_val_penal(a.pt_dans_val_penal)
  {  noeud = a.noeud;
     noeud_t = noeud;
     elfront = new Front(*(a.Elfront()));
     elfront_t=elfront;
     // on cree un nouvelle element generique frontiere identique
     // pour stocker les info specifique aux intersections et tangence
     // par contre les info generique (noeuds etc) sont identiques
// non, car normalement toutes les initialisation on déjà été faites dans "a"     Construction_TabNoeud(); // on construit le tableau de noeud global
  }; 
      // DESTRUCTEUR :
ElContact::~ElContact ()
 { Libere(); 
 };
    
    // METHODES PUBLIQUES :
// affectation de toute les donnees
ElContact& ElContact::operator = ( const ElContact& a)
  { tabNoeud=a.tabNoeud;tabNoeud_t=a.tabNoeud_t;
    Mtdt=a.Mtdt;Mt=a.Mt;M_noeud_tdt_avant_projection=a.M_noeud_tdt_avant_projection;
    phi_theta_0=a.phi_theta_0;
    force_contact=a.force_contact;tabForce_cont=a.tabForce_cont;
    force_contact_t=a.force_contact_t;tabForce_cont_t=a.tabForce_cont_t;
    F_N_max=a.F_N_max;F_T_max=a.F_T_max;
    F_N_max_t=a.F_N_max_t;F_T_max_t=a.F_T_max_t;
    ddlElement_assemblage=a.ddlElement_assemblage;
    raideur=a.raideur;residu=a.residu;
    energie_penalisation=a.energie_penalisation;
    energie_frottement=a.energie_frottement;
    num_zone_contact=a.num_zone_contact;normale_lisser=a.normale_lisser;
    nb_decol_t=a.nb_decol_t;nb_decol_tdt=a.nb_decol_tdt;
    gap_t=a.gap_t;gap_tdt=a.gap_tdt;nb_pene_t=a.nb_pene_t;nb_pene_tdt=a.nb_pene_tdt;
    dep_T_t=a.dep_T_t;dep_T_tdt=a.dep_T_tdt;
    mult_pene_t=a.mult_pene_t;mult_pene_tdt=a.mult_pene_tdt;
    mult_tang_t=a.mult_tang_t;mult_tang_tdt=a.mult_tang_tdt;
    actif=a.actif;type_trajectoire_t=a.type_trajectoire_t;
    type_trajectoire_tdt=a.type_trajectoire_tdt;
    cas_solide=a.cas_solide;cas_collant=a.cas_collant;nb_change_frontiere=a.nb_change_frontiere;
    tabNoeud_pour_assemblage=a.tabNoeud_pour_assemblage;tab_posi_esclave=a.tab_posi_esclave;
    nb_posi_esclave_stocker_t=a.nb_posi_esclave_stocker_t;nb_posi_esclave_stocker_tdt=a.nb_posi_esclave_stocker_tdt;
    indice_stockage_glissant_t=a.indice_stockage_glissant_t;indice_stockage_glissant_tdt=a.indice_stockage_glissant_t;
    penalisation=a.penalisation;penalisation_tangentielle=a.penalisation_tangentielle;
    penalisation_t=a.penalisation_t;penalisation_tangentielle_t=a.penalisation_tangentielle_t;
    // utilisation éventuelle de fonction nD
    fct_nD_contact=a.fct_nD_contact;
    N=a.N;dep_tangentiel=a.dep_tangentiel;
    // pour le contact 4, pour le calcul de la pénalisation avec une moyenne glissante
    val_penal=a.val_penal;pt_dans_val_penal=a.pt_dans_val_penal;
    noeud = a.noeud;noeud_t = noeud;
    elfront = new Front(*(a.Elfront()));
    elfront_t=elfront;
    return *this;
  };
          

// affichage à l'écran des informations liées au contact
// cas = 0 : affichage en fonction du niveau de commentaire voulu
// cas = 1 : affichage minimal nécessaire pour repérer l'élément
// cas = 2 : affichage minimal en fonction du niveau de commentaire pour les éléments de contact
void ElContact::Affiche(int cas) const
 {
   if (cas == 1)
    {cout << "el contact: noe "
          << noeud->Num_noeud() <<" mail. " << noeud->Num_Mail()
          << " avec front:";
     if (elfront->Angle_mort())
       cout << "(angle mort)";
     if (normale_lisser)
       cout << " (lissage N) ";
     if (!actif) cout << " (inactif) ";
     cout << elfront->Num_frontiere()
          << " de l'EF " << elfront->PtEI()->Num_elt() << " du mail. "
          << elfront->PtEI()->Num_maillage() <<   flush;;
     if (Permet_affichage() > 6)
       {const Tableau <Front*>* tabmit = elfront->TabMitoyen();
        if (tabmit != NULL)
          {int taille = tabmit->Taille();
           for (int i=1;i<= taille;i++)
             cout << " ele_mitoy: "<< (*tabmit)(i)->PtEI()->Num_elt()
                  << " du mail. "
                  << (*tabmit)(i)->PtEI()->Num_maillage() <<   flush;
          };
       };
    }
   else if (cas == 2)
    {if (Permet_affichage() > 5)
      {cout << "\n el contact: noe "
            << noeud->Num_noeud() <<" mail. " << noeud->Num_Mail()
            << " avec front:";
       if (elfront->Angle_mort())
         cout << "(angle mort)";
       if (normale_lisser)
         cout << " (lissage N) ";
       if (!actif) cout << " (inactif) ";
       cout << elfront->Num_frontiere()
            << " de l'EF " << elfront->PtEI()->Num_elt() << " du mail. "
            << elfront->PtEI()->Num_maillage() <<   flush;
       if (Permet_affichage() > 6)
         {const Tableau <Front*>* tabmit = elfront->TabMitoyen();
          if (tabmit != NULL)
            {int taille = tabmit->Taille();
             for (int i=1;i<= taille;i++)
               cout << " ele_mitoy: "<< (*tabmit)(i)->PtEI()->Num_elt()
                    << " du mail. "
                    << (*tabmit)(i)->PtEI()->Num_maillage() <<   flush;
            };
         };
      };
    // sinon pas d'affichage
    }
   else
    { cout << "\n element de contact ";
      if (!actif) cout << " Inactif ! ";
      if (normale_lisser) cout << " (lissage N) ";
      cout << " constitue du noeud " << noeud->Num_noeud() << " du maillage "
           << noeud->Num_Mail() << ", et de l'element de frontiere suivant: ";
      elfront->Affiche();
      cout << "\n F_contact_noeud: "<< force_contact;
      cout << "\n F_contact_facette: "<< tabForce_cont;
      cout << "\n norme_F_N: " << F_N_max
           << ", norme_F_T: " << F_T_max;
      cout << "\n normale: "<< N
           << " dep_tangentiel: " << dep_tangentiel;
      cout << "\n num_zone_contact " << num_zone_contact
           << ", Mtdt: "<< Mtdt
           << ", M_noeud_tdt_avant_projection " << M_noeud_tdt_avant_projection
           << ", energie_frottement= "<< energie_frottement
           << ", energie_penalisation= "<< energie_penalisation
           << ", cas_solide= " << cas_solide
           << ", cas_collant= " << cas_collant
           << ", nb_decol_tdt= " << nb_decol_tdt
           << ", gap_tdt= " << gap_tdt
           << ", dep_T_tdt= " << dep_T_tdt
           << ", nb_pene_tdt= " << nb_pene_tdt
           << ", mult_pene_tdt= " << mult_pene_tdt
           << ", mult_tang_tdt= " << mult_tang_tdt
           << ", penalisation= " << penalisation
           << ", penalisation_tangentielle= " << penalisation_tangentielle
           << ", type_trajectoire_tdt= " << type_trajectoire_tdt ;
     if (Permet_affichage() > 2)
      {cout << "\n fct_nD_penalisationPenetration: ";
       if (fct_nD_contact.fct_nD_penalisationPenetration != NULL )
            {cout << fct_nD_contact.fct_nD_penalisationPenetration->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_penetration_contact_maxi: ";
       if (fct_nD_contact.fct_nD_penetration_contact_maxi != NULL )
            {cout << fct_nD_contact.fct_nD_penetration_contact_maxi->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_penetration_borne_regularisation: ";
       if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL )
            {cout << fct_nD_contact.fct_nD_penetration_borne_regularisation->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_force_contact_noeud_maxi: ";
       if (fct_nD_contact.fct_nD_force_contact_noeud_maxi != NULL )
            {cout << fct_nD_contact.fct_nD_force_contact_noeud_maxi->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_penalisationTangentielle: ";
       if (fct_nD_contact.fct_nD_penalisationTangentielle != NULL )
            {cout << fct_nD_contact.fct_nD_penalisationTangentielle->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_tangentielle_contact_maxi: ";
       if (fct_nD_contact.fct_nD_tangentielle_contact_maxi != NULL )
            {cout << fct_nD_contact.fct_nD_tangentielle_contact_maxi->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_tangentielle_borne_regularisation: ";
       if (fct_nD_contact.fct_nD_tangentielle_borne_regularisation != NULL )
            {cout << fct_nD_contact.fct_nD_tangentielle_borne_regularisation->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_nD_force_tangentielle_noeud_maxi: ";
       if (fct_nD_contact.fct_nD_force_tangentielle_noeud_maxi != NULL )
            {cout << fct_nD_contact.fct_nD_force_tangentielle_noeud_maxi->NomFonction();}
       else { cout << " NULL "; };
       cout << "\n fct_niveau_commentaire: ";
       if (fct_nD_contact.fct_niveau_commentaire != NULL )
            {cout << fct_nD_contact.fct_niveau_commentaire->NomFonction();}
       else { cout << " NULL "; };
      };
   };

 };            

// calcul d'un contact eventuel entre le noeud esclave et la frontiere
// ramene true s'il y a contact
// si init = false, on recherche le contact a partir du precedent point sauvegarde
// sinon on commence a l'aide d'element de reference,
// et on calcule et sauvegarde les coordonnées
// initiale locales theta^i du point de contact

// si le contact existe  et si l'algo le demande (cf. ParaAlgoControle) :
// le noeud pourrait-être ramené sur la surface mais dans les faits:
// on ne fait pas de projection, sinon on ne peut pas tester plusieurs contacts pour
// choisir le meilleur, puisque les choses changent entre avant et après le test de contact
// donc ici la position du noeud esclave n'est pas modifiée
bool ElContact::Contact( bool init)
 {  int test = false; // pas de contact par défaut
    // on dimensionne éventuellement le tableau des positions successives
    if (tab_posi_esclave.Taille() != ParaGlob::param->ParaAlgoControleActifs().Nb_moy_glissant())
       tab_posi_esclave.Change_taille
             (ParaGlob::param->ParaAlgoControleActifs().Nb_moy_glissant(),Coordonnee(ParaGlob::Dimension()));
  
/*   // on exclue les cas non traité actuellement
//   if ((ParaGlob::Dimension() == 3) && (elfront->Eleme()->ElementGeometrique().Dimension()==1))
//   	{// le cas du contact avec des lignes en 3D n'est pas encore implanté, sauf si le ddl 3 est bloqué pour tous les noeuds
//		 int tail = tabNoeud.Taille();
//		 bool ok_apriori = true;
//		 for (int i=1;i<=tail;i++)
//		   if (!(tabNoeud(i)->Ddl_fixe(X3)))
//			   {ok_apriori=false;break;};
//		 
//   	 if ((ParaGlob::NiveauImpression() >= 6) && (!ok_apriori))
//   	 	{ cout << "\n *** warning : le contact avec des lignes n'est pas encore implante en 3D ";
//			  return test;
//			 }; 	
//   	};

//   // le cas de la trajectoire pose des pb donc on va tout d'abor regarder si le noeud est interne à l'élément
//   // si oui, on le déclare actif et dans ce cas la trajectoire sera la normale
//   {// récup de l'élément qui contient la frontière
//    Element& elem = *(elfront->PtEI());
//    if (actif ==0) // si actif est différent de 0, alors cela signifie qu'il y a déjà contact, on ne va pas plus loin
//      if (elem.In_boite_emcombrement_elem(noeud->Coord2()))
//       {actif=1;
//////debug
//if (noeud->Num_noeud() == 84)
// {cout << "\n debug ElContact::Contact( bool init)  temps " << ParaGlob::Variables_de_temps().TempsCourant();
//  noeud->Coord2().Affiche();
////  cout << " actif = " << actif;
////  elem.RecupBoite_encombre_element().Premier().Affiche();
////  elem.RecupBoite_encombre_element().Second().Affiche();
// };

//fin debug
//
//       
//       }; // cela va modifier la manière de calculer la trajectoire
//   }; */

    // on exclue un cas particulier qui ne peut pas être pris en compte a priori
    // c'est le cas 1D du contact de la frontière 1D d'un segment avec un point
    // le déplacement est forcément colinéaire avec la frontière
    // donc les vrais frontières pertinentes
    // sont les points d'extrémité mais pas le segment, mais c'est cependant la ligne que
    // l'on utilisera pour avoir la notion de in et out donc

////debug
//         cout << "\n ElContact::Contact: "; this->Affiche(1);
//// fin debug
    if ((ParaGlob::Dimension() == 1) // on est dans un espace 1D
        && (elfront->Eleme()->ElementGeometrique().Dimension()==1) // la frontière est une ligne ou un point
        && (elfront->Eleme_const()->TabNoeud_const().Taille() > 1 ) // + d'un seul noeud: c'est une ligne
       )
     {// le cas du contact avec des lignes en 1D n'est pas recevable a priori, ou plutôt il y a toujours intersection
      return false;
     };
    
// **** il y a des choses à faire au niveau de la droite tangente et du point d'impact
// lorsque l'on a une ligne d'angle mort, contrairement au cas du point, le point d'intersection
// doit-être calculé, et c'est important car il va servir pour RecupInfo et ensuite pour
// le calcul de la raideur et du second membre...
// bref il faut revoir la chose: c'est un peu moins simple que dans le cas d'un point_G

    // 1) === le cas des angles morts est particulier, on commence par le traiter
    if (elfront->Angle_mort())
      { //ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
        int dim = noeud->Dimension();
        Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
        N = Calcul_Normale(dim,pl,dr,indic);
        int permet_affiche = Permet_affichage();
        // la seule chose que l'on peut faire pour valider le fait qu'il y a contact ou pas
        // c'est de tester si le noeud est à l'intérieur de l'élément d'angle mort ou pas
        // c'est plus long qu'un test de surface mais il n'y a pas vraiment d'autre solution
        int cas = elfront->PtEI()->Interne_tdt(noeud->Coord2());
        if (permet_affiche > 5)
          {cout << "\n --  ElContact::Contact: ";
           cout << " Interne_tdt ? cas= " << cas <<", ";
           this->Affiche(1);
          };
        // retour : =0 le point est externe, =1 le point est interne ,
        //          = 2 le point est sur la frontière à la précision près
        if (cas == 0)
          return cas;
        // sinon on continue à tester
        // on déclare le noeud pas en contact si la distance entre le noeud et la frontière est trop grande:
        double dmax = ParaGlob::param->ParaAlgoControleActifs().DistanceMaxiAuPtProjete();
        //  pour cela il nous faut le point d'impact
        // calcul de la trajectoire a  utiliser
        Coordonnee V = Trajectoire(test);
        Coordonnee M1(dim);
        if (init)
          {M1 = Intersection(V,true);
          }
        else
          {M1 = Intersection(V,false);};
        // dans le cas ou M1 a une dimension  nulle, signifie qu'il n'y a pas d'intersection,
        // on peut alors dire qu'il ne peut pas y avoir contact donc retour
        if (M1.Dimension() == 0)
          return false;
        // on regarde la distance d'impact
        double d_MM1 = (M1-noeud->Coord2()).Norme();
        if (permet_affiche > 6)
         {cout << " intersection= ";M1.Affiche();
          cout << " d_MM1= "<<d_MM1<<", dmax= "<<dmax<<flush;
          this->Affiche(1);
         };

        if (d_MM1 > dmax)
          {return false;}
        else
          {return true;};
      };

   // 2) === suite pour le cas différent des angles morts

   // calcul de la trajectoire a  utiliser
   Coordonnee V = Trajectoire(test);
	
// -- le fait que l'on n'est pas en contact a priori est neutralisé, car on peut s'en accomoder (on essaie)
// s'il n'y a rien qui bouge: cas du tout début, ce sera la normale au point de référence qui sera la trajectoire
// dans le cas ou rien ne bouge -> pas de contact a priori  a l'initialisation
   // les solides doivent etre hors contacts, et ensuite le contact est a chaque fois regarde
//   if (!test) return false;

   // cas ou quelque chose a  bouge
	
// -- fin de la neutralisation

	
   // calcul de l'intersection de la trajectoire avec la frontiere
   int dim = noeud->Dimension();
   Coordonnee M1(dim);
   if (init)
     {M1 = Intersection(V,true);
     }
   else
     {M1 = Intersection(V,false);};
   // dans le cas ou M1 a une dimension  nulle, signifie qu'il n'y a pas d'intersection,
   // on peut alors dire qu'il ne peut pas y avoir contact donc retour
   if (M1.Dimension() == 0)
     return false;
     
     
   // pour la suite on suppose que l'on a trouvé une intersection
   // maintenant on regarde si l'intersection est acceptable
   // -- tout d'abord on arrete si l'intersection est vraiment trop loin du noeud
   //    donc 1) on teste le point projeter / à l'encombrement de l'élément frontière
   //    si le contact est collant, le point est concervé même en dehors de la boite
   //    c'est par exemple le cas point-facette, où le point est en dehors de la zone d'épaisseur
   if (!cas_collant)
    if (!(elfront->In_boite_emcombrement_front(M1)))
     return false;
   //         2) on regarde si la distance entre le point projeté et le noeud à tdt n'est pas
   //            anormalement grande: le test est également vrai pour le cas collant, pour
   //            éviter des contacts anormaux
   double d_MM1 = (M1-noeud->Coord2()).Norme();
  
//   if (d_MM1 > MaX(dep_max,ParaGlob::param->ParaAlgoControleActifs().DistanceMaxiAuPtProjete()))
//---- modif du 8 janvier: je ne comprend pas pourquoi j'ai pris en compte dep_max ici ??? donc je change
//     avec dep_max grand, cela supprime le paramétre de contrôle donc ... on ne contrôle plus rien ??
   if (d_MM1 > ParaGlob::param->ParaAlgoControleActifs().DistanceMaxiAuPtProjete())
     return false;
  // sinon c'est que le pt M1 est encore acceptable, on poursuit les tests	
//  cout << "\n debug 	ElContact::Contact distance =  " << d_MM1 << endl;

   // - cote surface
   ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
   // ici on recherche avec une précision donnée
   double extra_in_surface = ParaGlob::param->ParaAlgoControleActifs().PointInternePrecThetaiInterne();
   // dans le cas d'un contact avec pénalisation de type 7,
   // on considère le contact également pour un gap positif
   int type_penal = ParaGlob::param->ParaAlgoControleActifs().TypeCalculPenalisationPenetration();
   int contactType = ElContact::Recup_et_mise_a_jour_type_contact();
   if ((contactType == 2)|| (contactType == 41) || (contactType == 42))
    if ( abs(type_penal) == 7)
     {// recup des dernières différentes informations calculées au niveau de la cinématique de contact
      Coordonnee M_impact,N; // le point d'impact, la normale
      Vecteur phii;
      RecupInfo(N,M_impact,phii,false);
      
      // calcul de la pénétration normale: deltaX=de la surface au point pénétré
      Coordonnee deltaX = noeud->Coord2() - M_impact;
      // le N sort de la surface maître, donc rentre dans la surface esclave
      double gap= N * deltaX; // gap est négatif quand il y a pénétration
      double borne_regularisation;
      if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL)
           {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation
                  ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation
                  ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation
                  ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation);
           }
      else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();}
      // on calcul un facteur en fonction de la distance: = 0.25*exp(-gap/(0.25*borne))
      if (gap > borne_regularisation)
        // on considère qu'il n'y a pas contact dans ce cas
        {return false;}
      else if (elfro.InSurf(extra_in_surface))
        // on regarde si le point projeté est dans la surface de la frontière
       { Mtdt = M1; // sauvegarde
         return true;
       };
      // sinon on continue sur le cas normale
     };
  
   // on traite le cas particulier du contact collant
   if (cas_collant)
    { if (elfro.InSurf(extra_in_surface))
       // on ne s'occupe pas du coté de la projection pour le cas collant
       // on déclare le noeud pas en contact si la distance entre le noeud et la frontière est
       // trop grande
       { double dmax = ParaGlob::param->ParaAlgoControleActifs().DistanceMaxiAuPtProjete();
         double r=0;
         // on utilise BonCote_tdt uniquement pour calculer r
         elfront->BonCote_tdt(noeud->Coord2(),r);
         if (Dabs(r) > dmax)
          // pas de contact
          { return false;}
         else
          {Mtdt = M1; // sauvegarde
           return true;
          };
       };
    };
  
   // cas non collant: on regarde si le point projeté est dans la surface de la frontière
   double dismini = ParaGlob::param->ParaAlgoControleActifs().Precision_pt_sur_front();
   if (elfro.InSurf(extra_in_surface))
    { // on commence par traiter les cas classiques, en particulier on verra après le cas en 3D avec un point et une ligne
	     if ( !((ParaGlob::Dimension() == 3) && (elfront->Eleme()->ElementGeometrique().Dimension()==1)))
         // ok surface, -> cote droite ?
        {if (test == 1) // cas du noeud qui bouge
		        { double r=0;
            // le point a t est-il du bon cote (= hors matière) de la surface a t
			         if (elfront->BonCote_t(noeud->Coord1(),r))
				          // oui, deux cas selon que le point initial est dans l'element final ou pas
				          { if (elfront->BonCote_tdt(noeud->Coord1(),r))
					             // oui point initiale hors matière -> il faut maintenant regarder le point final
					             { bool bon_cote = elfront->BonCote_tdt(noeud->Coord2(),r);
                    if (bon_cote && (r > dismini))
						                // oui -> signifie que le noeud est sortie de l'element-> pas de contact
						                { return false;}
					               else // sinon on a le point initial hors matière finale et le point final dans la matière finale
						                { // donc contact sans original
						                  Mtdt = M1; // sauvegarde
						                  return true;
						                };
					             }
				            else // non le point initial est dans  la matiere final
					             { Coordonnee B = (noeud->Coord2() + noeud->Coord1()) / 2.; // point milieu entre t et tdt
					               double rl = (M1 - B).Norme();
					               if (rl <= V.Norme()/2.)
						                // le point est acceptable, on change les coordonnees a tdt du noeud
						                { // contact sans original
                        Mtdt = M1; // sauvegarde
							                 return true;
                						}
					               else
//             return false; // il n'y a pas contact 
// *** non la suite ne marche pas , je le laisse pour mémoire
// si on essaie
						                { // c'est un cas qui peut arriver en particulier dans le cas de la pénalisation 
						                  // on regarde si la position finale est dans l'élément fini associé à la frontière
							                 Element & elem = *(elfront->PtEI()); // l'element qui contiend la frontiere
							                 // on regarde si le noeud esclave est dans la boite d'encombrement de l'élément
							                 if (elem.In_boite_emcombrement_elem(noeud->Coord2()))
								                  // on teste alors plus précisemment, si la position finale est dans l'élément
								                  // il y a contact
								                  { int cas = elem.Interne_tdt(noeud->Coord2());
                            if (Permet_affichage() > 5)
                             {cout << "\n --  ElContact::Contact: ";
                              cout << " Interne_tdt (particulier 1 surface) ? cas= " << cas <<", ";
                              this->Affiche(1);
                             };
                            if (cas  > 0) // point interne ou sur la frontière
									                     { Mtdt = M1; // sauvegarde
									                      return true;
									                     }
								                    else // sinon il n'y a définitivement pas de contact
									                    { return false;};  
                          }
							                 else // si on n'est pas dans la boite d'encombrement  --> pas de contact   
								                  {return false;};
						                };
             
					             };
				          }
			         else // non point initial a t est du mauvais cote, donc dans la matière à t
      
// *** non la suite ne marche pas , je le laisse pour mémoire
// si on essaie
				         { // c'est un cas qui peut arriver en particulier dans le cas de la pénalisation 
				           // on regarde si la position finale est dans l'élément fini associé à la frontière
				           Element & elem = *(elfront->PtEI()); // l'element qui contiend la frontiere
				           // on regarde si le noeud esclave est dans la boite d'encombrement de l'élément
				           if (elem.In_boite_emcombrement_elem(noeud->Coord2()))
					            // on teste alors plus précisemment, si la position finale est dans l'élément
					            // il y a contact
					            { int cas = elem.Interne_tdt(noeud->Coord2());
                   if (Permet_affichage() > 5)
                    {cout << "\n --  ElContact::Contact: ";
                     cout << " Interne_tdt (particulier 2 surface) ? cas= " << cas <<", ";
                     this->Affiche(1);
                    };
                   if (cas > 0) // point interne ou sur la frontière
						              { Mtdt = M1; // sauvegarde
                      return true;
                    }
					              else // sinon il n'y a définitivement pas de contact
						              {return false;};  
					            }
				           else // si on n'est pas dans la boite d'encombrement  --> pas de contact   
					            {return false;};
				         };

            // en fait si on arrive ici, et qu'il y a vraiment contact, cela signifie qu'il y a une erreur
            // car le point à t ne devrait pas être dans la matière à t, car cela veut dire qu'il y avait déjà
            // contact à t, ce contact n'a pas été détecté !!
          
///              {return false;};
		        }// -- fin du test == 1
	        else if (test == 0) // cas rien ne bouge ni noeud ni frontière
          { double r=0;
            // le point final = initial, est-il dans la matiere finale
		          if (!elfront->BonCote_tdt(noeud->Coord2(),r))
              // oui -> contact
              { // contact sans original
                Mtdt = M1; // sauvegarde
                return true;
              }             
		          else // non -> pas de contact, 
              return false;         
          }
         else
          // cas du noeud fixe, et de la frontiere qui bouge  
          { double r=0;
            // le point a t est-il du bon cote de la surface a t
            if (elfront->BonCote_t(noeud->Coord1(),r))
			           // le point final = initial, est-il dans la matiere finale
			           { bool bon_cote = elfront->BonCote_tdt(noeud->Coord2(),r);
                if (!bon_cote || (r <= dismini))
                  // oui -> contact
                  { // contact sans original
                    Mtdt = M1; // sauvegarde
                    return true;
			               }             
                else // non -> pas de contact, 
                return false;
			           }          
            else // non
              return false; // on etait a l'arriere d'une facette
          }
        }
      else if ( (ParaGlob::Dimension() == 3) && (ParaGlob::AxiSymetrie()) && (elfront->Eleme()->ElementGeometrique().Dimension()==1))
       // cas en dimension 3, axisymétrique
       // on examine le contact en fonction de la trajectoire
        {if (test == 1) // cas du noeud qui bouge
          { double r = 0.;
            // le point a t est-il du bon cote (= hors matière) de la ligne a t
            if (elfront->BonCote_t(noeud->Coord1(),r))
              // oui, deux cas selon que le point initial est dans l'element final ou pas
              { if (elfront->BonCote_tdt(noeud->Coord1(),r))
                  // oui point initiale hors matière -> il faut maintenant regarder le point final
                  { bool bon_cote = elfront->BonCote_tdt(noeud->Coord2(),r);
                    if (bon_cote && (r > dismini))
                      // oui -> signifie que le noeud est sortie de l'element-> pas de contact
                      { return false;}
                    else // sinon on a le point initial hors matière finale et le point final dans la matière finale
                      { // donc contact sans original
                        Mtdt = M1; // sauvegarde
                        return true;
                      };
                  }
                else // non le point initial est dans  la matiere final
                  { // on regarde si le point final est aussi dans la matière
                    bool bon_cote = elfront->BonCote_tdt(noeud->Coord1(),r);
                    if (bon_cote && (r > dismini))
                      // oui -> signifie que le noeud est sortie de l'element-> pas de contact
                      { return false;}
                    else // sinon on a le point initial et le point final sont dans la matière finale
                      { // donc contact sans original
                        Mtdt = M1; // sauvegarde
                        return true;
                      };

                  // --- la suite je ne sais pas pourquoi il faut la garder ??
                  
                  
                  
                  
                  
                  Coordonnee B = (noeud->Coord2() + noeud->Coord1()) / 2.; // point milieu entre t et tdt
                    double rl = (M1 - B).Norme();
                    if (rl <= V.Norme()/2.)
                      // le point est acceptable, on change les coordonnees a tdt du noeud
                      { Mtdt = M1; // sauvegarde
                        return true;
                      }
                    else
//             return false; // il n'y a pas contact
// *** non la suite ne marche pas , je le laisse pour mémoire
// si on essaie
                      { // c'est un cas qui peut arriver en particulier dans le cas de la pénalisation
                      
                      
                      
                        // on regarde si la position finale est dans l'élément fini associé à la frontière
                        Element & elem = *(elfront->PtEI()); // l'element qui contiend la frontiere
                        // on regarde si le noeud esclave est dans la boite d'encombrement de l'élément
                        if (elem.In_boite_emcombrement_elem(noeud->Coord2()))
                          // on teste alors plus précisemment, si la position finale est dans l'élément
                          // il y a contact
                          { int cas = elem.Interne_tdt(noeud->Coord2());
                            if (Permet_affichage() > 5)
                             {cout << "\n --  ElContact::Contact: ";
                              cout << " Interne_tdt (particulier 3 surface) ? cas= " << cas <<", ";
                              this->Affiche(1);
                             };
                            if (cas > 0) // point interne ou sur la frontière
                              { Mtdt = M1; // sauvegarde
                               return true;
                              }
                            else // sinon il n'y a définitivement pas de contact
                             { return false;};
                          }
                        else // si on n'est pas dans la boite d'encombrement  --> pas de contact
                          {return false;};
                      };
             
                  };
              }
            else // non point initial a t est du mauvais cote, donc dans la matière à t
      
// *** non la suite ne marche pas , je le laisse pour mémoire
// si on essaie
             { // c'est un cas qui peut arriver en particulier dans le cas de la pénalisation
               // on regarde si la position finale est dans l'élément fini associé à la frontière
               Element & elem = *(elfront->PtEI()); // l'element qui contiend la frontiere
               // on regarde si le noeud esclave est dans la boite d'encombrement de l'élément
               if (elem.In_boite_emcombrement_elem(noeud->Coord2()))
                 // on teste alors plus précisemment, si la position finale est dans l'élément
                 // il y a contact
                 { int cas = elem.Interne_tdt(noeud->Coord2());
                   if (Permet_affichage() > 5)
                    {cout << "\n --  ElContact::Contact: ";
                     cout << " Interne_tdt (particulier 4 surface) ? cas= " << cas <<", ";
                     this->Affiche(1);
                    };
                   if (cas > 0) // point interne ou sur la frontière
                    { Mtdt = M1; // sauvegarde
                      return true;
                    }
                   else // sinon il n'y a définitivement pas de contact
                    {return false;};
                 }
               else // si on n'est pas dans la boite d'encombrement  --> pas de contact
                 {return false;};
             };

            // en fait si on arrive ici, et qu'il y a vraiment contact, cela signifie qu'il y a une erreur
            // car le point à t ne devrait pas être dans la matière à t, car cela veut dire qu'il y avait déjà
            // contact à t, ce contact n'a pas été détecté !!
          
///              {return false;};
          }// -- fin du test == 1
         else if (test == 0) // cas rien ne bouge ni noeud ni frontière
          { // le point final = initial, est-il dans la matiere finale
            double r=0.;
            bool bon_cote = elfront->BonCote_tdt(noeud->Coord2(),r);
            if (!bon_cote || (r <= dismini))
              // oui -> contact
              { // contact sans original
                Mtdt = M1; // sauvegarde
                return true;
              }
            else // non -> pas de contact,
              return false;
          }
         else
          // cas du noeud fixe, et de la frontiere qui bouge
          { // le point a t est-il du bon cote de la surface a t
            double r=0.;
            bool bon_cote = elfront->BonCote_t(noeud->Coord1(),r);
            if (!bon_cote || (r <= dismini))
              // le point final = initial, est-il dans la matiere finale
              { bool bon_cote2 = elfront->BonCote_tdt(noeud->Coord2(),r);
                if (!bon_cote2 || (r <= dismini))
                  // oui -> contact
                  { // contact sans original
                    Mtdt = M1; // sauvegarde
                    return true;
                  }
                else // non -> pas de contact,
                return false;
              }
            else // non
              return false; // on etait a l'arriere d'une facette
          }
        }
      //---------------------
      else // cas en 3D avec un point et une ligne
       { // on regarde si la position finale est dans l'élément fini associé à la frontière
	        Element & elem = *(elfront->PtEI()); // l'element qui contiend la frontiere
	        // on regarde si le noeud esclave est dans la boite d'encombrement de l'élément
         // ou de tous les éléments qui contiennent également les noeuds de la frontière
         // car on peut passer à travers l'élément de la frontière et être dans un autre élément
         // et c'est valide
         if (elem.In_boite_emcombrement_elem(noeud->Coord2()))
		        // on teste alors plus précisemment, si la position finale est dans l'élément
		        // il y a contact
		        { int cas = elem.Interne_tdt(noeud->Coord2());
            if (Permet_affichage() > 5)
             {cout << "\n --  ElContact::Contact: ";
              cout << " Interne_tdt (particulier 5 surface) ? cas= " << cas <<", ";
              this->Affiche(1);
             };
            if (elem.Interne_tdt(noeud->Coord2()) > 0) // point interne ou sur la frontière
			          { Mtdt = M1; // sauvegarde
				           return true;
			          }
            else // sinon il n'y a définitivement pas de contact
			          {return false;};  
          }
		       else // si on n'est pas dans la boite d'encombrement  --> pas de contact   
          {return false;};
       };
	   }
//=========	
   else // le point n'est pas sur la frontiere
	   {return  false; };	
 };
    
// calcul de la trajectoire a prendre en compte pour le contact
//
// ---- a) cas ou à t il n'y avait pas de contact, ou que l'on n'a pas fait de projection sur
//         la surface (cas du contact cinématique)
// == 4 cas :
//       1) le noeud bouge, dans ce cas la trajectoire est determinee
//          par la variation de la position du noeud
//       2) le noeud est immobile, mais la frontiere bouge, la trajectoire est determine
//          par une parallele a la variation moyenne de la frontiere (variation de G)
//       4) est une variante du 2), cas où l'on a une rotation autour de G, dans ce cas
//          on prend comme trajectoire le maxi des déplacements de noeud
//       3) rien ne bouge, on utilise la normale au point de reference de l'element de frontiere
//          pour calculer la trajectoire .
//  dans le cas 3 la variable test = 0 , sinon elle vaut 1 pour le cas 1 et 2 pour le cas 2 , 4 pour le cas 4
//
// ---- b) cas ou à t on était déjà en contact avec projection sur la surface 
//         la trajectoire est alors systématiquement la direction de la dernière normale,
//         retour : test=5

Coordonnee  ElContact::Trajectoire(int & test)
 { // ---> constitution de la droite de projection, et initialisation de V
   Coordonnee V = noeud->Coord2() - noeud->Coord1(); // trajectoire du noeud de t a tdt
   ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
//   //debug
//            cout << "\n debug ElContact::Trajectoire: "; this->Affiche(1);
//            cout << "\n Vinitial= "<< V
//                 << "\n type de frontiere: " << (elfront->Eleme())->TypeFrontiere();
//            if (actif)
//              cout << "\n actif " << actif;
//   // fin debug
   // ------- cas où le contact était déjà effectif à t et que l'on impose la projection sur la surface -------
//   if ((actif > 1) && (ParaGlob::param->ParaAlgoControleActifs().ContactType() == 1))
   if ((actif > 1) // en fait c'est vrai quelque soit le type de contact
       || cas_collant // dans le cas collant on projette systématiquement de manière normale
                       // sinon on peut avoir des différences importante entre l'intersection
                       // et le point collé qui ne sont pas bien prise en compte,
                       // même si le noeud est sur la facette
      )
     { // dans ce cas la trajectoire considérée est la normale, cela veut dire que
       // systématiquement, ce que l'on fait
	      // c'est de calculer selon la normale, le point sur la surface
       // recherche de la normale ,traitement suivant les different cas
       // dimensionnement des variables de tangence (plan ou droite)
       int dim = noeud->Dimension(); // dimension des points
       
       // ----- on traite les cas particuliers des éléments d'angle mort
       // ---- on regarde si on veut une normale lissée
       //     ou s'il s'agit d'angle mort
       if ((normale_lisser) || (elfront->Angle_mort()))
         {// dans le cas d'une normale lissée on va utiliser une interpolation des normales aux noeuds
          ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
          // recup et actualisation du dernier plan tangent, coordonnées locale etc.
          // dimensionnement des variables de tangence
          Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
//          //debug
//                   cout << "\n debug ElContact::Trajectoire: ";
//                   cout << "\n elfro.DernierTangent " << flush;
//          // fin debug
          elfro.DernierTangent(dr,pl,indic,false);
          // récup des fonctions phi
          const Vecteur& phii = elfro.Phi();
          // on parcours les noeuds de la frontière
          // retourne le tableau de noeud en lecture lecture/écriture
          Tableau <Noeud *>& tNfront = elfro.TabNoeud();
          int nbNfront = tNfront.Taille();
          Coordonnee Nnoe(dim); // variable inter
          for (int inoe = 1;inoe <= nbNfront;inoe++)
           { const TypeQuelconque& tiq = tNfront(inoe)->Grandeur_quelconque(N_FRONT_t);
             const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
             Nnoe = gr.ConteneurCoordonnee_const();
             N += phii(inoe)*Nnoe;
           };
          double norme = N.Norme();
//          //debug
//                   cout << "\n debug ElContact::Trajectoire: ";
//                   cout << "\n norme= "<< norme;
//                   cout << " elfro.TypeFrontiere()= "<< elfro.TypeFrontiere() << flush;
//          // fin debug

          // dans le cas improbable, mais qui arrive, que norme = 0:
          // cas par exemple d'un élément dont la somme des  normales aux == 0
          // on se rabat sur le calcul classique
          if (Dabs(norme) < ConstMath::petit)
           {// dans le cas où il s'agit d'une frontière ligne ou surface
            Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
            // constitution d'un plan tangent ou de la droite tangente au point de reference
            elfro.TangentRef(dr,pl,indic);
            V= Calcul_Normale(dim,pl,dr, indic);
           }
          else // sinon c'est ok
           {N /= norme;
            V=N;
           };

          test = type_trajectoire_tdt = 5;
         }
      // --- on continue avec les éléments qui ne sont pas d'angle mort et sans lissage
       else
        {
//        //debug
//                 cout << "\n debug ElContact::Trajectoire: ";
//                 cout << " elfro.TypeFrontiere()= "<< elfro.TypeFrontiere() << flush;
//        // fin debug
          if (elfro.TypeFrontiere() == "FrontPointF")
           {// cas où la frontière est un point on utilise l'élément linéaire
            // dont c'est la frontière,  pour calculer la direction de la normale
            // qui correspond à la tangent à l'élément
            Element * elemf = elfront->PtEI(); // récup de l'élément fini
            ElemGeomC0& elem_geom = elemf->ElementGeometrique(); // recup de l'élément géométrique correspondant
            // on peut répondre directement s'il s'agit d'un élément de type segment
            if (elem_geom.TypeGeometrie() == SEGMENT)
             { // on doit avoir un élément de mécanique
               ElemMeca & elMeca = (ElemMeca &) (*elemf);
               // il n'y a qu'un seul noeud d'ou le pointeur du noeud frontière
               Noeud * noe = elfro.TabNoeud()(1);
               // récupération de la base locales au noeud noe, pour le temps: temps
               const BaseB & baseB = elMeca.Gib_elemeca(TEMPS_tdt,noe);
               // le premier vecteur est le vecteur tangent que l'on considère normal à la frontière
               V = baseB(1).Coor();
               // on renseigne le point frontière du vecteur tangent qui lui est attribué, ceci pour l'utilisation
               // de la fonction: duboncote
             }
            else
             { cout << "\n *** cas actuellement non implante: element point, frontiere d'une geometrie autre que segment"
                    << "\n ElContact::Trajectoire(int & test)";
               Sortie(1);
             }
            }
          else
           { // dans le cas où il s'agit d'une frontière ligne ou surface
             Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
             // constitution d'un plan tangent ou de la droite tangente au point de reference
             elfro.TangentRef(dr,pl,indic);
             V = Calcul_Normale(dim,pl,dr, indic);
           };
          test = type_trajectoire_tdt = 5;
	   	   };
     }
   else
    { // ------- cas où soit le contact n'était pas effectif auparavant,
      // soit il n'y a pas de projection sur la facette -------
      // on cherche dans quel cas de contact on se trouve
      // 3 cas : le noeud bouge ou seul la frontiere bouge ou tout est immobile
      double L = V.Norme(); // la distance de t a tdt
      if (L <= ConstMath::trespetit)
      // cas ou le vecteur est nul c-est a dire, cas ou le noeud n'a pas bouge, on regarde si
      // l'element maitre a bougé
        { const Tableau <Noeud *>& tabnoeud = elfro.TabNoeud();
          int tail = tabnoeud.Taille();
          Coordonnee V11(V); // variation totale
          // en fait on peut regarder le mvt du centre de gravité  qui à NBE près est
          // le V11 qui suit
          for (int it=1;it<= tail;it++)
            V11 += tabnoeud(it)->Coord2() - tabnoeud(it)->Coord1();
          double L11 = V11.Norme();
          bool rien_ne_bouge = true;
          if (L11 <= ConstMath::trespetit)
           { // le centre de gravité ne bouge pas, donc soit la frontière est immobile,
             // soit il y a rotation autour de G,
             // dans ce cas on essaie le maxi du déplacement de chaque noeud
             Coordonnee V12(V);double L12=0.; V11=V; L11= V11.Norme();// init
             for (int itt=2;itt<=tail;itt++)
              { V12=tabnoeud(itt)->Coord2() - tabnoeud(itt)->Coord1(); L12=V12.Norme();
                if (L12 > L11) {V11=V12;L11=L12; };
              };
             // on reteste la nouvelle longueur L11 finale
             if (L11 > ConstMath::trespetit)
                // cas où il y a rotation autour de G, et V11 est le maxi des déplacements
                {  test = type_trajectoire_tdt = 4;   V=V11;rien_ne_bouge=false;};
              // sinon on sort des deux if (on sera à: // -<<>>-- ) et on en conclue que vraiment rien ne bouge !
           }
          else
           // cas où le centre de gravité bouge et V11 est le déplacement de G
           {   test = type_trajectoire_tdt = 2;   V = V11;rien_ne_bouge=false;};

          // -<<>>-- là vraiment il n'y a rien qui bouge sur la frontière
          // --on prendra comme direction la normale a l'element
          if (rien_ne_bouge)
          // recherche de la normale ,traitement suivant les different cas
           { // dimensionnement des variables de tangence (plan ou droite)
             int dim = noeud->Dimension(); // dimension des points
             
             // on traite les cas particuliers des éléments d'angle mort
             if (elfront->Angle_mort())
              { if (elfro.TypeFrontiere() == "FrontPointF")
                 {// s'il s'agit d'un point il n'y a pas d'existance de normale
   //               // on remplace la notion de normale par le vecteur qui joint le noeud à la position
   //               // du point frontière
   //               V = elfro.Ref()- noeud->Coord2();
                
                  // on utilise la normale au noeud
                  const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t);
                  const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
                  V = gr.ConteneurCoordonnee_const();
                 }
                else
                 { cout << "\n *** cas element Angle_mort actuellement non implante: , "
                        << elfro.TypeFrontiere()
                        << "\n ElContact::Trajectoire(int & test)";
                   Sortie(1);
                 };
              }
             else
              { if (elfro.TypeFrontiere() == "FrontPointF")
                 {// cas où la frontière est un point on utilise l'élément linéaire dont c'est la frontière
                  // pour calculer la direction de la normale qui correspond à la tangent à l'élément
                  Element * elemf = elfront->PtEI(); // récup de l'élément fini
                  ElemGeomC0& elem_geom = elemf->ElementGeometrique(); // recup de l'élément géométrique correspondant
                  // on peut répondre directement s'il s'agit d'un élément de type segment
                  if (elem_geom.TypeGeometrie() == SEGMENT)
                   { // on doit avoir un élément de mécanique
                     ElemMeca & elMeca = (ElemMeca &) (*elemf);
                     // il n'y a qu'un seul noeud d'ou le pointeur du noeud frontière
                     Noeud * noe = elfro.TabNoeud()(1);
                     // récupération de la base locales au noeud noe, pour le temps: temps
                     const BaseB & baseB = elMeca.Gib_elemeca(TEMPS_tdt,noe);
                     // le premier vecteur est le vecteur tangent que l'on considère normal à la frontière
                     V = baseB(1).Coor();
                     // on renseigne le point frontière du vecteur tangent qui lui est attribué, ceci pour l'utilisation
                     // de la fonction: duboncote
                    }
                  else
                   { cout << "\n **** cas actuellement non implante: element point, frontiere d'une geometrie autre que segment"
                          << "\n ElContact::Trajectoire(int & test)";
                     Sortie(1);
                    }
                  }
                else
                 { // dans le cas où il s'agit d'une frontière ligne ou surface
                   Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
                   // constitution d'un plan tangent ou de la droite tangente au point de reference
                   elfro.TangentRef(dr,pl,indic);
                   V = Calcul_Normale(dim,pl,dr, indic);
                 };
              };
             test = type_trajectoire_tdt = 0;
           };
        }
      else
      // lorsque l'on arrive ici, c'est donc que le noeud bouge, et c'est son déplacement
      //  qui est retenue comme trajectoire
       {test = type_trajectoire_tdt = 1;};
    };
//   //debug
//            cout << "\n debug ElContact::Trajectoire: V= " << V << flush;
//   // fin debug
   // retour
   return V;
 };
    
    
// calcul de l'Intersection de la trajectoire du noeud definit par le vecteur V
// avec l'element frontiere
// -> ramene les coordonnees du noeud projete
// dans le cas où il n'y a pas d'intersection, on ramène un point de dimension nulle
Coordonnee ElContact::Intersection( const Coordonnee& V,bool init)
 { ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
   int contactType = ElContact::Recup_et_mise_a_jour_type_contact();
   // dimensionnement des variables de tangence (plan ou droite)
   int dim = noeud->Dimension(); // dimension des points
   Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
   
   // on traite tout d'abord les cas particuliers des éléments d'angle mort
   if (elfront->Angle_mort())
    { if (elfro.TypeFrontiere() == "FrontPointF")
       {// s'il s'agit d'un point il n'y a pas d'existance de normale
        const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t);
        const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
        N = gr.ConteneurCoordonnee_const();
        // l'intersection est considérée toujours valide et est égale aux coordonnées
        // du point frontière
        return elfro.Ref();
       }
      else
       { cout << "\n *** cas element Angle_mort actuellement non implante: , "
              << elfro.TypeFrontiere()
              << "\n ElContact::Intersection(...";
         Sortie(1);
       };

    }
   else
    { // cas des éléments de contact autres que ceux qui décrivent les angles morts
      if (init)
        // constitution d'un plan tangent ou de la droite tangente au point de reference
        elfro.TangentRef(dr,pl,indic);
      else
        // on reactualise le dernier plan tangent
        elfro.DernierTangent(dr,pl,indic);
//   //---- debug
//   if (noeud->Num_noeud() == 29)
//    {cout << "\n debug : ElContact::Intersection(.. ";
//     // on vérifie la perpendicularité
//     double perpen = dr.VecDroite() * V;
//     cout << "\n diff normale : "<< perpen  ;
//     cout << " N= ";
//     cout <<endl;
//    };
//   //---- fin debug
      // --- def de la droite de trajectoire
      // dans le cas où le précédent point est libre (pas de projection systématique) on retient le point à t
      // dans le cas où il y a une projection systématique on retient le point à tdt, car il s'agit de la simple
      // projection du point final sur la facette
      Coordonnee pointDeBase = noeud->Coord1();
      if (((actif > 1)
           &&  (   (contactType == 1) || (contactType == 3)  )
          )
          // dans le cas du type 4 on fait toujours une projection normale
          || ((actif > 0) && (contactType == 4))
         )
         pointDeBase = noeud->Coord2();
//   //---- debug
//   if (noeud->Num_noeud() == 29)
//    {cout << "\n debug : ElContact::Intersection(.. ";
//     cout << "\n pointDeBase : "<< pointDeBase  ;
//    };
//   //---- fin debug
      Droite droite(pointDeBase,V); // la droite
      // dans le cas très particulier où l'élément de frontière maître est un point:
      if (indic == 0)
       { // il n'y a qu'un seul noeud d'ou le pointeur du noeud frontière
         Noeud * noe = elfro.TabNoeud()(1);
         Coordonnee point_frontiere=noe->Coord2();
         // on regarde si le point frontière est sur la trajectoire
         if (droite.Appartient_a_la_droite(point_frontiere))
          { // le point frontière est sur la trajectoire, l'intersection est donc au noeud frontière
            // les fonctions d'interpolations du premier point en contact
            if (init)
              {phi_theta_0 = elfro.Phi();};
            return point_frontiere;
          }
         else
          { // il n'y a pas d'intersection, on ramène un point de dimension nulle
            Coordonnee retour(0);
            return retour;
          }
       };
        
//   //---- debug
//   if (noeud->Num_noeud() == 29)
//    {cout << "\n debug : ElContact::Intersection(.. ";
//     // on vérifie la perpendicularité
//     double perpen = dr.VecDroite() * droite.VecDroite();
//     cout << " diff normale2 : "<< perpen  ;
//     cout <<endl;
//    };
//   //---- fin debug
      // cas général : algorithme de recherche d'intersection
      //  - initialisation
      int intp = 0; // indicateur de recherche de pt d'inter
      Coordonnee M(dim); // le point cherché d'intersection avec la surface
      Coordonnee M1; // le point courant de la surface
      int compteur = 0; double nm = 0.; // grandeurs de contrôle
      double max_compteur = ParaGlob::param->ParaAlgoControleActifs().Nb_boucle_newton_position_sur_frontiere();
      // distance maxi entre le point d'inter du plan et le point correspondant de la surface
      double dismini = ParaGlob::param->ParaAlgoControleActifs().Precision_pt_sur_front();
      //---------------------------------------
      // recherche du point d'intersection M
      //---------------------------------------
      do
       {
        // calcul de l'intersection
        if (indic == 1)
        // cas d'une droite
          { int ret = dr.Intersection(droite,M);
            if (( ret == 0) || ( ret == -1)) // pas d'intersection, on essaie deux autres points
             { for (int i=1;i<= 2; i++)
                 { elfro.AutreTangent(dr,pl,indic); // on suppose indic du meme type
                   ret = dr.Intersection(droite,M);
                   if (ret > 0 ) { intp = 1; break;};
                 };
               // dernier essai, on perturbe 2 fois la trajectoire
               for (int i=1;i<=2; i++)
                { Coordonnee newtra = droite.VecDroite() + 0.1 * droite.UneNormale();
   //// debug
   //{                  double d = newtra.Norme();
   //              //    if (d <= ConstMath::trespetit)
   //                  if (d <= ConstMath::petit) // *****pour test
   //                   { cout << "\n*** (1) debug ElContact::Intersection : la norme du vecteur est trop petite !";
   //                     cout <<"\nnorme = " << d << " newtra= " << newtra ;
   //                    }
   //}
   //// fin debug
                  droite.Change_vect(newtra);
                  for (int i=1;i<= 2; i++)
                    { elfro.AutreTangent(dr,pl,indic); // on suppose indic du meme type
                      ret = dr.Intersection(droite,M);
                      if (ret > 0) { intp = 1; break;};
                    };
                  if (intp) break;
                };
             }
            else
             intp = 1; // on a trouve un point
          }
        else
        // cas d'un plan
          { int ret = pl.Intersection(droite,M);
            if (( ret == 0) || ( ret == -1)) // pas d'intersection, on essaie trois autres points
             { for (int i=1;i<= 3; i++)
                 { elfro.AutreTangent(dr,pl,indic); // on suppose indic du meme type
                   ret = pl.Intersection(droite,M);
                   if (ret > 0) { intp = 1; break;};
                 };
               // dernier essai, on perturbe 2 fois la trajectoire
               for (int i=1;i<=2; i++)
                { Coordonnee newtra = droite.VecDroite() + 0.1 * droite.UneNormale();
                  droite.Change_vect(newtra);
   //  // debug
   //  {                  double d = newtra.Norme();
   //                //    if (d <= ConstMath::trespetit)
   //                    if (d <= ConstMath::petit) // *****pour test
   //                     { cout << "\n*** (2) debug ElContact::Intersection : la norme du vecteur est trop petite !";
   //                       cout <<"\nnorme = " << d << " newtra= " << newtra ;
   //                      }
   //  }
   //  // fin debug
                  for (int i=1;i<= 2; i++)
                    { elfro.AutreTangent(dr,pl,indic); // on suppose indic du meme type
                      ret = dr.Intersection(droite,M);
                      if (ret > 0) { intp = 1; break;};
                    };
                  if (intp) break;
                };
             }
            else
             intp = 1; // on a trouve un point
          };
        if (intp)
         // calcul du point M1 de la surface correspondant  a M
         { elfro.Tangent(M,M1,dr,pl,indic);
           nm = (M1 - M).Norme();
         }
        else
          { // on n'arrive pas à calculer un point d'intersection
            if (Permet_affichage() >= 7)
               {cout << "\n remarque , pas d'intersection "
                     << "  entre la trajectoire du noeud esclave "
                     << noeud->Num_noeud() <<" mail= " << noeud->Num_Mail()
                     << " avec frontiere:"
                     << elfront->Num_frontiere()
                     << " de l'element " << elfront->PtEI()->Num_elt() << " du maillage "
                     << elfront->PtEI()->Num_maillage() << " \n ";
               };
            if (Permet_affichage() >= 9)
             { noeud->Affiche();
               elfro.Affiche();
               cout << "ElContact::Intersection(bool init, bool forcee)"<< endl;
             };
            // on considère qu'il n'y a pas d'intersection (donc on n'en tient pas compte) et retour
            // évite au programme de planter globalement alors qu'il n'y a peut-être qu'un pb local
            Coordonnee retour(0);
            return retour;
          };
        // test d'arret si l'on est en dehors de la boite d'encombrement
        // sauf dans le cas d'un contact collant car pour le collant: l'intersection est toujours valide
        // ensuite il y aura éventuellement un choix entre différents éléments de contact
   //	  cout << "\n"; M1.Affiche();elfront->Encom_mini_FR().Affiche();elfront->Encom_maxi_FR().Affiche();cout << endl;
        if (!cas_collant)
         if (!(elfront->In_boite_emcombrement_front(M1)))
            break;
        compteur++;
       }
      while ((compteur <= max_compteur) && (nm >= dismini)) ;
      // retour du point d'intersection
   ////////debug
   // if (noeud->Num_noeud()==84)
   //     {
   //      cout << "\n debug ElContact::Intersection : noeud: " << noeud->Num_noeud() <<" mail= " << noeud->Num_Mail()
   //           << " avec frontiere:" << elfront->Num_frontiere()
   //                  << " de l'element " << elfront->PtEI()->Num_elt() << " du maillage "
   //					             << elfront->PtEI()->Num_maillage();
   //      cout <<"\n Noe_atdt= " << noeud->Coord2() << "\n point projete= " << M1 << endl;
   //     };
   ////// fin debug
      // les fonctions d'interpolations du premier point en contact
      if (init)
        {phi_theta_0 = elfro.Phi();};
      // retour du point d'intersection
      return M1;
   };
   // si on arrive ici on considère qu'il n'y a pas d'intersection
   Coordonnee a(0);
   return a;
 };
    
// construction de la condition lineaire de contact
// nb_assemb : indique le numéro d'assemblage correspondant    
Condilineaire ElContact::ConditionLi(int nb_assemb)
  { ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
    // recup de la dimension 
    int dim = noeud->Dimension();
    // --- def des tableaux pour la condition lineaire
    Tableau <Noeud *>& tabnoeudf = elfro.TabNoeud(); // le tableau de noeuds de l'élément frontière
    
    //     on cherche le nombre de ddl attaché à l'ensemble des noeuds concerné par le contact
    int nbnc = tabnoeudf.Taille(); int tail=0; // init
    // puis on boucle sur les noeuds en cumulant les ddl actif pour le cas d'assemblage donné
    // on considère que le cas d'assemblage est celui de ddl de déplacement uniquement
    if (cas_solide == 1)
       // on ajoute les ddl du noeud projectile, qui est ici le seul libre, le reste est solide
       {tail += noeud->NB_ddl_actif_casAssemb (nb_assemb);}
    else // pour les autres cas on considère pour l'instant tous les ddl
       {for (int ine=1;ine<=nbnc;ine++) tail += tabnoeudf(ine)->NB_ddl_actif_casAssemb (nb_assemb);
        //on rajoute les ddl du noeud projectile
        tail += noeud->NB_ddl_actif_casAssemb (nb_assemb);
       };
    Tableau<int> pt(tail); //définition du tableau de pointeurs des ddls concernées par la condition de contact
    Vecteur val(tail);     // coefs de la condition lineaire
    Tableau <Enum_ddl> t_enu(tail); // le tableau des énumérés correspondant aux ddl de la CLL
    int nb_noeud_CLL= 1; // init par défaut
    if (cas_solide != 1) nb_noeud_CLL += nbnc;
    Tableau <Noeud *> t_noeudCLL(nb_noeud_CLL); // le tableau de noeuds de la CLL
    // --- fin def des tableaux pour la condition lineaire
   
    // recup du dernier plan tangent (ou droite)
    // dimensionnement des variables de tangence
    Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
    elfro.DernierTangent(dr,pl,indic);
    // calcul de la normale en fonction de differente conditions
    Coordonnee N( Calcul_Normale(dim,pl,dr,indic)) ;
   

    // remplissage des tableaux
    // on commence par le noeud esclave
///    int numX = noeud->Existe(X1) -1;
///    int pointeur = noeud->PosiAssemb(nb_assemb) + numX;// pointeur d'assemblage
    int numX = noeud->Position_ddl(X1,nb_assemb)-1;
    #ifdef MISE_AU_POINT	 	 
	   if ( numX == -2 )
			  { cout << "\nErreur : ddl X1 inexistant  pour le cas de charge " << nb_assemb << '\n'
			         << "ElContact::ConditionLi(int nb_assemb) 1\n";
			    Sortie(1);
     };
    #endif  
              
    int pointeur = noeud->Pointeur_assemblage(X1,nb_assemb) - 1;// pointeur d'assemblage - 1
    #ifdef MISE_AU_POINT	 	 
	   if ( pointeur == -2 )
			  { cout << "\nErreur : ddl X1 inexistant  pour le cas de charge " << nb_assemb << '\n'
			         << "ElContact::ConditionLi(int nb_assemb) 2\n";
			    Sortie(1);
     };
    #endif  
    int itab ; // indice des tableaux  
    // le premier indice correspondra a la direction bloquee, il faut s'assurer qu'il ne 
    // soit pas nul et qu'il ne correspond pas a une direction deja bloquee par les conditions
    // limites de depart
    // on boucle sur la dim pour trouver le coef le plus grand correspondant a une direction
    // non bloquee
    int iti;double max = 0;int indi = 0;
    int posi = Id_nom_ddl("X1") -1; // position du ddl X1
    for (iti=1; iti<= dim; iti++)
        { if (!noeud->Ddl_fixe(Enum_ddl(posi+iti))) // cas d'un ddl libre
            if (Dabs(N(iti)) > max)
              { max = Dabs(N(iti)); 
                indi = iti;
              };
        };
    if (indi ==0)
        { // cas ou soit les trois ddl sont bloque, soit le vecteur normal = 0
          cout << "\n erreur dans les conditions lineaire de contact, " 
               << " soit les trois ddl de deplacement du noeud esclave  sont tous bloque, "
               << " soit le vecteur normal a la surface maitre = 0 \n";
          noeud->Affiche();
          cout << "\n et la normal : " ; N.Affiche();
          cout << "ElContact::ConditionLi()" << endl;			 
//			 cout << "\n fixite= " << noeud->Ddl_noeud_t(Enum_ddl(posi+1)).Fixe() 
//			      << "  " << noeud->Ddl_noeud_t(Enum_ddl(posi+1)).Fixe() << endl;
          Sortie (1);
        };
    int contactType = ElContact::Recup_et_mise_a_jour_type_contact();
    // dans le cas d'un contact de type 1
    if (    (contactType == 1)|| (contactType == 3) )
      {// on indique que la nouvelle direction "indi" est bloque, ceci permet
       // pour une autre condition de contact, de ne pas reprendre cette direction
       // --> ce marquage est effacer apres resolution, par la fonction efface blocage
       noeud->Change_fixe(Enum_ddl(posi+indi),true);
      };
    t_noeudCLL(1) = noeud;
    
    // constitution de la condition lineaire
    int indi_t = indi; // récup de l'indice du ddl du noeud esclave bloqué (de 1 à dim)
                       // dans une variable intermédiaire, car elle sera modifié par la boucle suivante
                       // et on en a besoin dans la constitution de la condition linéaire-> condi_retour      
    for (itab=1;itab<= dim;itab++)
       { val(itab) = N(indi_t); // val(1) correspond à N(indi), puis les autres suivent
         pt(itab) = pointeur + indi_t; // ici il s'agit du pointeur d'assemblage du noeud bloqué
         t_enu(itab) = Enum_ddl(posi+indi_t);
         indi_t ++;
         if (indi_t > dim) indi_t = 1;
       };
    if (cas_solide != 1) // le  cas == 1 correspond au cas d'un contact avec un solide
                         // dans ce cas, les noeuds maîtres ne font pas partie de la CLL
      {// cas des noeuds de l'element maitre
       Vecteur phi = elfro.Phi();
       for (int nn=1; nn<=tabnoeudf.Taille(); nn++)
          { Noeud & noe = *(tabnoeudf(nn));
            t_noeudCLL(nn+1)=&noe;
            // récup du pointeur d'assemblage des noeuds de l'élément maitre -1
            pointeur = noe.Pointeur_assemblage(X1,nb_assemb) -1;
            #ifdef MISE_AU_POINT
            if ( pointeur == -2 )
             { cout << "\nErreur : ddl X1 inexistant  pour le cas de charge " << nb_assemb
                    <<  '\n'
                    << "ElContact::ConditionLi(int nb_assemb) 3\n";
               Sortie(1);
             };
            #endif
            for (int i=1;i<= dim;i++,itab++)
                { val(itab) = - N(i) * phi(nn);
                  pt(itab) = pointeur + i ; // ici il s'agit du pointeur d'assemblage des noeuds maitres
                  t_enu(itab) = Enum_ddl(posi+i);
                };
          };
      };
    Condilineaire condi_retour(t_enu,pt,val,0.,numX+indi,t_noeudCLL);                   
    return condi_retour;
   };
    
// actualisation de  la projection du noeud esclave en fonction de la position de l'element
// maitre frontiere. Lorsque le noeud change d'element fontiere, on change l'element
// frontiere de l'element de contact en consequence
// retourne:
// 0 : dans le cas ou on ne trouve pas d'intersection, cas d'une intersection  qui n'est pas calculable
// 1  : contact ok, sur l'élément
// -1 : contact ok, mais hors de l'élément (là, on a parcouru tous les mitoyens et aucun ne convient
//      du coup on utilise l'élément initial pour la projection
// 2  : changement de frontière, et contact ok sur la nouvelle frontière
// NB: sauf pour le cas 0, dans tous les autres cas, les grandeurs locales sont correctement calculées
//     c-a-d : la projection, la normale
// en fonction de la méthode de contact, le noeud est ramené éventuellement sur la frontière
int ElContact::Actualisation()
{ // comme le point est déjà en contact, son déplacement est grosso modo // à la surface
  // donc on ne peut pas utiliser sa trajectoire et faire une intersection avec la surface, cela ne veut rien dire
  // donc on cherche la projection du point sur la surface. Pour cela on calcule l'intersection d'une trajectoire qui a
  // comme direction, la normale à la surface
  // dans le cas d'un point avec ligne c'est idem
  int nb_change_frontiere_max = ParaGlob::param->ParaAlgoControleActifs().NB_change_frontiere_max();
 #ifdef MISE_AU_POINT
 if (!actif)
     {cout << "\n ElContact::Actualisation(., actif=  "<< actif
           << " l'element n'est pas actif, ce n'est pas normal, Actualisation ne s'applique "
           << " qu'a un element actif ! ";
     };
 #endif

  if (Permet_affichage() > 4)
      {cout << "\n   -- ElContact::Actualisation ";
       this->Affiche(1);
      };
      
  int intersec = 0; // init: pour l'instant pas d'intersection
  int passage_par_angle_mort=0; // pour le passage d'angle mort à normal
  Coordonnee M1;
  Coordonnee N_sauve = N; // pour un retour si l'intersection ne peut-être calculé
  // cependant le cas des angles morts est particulier, on commence par le traiter
  if (elfront->Angle_mort())
   { ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
     if (elfro.TypeFrontiere() == "FrontPointF")
      {// s'il s'agit d'un point il n'y a pas d'existance de normale
//       // on remplace la notion de normale par le vecteur qui joint le noeud à la position
//       // du point frontière
//       Coordonnee V = elfro.Ref()- noeud->Coord2();
       const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t);
       const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
       N = gr.ConteneurCoordonnee_const();

       // l'intersection avec la frontière est systématiquement le point frontière
       // même si on change de frontière (c'est pour simplifier le calcul de M1
       //  qui sera mieux calculé au prochain appel)
       M1 = elfro.Ref();
       
       intersec=-1; // init à minima

       // la seule chose que l'on peut faire pour valider le fait qu'il y a contact ou pas
       // c'est de tester si le noeud est à l'intérieur de l'élément d'angle mort ou pas
       // c'est plus long qu'un test de surface mais il n'y a pas vraiment d'autre solution
       int cas = elfront->PtEI()->Interne_tdt(noeud->Coord2());
       if (Permet_affichage() > 5)
         {cout << " Frontiere: FrontPointF -> test Interne_tdt ? cas= " << cas;
         // retour : =0 le point est externe, =1 le point est interne ,
         //          = 2 le point est sur la frontière à la précision près
          switch (cas)
            {case 0 : cout << " (le point est externe a l'element) "; break;
             case 1 : cout << " (le point est interne a l'element) "; break;
             case 2 : cout << " (le point est sur la frontiere a la prec pres) "; break;
             default:
               cout << "\n **** erreur le retour cas = "<< cas
                    << " de la methode elfront->PtEI()->Interne_tdt(... "
                    << " n'est pas normal on ne peut pas continuer ! ";
               Sortie(1);
               break;
            };
         };
       
       // pour éviter les flip flop, si on a eu trop de changement de frontière
       // on valide de manière arbitraire le choix actuel
       if (nb_change_frontiere > abs(nb_change_frontiere_max))
         {cas = 1;
          if (Permet_affichage() > 5)
            {cout << "\n     >>  nb_change_frontiere: " << nb_change_frontiere
                  << " > nb_change_frontiere_max ("<<nb_change_frontiere_max << ") "
                  << " on impose cas = 1 (c-a-d le point est interne a l'element)" ;
            };
         };
       if (!cas)
         // le point n'est pas dans l'élément on essaie les éléments voisins
         // on itere sur les elements voisins
         {const Tableau <Front*>* ta = elfront->TabMitoyen();
          if (ta != NULL) // cas où il y a des éléments voisins !
            { int ta_taille_et1 = 1+ta->Taille();
              for (int i=1; i < ta_taille_et1; i++)
                 {
               #ifdef UTILISATION_MPI
                 // dans le cas // on ne continue que si la frontière est relative à un élément associé au cpu
                 Element & elem11 = *((*ta)(i)->PtEI()); // l'element qui contiend la frontiere
                 // on ne continue que si l'élément est concerné
                 if (ParaGlob::param->Element_concerner(elem11.Num_maillage(),elem11.Num_elt_const()) )
                  {
               #endif
                  if (Permet_affichage() > 5)
                    cout << "\n      >> test de l'element voisin "<<i <<" : " ;
                   cas = (*ta)(i)->PtEI()->Interne_tdt(noeud->Coord2());
                   if (Permet_affichage() > 5)
                     {cout << " Interne_tdt sur mitoyen d'angle mort ? cas= " << cas ;
                      switch (cas)
                        {case 0 : cout << " (le point est externe a l'element) "; break;
                         case 1 : cout << " (le point est interne a l'element) "; break;
                         case 2 : cout << " (le point est sur la frontiere a la prec pres) "; break;
                         default:
                           cout << "\n **** erreur le retour cas = "<< cas
                                << " de la methode elfront->PtEI()->Interne_tdt(... "
                                << " n'est pas normal on ne peut pas continuer ! ";
                           Sortie(1);
                           break;
                        };
                     };
                   if (cas)
                     // on a trouvé un élément mitoyen qui contient le point
                     {// --- changement d'element frontiere
                      if (Permet_affichage() > 2)
                        {cout << "\n       >> change front. du noeud  "<< noeud->Num_noeud() <<" mail= " << noeud->Num_Mail()
                              << " ==>  ";
                         elfront->Affiche(1);
                         if (Permet_affichage() > 4)
                           cout << " ( mitoyen: " << i << " )" ;
                        };

                      // on commence par sauvegarder (*ta)(i) en créant une copie
                      Front*  newelfront = new Front(*((*ta)(i)));
////----- debug
//if (newelfront->NumUnique()==0)
//  {cout << "\ndebug Elcontact::Actualisation , ***** NumUnique est nulle !! erreur ";
//   cout << "\n (*ta)(i)->Unique() = "<<(*ta)(i)->NumUnique();
//  }
////----- fin debug
//                      Libere(); // on supprime le elfront actuel, donc le tableau ta
                      // on ne peut pas supprimer le elfront, car il est également pointé par elfront_t
                      // du coup ce sera fait dans la méthode ElContact::TdtversT()
                      elfront = newelfront; // on récupére le nouvelle élément créé
                      if (Permet_affichage() > 2)
                        {cout << " --> nevez ";
                         elfront->Affiche(1);
                        };
                      Construction_TabNoeud(); // on reconstruit le tableau de noeud global
                      nb_change_frontiere++; // on incrémente
                      intersec = 2;
                      passage_par_angle_mort=2;
                     };
                   if (intersec == 2)
                     break;
               #ifdef UTILISATION_MPI
                  };
               #endif
                 };
            };
         }
       else // l'intersection est ok
         {intersec = 1; };
      }
     else
      { cout << "\n *** cas element Angle_mort actuellement non implante: , "
             << elfro.TypeFrontiere()
             << "\n ElContact::Actualisation()";
        Sortie(1);
      };
   }
   
  // ---- la suite concerne des éléments autres que ceux d'angle mort
  else
   {// recup des dernières différentes informations calculées au niveau de la cinématique de contact
    Coordonnee M_impact,N; // le point d'impact, la normale
    Vecteur phii;
    RecupInfo(N,M_impact,phii,false);
 
    // def de la pseudo-trajectoire a  utiliser: ici ce sera une trajectoire normale à la surface maître
    Coordonnee V = N;
    // calcul de l'intersection de la pseudo-trajectoire avec la frontiere
    M1 = Intersection(V,false);
   
    // dans le cas ou M1 a une dimension nulle, signifie qu'il n'y a pas d'intersection,
    // on peut alors dire qu'il ne peut pas y avoir contact donc retour
    if (M1.Dimension() == 0)
       {intersec = 0;}
    else // sinon on peut continuer
       {intersec=-1; // init à minima
        // maintenant on regarde si l'intersection est acceptable
        // - cote surface
        ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
        // ici on recherche avec une précision donnée
        double extra_in_surface = ParaGlob::param->ParaAlgoControleActifs().PointInternePrecThetaiInterne();
//        int contactType =
        ElContact::Recup_et_mise_a_jour_type_contact();
        if ((elfro.InSurf(extra_in_surface))
            || (nb_change_frontiere > abs(nb_change_frontiere_max)) // pour limiter les flip-flop
           )
         // l'intersection est ok pour l'element
         // cote droite  : contrairement a Intersection() , on oblige ici éventuellement en fonction du modèle de contact
         // le point a etre sur la  surface, car a priori il y a toujours contact
          { if (nb_change_frontiere_max > 0) intersec = 1; else { intersec = 0;};
          }
        else
          // l'intersection est en dehors de la surface ou courbe maitre, on essaie une projection sur
          //  l'element voisin
          // on itere sur les elements voisins
          {const Tableau <Front*>* ta = elfront->TabMitoyen();
           if (ta != NULL) // cas où il y a des éléments voisins !
             {for (int i=1; i<= ta->Taille(); i++)
               { bool mitoyen_valider = false;
               
               #ifdef UTILISATION_MPI
                 // dans le cas // on ne continue que si la frontière est relative à un élément associé au cpu
                 Element & elem11 = *((*ta)(i)->PtEI()); // l'element qui contiend la frontiere
                 // on ne continue que si l'élément est concerné
                 if (ParaGlob::param->Element_concerner(elem11.Num_maillage(),elem11.Num_elt_const()) )
                  {
               #endif
                 if (Permet_affichage() > 5)
                   cout << "\n      >> test de l'element voisin "<<i <<" : " ;
                 // le cas des angles morts est particulier, on commence par le traiter
                 if ((*ta)(i)->Angle_mort())
                  { ElFrontiere & elfroto = *((*ta)(i)->Eleme()); // pour commodite
                    if (elfroto.TypeFrontiere() == "FrontPointF")
                     {// s'il s'agit d'un point il n'y a pas d'existance de normale
                      const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t);
                      const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
                      V = gr.ConteneurCoordonnee_const();

                      // la seule chose que l'on peut faire pour valider le fait qu'il y a contact ou pas
                      // c'est de tester si le noeud est à l'intérieur de l'élément d'angle mort ou pas
                      // c'est plus long qu'un test de surface mais il n'y a pas vraiment d'autre solution
                      int cas = (*ta)(i)->PtEI()->Interne_tdt(noeud->Coord2());
                      if (Permet_affichage() > 5)
                        {cout << " Interne_tdt sur mitoyen d'angle mort ? cas= " << cas ;
                         switch (cas)
                           {case 0 : cout << " (le point est externe a l'element) "; break;
                            case 1 : cout << " (le point est interne a l'element) "; break;
                            case 2 : cout << " (le point est sur la frontiere a la prec pres) "; break;
                            default:
                              cout << "\n **** erreur le retour cas = "<< cas
                                   << " de la methode elfront->PtEI()->Interne_tdt(... "
                                   << " n'est pas normal on ne peut pas continuer ! ";
                              Sortie(1);
                              break;
                           };
                        };

                      if (cas) // si true, le point est interne
                        {M1 = elfroto.Ref();
                         N = V;
                         mitoyen_valider = true;
                        };
                     }
                    else
                     { cout << "\n *** cas element Angle_mort actuellement non implante: , "
                            << elfroto.TypeFrontiere()
                            << "\n ElContact::Actualisation()";
                       Sortie(1);
                     };
                  }
                 else // ---- la suite concerne des éléments autres que ceux d'angle mort
                  {ElContact elct((*ta)(i),noeud,fct_nD_contact); // on cree un element contact pour utiliser ses methodes

                   Coordonnee M1_new;
                   M1_new = elct.Intersection(V,true);
                   // ici on recherche avec une précision donnée
                   double extra_in_surface = ParaGlob::param->ParaAlgoControleActifs().PointInternePrecThetaiInterne();
                   // dans le cas ou M1 a une dimension nulle, signifie qu'il n'y a pas d'intersection,
                   // on n'examine que le cas où la dimension est différente de 0
                   if (M1_new.Dimension() != 0)
                    { if ((elct.Elfront()->Eleme()->InSurf(extra_in_surface))
                          || (nb_change_frontiere > nb_change_frontiere_max) // pour éviter les flip-flop
                         )
                       // on a trouve une bonne intersection
                       { mitoyen_valider = true;M1 = M1_new;
                       if (Permet_affichage() > 5)
                         cout << " intersection ok ";
                       };
                    };
                  };
               #ifdef UTILISATION_MPI
                  };
               #endif

                 // si le mitoyen est validé on change de Front pour l'élément de contact
                 if (mitoyen_valider)
                  { // --- changement d'element frontiere
                    if (Permet_affichage() > 3)
                      {cout << "\n   change front. du noeud  "<< noeud->Num_noeud() <<" mail= " << noeud->Num_Mail()
                            << " ==  ";
                       elfront->Affiche(1);
                       if (Permet_affichage() > 4)
                         cout << " ( mitoyen: " << i << " )" ;
                      };
                                           
                    // on commence par sauvegarder (*ta)(i)
                    Front*  newelfront = new Front(*((*ta)(i)));
                    
//                   //----- debug
//                   if (newelfront->NumUnique()==0)
//                     {cout << "\ndebug Elcontact::Actualisation , ***** NumUnique est nulle !! erreur 2 ";
//                      cout << "\n (*ta)(i)->Unique() = "<<(*ta)(i)->NumUnique();
//                     }
//                   //----- fin debug

                    // Libere(); // on supprime le elfront actuel, donc le tableau ta
                    // on ne peut pas supprimer le elfront, car il est également pointé par elfront_t
                    // du coup ce sera fait dans la méthode ElContact::TdtversT()
                    elfront = newelfront; // on récupére le nouvelle élément créé
                    if (Permet_affichage() > 2)
                      {cout << " --> nevez ";
                       elfront->Affiche(1);
                      };
                    Construction_TabNoeud(); // on reconstruit le tableau de noeud global
                    nb_change_frontiere++;
                    intersec = 2;
                    // si le nouvel élément de frontière est un angle mort, alors la normale
                    // a été calculée dans Intersection, donc pas la peine de la recalculer
                    break;
                };
               };
             };
          };
      };
  }; // fin du test sur les front d'angle mort

  if (intersec != 0)
    {Mtdt = M1; // sauvegarde du point d'intersection
     int contactType = ElContact::Recup_et_mise_a_jour_type_contact();
     // gestion du retour
     if ((contactType == 1)|| (contactType == 3)|| (contactType == 4))
       { M_noeud_tdt_avant_projection = noeud->Coord2();
         // on change les coordonnées s'ils sont libres sinon on laisse inchangé
         int dim = M_noeud_tdt_avant_projection.Dimension();
         switch (dim)
          { case 3: if (!(noeud->Ddl_fixe(X3))) noeud->Change_val_tdt(X3,M1(3));
            case 2: if (!(noeud->Ddl_fixe(X2))) noeud->Change_val_tdt(X2,M1(2));
            case 1: if (!(noeud->Ddl_fixe(X1))) noeud->Change_val_tdt(X1,M1(1));
            default:
                break;
          };
         if (Permet_affichage() > 6)
           {cout << "\n   actualisation du noeud esclave "
                 << noeud->Num_noeud() <<" mail= " << noeud->Num_Mail()
                 << " ancienne position " << M_noeud_tdt_avant_projection
                 << "  nouvelle position " << noeud->Coord2();
           };
       };
    }
  else // cas où l'intersection n'est pas calculable
     // on remet les choses à l'état initiale
    {Mtdt = Mt;
     N = N_sauve;
    }
    
  #ifdef MISE_AU_POINT
  if (Permet_affichage() > 4)
      {cout << ": noeud  "<< noeud->Num_noeud() <<" mail= " << noeud->Num_Mail()
            << ", intersec=  "<< intersec;
       if (intersec == 0) {cout << " (pas d'intersection trouvé)";}
       else if (intersec == 1) {cout << " (contact ok sur l'élément)";}
       else if (intersec == -1){cout << " (contact ok mais hors de l'élément)";}
       else if (intersec == 2){cout << " (changement de frontière, et contact ok sur la nouvelle frontière)";};
      };
  #endif
    
  // retour du type d'intersection
  return intersec;
};


// test et met à jour le compteur de décollage du noeud
// si le noeud decolle ou non en fonction de la force de reaction
// ramene 1:  s'il decolle
//        0: s'il ne décolle pas

bool ElContact::Decol()
{  // recup de la dimension
   int dim = noeud->Dimension();
   // recup du dernier plan tangent (ou droite)
   // dimensionnement des variables de tangence
   Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
   // calcul du dernier plan tangent
   (elfront->Eleme())->DernierTangent(dr,pl,indic);
   // recherche de la normale ,traitement suivant les different cas
   N = Calcul_Normale(dim,pl,dr, indic); // elle pointe vers l'extérieur de l'élément
	                                                // donc vers l'intérieur de l'élément esclave
   // def du decollement
   double r = - force_contact * N;
   
	 
//----------- debug--------------
/*{// a priori on suppose que le contact est collant --> calcul des efforts de contact
Coordonnee  Noe_atdt = noeud->Coord2(); // récup de la position actuelle du noeud projectile
 // recup du dernier des différentes informations
 Coordonnee M_impact,N; // le point d'impact, la normale
 Vecteur phii;
 RecupInfo(N,M_impact,phii);

// calcul de la pénétration normale
Coordonnee deltaX = Noe_atdt - M_impact;
double gap= N * deltaX;

cout << "\n noeud: " << noeud->Num_noeud() << ": gap= " << gap << ", F_N= " << r << endl;
}*/	 
//----------- fin debug--------------
	 
   double borne_regularisation;
   if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL)
        {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation
                                              ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation
                                              ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation
                                              ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation);
        }
   else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();}
	  // NB: nb_decol_tdt est également mis à 0 lors de l'inactivation ou de l'activation de l'élément: cf. Met_Inactif() et Met_actif()
	  if (r <= 0) 
	  // ici cas classique on dit qu'il n'y a pas décolement car la réaction est négative
      {nb_decol_tdt=0;return false;} // cas d'une compression
   else // cas d'une traction
	     {nb_decol_tdt++;  // cas d'une traction
		     int typededecolement = ParaGlob::param->ParaAlgoControleActifs().TypeDeDecolement();
		     if (!typededecolement) // dans le cas classique (==0), on regarde uniquement la réaction
         { if (nb_decol_tdt >= ParaGlob::param->ParaAlgoControleActifs().NbDecolAutorise())
                { if (Permet_affichage() > 4)
                    {cout << "\n cas d'une reaction negative:"<<-r<<",  nb decollement  ("
                          << nb_decol_tdt << "): >= le nombre de decollement ("
                          << ParaGlob::param->ParaAlgoControleActifs().NbDecolAutorise()
                          << ") tolere ";
                     Affiche(1);
                    };
                         
                  if (Permet_affichage() >= 7)
                   {cout << "\n force de contact: "<< force_contact;
                    cout << "\n normale: "<< N;
                    cout << "\n intensite du contact (r = - force_contact * N;) "
                         << r ;
                   };
                  return true;
                }
			        else {return false;};
			      } 
		     else // sinon on dit qu'il y a décollement que si la non pénétration, est supérieur à une certaine valeur
		        // cependant, le facteur de pénalisation sera quasi-nulle ou même strictement nulle donc finalement, la force de réaction 
			       // n'aura pas d'importance sur l'équilibre, par contre l'élément restera ce qui évite de le reconstruire 
			       // et peut-être évite les oscillation
		      {if (gap_tdt > Dabs(borne_regularisation) * typededecolement)
               {if (nb_decol_tdt >= ParaGlob::param->ParaAlgoControleActifs().NbDecolAutorise())
                     { if (Permet_affichage() > 4)
                        {cout << "\n cas d'une reaction negative:"<<-r<<", et "
                              << "d'un gap > Dabs(borne_regularisation) ("
                              << Dabs(borne_regularisation) << " * typededecolement ("
                              << typededecolement << "),  nb decollement ("
                              << nb_decol_tdt << "): >= le nombre de decollement ("
                              << ParaGlob::param->ParaAlgoControleActifs().NbDecolAutorise()
                              << ") tolere ";
                         Affiche(1);
                        };
                       if (Permet_affichage() >= 7)
                        {cout << "\n gap_tdt= " << gap_tdt
                              << ", force de contact: "<< force_contact;
                         cout << "\n normale: "<< N;
                         cout << "\n intensite du contact (r = - force_contact * N;) "
                              << r ;
                        };
                       return true;
                     }
			             else {return false;}; // pas encore de libération
               }
		       else {return false; }; // sinon on maintient le contact, donc pas de décollement	 
		      };
		    };
};

// change force: permet de changer la valeur de la force
// utile quand la force est calculée en dehors de l'élément de contact
void ElContact::Change_force(const Coordonnee& force)
 {force_contact = force;
  // -- on met à jour également la partie force tangentielle et normale individuellement
  // recup du dernier plan tangent (ou droite)
  // dimensionnement des variables de tangence
  int dim = noeud->Dimension();
  Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence
  ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite
  elfro.DernierTangent(dr,pl,indic);
  // calcul de la normale en fonction de differente conditions
  Coordonnee N( Calcul_Normale(dim,pl,dr,indic));
  
  F_N_max = (force_contact*N); // sauvegarde
  F_T_max = sqrt(force_contact.Norme()-F_N_max*F_N_max);

 
 };

// cas d'une méthode avec pénalisation: calcul éventuel d'un pas de temps idéal,
// permettant de limiter la pénétration
// si oui retour de la valeur delta_t proposé
// sinon dans tous les autres cas retour de 0.
// le calcul se fait en fonction du pas de temps courant et de la pénétration
// donc nécessite que le contact ait déjà été étudié
double ElContact::Pas_de_temps_ideal()const
{ double delta_optimal=ConstMath::tresgrand;
  // récupération du temps
  double delta_t = ParaAlgoControle::Variables_de_temps().IncreTempsCourant();
  // limitation au niveau du gap
  double borne_regularisation;
  if (fct_nD_contact.fct_nD_penetration_borne_regularisation != NULL)
       {borne_regularisation = Valeur_fct_nD(fct_nD_contact.fct_nD_penetration_borne_regularisation
                                             ,ElContact::Fct_nD_contact::tqi_fct_nD_penetration_borne_regularisation
                                             ,ElContact::Fct_nD_contact::tqi_const_fct_nD_penetration_borne_regularisation
                                             ,ElContact::Fct_nD_contact::t_num_ordre_fct_nD_penetration_borne_regularisation);
       }
  else {borne_regularisation = ParaGlob::param->ParaAlgoControleActifs().Penetration_borne_regularisation();}
  int contactType = ElContact::Recup_et_mise_a_jour_type_contact();
  if (actif  && (borne_regularisation < 0)
      && ((contactType == 2) || (contactType == 4)|| (contactType == 41)|| (contactType == 42))
     )
   {if (gap_tdt < - Dabs(borne_regularisation))
     delta_optimal = delta_t * Dabs(borne_regularisation / gap_tdt);
  };
  // retour du temps proposé
  if (delta_optimal==ConstMath::tresgrand)
    // cela veut dire qu'il n'a pas été modifié, on le met à 0 pour indiquer que l'on a rien
  // n'a proposer
  delta_optimal = 0.;
  return delta_optimal;
};


// permet de modifier le contact entre collant ou non suivant "change"
void ElContact::Change_contact_collant(bool change)
 { // dans une première mouture c'est assez simple mais on définit une méthode
   // interne car je pressent que cela pourrait se complexifier !
   if (change )
    {cas_collant = 1;}
   else
    {cas_collant = 0;};
 };


// récupération des ddl ou des grandeurs actives de tdt vers t
void ElContact::TdtversT()
 { Mt=Mtdt;gap_t = gap_tdt; nb_pene_t = nb_pene_tdt;
   dep_T_t = dep_T_tdt;
   mult_pene_t=mult_pene_tdt;
   mult_tang_t=mult_tang_tdt;
   type_trajectoire_t = type_trajectoire_tdt;
   nb_decol_t = nb_decol_tdt;
   nb_posi_esclave_stocker_t = nb_posi_esclave_stocker_tdt;
   indice_stockage_glissant_t = indice_stockage_glissant_tdt;
   force_contact_t=force_contact;
   tabForce_cont_t=tabForce_cont;
   F_N_max_t = F_N_max; F_T_max_t=F_T_max;
   penalisation_t = penalisation;
   penalisation_tangentielle_t = penalisation_tangentielle;
   // on regarde si on a changé de elfront
   if (elfront != elfront_t)
     // cas où on a changé de frontière
     { delete elfront_t; // suppression de la frontière ancienne
       elfront_t = elfront; // mise à jour
     };
   noeud_t = noeud;
   tabNoeud_t = tabNoeud;
   nb_change_frontiere = 0; // initialisé à chaque fois
 };
// actualisation des ddl et des grandeurs actives de t vers tdt
void ElContact::TversTdt()
 { if (Mt.Dimension() != 0)
      Mtdt=Mt;
   gap_tdt = gap_t; nb_pene_tdt = nb_pene_t;
   dep_T_tdt = dep_T_t;
   mult_pene_tdt=mult_pene_t;
   mult_tang_tdt=mult_tang_t;
   type_trajectoire_tdt = type_trajectoire_t;
   nb_decol_tdt = nb_decol_t;
   nb_posi_esclave_stocker_tdt = nb_posi_esclave_stocker_t;
   indice_stockage_glissant_tdt = indice_stockage_glissant_t;
   force_contact=force_contact_t;
   tabForce_cont=tabForce_cont_t;
   F_N_max = F_N_max_t; F_T_max=F_T_max_t;
   penalisation = penalisation_t;
   penalisation_tangentielle = penalisation_tangentielle_t;
   // on regarde si on a changé de elfront
   if (elfront != elfront_t)
     // cas où on a changé de frontière
     { delete elfront; // suppression de la nouvelle frontière
       elfront = elfront_t; // retour sur l'ancienne frontière
     };
   noeud = noeud_t;
   tabNoeud = tabNoeud_t;
   nb_change_frontiere = 0; // initialisé à chaque fois

 };

//----- lecture écriture base info -----
// lecture base info
void ElContact::Lec_base_info_ElContact(ifstream& ent)
 {  // tout d'abord les données propres à l'élément
    string toto;
////debug
// if (noeud->Num_noeud()==56)
//  {
//   cout << "\n debug : ElContact::Lec_base_info_ElContact: noeud 56"<< endl;
//  };
//// fin debug
    ent >> toto >> toto >> actif
        >> toto >> Mt
        >> toto >> energie_frottement
        >> toto >> energie_penalisation
        >> toto >> penalisation_t >> toto >> penalisation_tangentielle_t ;
    penalisation = penalisation_t; penalisation_tangentielle = penalisation_tangentielle_t;
    ent >> toto >> num_zone_contact >> toto >> cas_collant
        >> toto >> phi_theta_0;
    ent >> toto >> force_contact >> toto >>  tabForce_cont
        >> toto >> F_N_max >> toto >> F_T_max
        >> toto >> gap_t >> toto >> nb_decol_t >> toto >> nb_pene_t
        >> toto >> dep_T_t
        >> toto >> mult_pene_t
        >> toto >> mult_tang_t
        >> toto >> type_trajectoire_t ;
    ent >> toto >> N >> toto >> dep_tangentiel;
    // on identifie les grandeurs à tdt
    Mtdt = Mt;
    nb_decol_tdt = nb_decol_t;
    gap_tdt=gap_t;dep_T_tdt = dep_T_t;
    nb_pene_tdt = nb_pene_t;
    mult_pene_tdt = mult_pene_t;
    mult_tang_tdt = mult_tang_t;
    type_trajectoire_tdt = type_trajectoire_t;
  
    // puis les données frontière
    elfront->Lecture_base_info_front(ent);
    elfront_t=elfront;
    noeud_t = noeud;
    tabNoeud_t = tabNoeud;
 
  };
//   écriture base info
// ici comme les contacts peuvent apparaitre à tout moment, cela n'a pas de sens de faire des sauvegardes partielles
// donc on sauvegarde toujours tout
// on ne sauvegarde que ce qui est particulier à l'élément (en particulier ce qui est relatif au noeud esclave et à l'éléments finis
//    n'est évidemment pas sauvegardé ici)
void ElContact::Ecri_base_info_ElContact(ofstream& sort)
  { // écriture du type
    sort << "\n Elcontact: ";
    // les aspects liés au traitement du contact
    sort << " actif= " << actif << " Mtdt= " << Mt << " ener_frot= " << energie_frottement
         << " ener_penal= " << energie_penalisation
         << " penalisation= " << penalisation_t
         << " penalisation_tangentielle= " << penalisation_tangentielle_t;
    sort << "\n num_zone_contact= " << num_zone_contact << " cas_collant= " << cas_collant
         << " phi_theta_0= " << phi_theta_0;
    sort << "\n F_cont= " << force_contact << " F_facette= "<< tabForce_cont
         << " norme_F_N= " << F_N_max << " norme_F_T= " <<  F_T_max
         << "\n gap= " << gap_t
         << " nb_decol= " << nb_decol_t << " nb_pene= " << nb_pene_t
         << " dep_T_t= " << dep_T_t
         << " mult_pene_t= " << mult_pene_t
         << " mult_tang_t= " << mult_tang_t
         << " type_traj= " << type_trajectoire_t ;
    sort << "\n normale: " << N << " dep_tangentiel: " << dep_tangentiel << " ";
    // les données frontière
    elfront->Ecriture_base_info_front(sort);
    sort << " " ;
  };


void ElContact::Libere()
 { //delete (elfront->Eleme());
   delete elfront;
   elfront = elfront_t = NULL;
  };

    // construction du tableau de tous les noeuds, le premier est celui esclave
    // et mise à jour de ddlElement et de list_Ddl_global éventuellement
// et mise à jour des grandeurs liées à elfront: concerne le dimensionnement
// si on change d'elfront -> peut changer le nombre de noeud maître d'où
// il faut redimensionner: tabForce_cont,tabForce_cont_t;
void ElContact::Construction_TabNoeud()
  {  // def du tableau de l'element frontiere
     Tableau <Noeud *>& tabnoeuddd = (elfront->Eleme())->TabNoeud();
     int taile = tabnoeuddd.Taille() +1;
     // def du tableau tabNoeud
     tabNoeud.Change_taille(taile);
     tabNoeud(1) = noeud;
     for (int i=2; i<= taile;i++)
       tabNoeud(i) = tabnoeuddd(i-1);
     // --- on ajoute les ddl de réaction, si jamais cela n'est pas fait
     Enum_ddl enu_reac_1= R_X1; // le premier ddl
     int dim = ParaGlob::Dimension();
     // on crée les tableau de ddl de réaction à ajouter
     Tableau<Ddl> ta(dim);
     int posi = Id_nom_ddl("R_X1") -1;
     for (int i=1;i<=dim;i++)
       ta(i)= Ddl(Enum_ddl(i+posi),0.0,LISIBLE_FIXE);
     // on ajoute les ddl si ils ne sont pas présents
     Ddl_enum_etendu& ddl_reaction_normale=Ddl_enum_etendu::Tab_FN_FT()(1);
     Ddl_enum_etendu&  ddl_reaction_tangentielle=Ddl_enum_etendu::Tab_FN_FT()(2);
     for (int i=1; i<= taile;i++)
        { Noeud& noe = *tabNoeud(i); // pour simplifier
          if (!noe.Existe_ici(enu_reac_1))
             noe.PlusTabDdl(ta);
          // on introduit également spécifiquement la réaction normale et la réaction tangentielle
          if (!noe.Existe_ici_ddlEtendu(ddl_reaction_normale))
            noe.AjoutUnDdl_etendu(ddl_reaction_normale);
          if (!noe.Existe_ici_ddlEtendu(ddl_reaction_tangentielle))
            noe.AjoutUnDdl_etendu(ddl_reaction_tangentielle);
        };
        
     // mise à jour des grandeurs liées à elfront: concerne le dimensionnement
     // si on change d'elfront -> peut changer le nombre de noeud maître d'où
     // il faut redimensionner: tabForce_cont,tabForce_cont_t;
     tabForce_cont.Change_taille(taile-1);
     tabForce_cont_t.Change_taille(taile-1);
     // on init à 0
     for (int i=1; i< taile; i++)
      {tabForce_cont(i).Zero();
       tabForce_cont_t(i).Zero();
      };

     // --- cas du tableau des DdlElement
   
     // en fait tous les éléments du tableau sont identiques et il suffit qu'ils existent
     // on commence par parcourir la liste pour trouver un bon candidat
     list <DdlElement>::iterator ili,ilifin=list_Ddl_global.end();
     bool trouve = false;
     for (ili=list_Ddl_global.begin();ili !=ilifin; ili++)
       if ((*ili).NbNoeud() == taile) // on a trouvé un candidat
         {ddlElement_assemblage = &(*ili); trouve = true;};
     // dans le cas où on n'a pas trouvé, on en cré un
     if (!trouve)
      { int dima = ParaGlob::Dimension();
        // dans le cas où on est en axisymétrie il faut diminuer de 1
        // car on ne fera pas d'assemblage suivant la direction 3
        if (ParaGlob::AxiSymetrie())
          dima--;
        DdlElement  tab_ddl(taile,dima);
        int posi = Id_nom_ddl("X1") -1;
        for (int i =1; i<= dima; i++)
          for (int j=1; j<= taile; j++)
            tab_ddl.Change_Enum(j,i,Enum_ddl(i+posi));
        list_Ddl_global.push_front(tab_ddl);
        ddlElement_assemblage = &(*list_Ddl_global.begin());
      };
     tabNoeud_pour_assemblage = tabNoeud; // cas général bidéformable
     // on appel quand même la mise à jour car elle construit également cas_solide
     // et réduit éventuellement  ddlElement_assemblage
     Mise_a_jour_ddlelement_cas_solide_assemblage();
   };


// méthode statique de modification éventuelle du type de contact utilisé localement
// par exemple pour le type 4 en fonction des itérations
int ElContact::Recup_et_mise_a_jour_type_contact() 
  {int contactType = ParaGlob::param->ParaAlgoControleActifs().ContactType();
   // le contact 4 passe en 2 quand on dépasse un nombre donné d'itération
   if (contactType == 4)
     {if (fct_pilotage_contact4 == NULL)
        {// on récupère le pointeur correspondant à la grandeur COMPTEUR_ITERATION_ALGO_GLOBAL
         const void* pointe =  (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
         if (pointe == NULL)
           { cout << "\n *** pb ElContact::mise_a_jour_type_contact(... avec le cas contact 4  !! "
                  << " la variable globale  "<< Nom_GrandeurGlobale(COMPTEUR_ITERATION_ALGO_GLOBAL)
                  << ", n'est pas disponible, on ne peut pas continuer  "<<endl;
             Sortie(1);
           };
         // sinon c'est ok
         TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
         Grandeur_scalaire_entier& gr = *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
         int nb_changement = 1; // valeur par défaut
         if (*(gr.ConteneurEntier()) > nb_changement) // à partir du 2 ième itérations on passe en 41
           {contactType = 41;}
         else if  (*(gr.ConteneurEntier()) > 2 * nb_changement)  // on passe à 42 après une transition
           {contactType = 42;}
        }
      else // on a un pilotage via une fonction nD
        {Tableau <double> & tava = fct_pilotage_contact4->Valeur_pour_variables_globales();
//cout << "\n debug ElContact::mise_a_jour_type_contac: tava(1)= " << tava(1);
         if (tava(1) == 1)
           {contactType = 41;}
         else if (tava(1) == 2)
           {contactType = 42;}
        };
     };
   // retour
   return contactType;

  };