// FICHIER : Biel_axiQ.cc
// CLASSE : Biel_axiQ
//#include "Debug.h"


// 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 <iostream>
using namespace std;  //introduces namespace std
#include <stdlib.h>
#include "Sortie.h"


#include "Biel_axiQ.h"
#include "FrontPointF.h"
#include "FrontSegLine.h"
#include "TypeConsTens.h"
#include "TypeQuelconqueParticulier.h"

//----------------------------------------------------------------
// def des donnees commune a tous les elements
// la taille n'est pas defini ici car elle depend de la lecture 
//----------------------------------------------------------------
Biel_axiQ::DonneeCommune * Biel_axiQ::doCo = NULL; 
Biel_axiQ::UneFois  Biel_axiQ::unefois; 
Biel_axiQ::NombresConstruire Biel_axiQ::nombre_V; 
Biel_axiQ::ConstrucElementbiel Biel_axiQ::construcElementbiel;

// constructeur définissant les nombres (de noeud, de point d'integ ..)
// utilisé dans la construction des éléments
Biel_axiQ::NombresConstruire::NombresConstruire() 
 { nbne = 3; // le nombre de noeud de l'élément
   nbneA = 3;// le nombre de noeud des aretes
   nbi = 4;//10;  // le nombre de point d'intégration pour le calcul mécanique
   nbiEr = 4;// le nombre de point d'intégration pour le calcul d'erreur
   nbiA = 4;//10; // le nombre de point d'intégration pour le calcul de second membre linéique
   nbiMas = 4; // le nombre de point d'intégration pour le calcul de la matrice masse
   nbiHour = 0; // le nombre de point d'intégration un blocage d'hourglass
 };

// ---- definition du constructeur de la classe conteneur de donnees communes ------------	
Biel_axiQ::DonneeCommune::DonneeCommune (GeomSeg& seg,DdlElement& tab,DdlElement&  tabErr,DdlElement&  tab_Err1Sig,
	                       MetAxisymetrique2D&  met_bie,
	                       Tableau <Vecteur *> & resEr,Mat_pleine& raidEr,
	                       GeomSeg& seEr,Vecteur&  residu_int,Mat_pleine&  raideur_int,
	                       Tableau <Vecteur* > & residus_extN,Tableau <Mat_pleine* >&  raideurs_extN,
	                       Tableau <Vecteur* > & residus_extA,Tableau <Mat_pleine* >&  raideurs_extA,
	                       Mat_pleine&  mat_masse,GeomSeg& seMa,int nbi,GeomSeg* segHourg) :
	          segment(seg),tab_ddl(tab)
           ,segS(segment),point() // pour le second membre on utilise la même interpolation
           ,tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig)
	          ,met_biellette(met_bie)
	          ,matGeom(tab.NbDdl(),tab.NbDdl())
	          ,matInit(tab.NbDdl(),tab.NbDdl())
	          ,d_epsBB(tab.NbDdl()),d_sigHH(tab.NbDdl()),d2_epsBB(nbi)
	          ,resErr(resEr),raidErr(raidEr)
	          ,segmentEr(seEr)
	          ,residu_interne(residu_int),raideur_interne(raideur_int)
	          ,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN)
	          ,residus_externeA(residus_extA),raideurs_externeA(raideurs_extA)
	          ,matrice_masse(mat_masse)
	          ,segmentMas(seMa),segmentHourg(segHourg)
         {
           int nbddl = tab.NbDdl();      
           for (int ni=1;ni<=nbi;ni++)     
            {d2_epsBB(ni).Change_taille(nbddl);
             for (int i1=1; i1<= nbddl; i1++)
              for (int i2=1; i2<= nbddl; i2++)
               d2_epsBB(ni)(i1,i2) = NevezTenseurBB (2);
            };
           int tailledeps = d_epsBB.Taille();
           for (int i=1;i<= tailledeps; i++)
              d_epsBB(i) = NevezTenseurBB (2);
           int tailledsig = d_sigHH.Taille();
           for (int j=1;j<= tailledsig; j++)
              d_sigHH(j) = NevezTenseurHH (2);
         };

Biel_axiQ::DonneeCommune::DonneeCommune(DonneeCommune& a) :
	    segment(a.segment),tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr),tab_Err1Sig11(a.tab_Err1Sig11)
	    ,met_biellette(a.met_biellette),matGeom(a.matGeom),matInit(a.matInit)
	    ,d2_epsBB(a.d2_epsBB),resErr(a.resErr),raidErr(a.raidErr),segmentEr(a.segmentEr)
	    ,d_epsBB(a.d_epsBB),d_sigHH(a.d_sigHH)
	    ,residu_interne(a.residu_interne),raideur_interne(a.raideur_interne)
	    ,residus_externeN(a.residus_externeN),raideurs_externeN(a.raideurs_externeN)
	    ,residus_externeA(a.residus_externeA),raideurs_externeA(a.raideurs_externeA)
	    ,matrice_masse(a.matrice_masse),segmentMas(a.segmentMas),segmentHourg(a.segmentHourg)
	       { int nbddl = d_sigHH.Taille();     
	         int nbi=d2_epsBB.Taille();
	         for (int ni=1;ni<=nbi;ni++)     
             for (int i1=1; i1<= nbddl; i1++)
               for (int i2=1; i2<= nbddl; i2++)
                 d2_epsBB(ni)(i1,i2) = NevezTenseurBB (*(a.d2_epsBB(ni)(i1,i2)));
          int tailledeps = d_epsBB.Taille();
	         for (int i=1;i<= tailledeps; i++)
	            d_epsBB(i) = NevezTenseurBB (*(a.d_epsBB(i)));
          int tailledsig = d_sigHH.Taille();
	         for (int j=1;j<= tailledsig; j++)
	            d_sigHH(j) = NevezTenseurHH (*(d_sigHH(j)));
        };
Biel_axiQ::DonneeCommune::~DonneeCommune()
        { int nbddl = tab_ddl.NbDdl();
          int nbi=d2_epsBB.Taille();
          for (int ni=1;ni<=nbi;ni++)
            for (int i1=1; i1<= nbddl; i1++)
             for (int i2=1; i2<= nbddl; i2++)
               delete d2_epsBB(ni)(i1,i2);
          int tailledeps = d_epsBB.Taille();
          for (int i=1;i<= tailledeps; i++)
            delete d_epsBB(i);
          int tailledsig = d_sigHH.Taille();
          for (int j=1;j<= tailledsig; j++)
           delete d_sigHH(j);
        };

// ---------- fin definition de la  classe conteneur de donnees communes ------------	
// -+-+ definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
Biel_axiQ::UneFois::UneFois () : // constructeur par défaut
	doCoMemb(NULL),CalResPrem_t(0),CalResPrem_tdt(0),CalimpPrem(0),dualSortbiel(0)
	,CalSMlin_t(0),CalSMlin_tdt(0),CalSMRlin(0)
	,CalSMvol_t(0),CalSMvol_tdt(0),CalSMvol(0)
	,CalDynamique(0),CalPt_0_t_tdt(0)
	,nbelem_in_Prog(0)
       {};
Biel_axiQ::UneFois::~UneFois ()
       {	  delete doCoMemb;
       };

// -+-+ fin definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
	            
	             
// ---------- fin definition de la classe conteneur de donnees communes ------------	

Biel_axiQ::Biel_axiQ ()  :
  // Constructeur par defaut
  ElemMeca(),lesPtMecaInt(),donnee_specif()

      {// on intervient seulement à partir du deuxième élément,
       if (unefois.nbelem_in_Prog == 0)
         { unefois.nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
         }
       else 	// sinon on construit     
        { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
          id_interpol=BIE2;  //   donnees de la classe mere
          id_geom=SEG_AXI;      //
          tab_noeud.Change_taille(nombre_V.nbne);
          // vérif des dimensions
          if ( ParaGlob::Dimension() < 3)  // cas d'une dimension ou l'élément est impossible
           { if (ParaGlob::NiveauImpression() >= 7)
               cout << "\nATTENTION -> dimension " << 1 
                    <<", pas de definition d\'elements triangle Biel_axiQ "<< endl;
           }
          else
           {// initialisation par défaut
            doCo = Biel_axiQ::Init ();
            unefois.nbelem_in_Prog++;
           };
        };
      };  

Biel_axiQ::Biel_axiQ (double epai,int num_maill,int num_id):
// Constructeur utile si l'epaisseur de l'element et
// le numero de l'element  sont connus
ElemMeca(num_maill,num_id,BIE2,SEG_AXI),lesPtMecaInt(),donnee_specif(epai)
     {// on intervient seulement à partir du deuxième élément,
      if (unefois.nbelem_in_Prog == 0)
        { unefois.nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
        }
      else 	// sinon on construit     
       { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
         tab_noeud.Change_taille(nombre_V.nbne);
         // vérif des dimensions
          if ( ParaGlob::Dimension() < 3)  // cas d'une dimension ou l'élément est impossible
           { if (ParaGlob::NiveauImpression() >= 7)
               cout << "\nATTENTION -> dimension " << 1 
                    <<", pas de definition d\'elements triangle Biel_axiQ "<< endl;
           }
          else
           {// initialisation
            doCo = Biel_axiQ::Init (donnee_specif);
            unefois.nbelem_in_Prog++;
           };
       };
     };

// Constructeur fonction  d'un numero de maillage et d'identification  
Biel_axiQ::Biel_axiQ (int num_maill,int num_id) :
  ElemMeca(num_maill,num_id,BIE2,SEG_AXI),lesPtMecaInt(),donnee_specif()
    {// on intervient seulement à partir du deuxième élément,
     if (unefois.nbelem_in_Prog == 0)
       { unefois.nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
       }
     else 	// sinon on construit     
      { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
        tab_noeud.Change_taille(nombre_V.nbne);
        // vérif des dimensions
        if ( ParaGlob::Dimension() < 3)  // cas d'une dimension ou l'élément est impossible
           { if (ParaGlob::NiveauImpression() >= 7)
               cout << "\nATTENTION -> dimension " << 1 
                    <<", pas de definition d\'elements triangle Biel_axiQ "<< endl;
           }
          else
           {// initialisation par défaut
            doCo = Biel_axiQ::Init ();
            unefois.nbelem_in_Prog++;
           };
      };
    };

Biel_axiQ::Biel_axiQ (double epai,int num_maill,int num_id,const Tableau<Noeud *>& tab):
// Constructeur utile si l'épaisseur de l'element, le numero de l'element et
// le tableau des noeuds sont connus
ElemMeca(num_maill,num_id,tab,BIE2,SEG_AXI),lesPtMecaInt(),donnee_specif(epai)
    {// on intervient seulement à partir du deuxième élément,
     if (unefois.nbelem_in_Prog == 0)
       { unefois.nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
        }
     else 	// sinon on construit     
      { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
        if (tab_noeud.Taille() != nombre_V.nbne)
         {cout << "\n erreur de dimensionnement du tableau de noeud \n";
          cout << " Biel_axiQ::Biel_axiQ (double sect,int num_id,const Tableau<Noeud *>& tab)\n";
          Sortie (1);
         };
        // vérif des dimensions
        if ( ParaGlob::Dimension() < 3)  // cas d'une dimension ou l'élément est impossible
           { if (ParaGlob::NiveauImpression() >= 7)
               cout << "\nATTENTION -> dimension " << 1 
                    <<", pas de definition d\'elements triangle Biel_axiQ "<< endl;
           }
          else
           {// construction du tableau de ddl spécifique à l'élément pour ses
            ConstTabDdl();
            // initialisation
            bool sans_init_noeud = true;
            doCo = Biel_axiQ::Init (donnee_specif,sans_init_noeud);
            unefois.nbelem_in_Prog++;
           };
      };
    };

