// FICHIER : ElemPoint.cc
// CLASSE : ElemPoint

// 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 "Debug.h"


# include <iostream>
using namespace std;  //introduces namespace std
#include <stdlib.h>
#include "Sortie.h"


#include "ElemPoint.h"
#include "FrontPointF.h"
#include "TypeConsTens.h"
#include "Def_Umat.h"
#include "ExceptionsLoiComp.h"

//----------------------------------------------------------------
// def des donnees commune a tous les elements
// la taille n'est pas defini ici car elle depend de la lecture 
//----------------------------------------------------------------
ElemPoint::DonneeCommune * ElemPoint::doCo_Point = NULL;
ElemPoint::UneFois  ElemPoint::unefois_Point;
ElemPoint::NombresConstruire ElemPoint::nombre_V; 
ElemPoint::ConstrucElemPoint ElemPoint::construcElemPoint;

// constructeur définissant les nombres (de noeud, de point d'integ ..)
// utilisé dans la construction des éléments
ElemPoint::NombresConstruire::NombresConstruire() 
 { nbne = 1; // le nombre de noeud de l'élément
   nbiEr = 1;// le nombre de point d'intégration pour le calcul d'erreur
   nbiMas = 1; // le nombre de point d'intégration pour le calcul de la matrice masse
 }; 

// ---- definition du constructeur de la classe conteneur de donnees communes ------------	
ElemPoint::DonneeCommune::DonneeCommune (GeomPoint& pteg,DdlElement& tab,DdlElement&  tabErr,DdlElement&  tab_Err1Sig,
	                       Met_ElemPoint&  met_point, 
	                       Tableau <Vecteur *> & resEr,Mat_pleine& raidEr,
	                       GeomPoint& pteEr,Vecteur&  residu_int,Mat_pleine&  raideur_int,
	                       Tableau <Vecteur* > & residus_extN,Tableau <Mat_pleine* >&  raideurs_extN,
	                       Mat_pleine&  mat_masse,GeomPoint& pteMa,UmatAbaqus& umatAbaqusqus,int dim_tenseur,int nbi) :
	          ptpoint(pteg),tab_ddl(tab),tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig)
	          ,met_ElemPoint(met_point)
	          ,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)
	          ,pteEr(pteEr)
	          ,residu_interne(residu_int),raideur_interne(raideur_int)
	          ,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN)
	          ,matrice_masse(mat_masse)
	          ,pteMas(pteMa)
	          ,umatAbaqus(umatAbaqusqus),inne()
	{ 
   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 (dim_tenseur);
     };
   int tailledeps = d_epsBB.Taille();
	  for (int i=1;i<= tailledeps; i++)
      d_epsBB(i) = NevezTenseurBB (dim_tenseur);
   int tailledsig = d_sigHH.Taille();
	  for (int j=1;j<= tailledsig; j++)
      d_sigHH(j) = NevezTenseurHH (dim_tenseur);
	  // liaisons entre inne et   umatAbaqus
	  inne.incre = &umatAbaqus.nb_increment;
	  inne.step = &umatAbaqus.nb_step;
	  inne.nbe = &umatAbaqus.nb_elem;
	  inne.nbpti = &umatAbaqus.nb_pt_integ;
	  inne.temps_tdt = &umatAbaqus.temps_tdt;
   inne.delta_t = &umatAbaqus.delta_t;
   inne.nom_loi = &umatAbaqus.nom_materiau;
 };

ElemPoint::DonneeCommune::DonneeCommune(DonneeCommune& a) :
	    ptpoint(a.ptpoint),tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr),tab_Err1Sig11(a.tab_Err1Sig11)
	    ,met_ElemPoint(a.met_ElemPoint),matGeom(a.matGeom),matInit(a.matInit)
	    ,d2_epsBB(a.d2_epsBB),resErr(a.resErr),raidErr(a.raidErr),pteEr(a.pteEr)
	    ,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)
	    ,matrice_masse(a.matrice_masse),pteMas(a.pteMas)
	    ,umatAbaqus(a.umatAbaqus),inne()
	        { 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)));
	          // liaisons entre inne et   umatAbaqus
	          inne.incre = &umatAbaqus.nb_increment;
	          inne.step = &umatAbaqus.nb_step;
	          inne.nbe = &umatAbaqus.nb_elem;
	          inne.nbpti = &umatAbaqus.nb_pt_integ;
	          inne.temps_tdt = &umatAbaqus.temps_tdt;
           inne.delta_t = &umatAbaqus.delta_t;
           inne.nom_loi = &umatAbaqus.nom_materiau;
         }; 
ElemPoint::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 -+-+-+
ElemPoint::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)
	,CalDynamique(0),CalPt_0_t_tdt(0)
	,nbelem_in_Prog(0)
	        {};
ElemPoint::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 ------------	

// Constructeur par defaut
// les tenseurs on par défaut (dimension==-1) la dimension de l'espace
// sinon il ont la dimension donnée
ElemPoint::ElemPoint (int dimension)  :
   ElemMeca(),nbi(1),nb_appelsCalculUmat(0),lesPtMecaInt()

