// FICHIER : PoutSimple1.cc
// CLASSE : PoutSimple1

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


#include "PoutSimple1.h"
#include "FrontPointF.h"
#include "FrontSegLine.h"
#include "DeformationP2D.h"


// ---- definition du constructeur de la classe conteneur de donnees communes ------------	
PoutSimple1::DonnCommunPS1::DonnCommunPS1 (GeomSeg& segL,GeomSeg& segH,DdlElement& tab,
	                       Met_Pout2D&  met_bie,Tableau <Vecteur *> & resEr,Mat_pleine& raidEr,
	                       Mat_pleine&  mat_masse,int nbi) :
	          segmentL(segL),segmentH(segH),tab_ddl(tab),met_pout(met_bie)
	          ,matGeom(tab.NbDdl(),tab.NbDdl())
	          ,matInit(tab.NbDdl(),tab.NbDdl())
	          ,d2_epsBB(nbi)
	          ,resErr(resEr),raidErr(raidEr)
	          ,matrice_masse(mat_masse)
	{ 
      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 (1);
       };
	 };
// Il ne doit exister  qu'un exemplaire de la classe, donc au niveau du constructeur de copie les tableaux
// de pointeurs pointes sur les mêmes entités que la valeur passée en paramètre	              
PoutSimple1::DonnCommunPS1::DonnCommunPS1(DonnCommunPS1& a) :
	  segmentL(a.segmentL),segmentH(a.segmentH),tab_ddl(a.tab_ddl),met_pout(a.met_pout)
	  ,matGeom(a.matGeom),matInit(a.matInit),d2_epsBB(a.d2_epsBB)
	  ,resErr(a.resErr),raidErr(a.raidErr)
	  ,matrice_masse(a.matrice_masse)
	          {};
PoutSimple1::DonnCommunPS1::DonnCommunPS1::~DonnCommunPS1()
      { 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);
       };
	             
// ---------- fin definition de la classe conteneur de donnees communes ------------	

//----------------------------------------------------------------
// def des donnees commune a tous les elements
// la taille n'est pas defini ici car elle depend de la lecture 
//----------------------------------------------------------------
int PoutSimple1::nbNoeud = 3; // nombre de noeud indiqué en dur
int PoutSimple1::nbintL = 2;   // nombre de point d'intégration indiqué en dur sur l'axe
int PoutSimple1::nbintH = 2;   // nombre de point d'intégration indiqué en dur dans 
                              //l'épaisseur
Tableau <Vecteur>  PoutSimple1::tabD2phi; // par commodité
                              
Tableau <TenseurBB *> PoutSimple1::d_epsBB; // place pour la variation des def
Tableau <TenseurHH *> PoutSimple1::d_sigHH; // place pour la variation des sigma
PoutSimple1::DonnCommunPS1 * PoutSimple1::doCoPS1 = NULL;                                     
int PoutSimple1::CalculResidu_t_PoutSimple1_met_abstraite = 0;  
int PoutSimple1::CalculResidu_tdt_PoutSimple1_met_abstraite = 0;  
int  PoutSimple1::Calcul_implicit_PoutSimple1_met_abstraite = 0;
int  PoutSimple1::Calcul_VarDualSort = 0;
int  PoutSimple1::CalDynamique = 0;
int  PoutSimple1::CalPt_0_t_tdt = 0 ; // pour le calcul de point à 0 t et tdt