Biel_axiQ::Biel_axiQ (const Biel_axiQ& biel) :
 ElemMeca (biel),lesPtMecaInt(biel.lesPtMecaInt),donnee_specif(biel.donnee_specif) 
    // Constructeur de copie
    { if (unefois.nbelem_in_Prog == 1)
       { cout << "\n **** erreur pour l'element Biel_axiQ, le constructeur de copie ne doit pas etre utilise"
              << " pour le premier element !! " << endl;
         Sortie (1);
       }
     else 	    
      { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
        // initialisation
        unefois.nbelem_in_Prog++;
      };	 
    };	

Biel_axiQ::~Biel_axiQ ()
    // Destructeur
    { LibereTenseur();
      if (unefois.nbelem_in_Prog != 0) 
         unefois.nbelem_in_Prog--;
      // defArete pointe sur la même grandeur que def donc pour éviter
      // une destruction lors du destructeur générale dans elemméca on le met à null
      if (defArete.Taille() != 0)
         defArete(1) = NULL;
      Destruction();
   };

// Lecture des donnees de la classe sur fichier
void 
Biel_axiQ::LectureDonneesParticulieres
    (UtilLecture * entreePrinc,Tableau<Noeud  *> * tabMaillageNoeud)
  { int nb;
    tab_noeud.Change_taille(nombre_V.nbne);
    for (int i=1; i<= nombre_V.nbne; i++)
     { *(entreePrinc->entree) >> nb;
       if ((entreePrinc->entree)->rdstate() == 0) 
       // pour mémoire ici on a 
           /*       enum io_state
                    {  badbit   = 1<<0, // -> 1 dans rdstate()
                       eofbit   = 1<<1, // -> 2 
                       failbit  = 1<<2, // -> 4 
                       goodbit  = 0     // -> O 
                     };*/
          tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
     #ifdef ENLINUX
       else  if ((entreePrinc->entree)->fail())
               // on a atteind la fin de la ligne et on appelle un nouvel enregistrement 
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
               *(entreePrinc->entree) >> nb;
              tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
              }
     #else
   /*   #ifdef SYSTEM_MAC_OS_X_unix
       else  if ((entreePrinc->entree)->fail())
               // on a atteind la fin de la ligne et on appelle un nouvel enregistrement 
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
               *(entreePrinc->entree) >> nb;
              tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
              }
      #else*/
       else  if ((entreePrinc->entree)->eof())
        // la lecture est bonne mais on a atteind la fin de la ligne
        {  tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture
           // si ce n'est pas la fin de la lecture on appelle un nouvel enregistrement 
           if (i != nombre_V.nbne) entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
         }
    //  #endif  
     #endif  
       else // cas d'une erreur de lecture
        { cout << "\n erreur de lecture inconnue  ";
          entreePrinc->MessageBuffer("** lecture des données particulières **");
          cout << "Biel_axiQ::LectureDonneesParticulieres";
          Affiche();
          Sortie (1);
        }
     }
    // construction du tableau de ddl des noeuds de biellette
    ConstTabDdl(); 
  };
  		
// calcul d'un point dans l'élément réel en fonction des coordonnées dans l'élément de référence associé
// temps: indique si l'on veut les coordonnées à t = 0, ou t ou tdt
// 1) cas où l'on utilise la place passée en argument
Coordonnee & Biel_axiQ::Point_physique(const Coordonnee& c_int,Coordonnee & co,Enum_dure temps)
 { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture
   // a) on commence par définir les bonnes grandeurs dans la métrique
   if( !(unefois.CalPt_0_t_tdt ))
      { unefois.CalPt_0_t_tdt += 1;
        Tableau<Enum_variable_metrique> tab(3);
        tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt;
        doCo->met_biellette.PlusInitVariables(tab) ;
       };  
   // b) calcul de l'interpolation 
   const Vecteur& phi = doCo->segment.Phi_point(c_int);
   // c) calcul du point
   switch (temps)
    { case TEMPS_0 : co = doCo->met_biellette.PointM_0(tab_noeud,phi); break;
      case TEMPS_t : co = doCo->met_biellette.PointM_t(tab_noeud,phi); break;
      case TEMPS_tdt : co = doCo->met_biellette.PointM_tdt(tab_noeud,phi); break;
     }       
   // d) retour
   return co;  
  };
  
// 3) cas où l'on veut les coordonnées aux 1, 2 ou trois temps selon la taille du tableau t_co
void Biel_axiQ::Point_physique(const Coordonnee& c_int,Tableau <Coordonnee> & t_co)
 { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;  // pour simplifier l'écriture
   // a) on commence par définir les bonnes grandeurs dans la métrique
   if( !(unefois.CalPt_0_t_tdt ))
      { unefois.CalPt_0_t_tdt += 1;
        Tableau<Enum_variable_metrique> tab(3);
        tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt;
        doCo->met_biellette.PlusInitVariables(tab) ;
       };  
   // b) calcul de l'interpolation 
   const Vecteur& phi = doCo->segment.Phi_point(c_int);
   // c) calcul des point
   switch (t_co.Taille())
    { case 3 : t_co(3) = doCo->met_biellette.PointM_tdt(tab_noeud,phi); 
      case 2 : t_co(2) = doCo->met_biellette.PointM_t(tab_noeud,phi); 
      case 1 : t_co(1) = doCo->met_biellette.PointM_0(tab_noeud,phi); 
     }       
  };
  
// Calcul du residu local à t ou tdt en fonction du booleen atdt
 Vecteur* Biel_axiQ::CalculResidu (bool atdt,const ParaAlgoControle & pa)
    { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
      Tableau <TenseurBB *>& d_epsBB = (doCo->d_epsBB);//             "
      // dimensionnement de la metrique
      if (!atdt)
       {if( !(unefois.CalResPrem_t ))
        { unefois.CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         doCo->met_biellette.PlusInitVariables(tab) ;
         };}  
      else 
       {if( !(unefois.CalResPrem_tdt ))
        {unefois.CalResPrem_tdt += 1;
         Tableau<Enum_variable_metrique> tab(11);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
         tab(10) = igradVBB_tdt;tab(11) = iVtdt; 
         doCo->met_biellette.PlusInitVariables(tab) ;
         };};  
       // initialisation du résidu 
       residu->Zero();  
       Vecteur poids =(doCo->segment).TaWi();
       ElemMeca::Cal_explicit ( doCo->tab_ddl,d_epsBB,nombre_V.nbi,poids,pa,atdt);
       // prise en compte de l'épaisseur sauf dans le cas d'une loi rien
       double epaisseur_moyenne = 0.; // où on met tout à 0
       // et surtout on ne recalcule pas la section car on n'a pas de compressibilité avec une loi rien !!
       if (!Loi_rien(loiComp->Id_comport()))
            epaisseur_moyenne = CalEpaisseurMoyenne_et_vol_pti(atdt);
       (*residu) *= epaisseur_moyenne;
       energie_totale *= epaisseur_moyenne;
       E_Hourglass *= epaisseur_moyenne; // meme si l'énergie d'hourglass est nulle
       E_elem_bulk_tdt*= epaisseur_moyenne; // idem
       P_elem_bulk *= epaisseur_moyenne;  // idem
       volume *=  epaisseur_moyenne;

       return residu;  
     };

// Calcul du residu local et de la raideur locale,
//  pour le schema implicite
Element::ResRaid  Biel_axiQ::Calcul_implicit (const ParaAlgoControle & pa)
    { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
      Tableau <TenseurBB *>& d_epsBB = (doCo->d_epsBB);//             "
      Tableau <TenseurHH *>& d_sigHH = (doCo->d_sigHH);//             "
      bool cald_Dvirtuelle = false;
      if (unefois.CalimpPrem == 0)
       { unefois.CalimpPrem = 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
         // on ne calcul la dérivée de la déformation virtuelle qu'une fois
         // car elle est constante dans le temps et indépendante des coordonnées
         cald_Dvirtuelle=true; 
        };  

       // initialisation du résidu
       residu->Zero();   
       // initialisation de la raideur      
       raideur->Zero(); 
             
       Vecteur poids =(doCo->segment).TaWi();
       ElemMeca::Cal_implicit(doCo->tab_ddl, d_epsBB,(doCo->d2_epsBB),d_sigHH,nombre_V.nbi
                              ,poids,pa,cald_Dvirtuelle);

       // prise en compte de l'épaisseur  sauf dans le cas d'une loi rien
       double epaisseur_moyenne = 0.; // où on met tout à 0
       // et surtout on ne recalcule pas la section car on n'a pas de compressibilité avec une loi rien !!
       if (!Loi_rien(loiComp->Id_comport()))
           {const bool atdt=true;
            epaisseur_moyenne = CalEpaisseurMoyenne_et_vol_pti(atdt);
           };

       (*residu) *= epaisseur_moyenne;
       (*raideur) *= epaisseur_moyenne;
       energie_totale *= epaisseur_moyenne;
       E_Hourglass *= epaisseur_moyenne; // meme si l'énergie d'hourglass est nulle
       E_elem_bulk_tdt*= epaisseur_moyenne; // idem
       P_elem_bulk *= epaisseur_moyenne;  // idem
       volume *=  epaisseur_moyenne;
       Element::ResRaid el;
       el.res = residu; 
       el.raid = raideur;       
      return el;
    };
		
// Calcul de la matrice masse pour l'élément
Mat_pleine * Biel_axiQ::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse)
 {    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb; // pour simplifier l'écriture
      // dimensionement de la métrique si nécessaire
      if (!(unefois.CalDynamique ))
       { unefois.CalDynamique += 1;
         Tableau<Enum_variable_metrique> tab(5);
         tab(1) = igiB_0;  tab(2) = igiB_tdt; tab(3) = igijBB_0;
         tab(4) = igijBB_tdt;tab(5) = igradVmoyBB_t; 
         doCo->met_biellette.PlusInitVariables(tab) ;
         // on vérifie le bon dimensionnement de la matrice
         if (type_calcul_masse == MASSE_CONSISTANTE)
           // dans le cas où la masse est consistante il faut la redimensionner
          { int nbddl = doCo->tab_ddl.NbDdl();
            (doCo->matrice_masse).Initialise (nbddl,nbddl,0.);
           } 
        };
       Vecteur poids =(doCo->segmentMas).TaWi();
       // appel de la routine générale
       ElemMeca::Cal_Mat_masse (doCo->tab_ddl,type_calcul_masse,
                            nombre_V.nbiMas,(doCo->segmentMas).TaPhi(),nombre_V.nbne
                            ,poids);
       (*mat_masse) *= donnee_specif.epais.epaisseur0; // prise en compte de l'épaisseur
       return mat_masse;
  };
		