{ unefois=&unefois_Point; // le cas d'un vrai ElemPoint
  if (unefois->nbelem_in_Prog == 0)
   { unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
   }
  else  // sinon on construit
   {id_interpol=CONSTANT;  //   donnees de la classe mere
	   id_geom=POINT;
    Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE);     //mise à jour du type de pb
    tab_noeud.Change_taille(nombre_V.nbne);
	   // initialisation par défaut
	   int dim_tenseur = ParaGlob::Dimension();
	   if (dimension != -1) dim_tenseur = dimension;
    Init (dim_tenseur);
    unefois->nbelem_in_Prog++;
   };
};  

// Constructeur fonction  d'un numero de maillage et d'identification  
ElemPoint::ElemPoint (int num_mail,int num_id,int dimension) :
  ElemMeca(num_mail,num_id,CONSTANT,POINT),nbi(1),nb_appelsCalculUmat(0)
  ,lesPtMecaInt()
 {unefois=&unefois_Point; // le cas d'un vrai ElemPoint
  if (unefois->nbelem_in_Prog == 0)
   { unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
   }
  else  // sinon on construit
   {tab_noeud.Change_taille(nombre_V.nbne);
	   // initialisation par défaut
	   int dim_tenseur = ParaGlob::Dimension();
	   if (dimension != -1) dim_tenseur = dimension;
    unefois->nbelem_in_Prog++;
    Init (dim_tenseur);
    Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE);     //mise à jour du type de pb
   };
 };

  //------- pour des classes dérivées ------
  // CONSTRUCTEURS :
  // les tenseurs on par défaut (dimension==-1) la dimension de l'espace
  // sinon il ont la dimension donnée
ElemPoint::ElemPoint (ElemPoint::UneFois& unefoisfois, Enum_geom nouveau_id, int dimension)  :
   ElemMeca(),nbi(1),nb_appelsCalculUmat(0),lesPtMecaInt()

{ unefois=&unefoisfois; // le cas d'un vrai ElemPoint
  if (unefois->nbelem_in_Prog == 0)
   { unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
   }
  else  // sinon on construit
   {id_interpol=CONSTANT;  //   donnees de la classe mere
    id_geom=nouveau_id;
    Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE);     //mise à jour du type de pb
    tab_noeud.Change_taille(nombre_V.nbne);
    // initialisation par défaut
    int dim_tenseur = ParaGlob::Dimension();
    if (dimension != -1) dim_tenseur = dimension;
    unefois->nbelem_in_Prog++;
    Init (dim_tenseur);
   };
};
  // Constructeur fonction  d'un numero de maillage et d'identification
ElemPoint::ElemPoint (ElemPoint::UneFois& unefoisfois, Enum_geom nouveau_id, int num_mail,int num_id,int dimension) :
  ElemMeca(num_mail,num_id,CONSTANT,nouveau_id),nbi(1),nb_appelsCalculUmat(0)
  ,lesPtMecaInt()
 {unefois=&unefoisfois; // le cas d'un vrai ElemPoint
  if (unefois->nbelem_in_Prog == 0)
   { unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
   }
  else  // sinon on construit
   {tab_noeud.Change_taille(nombre_V.nbne);
    // initialisation par défaut
    int dim_tenseur = ParaGlob::Dimension();
    if (dimension != -1) dim_tenseur = dimension;
    unefois=&unefoisfois; // le cas d'un vrai ElemPoint
    unefois->nbelem_in_Prog++;
    Init (dim_tenseur);
    Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE);     //mise à jour du type de pb
   };
 };


ElemPoint::ElemPoint (const ElemPoint& a) :
 ElemMeca (a),nbi(a.nbi),nb_appelsCalculUmat(0)
 ,lesPtMecaInt(a.lesPtMecaInt) 
// Constructeur de copie
// a priori si on utilise le constructeur de copie, donc il y a déjà un élément
// par contre a priori on ne doit pas faire une copie du premier élément
{unefois= a.unefois; // le cas d'un vrai ElemPoint
 if (unefois->nbelem_in_Prog == 1)
   { cout << "\n **** erreur pour l'element ElemPoint, le constructeur de copie ne doit pas etre utilise"
          << " pour le premier element !! " << endl;
     Sortie (1);
   }
 else
   {// initialisation
    unefois = a.unefois;
	   int dim_tenseur=lesPtMecaInt(1).EpsBB()->Dimension();
    unefois->nbelem_in_Prog++;
    Init (dim_tenseur);
   };
};

ElemPoint::~ElemPoint ()
// Destructeur
{ LibereTenseur();
  unefois->nbelem_in_Prog--;
  Destruction();
};

// Lecture des donnees de la classe sur fichier
void 
ElemPoint::LectureDonneesParticulieres
    (UtilLecture * entreePrinc,Tableau<Noeud  *> * tabMaillageNoeud)
  { int nb;
    tab_noeud.Change_taille(1);	 
	   int i=1;
	   // on ne continue que si le tableau de noeud passé en paramètre a une taille non nulle
	   if ( (*tabMaillageNoeud).Taille() != 0)
     {*(entreePrinc->entree) >> nb;
	   if ( ((entreePrinc->entree)->rdstate() == 0) || ((entreePrinc->entree)->eof()) )
	   // ici il n'y a qu'une seule information à lire, donc dès le début on a
		// io_state == 1
 
       // 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 != 2) 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 << "ElemPoint::LectureDonneesParticulieres"; 
          Affiche();
          Sortie (1);
        };
     // construction du tableau de ddl des noeuds de ElemPoint
     ConstTabDdl(); 
	  };
  };
  		