PoutSimple1::ConstrucElementPoutSimple1 PoutSimple1::construcElementPoutSimple1;
// fonction privee
// dans cette fonction il ne doit y avoir que les données communes !!
void PoutSimple1::Def_DonnCommunPS1() 
 {  GeomSeg segmentL(nbintL,nbNoeud) ; // element geometrique correspondant de l'axe
    GeomSeg segmentH(nbintH) ; // element geometrique correspondant à l'épaisseur
    // par défaut le nombre de noeuds est 2
    // degre de liberte
    int dim = ParaGlob::Dimension();
    DdlElement  tab_ddl(nbNoeud,dim);
    int posi = Id_nom_ddl("X1") -1;
    for (int i =1; i<= ParaGlob::Dimension(); i++)
      for (int j =1; j<= nbNoeud; j++)
  	     tab_ddl.Change_Enum(j,i,Enum_ddl(i+posi));
    // cas des ddl éléments secondaires pour le calcul d'erreur
    // les tenseurs de contrainte sont de dimension 1 a priori
    // donc 1 ddl sig11
    DdlElement  tab_ddlErr(nbNoeud,1);
    posi = Id_nom_ddl("SIG11") -1;
    for (int j=1; j<= nbNoeud; j++)
	    tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi));
    // 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 = 1 ici,
	   // , tableau de ddl et la def de variables
	   Met_Pout2D  metri(ParaGlob::Dimension(),1,tab_ddl,tab,nbNoeud) ;
		  // ---- cas du calcul d'erreur sur sigma ou epsilon
    Tableau <Vecteur *>  resEr(nbNoeud);
    for (int i=1; i<= nbNoeud; i++)
        resEr(i)=new Vecteur (1); // tenseur à une dimension
    Mat_pleine  raidEr(nbNoeud,nbNoeud); // la raideur pour l'erreur
       // cas de la dynamique
    Mat_pleine  matmasse(1,nbNoeud);  // a priori on dimensionne en diagonale 

    // definition de la classe static contenant toute les variables communes aux PoutSimple1s
    doCoPS1 = new DonnCommunPS1(segmentL,segmentH,tab_ddl,metri,resEr,raidEr,matmasse,nbintL*nbintH);
    //dimensionnement des variations des deformations 
	   int nbddl = tab_ddl.NbDdl();
    d_epsBB.Change_taille(nbddl);
    d_sigHH.Change_taille(nbddl);
    for (int i=1; i<= nbddl; i++)
      { d_epsBB(i) = NevezTenseurBB (1);
        d_sigHH(i) = NevezTenseurHH (1);
      }
  };    
PoutSimple1::PoutSimple1 ()  :
// Constructeur par defaut
PiPoCo(),lesPtMecaInt(),donnee_specif()

  { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
    int dim = ParaGlob::Dimension();
    if (dim != 2) 
      { if (ParaGlob::NiveauImpression() >= 7)
          {cout << "\n attention !!, l'element PoutSimple1 n'est pas utilisable car \n" 
                "seul la dimension 2 est acceptable pour cet element \n";
           }     
        return;
       }
    id_interpol=CUBIQUE;  //   donnees de la classe mere
	   id_geom=PS1;      //
    // stockage des donnees particulieres de la loi de comportement mécanique au point d'integ
    tabSaveDon.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la loi de comportement thermo mécanique au point d'integ
    tabSaveTP.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la déformation mécanique au point d'integ
    tabSaveDefDon.Change_taille(nbintL*nbintH);
    tab_energ.Change_taille(nbintL*nbintH);
    tab_energ_t.Change_taille(nbintL*nbintH);
    tab_noeud.Change_taille(nbNoeud);
    // le fait de mettre les pointeurs a null permet
    // de savoir que l'element n'est pas complet
    for (int i=1; i<= nbNoeud; i++)
      tab_noeud(i) = NULL;
    // definition des donnees communes aux PoutSimple1s
    // a la premiere definition d'une PoutSimple1
    if (doCoPS1 == NULL) Def_DonnCommunPS1();
    met = &(doCoPS1->met_pout); // met est defini dans elemeca
        // --- cas de la dynamique -----
    mat_masse =  &(doCoPS1->matrice_masse);
    // pour les dérivées secondes, transformation de tableau en vecteur
    tabD2phi.Change_taille((doCoPS1->segmentL).taD2phi().Taille());
    for (int i=1;i<= (doCoPS1->segmentL).taD2phi().Taille(); i++)
       tabD2phi(i) = ((doCoPS1->segmentL).taD2phi())(i).Ligne(1);
    // def pointe sur la deformation specifique a l'element
    def = new DeformationP2D(*met,tab_noeud
                    ,(doCoPS1->segmentH).TaDphi(),(doCoPS1->segmentH).TaPhi()
                    ,(doCoPS1->segmentL).TaDphi(),(doCoPS1->segmentL).TaPhi()
                    ,tabD2phi);
    // idem pour la masse ici par simplicité
    defMas = new DeformationP2D(*met,tab_noeud
                    ,(doCoPS1->segmentH).TaDphi(),(doCoPS1->segmentH).TaPhi()
                    ,(doCoPS1->segmentL).TaDphi(),(doCoPS1->segmentL).TaPhi()
                    ,tabD2phi);

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

};  
	