//------- calcul d'erreur, remontée des contraintes -------------------
    // 1) calcul du résidu et de la matrice de raideur pour le calcul d'erreur
 Element::Er_ResRaid Biel_axiQ::ContrainteAuNoeud_ResRaid()
    { 
       if(!( unefois.CalResPrem_t ))
        {unefois.CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         doCo->met_biellette.PlusInitVariables(tab) ;
         };  
       // appel du programme général
       int tabn_taille = tab_noeud.Taille();
       ElemMeca::SigmaAuNoeud_ResRaid(tabn_taille
                  ,(doCo->segment).TaPhi()
                  ,(doCo->segment).TaWi()
                  ,doCo-> resErr,doCo->raidErr
                  ,(doCo->segmentEr).TaPhi()
                  ,(doCo->segmentEr).TaWi());
       // on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
       for (int i=1;i<= tabn_taille;i++)
        (*doCo-> resErr(i)) *= donnee_specif.epais.epaisseur_t;
       doCo->raidErr *= donnee_specif.epais.epaisseur_t;
       return  (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr)));
      };
      
    // 2) remontée aux erreurs aux noeuds
 Element::Er_ResRaid Biel_axiQ::ErreurAuNoeud_ResRaid()
    {  if(!( unefois.CalResPrem_t ))
       { unefois.CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         doCo->met_biellette.PlusInitVariables(tab) ;
         };  
       // appel du programme général
       int tabn_taille = tab_noeud.Taille();
       ElemMeca::Cal_ErrAuxNoeuds(tabn_taille, (doCo->segment).TaPhi(),
                (doCo->segment).TaWi(),doCo-> resErr );
       // on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
       for (int i=1;i<= tabn_taille;i++)
        (*doCo-> resErr(i)) *= donnee_specif.epais.epaisseur_t;
       doCo->raidErr *= donnee_specif.epais.epaisseur_t;
        return  (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr)));
      };
     
// actualisation des ddl et des grandeurs actives de t+dt vers t      
void Biel_axiQ::TdtversT()
 {   lesPtMecaInt.TdtversT(); // contrainte 
     // on boucle sur les sauvegardes
     for (int i=1; i<= nombre_V.nbi; i++)
       {if (tabSaveDon(i) != NULL) tabSaveDon(i)->TdtversT();
        if (tabSaveTP(i) != NULL) tabSaveTP(i)->TdtversT();
        if (tabSaveDefDon(i) != NULL) tabSaveDefDon(i)->TdtversT();
       };
	    donnee_specif.epais.epaisseur_t = donnee_specif.epais.epaisseur_tdt;
     ElemMeca::TdtversT_(); // appel de la procédure mère
 };
// actualisation des ddl et des grandeurs actives de t vers tdt      
void Biel_axiQ::TversTdt()
 {    lesPtMecaInt.TversTdt(); // contrainte
     // on boucle sur les sauvegardes
     for (int i=1; i<= nombre_V.nbi; i++)
       { if (tabSaveDon(i) != NULL) tabSaveDon(i)->TversTdt();
         if (tabSaveTP(i) != NULL) tabSaveTP(i)->TversTdt();
         if (tabSaveDefDon(i) != NULL)tabSaveDefDon(i)->TversTdt();
       };
	     donnee_specif.epais.epaisseur_tdt = donnee_specif.epais.epaisseur_t;
      ElemMeca::TversTdt_(); // appel de la procédure mère
 };

// calcul de l'erreur sur l'élément. Ce calcul n'est disponible
// qu'une fois la remontée aux contraintes effectuées sinon aucune
// action. En retour la valeur de l'erreur sur l'élément
// type indique le type de calcul d'erreur :
void Biel_axiQ::ErreurElement(int type,double& errElemRelative
                            ,double& numerateur, double& denominateur)
    {  if(!( unefois.CalResPrem_t ))
        { unefois.CalResPrem_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
          doCo->met_biellette.PlusInitVariables(tab) ;
        };  
       // appel du programme général
       ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur,
                            tab_noeud.Taille(),(doCo->segment).TaPhi(),
                (doCo->segment).TaWi(),
                (doCo->segmentEr).TaPhi(),(doCo->segmentEr).TaWi());
      };
// mise à jour de la boite d'encombrement de l'élément, suivant les axes I_a globales
// en retour coordonnées du point mini dans retour.Premier() et du point maxi dans .Second()
// la méthode est différente de la méthode générale car il faut prendre en compte l'épaisseur de l'élément
const DeuxCoordonnees& Biel_axiQ::Boite_encombre_element(Enum_dure temps)
  { // on commence par calculer la boite d'encombrement pour l'élément médian
    Element::Boite_encombre_element( temps);
    // ensuite on augmente sytématiquement dans toutes directions d'une valeur sqrt(s)/2 majorée
    double sSur2maj = sqrt(donnee_specif.epais.epaisseur_tdt) * 0.5
                         * ParaGlob::param->ParaAlgoControleActifs().Extra_boite_prelocalisation();
    // ajout d'un extra dans toutes les directions
    sSur2maj += ParaGlob::param->ParaAlgoControleActifs().Ajout_extra_boite_prelocalisation();
    // mise à jour
    boite_encombre.Premier().Ajout_meme_valeur(-sSur2maj);	// le min
	   boite_encombre.Second().Ajout_meme_valeur(sSur2maj);    // le max
    // retour
    return boite_encombre;
   };

//============= lecture écriture dans base info ==========
	
// cas donne le niveau de la récupération
// = 1 : on récupère tout
// = 2 : on récupère uniquement les données variables (supposées comme telles)
void Biel_axiQ::Lecture_base_info
    (ifstream& ent,const Tableau<Noeud  *> * tabMaillageNoeud,const int cas) 
{// tout d'abord appel de la lecture de la classe elem_meca
 ElemMeca::Lecture_bas_inf(ent,tabMaillageNoeud,cas);
 // traitement du cas particulier de la biellette
 switch (cas)
  { case 1 : // ------- on récupère tout -------------------------
     { // construction du tableau de ddl des noeuds
       ConstTabDdl(); 
       // récup contraintes et déformation
       lesPtMecaInt.Lecture_base_info(ent,cas);
       // les données spécifiques
		     string nom;
       ent >> nom;
       if (nom == "epaisStockeDansElement")
		         ent >>  nom >> donnee_specif.epais.epaisseur0
               >>  nom >> donnee_specif.epais.epaisseur_t
               >>  nom >> donnee_specif.epais.epaisseur_tdt;
	      break;
	    }
    case  2 : // ----------- lecture uniquement de se qui varie --------------------
     { // récup contraintes et déformation
       lesPtMecaInt.Lecture_base_info(ent,cas);
       // les données spécifiques
		     string nom;
       ent >> nom >> donnee_specif.epais.epaisseur_tdt;
       donnee_specif.epais.epaisseur_t = donnee_specif.epais.epaisseur_tdt;
	      break;
	      break;
	    }
   default :
    { cout << "\nErreur : valeur incorrecte du type de lecture !\n";
	     cout << "Biel_axiQ::Lecture_base_info(ofstream& sort,int cas)"
			        << " cas= " << cas << endl;
	     Sortie(1);
	   }
 };
};
  
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Biel_axiQ::Ecriture_base_info(ofstream& sort,const int cas)
{// tout d'abord appel de l'écriture de la classe elem_meca
 ElemMeca::Ecriture_bas_inf(sort,cas);
 // traitement du cas particulier de la biellette
 switch (cas)
  { case 1 : // ------- on sauvegarde tout -------------------------
     { 
       // des tenseurs déformation et contrainte,
       lesPtMecaInt.Ecriture_base_info(sort,cas);
       // les données spécifiques
       sort << "\n epaisStockeDansElement " << " epaisseur0= " << donnee_specif.epais.epaisseur0
                << " epaisseur_t= " << donnee_specif.epais.epaisseur_t
                << " epaisseur_tdt= " << donnee_specif.epais.epaisseur_tdt;
	      break;
	    }
    case  2 : // ----------- sauvegarde uniquement de se qui varie --------------------
     { // des tenseurs déformation et contrainte, 
       lesPtMecaInt.Ecriture_base_info(sort,cas);
		     sort << "\n epaisseur_tdt= " << donnee_specif.epais.epaisseur_tdt;
	      break;
	    }
    default :
	    { cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
	      cout << "Biel_axiQ::Ecriture_base_info(ofstream& sort,int cas)"
			         << " cas= " << cas << endl;
	      Sortie(1);
	    }
   };
};
 
      
// Calcul de la matrice géométrique et initiale
ElemMeca::MatGeomInit Biel_axiQ::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & pa)
    { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
      Tableau <TenseurBB *>& d_epsBB = (doCo->d_epsBB);//             "
      Tableau <TenseurHH *>& d_sigHH = (doCo->d_sigHH);//             "
      bool cald_Dvirtuelle = false;
      if (unefois.CalimpPrem == 0)
       { unefois.CalimpPrem = 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
         // on ne calcul la dérivée de la déformation virtuelle qu'une fois
         // car elle est constante dans le temps et indépendante des coordonnées
         cald_Dvirtuelle=true; 
        };
    // Par simplicité
    Mat_pleine & matGeom = doCo->matGeom;    
    Mat_pleine & matInit = doCo->matInit;    
    // mise à zéro de la matrice géométrique
    matGeom.Initialise();
    Vecteur poids =(doCo->segment).TaWi();
    // prise en compte de l'épaisseur  sauf dans le cas d'une loi rien
    double epaisseur_moyenne = 0.;
    // et surtout on ne recalcule pas la section car on n'a pas de compressibilité avec une loi rien !!
    if (!Loi_rien(loiComp->Id_comport()))
        {const bool atdt=true;
         epaisseur_moyenne = CalEpaisseurMoyenne_et_vol_pti(atdt);
        };
    poids *=  epaisseur_moyenne;
    ElemMeca::Cal_matGeom_Init 
          (matGeom,matInit,doCo->tab_ddl, d_epsBB,
          doCo->d2_epsBB,d_sigHH,nombre_V.nbi,poids,pa,cald_Dvirtuelle);
    return  MatGeomInit(&matGeom,&matInit);     
  } ;
	
// retourne les tableaux de ddl associés aux noeuds, gere par l'element
// ce tableau et specifique a l'element
const DdlElement & Biel_axiQ::TableauDdl() const
     { return doCo->tab_ddl; };
// liberation de la place pointee     
void Biel_axiQ::Libere ()
  {Element::Libere (); // liberation de residu et raideur
   LibereTenseur() ; // liberation des tenseur intermediaires
  };