// on associe un noeud à l'élément, remplace le noeud existant, s'il existe déjà
void ElemPoint::Associer_noeud (Noeud  *  noeu)
{ // on dimensionne le tableau de noeud éventuellement
  tab_noeud.Change_taille(1);
  // on réaffecte
  tab_noeud(1)=noeu;
  // construction du tableau de ddl des noeuds de ElemPoint
  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 & ElemPoint::Point_physique(const Coordonnee& c_int,Coordonnee & co,Enum_dure temps)
 { ElemPoint::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_ElemPoint.PlusInitVariables(tab) ;
      };  
   // b) calcul de l'interpolation 
   const Vecteur& phi = doCo->ptpoint.Phi_point(c_int);
   // c) calcul du point
   switch (temps)
    { case TEMPS_0 : co = doCo->met_ElemPoint.PointM_0(tab_noeud,phi); break;
      case TEMPS_t : co = doCo->met_ElemPoint.PointM_t(tab_noeud,phi); break;
      case TEMPS_tdt : co = doCo->met_ElemPoint.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 ElemPoint::Point_physique(const Coordonnee& c_int,Tableau <Coordonnee> & t_co)
 { ElemPoint::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_ElemPoint.PlusInitVariables(tab) ;
      };  
   // b) calcul de l'interpolation 
   const Vecteur& phi = doCo->ptpoint.Phi_point(c_int);
   // c) calcul des point
   switch (t_co.Taille())
    { case 3 : t_co(3) = doCo->met_ElemPoint.PointM_tdt(tab_noeud,phi); 
      case 2 : t_co(2) = doCo->met_ElemPoint.PointM_t(tab_noeud,phi); 
      case 1 : t_co(1) = doCo->met_ElemPoint.PointM_0(tab_noeud,phi); 
    };
  };
  
// Calcul du residu local à t ou tdt en fonction du booleen atdt
 Vecteur* ElemPoint::CalculResidu (bool atdt,const ParaAlgoControle & pa)
    { ElemPoint::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_ElemPoint.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_ElemPoint.PlusInitVariables(tab) ;
         };
        };
      // initialisation du résidu 
      residu->Zero();  
      Vecteur poids =(doCo->ptpoint).TaWi();   // poids d'interpolation  = 2
      ElemMeca::Cal_explicit ( doCo->tab_ddl,d_epsBB,1,poids,pa,atdt);
      return residu;  
    };

// Calcul du residu local et de la raideur locale,
//  pour le schema implicite
Element::ResRaid  ElemPoint::Calcul_implicit (const ParaAlgoControle & pa)
    { ElemPoint::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_ElemPoint.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->ptpoint).TaWi();   
// --- pour l'instant il n'y a pas de calcul, on retourne les matrices vides		 
//       ElemMeca::Cal_implicit 
//           (doCo->tab_ddl, d_epsBB,(doCo->d2_epsBB),d_sigHH,nbi,poids,pa,cald_Dvirtuelle);
      Element::ResRaid el;
      el.res = residu; 
      el.raid = raideur;
      return el;
    };

// calcul de l'UMat pour abaqus
void ElemPoint::CalculUmatAbaqus(ParaAlgoControle & pa)
	{ // -- récupération des variables Umat abaqus
	  UmatAbaqus& umatA = unefois->doCoMemb->umatAbaqus; // pour simplifier
	  // cas du premier passage
	  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; 
       unefois->doCoMemb->met_ElemPoint.PlusInitVariables(tab) ;
      };
	  // -- appel de la loi de comportement
	  int ni=umatA.nb_pt_integ; // num du pt d'integ (pour simplifier)
   // on change les coordonnées du point associé au premier pti
   // (a priori ce n'est pas utilisé ici mais on le fait quand même pour le futur)
   tab_noeud(1)->Change_coord2(umatA.coor_pt);
   // on change la température au noeud en fonction de ce qui existe à l'umat
   // car ce sera utilisé pour les grandeurs interpolées aux noeuds
   if (tab_noeud(1)->Existe_ici(TEMP))
     {tab_noeud(1)->Change_val_t(TEMP,umatA.temper_t);
      tab_noeud(1)->Change_val_tdt(TEMP,umatA.temper_t + umatA.delta_temper);
     };  

   // on change le numéro d'intégration pour les calculs d'interpolation c-a-d des acces
   // aux tableaux phi et dphi : pour le point il n'y a qu'un pti !!
	  def->ChangeNumInteg(1); //ni);
  
   // on associe la déformation à la bonne sauvegarde
	  def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
	  // mise à jour de la métrique avec les infos de umat
	  ((Def_Umat*)def)->Mise_a_jourUmat(umatA);
	  // valeurs aux pt d'integ 
   PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);

   // on regarde si c'est le premier incrément ou pas (d'où le  calcul du premier pas)
   premier_calcul_meca_impli_expli=false; 
   if (umatA.nb_increment <= 1) premier_calcul_meca_impli_expli  = true;
   CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL;  // les données spécifique thermo physiques
   if (loiTP != NULL) {sTP = tabSaveTP(ni);};         // au pt d'integ si elles existes 
   #ifdef MISE_AU_POINT
   if ( ParaGlob::NiveauImpression() > 5)
     {cout << "\n ElemPoint::CalculUmatAbaqus(... "
           << " dilatation= "<< dilatation
           << flush;
     };
   #endif
   // on gère les exceptions éventuelles en mettant le bloc sous surveillance
   try 
    { // appel du comportement Umat
	     loiComp->ComportementUmat(tabSaveDon(ni),*def,ptIntegMeca,pa,sTP,loiTP,dilatation,umatA
	                              ,premier_calcul_meca_impli_expli);
    }     
   catch (ErrNonConvergence_loiDeComportement excep) 
         // cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
         // d'une loi de comportement 
    { // l'idée est que dans le cas où il y a un pb avec la loi on signale qu'il faut diminuer le pas de temps ??
	     // c'est surtout pour une utilisation avec Abaqus 
	     umatA.pnewdt = sqrt(2.)/2.;
      if (ParaGlob::NiveauImpression() >= 1) 
           cout << "\n warning: CalculUmatAbaqus: modification du pas de temps -> du parametre pnewt";
    }
   catch (ErrSortieFinale)
        // cas d'une direction voulue vers la sortie
        // on relance l'interuption pour le niveau supérieur
      { ErrSortieFinale toto;
        throw (toto);
      }
   catch ( ... )  
    { if (ParaGlob::NiveauImpression() >= 1) 
        {cout << "\n warning: exception generee par la loi de comportement sous Umat mais dont la prise en compte "
              << " n'est pas prevu !, on demande une diminution du pas de temps de 0.5 ";
//              << " n'est pas prevu !, on ne fait rien et on continue le calcul";
	        umatA.pnewdt = 0.5;
         if (ParaGlob::NiveauImpression() >= 4) cout <<  "\n ElemPoint::CalculUmatAbaqus(..";
        };
     };                	                             
	                              
	  // comptabilisation des appels si c'est le premier point d'intégration
	  if (ni == 1)
	    nb_appelsCalculUmat++;                             						
	};
	