PoutSimple1::PoutSimple1 (double epais,double larg,int num_mail,int num_id):
// Constructeur utile si la hauteur, la largeur de l'element et
// le numero de l'element  sont connus
PiPoCo(num_mail,num_id,CUBIQUE,PS1),lesPtMecaInt()
,donnee_specif(epais,larg)
  { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
    int dim = ParaGlob::Dimension();
    if (dim != 2) 
      { if (ParaGlob::NiveauImpression() >= 7)
          {cout << "\n attention !!, l'''\'element PoutSimple1 n''' est pas utilisable car \n" 
                "seul la dimension 2 est acceptable pour cet élément \n";
           }     
        return;
       }
    // stockage des donnees particulieres de la loi de comportement mécanique au point d'integ
    tabSaveDon.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la loi de comportement thermo mécanique au point d'integ
    tabSaveTP.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la déformation mécanique au point d'integ
    tabSaveDefDon.Change_taille(nbintL*nbintH);
    tab_energ.Change_taille(nbintL*nbintH);
    tab_energ_t.Change_taille(nbintL*nbintH);
    tab_noeud.Change_taille(nbNoeud);
    // le fait de mettre les pointeurs a null permet
    // de savoir que l'element n'est pas complet
    for (int i=1; i<= nbNoeud; i++)
      tab_noeud(i) = NULL;
    // definition des donnees communes aux PoutSimple1s
    // a la premiere definition d'une PoutSimple1
    if (doCoPS1 == NULL) Def_DonnCommunPS1();
    met = &(doCoPS1->met_pout);// met est defini dans elemeca
        // --- cas de la dynamique -----
    mat_masse =  &(doCoPS1->matrice_masse);
    // pour les dérivées secondes, transformation de tableau en vecteur
    tabD2phi.Change_taille((doCoPS1->segmentL).taD2phi().Taille());
    for (int i=1;i<= (doCoPS1->segmentL).taD2phi().Taille(); i++)
       tabD2phi(i) = ((doCoPS1->segmentL).taD2phi())(i).Ligne(1);
    // def pointe sur la deformation specifique a l'element
    def = new DeformationP2D(*met,tab_noeud
                    ,(doCoPS1->segmentH).TaDphi(),(doCoPS1->segmentH).TaPhi()
                    ,(doCoPS1->segmentL).TaDphi(),(doCoPS1->segmentL).TaPhi()
                    ,tabD2phi);

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

};

// Constructeur fonction  d'un numero de maillage et d'identification  
PoutSimple1::PoutSimple1 (int num_mail,int num_id) :
  PiPoCo(num_mail,num_id,CUBIQUE,PS1),lesPtMecaInt(),donnee_specif()
  { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
    int dim = ParaGlob::Dimension();
    if (dim != 2) 
      { if (ParaGlob::NiveauImpression() >= 7)
          {cout << "\n attention !!, l'''\'element PoutSimple1 n''' est pas utilisable car \n" 
                "seul la dimension 2 est acceptable pour cet élément \n";
           }     
        return;
       }
    // stockage des donnees particulieres de la loi de comportement mécanique au point d'integ
    tabSaveDon.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la loi de comportement thermo mécanique au point d'integ
    tabSaveTP.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la déformation mécanique au point d'integ
    tabSaveDefDon.Change_taille(nbintL*nbintH);
    tab_energ.Change_taille(nbintL*nbintH);
    tab_energ_t.Change_taille(nbintL*nbintH);
    tab_noeud.Change_taille(nbNoeud);
    // le fait de mettre les pointeurs a null permet
    // de savoir que l'element n'est pas complet
    for (int i=1; i<= nbNoeud; i++)
      tab_noeud(i) = NULL;
    // definition des donnees communes aux PoutSimple1s
    // a la premiere definition d'une PoutSimple1
    if (doCoPS1 == NULL) Def_DonnCommunPS1();
    met = &(doCoPS1->met_pout);// met est defini dans elemeca
        // --- cas de la dynamique -----
    mat_masse =  &(doCoPS1->matrice_masse);
    // pour les dérivées secondes, transformation de tableau en vecteur
    tabD2phi.Change_taille((doCoPS1->segmentL).taD2phi().Taille());
    for (int i=1;i<= (doCoPS1->segmentL).taD2phi().Taille(); i++)
       tabD2phi(i) = ((doCoPS1->segmentL).taD2phi())(i).Ligne(1);
    // def pointe sur la deformation specifique a l'element
    def = new DeformationP2D(*met,tab_noeud
                    ,(doCoPS1->segmentH).TaDphi(),(doCoPS1->segmentH).TaPhi()
                    ,(doCoPS1->segmentL).TaDphi(),(doCoPS1->segmentL).TaPhi()
                    ,tabD2phi);

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

 };