// acquisition ou modification d'une loi de comportement
void Biel_axiQ::DefLoi (LoiAbstraiteGeneral * NouvelleLoi)
    { // verification du type de loi
      if ((NouvelleLoi->Dimension_loi() != 2) && (NouvelleLoi->Dimension_loi() != 4))
        { cout << "\n Erreur, la loi de comportement a utiliser avec des biel axi quadratique ";
          cout << " doit etre de type 2D (contrainte plane ou def plane), \n ici est de type = "
               << (NouvelleLoi->Dimension_loi()) << " !!! " << endl;
          Sortie(1);     
        }
      // cas d'une loi mécanique
      if (GroupeMecanique(NouvelleLoi->Id_categorie()))
       {loiComp = (Loi_comp_abstraite *) NouvelleLoi;
        // on boucle sur les points d'intégration
        for (int i=1; i<= nombre_V.nbi; i++)
         {tabSaveDon(i) = loiComp->New_et_Initialise();
          // idem pour le type de déformation mécanique associé
          tabSaveDefDon(i) = def->New_et_Initialise();
         };
        // définition du type de déformation associé à la loi
        loiComp->Def_type_deformation(*def);
        // on active les données particulières nécessaires au fonctionnement de la loi de comp
        loiComp->Activation_donnees(tab_noeud,dilatation,lesPtMecaInt); 
        };
      // cas d'une loi thermo physique
      if (GroupeThermique(NouvelleLoi->Id_categorie()))
       {loiTP = (CompThermoPhysiqueAbstraite *) NouvelleLoi;
        // on boucle sur les points d'intégration
        for (int i=1; i<= nombre_V.nbi; i++)
         { tabSaveTP(i) = loiTP->New_et_Initialise();};
        // on active les données particulières nécessaires au fonctionnement de la loi de comp
        loiTP->Activation_donnees(tab_noeud); 
        };
      // cas d'une loi de frottement 
      if (GroupeFrottement(NouvelleLoi->Id_categorie()))
        loiFrot = (CompFrotAbstraite *) NouvelleLoi;
    };
  
// test si l'element est complet
int Biel_axiQ::TestComplet()
  { int res = ElemMeca::TestComplet(); // test dans la fonction mere
    if ( donnee_specif.epais.epaisseur0 == epaisseur_defaut)
     { cout << "\n l'épaisseur de la biellette n'est pas defini  \n";
       res = 0; }
    if ( tab_noeud(1) == NULL)
     { cout << "\n les noeuds de la biellette ne sont pas defini  \n";
       res = 0; }
    else 
     { int testi =1;
       int posi = Id_nom_ddl("X1") -1;
       int dim = ParaGlob::Dimension();    
       for (int i =1; i<= dim; i++)
          for (int j=1;j<=nombre_V.nbne;j++)
           if(!(tab_noeud(j)->Existe_ici(Enum_ddl(posi+i))))
              testi = 0;
       if(testi == 0)
         { cout << "\n les ddls X1,X2 etc des noeuds de la biellette ne sont pas defini  \n";
           cout << " \n utilisez Biel_axiQ::ConstTabDdl() pour completer " ;
           res = 0; }            
      }          
    return res;   
  };
    
// ajout du tableau de ddl des noeuds de biellette
 void Biel_axiQ::ConstTabDdl()
    {  
    Tableau <Ddl> ta(ParaGlob::Dimension());
    int posi = Id_nom_ddl("X1") -1;    
    int dim = ParaGlob::Dimension();    
    for (int i =1; i<= dim-1; i++)
	    {Ddl inter((Enum_ddl(i+posi)),0.,LIBRE);
	     ta(i) = inter;
     };
    // le dernier ddl z est mis en HSLIBRE, car on ne le prend pas en compte dans le calcul
    // axisymétrique
    ta(3)=Ddl ((Enum_ddl(3+posi)),0.,HSLIBRE);

    // attribution des ddls aux noeuds
    for (int i=1; i<= nombre_V.nbne;i++)
        tab_noeud(i)->PlusTabDdl(ta);
    };
		
// procesure permettant de completer l'element apres
// sa creation avec les donnees du bloc transmis
// peut etre appeler plusieurs fois
 Element* Biel_axiQ::Complete(BlocGen & bloc,LesFonctions_nD*  lesFonctionsnD)
    { // complétion avec bloc
      if (bloc.Nom(1) == "epaisseurs")
       { donnee_specif.epais.epaisseur0 = bloc.Val(1);
         // on initialise aussi les grandeurs à t et tdt
         donnee_specif.epais.epaisseur_tdt = donnee_specif.epais.epaisseur_t = bloc.Val(1);
         return this;
        }
       else
         return ElemMeca::Complete_ElemMeca(bloc,lesFonctionsnD);
    };       

// Compléter pour la mise en place de la gestion de l'hourglass
Element* Biel_axiQ::Complet_Hourglass(LoiAbstraiteGeneral * loiHourglass, const BlocGen & bloc)
{ // on initialise le traitement de l'hourglass
  string str_precision; // string vide indique que l'on veut utiliser un élément normal
  Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;
  ElemMeca::Init_hourglass_comp(*(doCo->segmentHourg),str_precision,loiHourglass,bloc);
  // dans le cas où l'hourglass a été activé mais que l'élément n'a pas
  // de traitement particulier associé, alors on désactive l'hourglass
  if ( ((type_stabHourglass == STABHOURGLASS_PAR_COMPORTEMENT) || (type_stabHourglass == STABHOURGLASS_PAR_COMPORTEMENT_REDUIT))
     &&(doCo->segmentHourg == NULL))
    type_stabHourglass = STABHOURGLASS_NON_DEFINIE;
  return this;
};

// affichage dans la sortie transmise, des variables duales "nom"
// dans le cas ou nom est vide, affichage de "toute" les variables
 void Biel_axiQ::AfficheVarDual(ofstream& sort, Tableau<string>& nom)
     {// affichage de l'entête de l'element
      sort << "\n******************************************************************";
      sort << "\n Element bielette ("<<nombre_V.nbne<<" noeuds"<<nombre_V.nbi <<" point(s) d'integration) ";
      sort << "\n******************************************************************";
      // appel de la procedure de elem meca
      if (!(unefois.dualSortbiel) && (unefois.CalimpPrem))
        { VarDualSort(sort,nom,1,1);
          unefois.dualSortbiel += 1;
         } 
      else if ((unefois.dualSortbiel) && (unefois.CalimpPrem))       
         VarDualSort(sort,nom,1,11);
      else if (!(unefois.dualSortbiel) && (unefois.CalResPrem_tdt))       
        { VarDualSort(sort,nom,1,2);
          unefois.dualSortbiel += 1;
         }         
      else if ((unefois.dualSortbiel) && (unefois.CalResPrem_tdt))       
          VarDualSort(sort,nom,1,12); 
      // sinon on ne fait rien     
      };
		            
// récupération des  valeurs au numéro d'ordre  = iteg pour
// les grandeur enu
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
Tableau <double> Biel_axiQ::Valeur_a_diff_temps(bool absolue,Enum_dure enu_t,const List_io<Ddl_enum_etendu>& enu,int iteg)
    { // appel de la procedure de elem meca
	  int cas;       
      if (!(unefois.dualSortbiel) && (unefois.CalimpPrem))
        { cas=1;unefois.dualSortbiel = 1;
         } 
      else if ((unefois.dualSortbiel) && (unefois.CalimpPrem))       
         { cas = 11;}
      else if (!(unefois.dualSortbiel) && (unefois.CalResPrem_tdt))       
        { cas=2;unefois.dualSortbiel += 1;
         }         
      else if ((unefois.dualSortbiel) && (unefois.CalResPrem_tdt))       
         { cas = 12;}
      // sinon pour l'instant pb, car il faut définir des variable dans la métrique
      else 
        { cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
               << " pour inialiser les conteneurs des tenseurs  resultats ";
          if (ParaGlob::NiveauImpression() >= 4)     
            cout << "\n cas non prévu, unefois.dualSortbiel= " << unefois.dualSortbiel
                 << " unefois.CalimpPrem= " << unefois.CalimpPrem
                 << "\n Biel_axiQ::Valeur_a_diff_temps(Enum_dure enu_t...";
          Sortie(1);
         };       
      return ElemMeca::Valeur_multi(absolue,enu_t,enu,iteg,cas);
    };
    
// récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu
// ici il s'agit de grandeurs tensorielles, le retour s'effectue dans la liste
// de conteneurs quelconque associée
void Biel_axiQ::ValTensorielle_a_diff_temps(bool absolue,Enum_dure enu_t,List_io<TypeQuelconque>& enu,int iteg)
    { // appel de la procedure de elem meca
	  int cas;       
      if (!(unefois.dualSortbiel ) && (unefois.CalimpPrem ))
        { cas=1;unefois.dualSortbiel += 1;
         } 
      else if ((unefois.dualSortbiel ) && (unefois.CalimpPrem ))       
         { cas = 11;}
      else if (!(unefois.dualSortbiel ) && (unefois.CalResPrem_tdt ))       
        { cas=2;unefois.dualSortbiel += 1;
         }         
      else if ((unefois.dualSortbiel ) && (unefois.CalResPrem_tdt ))       
         { cas = 12;}
      // sinon pour l'instant pb, car il faut définir des variable dans la métrique
      else 
        { cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
               << " pour inialiser les conteneurs des tenseurs  resultats ";
          if (ParaGlob::NiveauImpression() >= 4)     
            cout << "\n cas non prévu, unefois.dualSortbiel= " << unefois.dualSortbiel
               << " unefois.CalimpPrem= " << unefois.CalimpPrem
               << "\n Biel_axiQ::ValTensorielle_a_diff_temps(Enum_dure enu_t...";
          Sortie(1);
         };       
      ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas);
    };
			
// cas d'un chargement volumique, 
// force indique la force volumique appliquée
// retourne  le second membre résultant
// ici on considère l'épaisseur de la biellette pour constituer le volume
// -> explicite à t ou tdt en fonction de la variable booléenne atdt
Vecteur Biel_axiQ::SM_charge_volumique_E
      (const Coordonnee& force,Fonction_nD* pt_fonct,bool atdt,const ParaAlgoControle & pa,bool sur_volume_finale_)
    { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
      // dimensionnement de la metrique
      if (!atdt)
      {if(!(unefois.CalSMvol_t ))
       { unefois.CalSMvol_t += 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };} 
      else
      {if(!(unefois.CalSMvol_tdt ))
       { unefois.CalSMvol_tdt += 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };}; 
      // initialisation du résidu
      residu->Zero();   
    // appel du programme général d'elemmeca et retour du vecteur second membre
    // multiplié par l'épaisseur pour avoir le volume, on considère que l'épaisseur est à jour
    double epaisseur = donnee_specif.epais.epaisseur_tdt;
    if (!atdt) epaisseur = donnee_specif.epais.epaisseur_t;
      return (epaisseur * ElemMeca::SM_charge_vol_E (doCo->tab_ddl,(doCo->segment).TaPhi()
               ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa,sur_volume_finale_));
   };

// calcul des seconds membres suivant les chargements 
// cas d'un chargement volumique, 
// force indique la force volumique appliquée
// retourne  le second membre résultant
// ici on l'épaisseur de l'élément pour constituer le volume -> implicite,
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur  
// retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_volumique_I
                 (const Coordonnee& force,Fonction_nD* pt_fonct,const ParaAlgoControle & pa,bool sur_volume_finale_)
    { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture

      // initialisation du résidu
      residu->Zero();   
      // initialisation de la raideur      
      raideur->Zero();       

      // -- définition des constantes de la métrique si nécessaire
      // en fait on fait appel aux même éléments que pour le calcul implicite
      if (!(unefois.CalSMvol_tdt ))
       { unefois.CalSMvol_tdt += 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };

      // appel du programme général d'elemmeca 
      ElemMeca::SMR_charge_vol_I (doCo->tab_ddl,(doCo->segment).TaPhi()
               ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa,sur_volume_finale_);
      // prise en compte de l'épaisseur
      (*residu) *= donnee_specif.epais.epaisseur_t;
      (*raideur) *= donnee_specif.epais.epaisseur_t;
      Element::ResRaid el;
      el.res = residu; 
      el.raid = raideur;       
      return el;
   };