// initialisation éventuelle: ajout de point d'intégration si nécessaire
// les tenseurs on par défaut (dim_tens==-1) la dimension de l'espace
void ElemPoint::InitialisationUmatAbaqus_interne(int dim_tens)
{ // dans le cas où le nombre de point d'intégration lue est supérieur à celui de l'élément
  // on agrandi les tableaux (qui sont des tableaux de pointeurs
  if (unefois->doCoMemb->umatAbaqus.nb_pt_integ > nbi)
   { 	int dim_tenseur = ParaGlob::Dimension();
	     if (dim_tens != -1) dim_tenseur = dim_tens;
      ChangeNombrePtinteg(unefois->doCoMemb->umatAbaqus.nb_pt_integ,dim_tenseur);
   };
};
		
// Calcul de la matrice masse pour l'élément
Mat_pleine * ElemPoint::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse) 
 {    ElemPoint::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_ElemPoint.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->pteMas).TaWi();   // poids d'intégration  = 2
      // appel de la routine générale 
      ElemMeca::Cal_Mat_masse (doCo->tab_ddl,type_calcul_masse,
                           nombre_V.nbiMas,(doCo->pteMas).TaPhi(),nombre_V.nbne
                           ,poids);
      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 ElemPoint::ContrainteAuNoeud_ResRaid()
    {  DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
       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_ElemPoint.PlusInitVariables(tab) ;
         };  
       // appel du programme général
       ElemMeca:: SigmaAuNoeud_ResRaid(tab_noeud.Taille()
                 ,(doCo->ptpoint).TaPhi()
                 ,(doCo->ptpoint).TaWi()
                 ,doCo-> resErr,doCo->raidErr
                 ,(doCo->pteEr).TaPhi()
                 ,(doCo->pteEr).TaWi());
       return  (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr)));       
      };
      
    // 2) remontée aux erreurs aux noeuds
 Element::Er_ResRaid ElemPoint::ErreurAuNoeud_ResRaid()
    {  DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
       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_ElemPoint.PlusInitVariables(tab) ;
       };  
       // appel du programme général
       ElemMeca::Cal_ErrAuxNoeuds(tab_noeud.Taille(), (doCo->ptpoint).TaPhi(), 
               (doCo->ptpoint).TaWi(),doCo-> resErr );
       return  (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr)));       
    };
     
// actualisation des ddl et des grandeurs actives de t+dt vers t      
void ElemPoint::TdtversT()
 { for (int ni=1;ni<= nbi; ni++)
     { lesPtMecaInt.TdtversT(); 
       if (tabSaveDon(ni) != NULL)  tabSaveDon(ni)->TdtversT();
       if (tabSaveTP(ni) != NULL)  tabSaveTP(ni)->TdtversT();
       if (tabSaveDefDon(ni) != NULL)  tabSaveDefDon(ni)->TdtversT();
     };
   ElemMeca::TdtversT_(); // appel de la procédure mère
   // normalement ça correspond à la convergence, on met à -1 le nombre d'appel
   // ainsi au premier appel on aura l'itération 0
   nb_appelsCalculUmat = -1;
  };