PoutSimple1::PoutSimple1 (double epais,double larg,int num_mail,int num_id,const Tableau<Noeud *>& tab):
// Constructeur utile si la hauteur, la largeur de l'element, le numero de l'element et
// le tableau des noeuds sont connus
PiPoCo(num_mail,num_id,tab,CUBIQUE,PS1),lesPtMecaInt()
,donnee_specif(epais,larg)
  { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
    int dim = ParaGlob::Dimension();
    if (dim != 2) 
      { if (ParaGlob::NiveauImpression() >= 7)
          {cout << "\n attention !!, l'''\'element PoutSimple1 n''' est pas utilisable car \n" 
                "seul la dimension 2 est acceptable pour cet élément \n";
           }     
        return;
       }
    // stockage des donnees particulieres de la loi de comportement mécanique au point d'integ
    tabSaveDon.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la loi de comportement thermo mécanique au point d'integ
    tabSaveTP.Change_taille(nbintL*nbintH);
    // stockage des donnees particulieres de la déformation mécanique au point d'integ
    tabSaveDefDon.Change_taille(nbintL*nbintH);
    tab_energ.Change_taille(nbintL*nbintH);
    tab_energ_t.Change_taille(nbintL*nbintH);
    if (tab_noeud.Taille() != nbNoeud)
      { cout << "\n erreur de dimensionnement du tableau de noeud \n";
       cout << " PoutSimple1::PoutSimple1 (double sect,int num_id,const Tableau<Noeud *>& tab)\n";
       Sortie (1); }
    // definition des donnees communes aux PoutSimple1s
    // a la premiere definition d'une PoutSimple1
    if (doCoPS1 == NULL) Def_DonnCommunPS1();
    met = &(doCoPS1->met_pout);// met est defini dans elemeca
        // --- cas de la dynamique -----
    mat_masse =  &(doCoPS1->matrice_masse);
    // pour les dérivées secondes, transformation de tableau en vecteur
    tabD2phi.Change_taille((doCoPS1->segmentL).taD2phi().Taille());
    for (int i=1;i<= (doCoPS1->segmentL).taD2phi().Taille(); i++)
       tabD2phi(i) = ((doCoPS1->segmentL).taD2phi())(i).Ligne(1);
    // def pointe sur la deformation specifique a l'element
    def = new DeformationP2D(*met,tab_noeud
                    ,(doCoPS1->segmentH).TaDphi(),(doCoPS1->segmentH).TaPhi()
                    ,(doCoPS1->segmentL).TaDphi(),(doCoPS1->segmentL).TaPhi()
                    ,tabD2phi);

    //dimensionnement des deformations et contraintes etc..
    int dimtens = 1;
    lesPtMecaInt.Change_taille_PtIntegMeca(nbintL*nbintH,dimtens);
    // attribution des numéros de référencement dans le conteneur
    for (int ni = 1; ni<= nbintL*nbintH; ni++)
       {lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage);
        lesPtMecaInt(ni).Change_Nb_ele(this->num_elt);
        lesPtMecaInt(ni).Change_Nb_pti(ni);
       };
    // construction du tableau de ddl des noeuds de PoutSimple1
    ConstTabDdl(); 
};

PoutSimple1::PoutSimple1 (const PoutSimple1& pout) :
 PiPoCo (pout),lesPtMecaInt(pout.lesPtMecaInt),donnee_specif(pout.donnee_specif)
 
// Constructeur de copie
  { lesPtIntegMecaInterne = &lesPtMecaInt; //  association avec le pointeur d'ElemMeca  
    // stockage des donnees particulieres de la loi de comportement au point d'integ
    //   tabSaveDon.Change_taille(1);
    *(residu) = *(pout.residu);
    *(raideur) = *(pout.raideur);
  };