// calcul des seconds membres suivant les chargements 
// cas d'un chargement surfacique, sur les frontières des éléments
// force indique la force surfacique appliquée
// numface indique le numéro de la face chargée
// retourne  le second membre résultant
// -> explicite à t ou tdt en fonction de la variable booléenne atdt
Vecteur Biel_axiQ::SM_charge_surfacique_E
                        (const Coordonnee& force,Fonction_nD* pt_fonct,int ,bool atdt,const ParaAlgoControle & pa) 
  {  // initialisation du vecteur résidu
     ((*res_extS)(1))->Zero();
    // on récupère ou on crée la frontière surfacique
    Frontiere_surfacique(1,true);
    // on pourrait utiliser la métrique des éléments de frontière
    // avec l'instance déformation dédiée pour
    // mais on utilise celle de l'élément   
    
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if (!atdt)
    {if ( unefois.CalSMsurf_t == 0)
       { unefois.CalSMsurf_t = 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };}  
    else
     {if ( unefois.CalSMsurf_tdt == 0)
       { unefois.CalSMsurf_tdt = 1;
         Tableau<Enum_variable_metrique> tab(11);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
         tab(10) = igradVBB_tdt;tab(11) = iVtdt; 
         doCo->met_biellette.PlusInitVariables(tab) ;
        };};  
    // on définit la déformation a doc si elle n'existe pas déjà
    if (defSurf(1) == NULL) 
      defSurf(1) = new Deformation
            (*met,tabb(1)->TabNoeud(),(doCo->segS).TaDphi(),(doCo->segS).TaPhi());
    // appel du programme général d'elemmeca et retour du vecteur second membre
    // 1 pour dire que c'est la première surface externe
    return ElemMeca::SM_charge_surf_E (tabb(1)->DdlElem(),1
               ,(doCo->segS).TaPhi(),tab_noeud.Taille()
               ,(doCo->segS).TaWi(),force,pt_fonct,pa);
   }; 
	   
// cas d'un chargement surfacique, sur les frontières des éléments
// force indique la force surfacique appliquée
// numface indique le numéro de la face chargée
// retourne  le second membre résultant
// -> implicite, 
// pa : permet de déterminer si oui ou non on  calcul la contribution à la raideur 
// retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_surfacique_I
        (const Coordonnee& force,Fonction_nD* pt_fonct,int ,const ParaAlgoControle & pa) 
  {  // initialisation du vecteur résidu et de la raideur
     // normalement numface = 1
     ((*res_extS)(1))->Zero();
     ((*raid_extS)(1))->Zero();
    // on récupère ou on crée la frontière surfacique
    Frontiere_surfacique(1,true);
    // on pourrait utiliser la métrique des éléments de frontière
    // avec l'instance déformation dédiée pour
    // mais on utilise celle de l'élément   
    
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if ( unefois.CalSMRsurf == 0)
       { unefois.CalSMRsurf = 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         met->PlusInitVariables(tab) ;
        }; 
    // on définit la déformation a doc si elle n'existe pas déjà
    // on utilise la même déformation pour toutes les arrêtes.
    if (defSurf(1) == NULL)
      defSurf(1) = new Deformation(*met,tabb(1)->TabNoeud(),
            (doCo->segS).TaDphi(),(doCo->segS).TaPhi());
   
    // appel du programme général d'elemmeca et retour du vecteur second membre
    return ElemMeca::SMR_charge_surf_I (tabb(1)->DdlElem(),1
               ,(doCo->segS).TaPhi(),(tabb(1)->TabNoeud()).Taille()
               ,(doCo->segS).TaWi(),force,pt_fonct,pa);
   };

// cas d'un chargement de type pression, sur les frontières des éléments
// pression indique la pression appliquée
// numface indique le numéro de la face chargée
// retourne  le second membre résultant
// NB: il y a une définition par défaut pour les éléments qui n'ont pas de 
// surface externe -> message d'erreur  d'où le virtuel et non virtuel pur
// -> explicite à t ou tdt en fonction de la variable booléenne atdt
Vecteur Biel_axiQ::SM_charge_pression_E(double pression,Fonction_nD* pt_fonct,int numArete,bool atdt,const ParaAlgoControle & pa)
  {  // initialisation du vecteur résidu
     ((*res_extA)(1))->Zero();
     //dans le cas des éléments axisymétriques on utilise les déformations d'arêtes pour les pressions
     // le surface, n'est pas utilisable pour la pression car en rotation cela donne un volume !!
     ElFrontiere* elf =Frontiere_lineique(numArete,true);
     // on récupère ou on crée la frontière , si elle n'existe pas il faut la creer d'ou le true

     // on  utilise la métrique des éléments de frontière
     // avec l'instance déformation dédiée pour
     // récupération de la métrique associée à l'élément qui est commune a tous les triangles
     // du même type
     Met_abstraite * meta= elf->Metrique();
    
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // on pourrait utiliser la métrique des éléments de frontière
    // avec l'instance déformation dédiée pour
    // mais on utilise celle de l'élément
    // dimensionnement de la metrique
    if (!atdt)
     {if ( unefois.CalSMsurf_t == 0)
       { unefois.CalSMsurf_t = 1;
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         meta->PlusInitVariables(tab) ;
        };}  
    else 
     {if ( unefois.CalSMsurf_tdt == 0)
       { unefois.CalSMsurf_tdt = 1;
         Tableau<Enum_variable_metrique> tab(11);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
         tab(10) = igradVBB_tdt;tab(11) = iVtdt; 
         meta->PlusInitVariables(tab) ;
        };};  
    // on définit la déformation ad hoc si elle n'existe pas déjà
    if (defArete(numArete) == NULL)
      defArete(numArete) = new Deformation(*meta,elf->TabNoeud(),
            (doCo->segS).TaDphi(),(doCo->segS).TaPhi());
    // appel du programme général d'elemmeca et retour du vecteur second membre
    // 1 pour dire que c'est la première surface externe
    return ElemMeca::SM_charge_pres_E (elf->DdlElem(),numArete
               ,(doCo->segS).TaPhi(),(elf->TabNoeud()).Taille()
               ,(doCo->segS).TaWi(),pression,pt_fonct,pa);

   }; 
	   
// cas d'un chargement de type pression, sur les frontières des éléments
// pression indique la pression appliquée
// numface indique le numéro de la face chargée
// retourne  le second membre résultant
// NB: il y a une définition par défaut pour les éléments qui n'ont pas de 
// surface externe -> message d'erreur  d'où le virtuel et non virtuel pur -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
// retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_pression_I
       (double pression,Fonction_nD* pt_fonct,int numArete,const ParaAlgoControle & pa)
  {  // initialisation du vecteur résidu et de la raideur ici d'arête vu l'axisymétrie
    ((*res_extA)(numArete))->Zero();
    ((*raid_extA)(numArete))->Zero();
     //dans le cas des éléments axisymétriques on utilise les déformations d'arêtes pour les pressions
     // le surface, n'est pas utilisable pour la pression car en rotation cela donne un volume !!
     ElFrontiere* elf =Frontiere_lineique(numArete,true);
     // on  utilise la métrique des éléments de frontière
     // avec l'instance déformation dédiée pour
     // récupération de la métrique associée à l'élément qui est commune a tous les triangles
     // du même type
     Met_abstraite * meta= elf->Metrique();
    
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if( unefois.CalSMRlin == 0)
      {  unefois.CalSMRlin = 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         meta->PlusInitVariables(tab) ;
        }; 
    // on définit la déformation ad hoc si elle n'existe pas déjà
    if (defArete(numArete) == NULL)
      defArete(numArete) = new Deformation(*meta,elf->TabNoeud(),
            (doCo->segS).TaDphi(),(doCo->segS).TaPhi());
           
    // appel du programme général d'elemmeca et retour du vecteur second membre
    return ElemMeca::SMR_charge_pres_I (elf->DdlElem(),numArete
               ,(doCo->segS).TaPhi(),(elf->TabNoeud()).Taille()
               ,(doCo->segS).TaWi(),pression,pt_fonct,pa);

   };


// cas d'un chargement de type pression hydrostatique, sur les frontières des éléments
// la charge dépend de la hauteur à la surface libre du liquide déterminée par un point
// et une  direction  normale à la surface libre:
// nSurf : le numéro de la surface externe
// poidvol: indique le poids volumique du liquide
// M_liquide : un point de la surface libre
// dir_normal_liquide : direction normale à la surface libre
// retourne  le second membre résultant
// calcul à l'instant tdt ou t en fonction de la variable atdt
// retourne  le second membre résultant
// NB: il y a une définition par défaut pour les éléments qui n'ont pas de
// surface externe -> message d'erreur  d'où le virtuel et non virtuel pur
// -> explicite à t ou à tdt en fonction de la variable booléenne atdt
Vecteur Biel_axiQ::SM_charge_hydrostatique_E(const Coordonnee& dir_normal_liquide,const double& poidvol
                              ,int ,const Coordonnee& M_liquide,bool atdt
                              ,const ParaAlgoControle & pa
                              ,bool sans_limitation)
  {  int numArete = 1; // ici c'est la seule arête qui nous intéresse
     // initialisation du vecteur résidu
     ((*res_extA)(numArete))->Zero();
     //dans le cas des éléments axisymétriques on utilise les déformations d'arêtes pour les pressions
     // en rotation cela donne une surface
     ElFrontiere* elf =Frontiere_lineique(numArete,true);

     //        // récupération de la métrique associée à l'élément
     Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
      // on pourrait utiliser la métrique des éléments de frontière
      // avec l'instance déformation dédiée pour
      // mais on utilise celle de l'élément

      // dimensionnement de la metrique
      if (!atdt)
       {if ( unefois.CalSMsurf_t == 0)
         { unefois.CalSMsurf_t = 1;
           Tableau<Enum_variable_metrique> tab(10);
           tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
           tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
           tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
           doCo->met_biellette.PlusInitVariables(tab) ;
          };}
      else
       {if ( unefois.CalSMsurf_tdt == 0)
         { unefois.CalSMsurf_tdt = 1;
           Tableau<Enum_variable_metrique> tab(11);
           tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
           tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
           tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
           tab(10) = igradVBB_tdt;tab(11) = iVtdt;
           doCo->met_biellette.PlusInitVariables(tab) ;
          };
       };
      // on définit la déformation ad hoc si elle n'existe pas déjà
      if (defArete(numArete) == NULL)
        defArete(numArete) = new Deformation(doCo->met_biellette,elf->TabNoeud(),
              (doCo->segS).TaDphi(),(doCo->segS).TaPhi());
      // appel du programme général d'elemmeca et retour du vecteur second membre
      return ElemMeca::SM_charge_hydro_E (tabb(1)->DdlElem(),numArete
               ,(doCo->segS).TaPhi(),(elf->TabNoeud()).Taille()
               ,(doCo->segS).TaWi(),dir_normal_liquide,poidvol,M_liquide,sans_limitation,pa,atdt);
   };