// actualisation des ddl et des grandeurs actives de t vers tdt      
void ElemPoint::TversTdt()
 { for (int ni=1;ni<= nbi; ni++)
    { lesPtMecaInt.TversTdt(); 
      if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->TversTdt();
      if (tabSaveTP(ni) != NULL) tabSaveTP(ni)->TversTdt();
      if (tabSaveDefDon(ni) != NULL) tabSaveDefDon(ni)->TversTdt();
    };
   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 ElemPoint::ErreurElement(int type,double& errElemRelative
                            ,double& numerateur, double& denominateur)
    {  DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
       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_ElemPoint.PlusInitVariables(tab) ;
        };  
       // appel du programme général
       ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur,
                            tab_noeud.Taille(),(doCo->ptpoint).TaPhi(),
                (doCo->ptpoint).TaWi(),
                (doCo->pteEr).TaPhi(),(doCo->pteEr).TaWi());
    };
	
//============= 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 ElemPoint::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 ElemPoint
  // car on n'a pas dimensionné les tableaux qui suivent !!
  switch (cas)
  { case 1 : // ------- on récupère tout -------------------------
     { string toto;
	      int nbi_old = nbi;
       ent >> toto >>  nbi ; // récup du nombre de point d'intégration
       // construction du tableau de ddl des noeuds de ElemPoint
       ConstTabDdl(); 
       // des déformations et contraintes éventuelles
       int dim_tens;
       ent >> dim_tens; // lecture de l'indicateur pour savoir s'il y a des tenseurs ou pas
		     if (nbi != nbi_old)
		     ChangeNombrePtinteg(nbi,dim_tens);
		     // récupération
		     lesPtMecaInt.Lecture_base_info(ent,cas);
	      break;
	    }
    case  2 : // ----------- lecture uniquement de se qui varie --------------------
     { // récup info
       // des déformations et contraintes éventuelles
       int dim_tens;string toto;
       ent >> toto >> dim_tens; // lecture de l'indicateur pour savoir s'il y a des tenseurs ou pas
		     // récupération
		     lesPtMecaInt.Lecture_base_info(ent, cas);
	      break;
	    }
    default :
	    { cout << "\nErreur : valeur incorrecte du type de lecture !\n";
	      cout << "ElemPoint::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 ElemPoint::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 ElemPoint
 // sinon on essaie de sauvegarder
 switch (cas)
  { case 1 : // ------- on sauvegarde tout -------------------------
     {sort << "\n nbpti= " << nbi << " ";
		    sort << lesPtMecaInt.DimTens() << " ";
		    // des tenseurs déformation et contrainte,
		    lesPtMecaInt.Ecriture_base_info(sort,cas);
	     break;
	    }
    case  2 : // ----------- sauvegarde uniquement de se qui varie --------------------
     { //  la déformation et  la contrainte
       sort << "\n dim_tenseur: " << lesPtMecaInt.DimTens() << " ";
		     // des tenseurs déformation et contrainte,
		     lesPtMecaInt.Ecriture_base_info(sort,cas);
       break;
	    }
    default :
	    { cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
	      cout << "ElemPoint::Ecriture_base_info(ofstream& sort,int cas)"
			         << " cas= " << cas << endl;
	      Sortie(1);
	    };
 };
};
 
      
// Calcul de la matrice géométrique et initiale
ElemMeca::MatGeomInit ElemPoint::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & pa) 
    { ElemPoint::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_ElemPoint.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->ptpoint).TaWi();   // poids d'interpolation  = 2
    ElemMeca::Cal_matGeom_Init (matGeom,matInit,doCo->tab_ddl, d_epsBB,
                              doCo->d2_epsBB,d_sigHH,1,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 & ElemPoint::TableauDdl() const
     { DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
       return doCo->tab_ddl; };
// liberation de la place pointee     
void ElemPoint::Libere ()
  {Element::Libere (); // liberation de residu et raideur
   LibereTenseur() ; // liberation des tenseur intermediaires
  };

// acquisition ou modification d'une loi de comportement
void ElemPoint::DefLoi (LoiAbstraiteGeneral * NouvelleLoi)
    { // verification du type de loi
      int dim=ParaGlob::Dimension(); // pour l'instant par défaut
     
// non avec les CP on passe en 2D au niveau de la loi, donc a priori on peut avoir plusieurs dimensions
//      if (NouvelleLoi->Dimension_loi() != dim)
//        { cout << "\n Erreur, la loi de comportement a utiliser avec des ElemPoints";
//          cout << " est de type = " << (NouvelleLoi->Dimension_loi()) << " au lieu de " << dim << "!!! " << endl;
//          Sortie(1);
//        };

      // cas d'une loi mécanique
      if (GroupeMecanique(NouvelleLoi->Id_categorie()))
       {loiComp = (Loi_comp_abstraite *) NouvelleLoi;
        // initialisation du stockage particulier, ici 1 pt d'integ
        tabSaveDon(1) = loiComp->New_et_Initialise();  
        // idem pour le type de déformation mécanique associé            
        tabSaveDefDon(1) = 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;
        // initialisation du stockage particulier, ici 1 pt d'integ
        tabSaveTP(1) = 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 ElemPoint::TestComplet()
  { int res = ElemMeca::TestComplet(); // test dans la fonction mere
    return res;   
  };
    
// ajout du tableau de ddl des noeuds de ElemPoint
 void ElemPoint::ConstTabDdl()
  {  
    Tableau <Ddl> ta(ParaGlob::Dimension());
    int posi = Id_nom_ddl("X1") -1;    
    int dim = ParaGlob::Dimension();    
    for (int i =1; i<= dim; i++)
	    {Ddl inter((Enum_ddl(i+posi)),0.,LIBRE);
	     ta(i) = inter;
     };
    // attribution des ddls aux noeuds    
    tab_noeud(1)->PlusTabDdl(ta); 
  };
		
// procesure permettant de completer l'element apres
// sa creation avec les donnees du bloc transmis
// peut etre appeler plusieurs fois
 Element* ElemPoint::Complete(BlocGen & bloc,LesFonctions_nD*  lesFonctionsnD)
    { // complétion avec bloc
      return ElemMeca::Complete_ElemMeca(bloc,lesFonctionsnD);
    };       

// lecture de données diverses sur le flot d'entrée
void ElemPoint::LectureContraintes(UtilLecture * entreePrinc) 
 { DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
   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_ElemPoint.PlusInitVariables(tab) ;
        };
	  ElemMeca::LectureDesContraintes (true,entreePrinc,lesPtMecaInt.TabSigHH_t());
 };
        
// retour des contraintes en absolu retour true si elle existe sinon false
bool ElemPoint::ContraintesAbsolues(Tableau <Vecteur>& tabSig)
 { DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
   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_ElemPoint.PlusInitVariables(tab) ;
         };  
	  ElemMeca::ContraintesEnAbsolues(true,lesPtMecaInt.TabSigHH_t(),tabSig);
	  return true; 
 };
               
// affichage dans la sortie transmise, des variables duales "nom"
// dans le cas ou nom est vide, affichage de "toute" les variables
 void ElemPoint::AfficheVarDual(ofstream& sort, Tableau<string>& nom)
   {// affichage de l'entête de l'element
    sort << "\n******************************************************************";
    sort << "\n Element bielette (2 noeuds 1 point 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     
	 };
          
  // recuperation des coordonnées du point de numéro d'ordre iteg pour 
  // la grandeur enu
// temps: dit si c'est à 0 ou t ou tdt
  // si erreur retourne erreur à true
  Coordonnee ElemPoint::CoordPtInteg(Enum_dure temps,Enum_ddl enu,int iteg,bool& erreur)
    { // ici les coordonnées du point d'intégration et celles du noeud sont identique
	     // donc on ramène les coordonnées du noeud associé
		    switch (temps) {
			    case TEMPS_0: {return (tab_noeud(1)->Coord0()); break;};
       case TEMPS_t: {return (tab_noeud(1)->Coord1()); break;};
			    case TEMPS_tdt:{ return (tab_noeud(1)->Coord2()); break;};
      };
    };
 
 // 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> ElemPoint::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 prevu, unefois->dualSortbiel= " << unefois->dualSortbiel
               << " unefois->CalimpPrem= " << unefois->CalimpPrem
               << "\n ElemPoint::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 ElemPoint::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 ElemPoint::ValTensorielle_a_diff_temps(Enum_dure enu_t...";
          Sortie(1);
        };       
		    ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas);
    };