PoutSimple1::~PoutSimple1 ()
// Destructeur
{  // dans le cas ou on est en dimension 3, il n'y a rien de créé
   // donc pas de delete
   if (ParaGlob::Dimension() != 2) return;
   // cas normal
   LibereTenseur();
 };

// Lecture des donnees de la classe sur fichier
void 
PoutSimple1::LectureDonneesParticulieres
    (UtilLecture * entreePrinc,Tableau<Noeud  *> * tabMaillageNoeud)
  { int nb;
    tab_noeud.Change_taille(nbNoeud);
    for (int i=1; i<= nbNoeud; 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
       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 != nbNoeud) entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
         }
     #endif  
       else // cas d'une erreur de lecture
        { cout << "\n erreur de lecture inconnue  ";
          entreePrinc->MessageBuffer("** lecture des données particulières **");
          cout << "PoutSimple1::LectureDonneesParticulieres"; 
          Affiche();
          Sortie (1);
        }
     }
    // construction du tableau de ddl des noeuds de PoutSimple1
    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 & PoutSimple1::Point_physique(const Coordonnee& ,Coordonnee & co ,Enum_dure )
 { /*// a) on commence par définir les bonnes grandeurs dans la métrique
   if( PoutSimple1::CalPt_0_t_tdt == 0)
      { PoutSimple1::CalPt_0_t_tdt = 1;
        Tableau<Enum_variable_metrique> tab(3);
        tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt;
        doCoPS1->met_pout.PlusInitVariables(tab) ;
       };  
   // b) calcul de l'interpolation 
   const Vecteur& phix = doCoPS1->segmentL->Phi(c_int(1)); // en x
   const Vecteur& phiy = doCoPS1->segmentH->Phi(c_int(2)); // en hauteur
   // c) calcul du point
   switch (temps)
    { case TEMPS_0 : co = met->PointM_0(tab_noeud,phi); break;
      case TEMPS_t : co = met->PointM_t(tab_noeud,phi); break;
      case TEMPS_tdt : co = met->PointM_tdt(tab_noeud,phi); break;
     }       
   // d) retour*/
   cout << "\n methode non implantee pour l'instant !"
        << "\n Coordonnee & PoutSimple1::Point_physique(...";
   Sortie(1);     
   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 PoutSimple1::Point_physique(const Coordonnee& ,Tableau <Coordonnee> & )
 { /*t_co.Change_taille(3);
   // a) on commence par définir les bonnes grandeurs dans la métrique
   if( PoutSimple1::CalPt_0_t_tdt == 0)
      { PoutSimple1::CalPt_0_t_tdt = 1;
        Tableau<Enum_variable_metrique> tab(3);
        tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt;
        doCoPS1->met_pout.PlusInitVariables(tab) ;
       };  
   // b) calcul de l'interpolation 
   const Vecteur& phi = doCoPS1->hexaed->Phi(c_int);
   // c) calcul des point
   switch (t_co.Taille())
    { case 3 : t_co(3) = met->PointM_tdt(tab_noeud,phi); 
      case 2 : t_co(2) = met->PointM_t(tab_noeud,phi); 
      case 1 : t_co(1) = met->PointM_0(tab_noeud,phi); 
     } */
   cout << "\n methode non implantee pour l'instant !"
        << "\n void PoutSimple1::Point_physique(...";
   Sortie(1);     
  };
  
// Calcul du residu local à t ou tdt en fonction du booleen atdt
 Vecteur* PoutSimple1::CalculResidu (bool atdt,const ParaAlgoControle & pa)
     { // dimensionnement de la metrique
      if (!atdt)
       {if( CalculResidu_t_PoutSimple1_met_abstraite == 0)
        { CalculResidu_t_PoutSimple1_met_abstraite = 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;
          doCoPS1->met_pout.PlusInitVariables(tab) ;
         };}  
      else 
       {if( CalculResidu_tdt_PoutSimple1_met_abstraite == 0)
        {CalculResidu_tdt_PoutSimple1_met_abstraite = 1;
         Tableau<Enum_variable_metrique> tab(10);
         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) = igradVBB_tdt;
         doCoPS1->met_pout.PlusInitVariables(tab) ;
         };};  
       // dimensionnement du residu
       int nbddl = doCoPS1->tab_ddl.NbDdl();
       if ( residu == NULL)
          residu = new Vecteur(nbddl); // cas du premier passage
       else
          for (int i =1; i<= nbddl; i++)  // cas des autres passages
             (*residu)(i) = 0.;             // on initialise a zero
       Vecteur poidsL =(doCoPS1->segmentL).TaWi();   // poids d'intergration sur l'axe      
       Vecteur poidsH =(doCoPS1->segmentH).TaWi();   // poids d'intergration dans l'épaisseur      
       poidsL *=  donnee_specif.largeur; // en fait il n'y a pas d'intégration suivant la largeur
       poidsH *=  0.5 * donnee_specif.epaisseur; //  car le domaine d'intégration dans l'épaisseur
       // suivant l'élément de référence est 2
       PiPoCo::Cal_explicitPiPoCo (doCoPS1->tab_ddl,d_epsBB,nbintL,poidsL,nbintH,poidsH,pa,atdt);
       return residu;  
     };