// cas d'un chargement de type pression hydrostatique, sur les frontières des éléments
// la charge dépend de la hauteur à la surface libre du liquide déterminée par un point
// et une  direction  normale à la surface libre:
// nSurf : le numéro de la surface externe
// poidvol: indique le poids volumique du liquide
// M_liquide : un point de la surface libre
// dir_normal_liquide : direction normale à la surface libre
// retourne  le second membre résultant
// calcul à l'instant tdt ou t en fonction de la variable atdt
// retourne  le second membre résultant
// NB: il y a une définition par défaut pour les éléments qui n'ont pas de
// surface externe -> message d'erreur  d'où le virtuel et non virtuel pur
// -> implicite,
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
// retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_hydrostatique_I
                        (const Coordonnee& dir_normal_liquide,const double& poidvol
                         ,int ,const Coordonnee& M_liquide
                         ,const ParaAlgoControle & pa
                         ,bool sans_limitation
                        )
  {  int numArete = 1; // ici c'est la seule arête qui nous intéresse
     // initialisation du vecteur résidu et de la raideur
     ((*res_extA)(numArete))->Zero();
     ((*raid_extA)(numArete))->Zero();
     //dans le cas des éléments axisymétriques on utilise les déformations d'arêtes pour les pressions
     // le surface, n'est pas utilisable pour la pression car en rotation cela donne un volume !!
     ElFrontiere* elf =Frontiere_lineique(numArete,true);
     // on  utilise la métrique des éléments de frontière
     // avec l'instance déformation dédiée pour
     // récupération de la métrique associée à l'élément qui est commune a tous les triangles
     // du même type
     Met_abstraite * meta= elf->Metrique();
    
     Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
     // dimensionnement de la metrique
     if( unefois.CalSMRlin == 0)
      {  unefois.CalSMRlin = 1;
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt;
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         meta->PlusInitVariables(tab) ;
        };
     // on définit la déformation ad hoc si elle n'existe pas déjà
     if (defArete(numArete) == NULL)
      defArete(numArete) = new Deformation(*meta,elf->TabNoeud(),
            (doCo->segS).TaDphi(),(doCo->segS).TaPhi());

     // appel du programme général d'elemmeca et retour du vecteur second membre
     return ElemMeca::SMR_charge_hydro_I (tabb(1)->DdlElem(),1
               ,(doCo->segS).TaPhi(),(elf->TabNoeud()).Taille()
               ,(doCo->segS).TaWi(),dir_normal_liquide,poidvol,M_liquide,sans_limitation,pa);
   };

// cas d'un chargement lineique, sur l'aretes de la biellette
// force indique la force lineique appliquée
// numarete indique le numéro de l'arete chargée ici 1
// retourne  le second membre résultant
// -> explicite à t ou tdt en fonction de la variable booléenne atdt
Vecteur Biel_axiQ::SM_charge_lineique_E(const Coordonnee& force,Fonction_nD* pt_fonct,int ,bool atdt,const ParaAlgoControle & pa)
  { // initialisation du vecteur résidu 
    ((*res_extA)(1))->Zero();
    // on récupère ou on crée la frontière arrête
    Frontiere_lineique(1,true);
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if (!atdt)
    {if( !(unefois.CalSMlin_t ))
      {  unefois.CalSMlin_t += 1;
         Tableau<Enum_variable_metrique> tab(8);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = igradVBB_t;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };} 
    else
    {if( !(unefois.CalSMlin_tdt ))
      {  unefois.CalSMlin_tdt += 1;
         Tableau<Enum_variable_metrique> tab(8);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = igradVBB_tdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };}; 
    // on définit la déformation a doc si elle n'existe pas déjà
    if (defArete(1) == NULL) 
      defArete(1) = def; // a priori idem que la biellette

    // appel du programme général d'elemmeca et retour du vecteur second membre
    return ElemMeca::SM_charge_line_E (doCo->tab_ddl,1,(doCo->segment).TaPhi()
               ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa);
   };
    
// cas d'un chargement lineique, sur les aretes frontières des éléments
// force indique la force lineique appliquée
// numarete indique le numéro de l'arete chargée ici 1 par défaut
// retourne  le second membre résultant
// -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
// retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_lineique_I
                (const Coordonnee& force,Fonction_nD* pt_fonct,int ,const ParaAlgoControle & pa)
  { // initialisation du vecteur résidu et de la raideur
    ((*res_extA)(1))->Zero();
    ((*raid_extA)(1))->Zero();
    // on récupère ou on crée la frontière arrête
    Frontiere_lineique(1,true);
    
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if( !(unefois.CalSMRlin ))
      {  unefois.CalSMRlin += 1;
         Tableau<Enum_variable_metrique> tab(15);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = igradVBB_tdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        }; 
    // on définit la déformation a doc si elle n'existe pas déjà
    if (defArete(1) == NULL) 
      defArete(1) = def; // a priori idem que la biellette
           
    // appel du programme général d'elemmeca et retour du vecteur second membre
    return ElemMeca::SMR_charge_line_I (doCo->tab_ddl,1
               ,(doCo->segment).TaPhi(),tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa);
   };


// cas d'un chargement lineique suiveuse, sur l'aretes frontière de la biellette (2D uniquement)
// force indique la force lineique appliquée
// numarete indique le numéro de l'arete chargée
// retourne  le second membre résultant
// -> explicite à t ou tdt en fonction de la variable booléenne atdt
Vecteur Biel_axiQ::SM_charge_lineique_Suiv_E(const Coordonnee& force,Fonction_nD* pt_fonct,int ,bool atdt,const ParaAlgoControle & pa)
  { // initialisation du vecteur résidu 
    ((*res_extA)(1))->Zero();
    // on récupère ou on crée la frontière arrête
    Frontiere_lineique(1,true);
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if (!atdt)
    {if( !(unefois.CalSMlin_t ))
      {  unefois.CalSMlin_t += 1;
         Tableau<Enum_variable_metrique> tab(8);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = igradVBB_t;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };} 
    else
    {if( !(unefois.CalSMlin_tdt ))
      {  unefois.CalSMlin_tdt += 1;
         Tableau<Enum_variable_metrique> tab(8);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = igradVBB_tdt;
         doCo->met_biellette.PlusInitVariables(tab) ;
        };}; 
    // on définit la déformation a doc si elle n'existe pas déjà
    if (defArete(1) == NULL) 
      defArete(1) = def; // a priori idem que la biellette

    // appel du programme général d'elemmeca et retour du vecteur second membre
    return ElemMeca::SM_charge_line_Suiv_E (doCo->tab_ddl,1,(doCo->segment).TaPhi()
               ,tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa);
   };
    
// cas d'un chargement lineique suiveuse, sur l'aretes frontière de la biellette (2D uniquement)
// force indique la force lineique appliquée
// numarete indique le numéro de l'arete chargée
// retourne  le second membre résultant
// -> implicite, 
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
// retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_lineique_Suiv_I
                (const Coordonnee& force,Fonction_nD* pt_fonct,int ,const ParaAlgoControle & pa)
  { // initialisation du vecteur résidu et de la raideur
    ((*res_extA)(1))->Zero();
    ((*raid_extA)(1))->Zero();
    // on récupère ou on crée la frontière arrête
    Frontiere_lineique(1,true);
    
    Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
    // dimensionnement de la metrique
    if( !(unefois.CalSMRlin ))
      {  unefois.CalSMRlin += 1;
         Tableau<Enum_variable_metrique> tab(15);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt;
         tab(15) = igradVBB_tdt; 
         doCo->met_biellette.PlusInitVariables(tab) ;
        }; 
    // on définit la déformation a doc si elle n'existe pas déjà
    if (defArete(1) == NULL) 
      defArete(1) = def; // a priori idem que la biellette
           
    // appel du programme général d'elemmeca et retour du vecteur second membre
    return ElemMeca::SMR_charge_line_Suiv_I (doCo->tab_ddl,1
               ,(doCo->segment).TaPhi(),tab_noeud.Taille(),(doCo->segment).TaWi(),force,pt_fonct,pa);
   };

	   // cas d'un chargement surfacique hydro-dynamique, 
       // Il y a trois forces: une suivant la direction de la vitesse: de type traînée aerodynamique
       // Fn = poids_volu * fn(V) * S * (normale*u) * u, u étant le vecteur directeur de V (donc unitaire)
       // une suivant la direction normale à la vitesse de type portance
       // Ft = poids_volu * ft(V) * S * (normale*u) * w, w unitaire, normal à V, et dans le plan n et V
       // une suivant la vitesse tangente de type frottement visqueux
       // T = to(Vt) * S * ut, Vt étant la vitesse tangentielle et ut étant le vecteur directeur de Vt
       // coef_mul: est un coefficient multiplicateur global (de tout)
	   // retourne  le second membre résultant
	   // // -> explicite à t ou à tdt en fonction de la variable booléenne atdt	   
Vecteur Biel_axiQ::SM_charge_hydrodynamique_E( Courbe1D* frot_fluid,const double& poidvol
	                                                 , Courbe1D* coef_aero_n,int num,const double& coef_mul
	                                                 , Courbe1D* coef_aero_t,bool atdt,const ParaAlgoControle & pa) 
  {  Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
     int dime = ParaGlob::Dimension();
     Met_abstraite * meta; ElemGeomC0* elemGeom; // définit dans le choix suivant
     // on  utilise la métrique des éléments de frontière
     // avec l'instance déformation dédiée pour
     // récupération de la métrique associée à l'élément qui est commune a tous les bel_axi
     // du même type
     if (dime == 3) // il s'agit de surface
       // initialisation du vecteur résidu  pour la surface
      {((*res_extA)(num))->Zero();
       ElFrontiere* elf = Frontiere_lineique(num,true);// on récupère ou on crée la frontière linéique
       meta= elf->Metrique();
       elemGeom = &(doCo->segS); // récup de l'élément géométrique
      }
     else // dime 2: il s'agit d'un point
       // initialisation du vecteur résidu  pour un point
       // initialisation du vecteur résidu  pour une arrête
      {((*res_extN)(num))->Zero();
       ElFrontiere* elf = Frontiere_points(num,true);// on récupère ou on crée la frontière point
       meta= elf->Metrique();
       elemGeom = &(doCo->point); // récup de l'élément géométrique
      };

    // dimensionnement de la metrique axi symétrique
    // définition des constantes de la métrique si nécessaire
    if (!atdt)
     {if (((dime == 3)&&(unefois.CalSMsurf_t == 0)) || ((dime == 2)&&(unefois.CalSMlin_t == 0)))
       { if (dime == 3) {unefois.CalSMsurf_t = 1;} else {unefois.CalSMlin_t=1;};
         Tableau<Enum_variable_metrique> tab(10);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
         tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
         tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
         meta->PlusInitVariables(tab) ;
        };}  
    else 
     {if (((dime == 3)&&(unefois.CalSMsurf_tdt == 0)) || ((dime == 2)&&(unefois.CalSMlin_tdt == 0)))
       { if (dime == 3) {unefois.CalSMsurf_tdt = 1;} else {unefois.CalSMlin_tdt=1;};
         Tableau<Enum_variable_metrique> tab(11);
         tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
         tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
         tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
         tab(10) = igradVBB_tdt;tab(11) = iVtdt; 
         meta->PlusInitVariables(tab) ;
        };};  
    // on définit la déformation a doc si elle n'existe pas déjà
    if (dime == 3)
     { if (defSurf(num) == NULL)
        defSurf(num) = new Deformation(*meta,tabb(num)->TabNoeud(),
            (doCo->segS).TaDphi(),(doCo->segS).TaPhi());
     }
   else 
     { if (defArete(num) == NULL)
        defArete(num) = new Deformation(*meta,tabb(num)->TabNoeud(),
            (doCo->point).TaDphi(),(doCo->point).TaPhi());
     };
    // appel du programme général d'ElemMeca et retour du vecteur second membre
    // tous les infos ne sont pas utilisées car pour les frontières points, pas de points d'integ
    return  ElemMeca::SM_charge_hydrodyn_E (poidvol,elemGeom->TaPhi(),(tabb(num)->TabNoeud()).Taille()
                                         ,frot_fluid,elemGeom->TaWi()
                                         ,coef_aero_n,num,coef_mul,coef_aero_t,pa,atdt);
   
   };

	   // -> implicite, 
	   // pa: permet de déterminer si oui ou non on calcul la contribution à la raideur 
	   // retourne le second membre et la matrice de raideur correspondant