// 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 & ElemPoint::Frontiere(bool force)
  { int cas = 3; // on veut des points
    return Frontiere_elemeca(cas,force);
//  { // le calcul et la création ne sont effectués qu'au premier appel
//    // ou lorsque l'on veut forcer une recréation
//	if (((ind_front_lin == 0) && (ind_front_surf == 0) && (ind_front_point == 0))
//	    || force )
////     if ((ind_front_point == 0) || force || (ind_front_point == 2)) 
//     { // dimensionnement des tableaux intermediaires
//       Tableau <Noeud *> tab(1); // les noeuds des points frontieres
//       DdlElement ddelem(1);  // les ddlelements des points frontieres
//       int tail=1;
//       tabb.Change_taille(tail); // le tableau total de frontières
//
//       // le  point  
////       tab(1) = tab_noeud(1);
////       ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(1));
////       if (tabb(1+posi_tab_front_point) == NULL)
////         tabb(1+posi_tab_front_point) = new  FrontPointF (tab,ddelem);
//       // mise à jour des indicateurs
////       ind_front_point = 1;    
//      }  
//    
//    return tabb;          
   };
		
//// ramène la frontière point
//// éventuellement création des frontieres points de l'element et stockage dans l'element 
//// si c'est la première fois  sinon il y a seulement retour de l'elements
//// a moins que le paramètre force est mis a true
//// dans ce dernier cas la frontière effacéee est recréée
//// num indique le numéro du point à créer (numérotation EF)
//ElFrontiere* const  ElemPoint::Frontiere_points(int num,bool force) 
// { // le calcul et la création ne sont effectués qu'au premier appel
//   // ou lorsque l'on veut forcer une recréation
//   #ifdef MISE_AU_POINT
//      if (num != 1)
//        { cout << "\n *** erreur, pour les ElemPoints il n'y a qu'une frontieres point ! "
//               << "\n Frontiere_points(int num,bool force)";
//          Sortie(1);
//         }       
//   #endif
//
//   if ((ind_front_point == 0) || force || (ind_front_point == 2)) 
//     {Tableau <Noeud *> tab(1); // les noeuds des points frontieres
//      DdlElement ddelem(1);  // les ddlelements des points frontieres
//      //   on regarde si la frontière point existe sinon on la crée
//      if (ind_front_point == 1)
//        return (ElFrontiere*)tabb(posi_tab_front_point+num); 
//      else if ( ind_front_point == 2)
//        // cas où certaines frontières existent
//        if (tabb(posi_tab_front_point+num) != NULL)            
//          return (ElFrontiere*)tabb(posi_tab_front_point+num);
//      // dans tous les autres cas on construit la frontière point
//      // on commence par dimensionner le tableau de frontière
//      if  (ind_front_point==0)
//        // cas où aucune frontière existe, on crée pour les points
//        { tabb.Change_taille(1); posi_tab_front_point = 0;}; 
//      // dans les deux autres cas  ((ind_front_lin == 0) && (ind_front_point>0)) et
//      //  ((ind_front_lin > 0) && (ind_front_point>0)), les points existant déjà on n'a rien n'a faire
//      // on définit le point
//      // premier point  
////      tab(1) = tab_noeud(1);
////      ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(1));
////      if (tabb(1+posi_tab_front_point) == NULL)
////         tabb(1+posi_tab_front_point) = new  FrontPointF (tab,ddelem);
//      ind_front_point=1; // mise à jour de l'indicateur
//     };
//    return (ElFrontiere*)tabb(num+posi_tab_front_point); 
//   };
//		
//// ramène la frontière linéique
//// éventuellement création des frontieres linéique de l'element et stockage dans l'element 
//// si c'est la première fois et en 3D sinon il y a seulement retour de l'elements
//// a moins que le paramètre force est mis a true
//// dans ce dernier cas la frontière effacéee est recréée
//// num indique le numéro de l'arête à créer (numérotation EF)
//// ici normalement la fonction ne doit pas être appelée
//ElFrontiere* const  ElemPoint::Frontiere_lineique(int ,bool )
//	{ cout << "\n *** erreur, pour les ElemPoints il n'y a pas de frontiere lineique ! "
//           << "\n Frontiere_lineique(int ,bool force = false)";
//      Sortie(1);
//	  return NULL;
//	 };
//		
//// ramène la frontière surfacique
//// éventuellement création des frontieres surfacique de l'element et stockage dans l'element 
//// si c'est la première fois sinon il y a seulement retour de l'elements
//// a moins que le paramètre force est mis a true
//// dans ce dernier cas la frontière effacéee est recréée
//// num indique le numéro de la surface à créer (numérotation EF) 
//// ici normalement la fonction ne doit pas être appelée
//ElFrontiere* const  ElemPoint::Frontiere_surfacique(int ,bool )
//	{ cout << "\n *** erreur, pour les ElemPoints il n'y a pas de frontiere surface ! "
//           << "\n Frontiere_surfacique(int ,bool force = false)";
//      Sortie(1);
//	  return NULL;
//	 };

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