// Calcul du residu local et de la raideur locale,
//  pour le schema implicite
Element::ResRaid  PoutSimple1::Calcul_implicit (const ParaAlgoControle & pa)
    {  
      bool cald_Dvirtuelle = false;
      if( Calcul_implicit_PoutSimple1_met_abstraite == 0)
        { 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; 
         doCoPS1->met_pout.PlusInitVariables(tab) ;
         Calcul_implicit_PoutSimple1_met_abstraite = 1;
         // 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; 
        };

      // dimensionnement du residu
       int nbddl = doCoPS1->tab_ddl.NbDdl();
       if ( residu == NULL)
          residu = new Vecteur(nbddl); // cas du premier passage
       else
          for (int i =1; i<= nbddl; i++)  // cas des autres passages
             (*residu)(i) = 0.;             // on initialise a zero
      // dimensionnement de la raideur
       if ( raideur == NULL)
          raideur = new Mat_pleine(nbddl,nbddl); // cas du premier passage
       else
          for (int i =1; i<= nbddl; i++)  // cas des autres passages
            for (int j=1; j<= nbddl; j++) //
              (*raideur)(i,j) = 0.;       // on initialise a zero
       Vecteur poidsL =(doCoPS1->segmentL).TaWi();   // poids d'intergration sur l'axe      
       Vecteur poidsH =(doCoPS1->segmentH).TaWi();   // poids d'intergration dans l'épaisseur      
       poidsL *=  donnee_specif.largeur; // en fait il n'y a pas d'intégration suivant la largeur
       poidsH *=   0.5 * donnee_specif.epaisseur; // 0.5 car le domaine d'intégration dans l'épaisseur
       // suivant l'élément de référence est 2
       PiPoCo::Cal_implicitPiPoCo 
                 (doCoPS1->tab_ddl,d_epsBB,doCoPS1->d2_epsBB,d_sigHH
                  ,nbintL,poidsL,nbintH,poidsH,pa,cald_Dvirtuelle);
       Element::ResRaid el;
       el.res = residu; 
       el.raid = raideur;
       return el;
    };
		
// Calcul de la matrice masse pour l'élément
Mat_pleine * PoutSimple1::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse) 
 {    // dimensionement de la métrique si nécessaire
      if (CalDynamique == 0)
       { 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; 
         doCoPS1->met_pout.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 = doCoPS1->tab_ddl.NbDdl();
            (doCoPS1->matrice_masse).Initialise (nbddl,nbddl,0.);
           } 
        };
       Vecteur poidsL =(doCoPS1->segmentL).TaWi();   // poids d'intergration sur l'axe
       // prise en compte de l'épaisseur et de la largeur
       poidsL *= (donnee_specif.largeur * donnee_specif.epaisseur);      
       // appel de la routine générale 
       ElemMeca::Cal_Mat_masse (doCoPS1->tab_ddl,type_calcul_masse,
                            nbintL,(doCoPS1->segmentL).TaPhi(),nbNoeud
                            ,poidsL);
       return mat_masse;                     
  };
		