Element::ResRaid Biel_axiQ::SMR_charge_hydrodynamique_I(  Courbe1D* frot_fluid,const double& poidvol
	                                                ,  Courbe1D* coef_aero_n,int num,const double& coef_mul
	                                                ,  Courbe1D* coef_aero_t,const ParaAlgoControle & pa) 
  {  Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
     int dime = ParaGlob::Dimension();
     Met_abstraite * meta; ElemGeomC0* elemGeom; // définit dans le choix suivant
     // on  utilise la métrique des éléments de frontière
     // avec l'instance déformation dédiée pour
     // récupération de la métrique associée à l'élément qui est commune a tous les bel_axi
     // du même type
     if (dime == 3) // il s'agit de surface
       // initialisation du vecteur résidu  pour la surface
      {((*res_extA)(num))->Zero();((*raid_extA)(num))->Zero();
       ElFrontiere* elf = Frontiere_lineique(num,true);// on récupère ou on crée la frontière linéique
       meta= elf->Metrique();
       elemGeom = &(doCo->segS); // récup de l'élément géométrique
      }
     else // dime 2: il s'agit d'un point
       // initialisation du vecteur résidu  pour un point
       // initialisation du vecteur résidu  pour une arrête
      {((*res_extN)(num))->Zero();((*raid_extN)(num))->Zero();
       ElFrontiere* elf = Frontiere_points(num,true);// on récupère ou on crée la frontière point
       meta= elf->Metrique();
       elemGeom = &(doCo->point); // récup de l'élément géométrique
      };

    // dimensionnement de la metrique
    // définition des constantes de la métrique si nécessaire
    if (((dime == 3)&&(unefois.CalSMRsurf == 0)) || ((dime == 2)&&(unefois.CalSMRlin == 0)))
       { if (dime == 3) {unefois.CalSMRsurf = 1;} else {unefois.CalSMRlin=1;};
         Tableau<Enum_variable_metrique> tab(20);
         tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
         tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
         tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) =  id_giH_tdt;
         tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt; 
         tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
         tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
         meta->PlusInitVariables(tab) ;
        }; 

    // on définit la déformation a doc si elle n'existe pas déjà
    if (dime == 3)
     { if (defSurf(num) == NULL)
        defSurf(num) = new Deformation(*meta,tabb(num)->TabNoeud(),
            (doCo->segS).TaDphi(),(doCo->segS).TaPhi());
     }
   else 
     { if (defArete(num) == NULL)
        defArete(num) = new Deformation(*meta,tabb(num)->TabNoeud(),
            (doCo->point).TaDphi(),(doCo->point).TaPhi());
     };
    // appel du programme général d'ElemMeca et retour du vecteur second membre
    // tous les infos ne sont pas utilisées car pour les frontières points, pas de points d'integ
    Element::ResRaid el(ElemMeca::SM_charge_hydrodyn_I (poidvol,doCo->segment.TaPhi(),1
                                         ,frot_fluid,doCo->segment.TaWi(),tabb(posi_tab_front_point+num)->DdlElem()
                                         ,coef_aero_n,num,coef_mul,coef_aero_t,pa));
    return el;
   };

	
// calcul de la nouvelle épaisseur moyenne finale (sans raideur) 
// mise à jour des volumes aux pti
// ramène l'épaisseur moyenne calculée à atdt
const double& Biel_axiQ::CalEpaisseurMoyenne_et_vol_pti(bool atdt)
  {  //     ..  en bouclant sur les pt d'integ enregistré ..
	  // -- on récupère et on calcule les jacobiens moyens à 0 et final	  
	  double jacobien_moy_fin = 0.; // init
   double jacobien_moy_ini = 0.;
	  // -- de même on récupère et on calcul la trace moyenne de la contrainte
	  double traceSig_moy = 0.;double traceSig_moy_ini = 0.;
	  // -- de même on récupère et on calcul le module de compressibilité moyen
	  double troisK_moy = 0.;
	  Epai&  epais = donnee_specif.epais; // pour simplifier
	  for (int i=1;i<= nombre_V.nbi;i++)
	    { // cas de la compressibilité
       const double& troisK = 3. * (*lesPtIntegMecaInterne)(i).ModuleCompressibilite_const();
			    troisK_moy += troisK;
       // cas des jacobiens
//       const double& jacobien_0 = *(tabSaveDefDon(i)->Meti_00().jacobien_);
			    if (atdt)
	        { const double& jacobien_ini = *(tabSaveDefDon(i)->Meti_t().jacobien_);
           jacobien_moy_ini += jacobien_ini;
           const double& jacobien_fin = *(tabSaveDefDon(i)->Meti_tdt().jacobien_);
           jacobien_moy_fin += jacobien_fin;
			        // cas de la trace de sigma
			        const double traceSig = (*lesPtIntegMecaInterne)(i).SigHH_const() && (*(tabSaveDefDon(i)->Meti_tdt().gijBB_));
			        traceSig_moy += traceSig;
           const double traceSig_ini = (*lesPtIntegMecaInterne)(i).SigHH_t_const() && (*(tabSaveDefDon(i)->Meti_t().gijBB_));
           traceSig_moy_ini += traceSig_ini;
          
           (*lesPtIntegMecaInterne)(i).Volume_pti() *= (epais.epaisseur_t * jacobien_ini) / jacobien_fin
		                           * troisK / (troisK - traceSig+traceSig_ini); // la pression = - traceSig/3 !!!
			      }
			    else
	        { const double& jacobien_ini = *(tabSaveDefDon(i)->Meti_00().jacobien_);
           jacobien_moy_ini += jacobien_ini;
           const double& jacobien_fin = *(tabSaveDefDon(i)->Meti_t().jacobien_);
           jacobien_moy_fin += jacobien_fin;
			        // cas de la trace de sigma
			        double traceSig = (*lesPtIntegMecaInterne)(i).SigHH_const() && (*(tabSaveDefDon(i)->Meti_t().gijBB_));
           traceSig_moy += traceSig;
           (*lesPtIntegMecaInterne)(i).Volume_pti() *= (epais.epaisseur0 * jacobien_ini) / jacobien_fin
                               * troisK / (troisK - traceSig);
			      };
			  };
	  jacobien_moy_ini /= nombre_V.nbi;
	  jacobien_moy_fin /= nombre_V.nbi;
	  traceSig_moy_ini /= nombre_V.nbi;
	  traceSig_moy /= nombre_V.nbi;
	  troisK_moy /= nombre_V.nbi;
     // d'où le calcul de la nouvelle épaisseur en utilisant la relation:
	  // (V-V_0)/V = trace(sigma)/3 /K_moy
	  if (atdt)
      {epais.epaisseur_tdt = (epais.epaisseur_t * jacobien_moy_ini) / jacobien_moy_fin
                               * troisK_moy / (troisK_moy - traceSig_moy+traceSig_moy_ini);
      
      //epais.epaisseur_tdt = (epais.epaisseur0 * jacobien_moy_0) / jacobien_moy_fin
		    //                       * troisK_moy / (troisK_moy + traceSig_moy);
//--debug
//if (num_elt==1)
// { Noeud* noe = tab_noeud(1);
//   double R_0 = noe->Coord0().Norme(); double R = noe->Coord2().Norme();
//cout << "\n e0= " << epais.epaisseur_tdt<< " troisK_moy=" << troisK_moy << " traceSig_moy=" << traceSig_moy
//     << " J0= " << jacobien_moy_0 << " J= " << jacobien_moy_fin << " R_0 " << R_0 << " R= " << R;
//	};
//-- fin debug									
		     return epais.epaisseur_tdt;
		    }
	  else  
      {epais.epaisseur_t = (epais.epaisseur0 * jacobien_moy_ini) / jacobien_moy_fin
                           * troisK_moy / (troisK_moy - traceSig_moy);
       return epais.epaisseur_t;
      };
	};

  // retourne la liste abondée de tous les  données particulières interne actuellement utilisées
// par l'élément (actif ou non), sont exclu de cette liste les données particulières des noeuds
// reliées à l'élément
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
List_io <TypeQuelconque> Biel_axiQ::Les_types_particuliers_internes(bool absolue) const
  { // on commence par récupérer la liste général provenant d'ElemMeca
    List_io <TypeQuelconque> ret = ElemMeca::Les_types_particuliers_internes(absolue);
    // ensuite on va ajouter les données particulières aux sfe
    Grandeur_scalaire_double grand_courant; // def d'une grandeur courante 
    //  $$$ cas de l'épaisseur initiale
    TypeQuelconque typQ1(EPAISSEUR_MOY_INITIALE,SIG11,grand_courant);
    ret.push_back(typQ1);
    //  $$$ cas de l'épaisseur finale 
    TypeQuelconque typQ2(EPAISSEUR_MOY_FINALE,SIG11,grand_courant);
    ret.push_back(typQ2);
	   return ret;
  }; 

// récupération de grandeurs particulières au numéro d'ordre  = iteg  
// celles-ci peuvent être quelconques
// en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void Biel_axiQ::Grandeur_particuliere (bool absolue,List_io<TypeQuelconque>& liTQ,int iteg)
   { 
     // on balaie la liste transmise pour les grandeurs propres
     List_io<TypeQuelconque>::iterator il,ilfin = liTQ.end();
     // on commence par appeler la fonction de la classe m่re
     //  il n'y aura pas de calcul des grandeurs inactivées
     ElemMeca::Grandeur_particuliere (absolue,liTQ,iteg);
     
     // puis les grandeurs sp้cifiques 
     for (il=liTQ.begin();il!=ilfin;il++)
      {TypeQuelconque& tipParticu = (*il); // pour simplifier
       if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
        switch (tipParticu.EnuTypeQuelconque().EnumTQ())
         {
          // 1) -----cas de l'้paisseur moyenne initiale, ici elle ne d้pend pas du point d'int้gration
          case EPAISSEUR_MOY_INITIALE: 
           { *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=donnee_specif.epais.epaisseur0;
             (*il).Active();
             break;
           }
          // 2) -----cas de l'้paisseur moyenne finale, ici elle ne d้pend pas du point d'int้gration
          case EPAISSEUR_MOY_FINALE: // on inactive la grandeur quelconque
           { *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=donnee_specif.epais.epaisseur_tdt;
             (*il).Active();
             break;
           }
           default: // on ne fait rien
             break;
         };
       };
    };                             	