// fonction d'initialisation servant  au niveau du constructeur
// dim_tenseur: = dimension des tenseurs
void ElemPoint::Init(int dim_tenseur) 
{ // 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);
	 for (int i =1;i<= nombre_V.nbne;i++) tab_noeud(i) = NULL;
	 // definition des donnees communes aux ElemPointxxx 
	 // a la premiere definition d'une instance
	 if (unefois->doCoMemb == NULL)
	   unefois->doCoMemb = ElemPoint::Def_DonneeCommune(dim_tenseur);
//	 unefois->doCoMemb = doCo  ;
	 met = &(unefois->doCoMemb->met_ElemPoint); // met est defini dans elemeca
	 DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
  // def pointe sur la deformation specifique a l'element pour le calcul mecanique
	 def = new Def_Umat(*met,tab_noeud,(doCo->ptpoint).TaDphi(),(doCo->ptpoint).TaPhi());
	 // idem pour la remontee aux contraintes et le calcul d'erreur   		
	 defEr = new Def_Umat(*met,tab_noeud,(doCo->pteEr).TaDphi(),(doCo->pteEr).TaPhi());
	 // idem pour la remontee aux contraintes et le calcul d'erreur   		
	 defMas = new Def_Umat(*met,tab_noeud,(doCo->pteMas).TaDphi(),(doCo->pteMas).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
	 lesPtMecaInt.Change_taille_PtIntegMeca(nbi,dim_tenseur);
  // attribution des numéros de référencement dans le conteneur
  for (int ni = 1; ni<= nbi; ni++)
     {lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage);
      lesPtMecaInt(ni).Change_Nb_ele(this->num_elt);
      lesPtMecaInt(ni).Change_Nb_pti(ni);
     };
 
	 lesPtIntegMecaInterne = &lesPtMecaInt;	//  association avec le pointeur d'ElemMeca

  // stockage des donnees particulieres de la loi de comportement au point d'integ
  tabSaveDon.Change_taille(nbi);
  tabSaveTP.Change_taille(nbi);
  tabSaveDefDon.Change_taille(nbi);
  tab_energ.Change_taille(nbi);
  tab_energ_t.Change_taille(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
};


// fonction privee
// dans cette fonction il ne doit y avoir que les données communes !!
ElemPoint::DonneeCommune* ElemPoint::Def_DonneeCommune(int dim_tenseur)
  {	 int nbn = nombre_V.nbne;
     // interpollation : element geometrique de base correspondant: 1 pt integ, 1 noeuds
     GeomPoint pteg ; 
     // degre de liberte
     int dim = ParaGlob::Dimension();     
     DdlElement  tab_ddl(nbn,dim);
     int posi = Id_nom_ddl("X1") -1;
     for (int i =1; i<= dim; 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
     int nbcomposante = 0; // init ParaGlob::NbCompTens ();
     int nb_vect=0; // init: nb de vecteur de la base
     switch (dim_tenseur)
      { case 3: nbcomposante=6;nb_vect=3; break;
        case 2: nbcomposante=3;nb_vect=2; break;
        case 1: nbcomposante=1;nb_vect=1; break;
        default: break;
      };
     DdlElement  tab_ddlErr(nbn,nbcomposante);
     posi = Id_nom_ddl("SIG11") -1;
     for (int j=1; j<= nbn; j++)
      {// on definit le nombre de composante de sigma en absolu
       if (nbcomposante == 3) 
        // dans l'énumération les composantes ne sont pas a suivre
        {//tab_ddlErr (j,1) = Enum_ddl(1+posi);  // cas de SIG11
         //tab_ddlErr (j,2) = Enum_ddl(2+posi);  // cas de SIG22
         //tab_ddlErr (j,3) = Enum_ddl(4+posi);  // cas de SIG12
	        tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi)); // cas de SIG11
	        tab_ddlErr.Change_Enum(j,2,Enum_ddl(2+posi)); // cas de SIG22
	        tab_ddlErr.Change_Enum(j,3,Enum_ddl(4+posi)); // cas de SIG12
        }
       else if (nbcomposante == 6)
         // les composantes sont a suivre dans l'enumération   
        {for (int i= 1;i<= nbcomposante; i++)
  //       tab_ddlErr (j,i) = Enum_ddl(i+posi);
	        tab_ddlErr.Change_Enum(j,i,Enum_ddl(i+posi));
        }
       else
       // tab_ddlErr (j,1) = Enum_ddl(1+posi); // uniquement SIG11  
        {tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi));}; // uniquement SIG11
      };
     // 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 : 
     // 1 points d'integration et 1 noeuds
     GeomPoint ptegEr ; 
     // pour le calcul de la matrice masse  definition des fonctions d'interpolation
     // 1 points d'integration et 1 noeuds, en particulier pour le calcul de la masse consistante
     GeomPoint ptegMa ; 
  // 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
	    Met_ElemPoint  metri(ParaGlob::Dimension(),nb_vect,tab_ddl,tab,nbn) ;
		// ---- cas du calcul d'erreur sur sigma ou epsilon
		// les tenseurs sont exprimees en absolu donc nombre de composante fonction 
		// de la dimension absolue
     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 une extrémité a priori
     Tableau <Vecteur* >  residus_extN(1); residus_extN(1) = new Vecteur(dim); 
     // pas d'arêtes, un noeud
     Tableau <Mat_pleine* >  raideurs_extN(1);raideurs_extN(1) = new Mat_pleine(dim,dim);
     // -- particularités pour la routine Umat pour Abaqus
     UmatAbaqus umatAbaqus(ParaGlob::Dimension(),dim_tenseur);
     // definition de la classe static contenant toute les variables communes aux ElemPoints
     int nbi=1; // par défaut on met un pt d'integ, ensuite il faudra remettre à jour 
                // les données communes si le nb de pt augmente
     DonneeCommune* doCoco = new DonneeCommune(pteg,tab_ddl,tab_ddlErr,tab_Err1Sig11,metri,resEr,raidEr,ptegEr
                            ,residu_int,raideur_int,residus_extN,raideurs_extN
                            ,matmasse,ptegMa,umatAbaqus,dim_tenseur,nbi);
     return doCoco;
 };    