// Calcul de la matrice géométrique et initiale
ElemMeca::MatGeomInit PoutSimple1::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & pa) 
 {  
     bool cald_Dvirtuelle = false;
    if( Calcul_implicit_PoutSimple1_met_abstraite == 0)
        { 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; 
         doCoPS1->met_pout.PlusInitVariables(tab) ;
         Calcul_implicit_PoutSimple1_met_abstraite = 1;
         // 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 = doCoPS1->matGeom;    
    Mat_pleine & matInit = doCoPS1->matInit;    
    // mise à zéro de la matrice géométrique
    matGeom.Initialise();
    Vecteur poidsL =(doCoPS1->segmentL).TaWi();   // poids d'intergration sur l'axe      
    Vecteur poidsH =(doCoPS1->segmentH).TaWi();   // poids d'intergration dans l'épaisseur      
    poidsL *=  donnee_specif.largeur; // en fait il n'y a pas d'intégration suivant la largeur
    poidsH *=   0.5 * donnee_specif.epaisseur; // 0.5 car le domaine d'intégration dans l'épaisseur
    // suivant l'élément de référence est 2
    PiPoCo::Cal_matGeom_InitPiPoCo (matGeom,matInit,doCoPS1->tab_ddl, d_epsBB,
                doCoPS1->d2_epsBB,d_sigHH,nbintL,poidsL,nbintH
                ,poidsH,pa,cald_Dvirtuelle);
    return  MatGeomInit(&matGeom,&matInit);     
  } ;
     
// actualisation des ddl et des grandeurs actives de t+dt vers t      
void PoutSimple1::TdtversT()
 { lesPtMecaInt.TdtversT(); // contrainte 
   int taille = lesPtMecaInt.NbPti();
   for (int ni=1;ni<= taille; ni++)
     { 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
  };
// actualisation des ddl et des grandeurs actives de t vers tdt      
void PoutSimple1::TversTdt()
 { lesPtMecaInt.TversTdt(); // contrainte
   int taille = lesPtMecaInt.NbPti();
   for (int ni=1;ni<= taille; ni++)
    { 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
  };

//============= 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 PoutSimple1::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 poutre
 switch (cas)
  { case 1 : // ------- on récupère tout -------------------------
     { // construction du tableau de ddl des noeuds du triangle
       ConstTabDdl(); 
       // récup contraintes et déformation
       lesPtMecaInt.Lecture_base_info(ent,cas);
	   break;
	 }
    case  2 : // ----------- lecture uniquement de se qui varie --------------------
     { // récup contraintes et déformation
       lesPtMecaInt.Lecture_base_info(ent,cas);
	   break;
	 }
    default :
	 { cout << "\nErreur : valeur incorrecte du type de lecture !\n";
	   cout << "Lecture_base_info(ifstream& ent,const Tableau<Noeud "
	        << " *> * tabMaillageNoeud,const 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 PoutSimple1::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 poutre
 switch (cas)
  { case 1 : // ------- on sauvegarde tout -------------------------
     { 
       // des tenseurs déformation et contrainte,
       lesPtMecaInt.Ecriture_base_info(sort,cas);
	   break;
	 }
    case  2 : // ----------- sauvegarde uniquement de se qui varie --------------------
     { // des tenseurs déformation et contrainte, 
       lesPtMecaInt.Ecriture_base_info(sort,cas);
	   break;
	 }
    default :
	 { cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
	   cout << "PoutSimple1::Ecriture_base_info(ofstream& sort,const int cas)"
			<< " cas= " << cas << endl;
	   Sortie(1);
	  }
   }	   
  };
 
// retourne les tableaux de ddl associés aux noeuds, gere par l'element
// ce tableau et specifique a l'element
const DdlElement & PoutSimple1::TableauDdl() const
     { return doCoPS1->tab_ddl; };
// liberation de la place pointee     
void PoutSimple1::Libere ()
  {Element::Libere (); // liberation de residu et raideur
   LibereTenseur() ; // liberation des tenseur intermediaires
  };
// acquisition ou modification d'une loi de comportement
void PoutSimple1::DefLoi (LoiAbstraiteGeneral * NouvelleLoi)
    { // verification du type de loi
      if ((NouvelleLoi->Dimension_loi() != 1) && (NouvelleLoi->Dimension_loi() != 4))
        { cout << "\n Erreur, la loi de comportement a utiliser avec des PoutSimple1s";
          cout << " doit etre de type 1D, \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;
        // initialisation du stockage particulier
        for (int i=1;i<=nbintL*nbintH;i++) tabSaveDon(i) = loiComp->New_et_Initialise();  
        // idem pour le type de déformation mécanique associé            
        int iDefmax = tabSaveDefDon.Taille();
        for (int i=1;i<= iDefmax;i++) 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;
        // initialisation du stockage particulier thermo physique, 
        int imax = tabSaveTP.Taille();
        for (int i=1;i<=nbintL*nbintH;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 PoutSimple1::TestComplet()
  { int res = ElemMeca::TestComplet(); // test dans la fonction mere
    if (( donnee_specif.epaisseur == epaisseur_defaut) 
          || (donnee_specif.largeur == largeur_defaut))
     { cout << "\n l\'epaisseur ou la largeur de la PoutSimple1 ne sont pas defini  \n";
       res = 0; }
    if ( tab_noeud(1) == NULL)
     { cout << "\n les noeuds de la PoutSimple1 ne sont pas defini  \n";
       res = 0; }
    else 
     { int testi =1;
       int posi = Id_nom_ddl("X1") -1;    
       for (int i =1; i<= ParaGlob::Dimension(); i++)
          for (int j=1;j<=2;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 PoutSimple1 ne sont pas defini  \n";
           cout << " \n utilisez PoutSimple1::ConstTabDdl() pour completer " ;
           res = 0; }            
      }          
    return res;   
  };
    
// ajout du tableau de ddl des noeuds de PoutSimple1
 void PoutSimple1::ConstTabDdl()
    {  
    Tableau <Ddl> ta(ParaGlob::Dimension());
    int posi = Id_nom_ddl("X1") -1;    
    for (int i =1; i<= ParaGlob::Dimension(); i++)
	    {Ddl inter((Enum_ddl(i+posi)),0.,LIBRE);
	     ta(i) = inter;
	     }
    // attribution des ddls aux noeuds  
    for (int i1=1;i1<=nbNoeud;i1++)  
     tab_noeud(i1)->PlusTabDdl(ta); 
    };
		
// procesure permettant de completer l'element apres
// sa creation avec les donnees du bloc transmis
// peut etre appeler plusieurs fois
 Element* PoutSimple1::Complete(BlocGen & bloc,LesFonctions_nD*  lesFonctionsnD)
    { // complétion avec bloc
      if (bloc.Nom(1) == "largeurs")
       { donnee_specif.largeur = bloc.Val(1);
         return this;
        }
       else if (bloc.Nom(1) == "epaisseurs")
       { donnee_specif.epaisseur = bloc.Val(1);
         return this;
        }
       else  
         return ElemMeca::Complete_ElemMeca(bloc,lesFonctionsnD);
    };       
                         
// affichage dans la sortie transmise, des variables duales "nom"
// dans le cas ou nom est vide, affichage de "toute" les variables
 void PoutSimple1::AfficheVarDual(ofstream& sort, Tableau<string>& nom)
     {// affichage de l'entête de l'element
      sort << "\n******************************************************************";
      sort << "\n Element PoutSimple1 (" << nbNoeud << " noeuds " 
           << nbintL<<"*"<<nbintH << " points d'integration L*H) ";
      sort << "\n******************************************************************";
      // appel de la procedure de elem meca
      if (!(Calcul_VarDualSort) && (Calcul_implicit_PoutSimple1_met_abstraite))
        { VarDualSort(sort,nom,1,1);
          Calcul_VarDualSort += 1;
         } 
      else if ((Calcul_VarDualSort) && (Calcul_implicit_PoutSimple1_met_abstraite))       
         VarDualSort(sort,nom,1,11);
      else if (!(Calcul_VarDualSort) && 
               ((CalculResidu_t_PoutSimple1_met_abstraite)||
               (CalculResidu_tdt_PoutSimple1_met_abstraite)))       
        { VarDualSort(sort,nom,1,2);
          Calcul_VarDualSort += 1;
         }         
      else if ((Calcul_VarDualSort) && 
               ((CalculResidu_t_PoutSimple1_met_abstraite)||
               (CalculResidu_tdt_PoutSimple1_met_abstraite)))       
         VarDualSort(sort,nom,1,12);
      // sinon on ne fait rien     
      };
      
// 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 & PoutSimple1::Frontiere(bool force)
  { int cas = 6; // on veut des lignes et des points
    return Frontiere_elemeca(cas,force);
  };