// Calcul des frontieres de l'element
//  creation des elements frontieres et retour du tableau de ces elements
// la création n'a lieu qu'au premier appel
// ou lorsque l'on force le paramètre force a true
// dans ce dernier cas seul les frontière effacées sont recréée
Tableau <ElFrontiere*> const & Biel_axiQ::Frontiere(bool force)
  { int cas = 6; // on veut des lignes et des points
    return Frontiere_elemeca(cas,force);
   };
		

// =====>>>> methodes privées appelees par les classes dérivees <<<<=====

// fonction d'initialisation servant  au niveau du constructeur
Biel_axiQ::DonneeCommune * Biel_axiQ::Init
     (Donnee_specif donnee_spec,bool sans_init_noeud) 
{   // bien que la grandeur donnee_specif est défini dans la classe generique
    // le fait de le passer en paramètre permet de tout initialiser dans Init
    // et ceci soit avec les valeurs par défaut soit avec les bonnes valeurs
    donnee_specif =donnee_spec;
    // le fait de mettre les pointeurs a null permet
    // de savoir que l'element n'est pas complet
       tab_noeud.Change_taille(nombre_V.nbne);
    // dans le cas d'un constructeur avec tableau de noeud, il ne faut pas mettre
    // les pointeurs à nuls d'où le test
    if (!sans_init_noeud)
        for (int i =1;i<= nombre_V.nbne;i++) tab_noeud(i) = NULL;
    // definition des donnees communes aux Biel_axiQxxx
    // a la premiere definition d'une instance
    if (unefois.doCoMemb == NULL) 
       Biel_axiQ::Def_DonneeCommune();
    unefois.doCoMemb = doCo  ; 
    met = &(doCo->met_biellette); // met est defini dans elemeca
    // def pointe sur la deformation specifique a l'element pour le calcul mecanique   		
    def = new Deformation(*met,tab_noeud,(doCo->segment).TaDphi(),(doCo->segment).TaPhi());
    // idem pour la remontee aux contraintes et le calcul d'erreur   		
    defEr = new Deformation(*met,tab_noeud,(doCo->segmentEr).TaDphi(),(doCo->segmentEr).TaPhi());
    // idem pour la remontee aux contraintes et le calcul d'erreur   		
    defMas = new Deformation(*met,tab_noeud,(doCo->segmentMas).TaDphi(),(doCo->segmentMas).TaPhi());
    // idem pour le calcul de second membre 
    defArete.Change_taille(1); // 1 arrête utilisée pour  le second membre
    // la déformation sera construite si nécessaire au moment du calcul de second membre
    defArete(1) = NULL; 

    //dimensionnement des deformations et contraintes etc..
    int dimtens = 2;
    lesPtMecaInt.Change_taille_PtIntegMeca(nombre_V.nbi,dimtens);
    // attribution des numéros de référencement dans le conteneur
    for (int ni = 1; ni<= nombre_V.nbi; ni++)
       {lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage);
        lesPtMecaInt(ni).Change_Nb_ele(this->num_elt);
        lesPtMecaInt(ni).Change_Nb_pti(ni);
       };

    // stockage des donnees particulieres de la loi de comportement au point d'integ
    tabSaveDon.Change_taille(nombre_V.nbi);
    tabSaveTP.Change_taille(nombre_V.nbi);
    tabSaveDefDon.Change_taille(nombre_V.nbi);
    tab_energ.Change_taille(nombre_V.nbi);
    tab_energ_t.Change_taille(nombre_V.nbi);
    // initialisation des pointeurs définis dans la classe Element concernant les résidus et 
    // raideur
    		// --- cas de la puissance interne ---
    residu = &(doCo->residu_interne); // residu local
    raideur = &(doCo->raideur_interne); // raideur locale 
        // --- cas de la dynamique -----
    mat_masse =  &(doCo->matrice_masse);  
     // --- cas des efforts externes concernant les noeuds ------
    res_extN = &(doCo->residus_externeN); // pour les résidus et second membres
    raid_extN= &(doCo->raideurs_externeN);// pour les raideurs
     // --- cas des efforts externes concernant les aretes ------
    res_extA = &(doCo->residus_externeA); // pour les résidus et second membres
    raid_extA= &(doCo->raideurs_externeA);// pour les raideurs
	
	   return doCo;
};  


// fonction privee
// dans cette fonction il ne doit y avoir que les données communes !!
void Biel_axiQ::Def_DonneeCommune()
  {	 int nbn = nombre_V.nbne;
     // interpollation : element geometrique correspondant: nombre_V.nbi pt integ, nombre_V.nbne noeuds
     GeomSeg seg(nombre_V.nbi,nbn) ; 
     // degre de liberte
     int dim = ParaGlob::Dimension();     
    // cas des ddl éléments primaires
    // ici c'est dim-1 car seules les ddl en x1 et x2 sont gérés
     DdlElement  tab_ddl(nbn,dim-1);
     int posi = Id_nom_ddl("X1") -1;
     for (int i =1; i<= dim-1; i++)
       for (int j=1; j<= nbn; j++)
	//    tab_ddl (j,i) = Enum_ddl(i+posi);
	    tab_ddl.Change_Enum(j,i,Enum_ddl(i+posi));
     // cas des ddl éléments secondaires pour le calcul d'erreur
     // def du nombre de composantes du tenseur de contrainte en absolu
     // en fait 6 ici car les tenseurs sont 3D
     int nbcomposante = 6;
     DdlElement  tab_ddlErr(nbn,nbcomposante);
     posi = Id_nom_ddl("SIG11") -1;
     for (int j=1; j<= nbn; j++)
     // les composantes sont a suivre dans l'enumération
        for (int i= 1;i<= nbcomposante; i++)
	          tab_ddlErr.Change_Enum(j,i,Enum_ddl(i+posi));
    // egalement pour tab_Err1Sig11, def d'un tableau de un ddl : enum SIG11
    // par noeud
    DdlElement tab_Err1Sig11(nbn,DdlNoeudElement(SIG11));
    // toujours pour le calcul d'erreur definition des fonctions d'interpolation
    // pour le calcul du hession de la fonctionnelle : 
    //  points d'integration et nombre_V.nbne noeuds
    GeomSeg segEr(nombre_V.nbiEr,nbn) ; 
    // pour le calcul de la matrice masse  definition des fonctions d'interpolation
    //  en particulier pour le calcul de la masse consistante
    GeomSeg segMa(nombre_V.nbiMas,nbn) ; 
    // pour le calcul relatifs à la stabilisation d'hourglass
    GeomSeg* segmentHourg = NULL;
    if (nombre_V.nbiHour > 0)
      segmentHourg = new GeomSeg(nombre_V.nbiHour,nombre_V.nbne);
	   // def metrique
	   // on definit les variables a priori toujours utiles
    Tableau<Enum_variable_metrique> tab(24);
    tab(1) = iM0;          tab(2) = iMt;            tab(3) = iMtdt ;
    tab(4) = igiB_0;       tab(5) = igiB_t;         tab(6) = igiB_tdt;
    tab(7) = igiH_0;       tab(8) = igiH_t;         tab(9) = igiH_tdt ;
    tab(10)= igijBB_0;     tab(11)= igijBB_t;       tab(12)= igijBB_tdt;
    tab(13)= igijHH_0;     tab(14)= igijHH_t;       tab(15)= igijHH_tdt ;
    tab(16)= id_gijBB_tdt; tab(17)= id_giH_tdt;     tab(18)= id_gijHH_tdt;
    tab(19)= idMtdt ;      tab(20)= id_jacobien_tdt;tab(21)= id2_gijBB_tdt; 
    tab(22)= igradVBB_tdt; tab(23) = iVtdt;         tab(24)= idVtdt;
	   // dim du pb , nb de vecteur de la base , tableau de ddl et la def de variables
	   MetAxisymetrique2D  metri(ParaGlob::Dimension(),tab_ddl,tab,nbn) ;
    // ---- cas du calcul d'erreur sur sigma ou epsilon
    // les tenseurs sont exprimees en absolu donc nombre de composante 6
    Tableau <Vecteur *>  resEr(nbcomposante);
    for (int  i = 1;i<= nbcomposante; i++)
       resEr(i)=new Vecteur (nbn); // une composante par noeud
    Mat_pleine  raidEr(nbn,nbn); // la raideur pour l'erreur
    // dimensionnement des différents résidus et raideurs pour le calcul mécanique
    int nbddl = tab_ddl.NbDdl();
    Vecteur  residu_int(nbddl); Mat_pleine  raideur_int(nbddl,nbddl);
       // cas de la dynamique
    Mat_pleine  matmasse(1,nbddl);  // a priori on dimensionne en diagonale 
    // il y a deux extrémités mais identiques
    Tableau <Vecteur* >  residus_extN(2); residus_extN(1) = new Vecteur(dim); 
    residus_extN(2) = residus_extN(1);
    int nbddlA = nombre_V.nbneA * dim; int nbA = 1; // 1 arêtes
    Tableau <Vecteur* >  residus_extA(nbA); residus_extA(1) = new Vecteur(nbddlA); 
    Tableau <Mat_pleine* >  raideurs_extA(nbA); raideurs_extA(1) = new Mat_pleine(nbddlA,nbddlA);
    Tableau <Mat_pleine* >  raideurs_extN(2);raideurs_extN(1) = new Mat_pleine(dim,dim);
    raideurs_extN(2) = raideurs_extN(1); 
    // definition de la classe static contenant toute les variables communes aux biellettes
    doCo = new DonneeCommune(seg,tab_ddl,tab_ddlErr,tab_Err1Sig11,metri,resEr,raidEr,segEr,
                           residu_int,raideur_int,residus_extN,raideurs_extN,residus_extA,raideurs_extA
                           ,matmasse,segMa,nombre_V.nbi,segmentHourg);
 };    

// destructions de certaines grandeurs pointées, créées  au niveau de l'initialisation
void Biel_axiQ::Destruction()
  { // tout d'abord l'idée est de détruire certaines grandeurs pointées que pour le dernier élément
    if ((unefois.nbelem_in_Prog == 0)&& (unefois.doCoMemb != NULL))
      // cas de la destruction du dernier élément
      { Biel_axiQ::DonneeCommune* doCo = unefois.doCoMemb;         // pour simplifier l'écriture
        int  resErrTaille = doCo->resErr.Taille();
        for (int i=1;i<= resErrTaille;i++)
          delete doCo->resErr(i);
        delete doCo->residus_externeN(1);
        delete doCo->raideurs_externeN(1);
        delete doCo->residus_externeA(1);
        delete doCo->raideurs_externeA(1);
        if (doCo->segmentHourg != NULL) delete doCo->segmentHourg;
       }
   };