// destructions de certaines grandeurs pointées, créées  au niveau de l'initialisation
void ElemPoint::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))
//    if (unefois->unefois->nbelem_in_Prog == 0)
      // cas de la destruction du dernier élément
      { ElemPoint::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);
      };
   };
   
// changement du nombre de point d'intégration
void ElemPoint::ChangeNombrePtinteg(int nevez_nbi, int dim_tens)
{ // les numéroq
  int old_nbi = nbi;
  nbi = nevez_nbi;
  // stockage des donnees particulieres de la loi de comportement au point d'integ
  tabSaveDon.Change_taille(nbi);
  tabSaveTP.Change_taille(nbi);
  tabSaveDefDon.Change_taille(nbi);
  // puis les tenseurs internes
  lesPtMecaInt.Change_taille_PtIntegMeca(nbi,dim_tens);
  // attribution des numéros de référencement dans le conteneur
  for (int ni = 1; ni<= nbi; ni++)
     {lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage);
      lesPtMecaInt(ni).Change_Nb_ele(this->num_elt);
      lesPtMecaInt(ni).Change_Nb_pti(ni);
     };
     
  for (int ni = old_nbi+1; ni <= nbi; ni++)
   {  // cas d'une loi mécanique
     // initialisation du stockage particulier pour le nouveau pt d'integ
     tabSaveDon(ni) = loiComp->New_et_Initialise();  
     // idem pour le type de déformation mécanique associé            
     tabSaveDefDon(ni) = def->New_et_Initialise(); 
     // cas éventuelle d'une loi thermo physique
     if (loiTP != NULL)
       {// initialisation du stockage particulier, ici 1 pt d'integ
        tabSaveTP(ni) = loiTP->New_et_Initialise();  
       };
   };
};