//#include "Debug.h"


// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

#include "GeomHexaQuad.h"
#include <math.h>
#include "GeomSeg.h"
#include "GeomQuadrangle.h"
#include "MathUtil.h"
#include "GeomHexalin.h"
#include "GeomHexaQuadComp.h"

// constructeur
//  la dimension est 3, on a 8 pt d'integration par défaut, 
// 20 noeuds et 6 faces, 12 aretes
GeomHexaQuad::GeomHexaQuad(int nbi) :
    GeomHexaCom(nbi,20,QUADRATIQUE),phi_M(),dphi_M()
   {  // coordonnees dans l'élément de référence des noeuds 
      ptelem(1) = Coordonnee(-1.,-1.,-1.); ptelem(2) = Coordonnee(1.,-1.,-1.); 
      ptelem(3) = Coordonnee(1.,1.,-1.); ptelem(4) = Coordonnee(-1.,1.,-1.);
      ptelem(5) = Coordonnee(-1.,-1.,1.); ptelem(6) = Coordonnee(1.,-1.,1.);
      ptelem(7) = Coordonnee(1.,1.,1.); ptelem(8) = Coordonnee(-1.,1.,1.);
      ptelem(9) = Coordonnee(0.,-1.,-1.); ptelem(10) = Coordonnee(1.,0.,-1.);
      ptelem(11) = Coordonnee(0.,1.,-1.); ptelem(12) = Coordonnee(-1.,0.,-1.); 
      ptelem(13) = Coordonnee(-1.,-1.,0.); ptelem(14) = Coordonnee(1.,-1.,0.);
      ptelem(15) = Coordonnee(1.,1.,0.); ptelem(16) = Coordonnee(-1.,1.,0.);
      ptelem(17) = Coordonnee(0.,-1.,1.); ptelem(18) = Coordonnee(1.,0.,1.);
      ptelem(19) = Coordonnee(0.,1.,1.); ptelem(20) = Coordonnee(-1.,0.,1.);
      // définition de la numérotation locale de l'élément de direction inverse
      INVCONNEC(1) = 1;INVCONNEC(2) = 4; INVCONNEC(3) = 3;INVCONNEC(4) = 2;
      INVCONNEC(5) = 5;INVCONNEC(6) = 8; INVCONNEC(7) = 7;INVCONNEC(8) = 6;
      INVCONNEC(9) = 12;INVCONNEC(10) = 11; INVCONNEC(11) = 10;INVCONNEC(12) = 9;
      INVCONNEC(13) = 13;INVCONNEC(14) = 16; INVCONNEC(15) = 15;INVCONNEC(16) = 14;
      INVCONNEC(17) = 20;INVCONNEC(18) = 19; INVCONNEC(19) = 18;INVCONNEC(20) = 17;
      // le tableau des tranches
      IND.Change_taille(4);
      IND(1)=8; // les sommets
      IND(2)=4; // les 4 noeuds quadratiques de la face du dessous
      IND(3)=4;	 // les 4 noeuds quadratiques verticaux
      IND(4)=4;	 // les 4 noeuds quadratiques de la face du dessus
      //--------------------------------
      //def des arretes
      //--------------------------------
      int nbil =2;  // nb de pt d'integ par ligne
      int nbnel =3;   // nb de noeud du segment
      seg(1) = new GeomSeg(nbil,nbnel);
      for (int il=2;il<= NBSE; il++)  // ici NBSE = 12
           seg(il) = seg(1);
      // def des tableaux de connection des noeuds  des aretes
      for (int i =1;i<=NBSE;i++) NONS(i).Change_taille(nbnel);
      // la description est fait selon le fichier EIMail
      NONS(1)(1) = 1;NONS(1)(2) = 9;NONS(1)(3) = 2;
      NONS(2)(1) = 2;NONS(2)(2) = 10;NONS(2)(3) = 3;
      NONS(3)(1) = 3;NONS(3)(2) = 11;NONS(3)(3) = 4;
      NONS(4)(1) = 4;NONS(4)(2) = 12;NONS(4)(3) = 1;
      NONS(5)(1) = 1;NONS(5)(2) = 13;NONS(5)(3) = 5;
      NONS(6)(1) = 2;NONS(6)(2) = 14;NONS(6)(3) = 6;
      NONS(7)(1) = 3;NONS(7)(2) = 15;NONS(7)(3) = 7;
      NONS(8)(1) = 4;NONS(8)(2) = 16;NONS(8)(3) = 8;
      NONS(9)(1) = 5;NONS(9)(2) = 17;NONS(9)(3) = 6;
      NONS(10)(1) = 6;NONS(10)(2) = 18;NONS(10)(3) = 7;
      NONS(11)(1) = 7;NONS(11)(2) = 19;NONS(11)(3) = 8;
      NONS(12)(1) = 8;NONS(12)(2) = 20;NONS(12)(3) = 5;
      //--------------------------------
      //def des faces
      //--------------------------------
      int nbis =4;  // nb de pt d'integ par facee
      int nbnes =8;   // nb de noeud de la face
      face(1) = new GeomQuadrangle(nbis,nbnes);
      for (int is=2;is<= NBFE; is++)  // ici NBFE = 6
           face(is) = face(1);
      // def des tableaux de connection des noeuds  des faces
      for (int i =1;i<=NBFE;i++) NONF(i).Change_taille(nbnes); 
      // connection entre les noeuds des faces et les noeuds des elements
      NONF(1)(1)= 1; NONF(1)(2)= 4; NONF(1)(3)= 3; NONF(1)(4)= 2;
      NONF(1)(5)= 12; NONF(1)(6)= 11; NONF(1)(7)= 10; NONF(1)(8)= 9;

      NONF(2)(1)= 1;NONF(2)(2)= 5; NONF(2)(3)= 8;NONF(2)(4)= 4;
      NONF(2)(5)= 13;NONF(2)(6)= 20; NONF(2)(7)= 16;NONF(2)(8)= 12;

      NONF(3)(1)= 1; NONF(3)(2)= 2;NONF(3)(3)= 6;NONF(3)(4)= 5;
      NONF(3)(5)= 9; NONF(3)(6)= 14;NONF(3)(7)= 17;NONF(3)(8)= 13;

      NONF(4)(1)= 5; NONF(4)(2)= 6; NONF(4)(3)= 7; NONF(4)(4)= 8;
      NONF(4)(5)= 17; NONF(4)(6)= 18; NONF(4)(7)= 19; NONF(4)(8)= 20;

      NONF(5)(1)= 2; NONF(5)(2)= 3;NONF(5)(3)= 7;NONF(5)(4)= 6;
      NONF(5)(5)= 10; NONF(5)(6)= 15;NONF(5)(7)= 18;NONF(5)(8)= 14;

      NONF(6)(1)= 3; NONF(6)(2)= 4; NONF(6)(3)= 8; NONF(6)(4)= 7;
      NONF(6)(5)= 11; NONF(6)(6)= 16; NONF(6)(7)= 19; NONF(6)(8)= 15;
      
      // triangulation des différentes  faces
      // on se sert d'une part de l'élément de référence de chaque face
      // puis de la connection les faces par rapport à celle de l'élément
      // ici c'est le même élément pour toutes les faces
       // 1) récup du tableau de l'élément de référence de la face
      const Tableau<Tableau<Tableau<int> > > & tabi = face(1)->Trian_lin();
      int nbtria = tabi(1).Taille(); // nombre de triangle par face 
      // on est obligé de boucler sur tous les indices et de faire
      // de l'adressage indirecte
      for (int isf=1;isf<= NBFE; isf++)  // boucle sur les faces
       { NONFt(isf).Change_taille(nbtria);
         for (int if1=1;if1<= nbtria; if1++)  // boucle sur les triangles de la face
          { NONFt(isf)(if1).Change_taille(3);
            for (int in1=1;in1<= 3; in1++)  // boucle sur les noeuds du triangle
              NONFt(isf)(if1)(in1) = NONF(isf)(tabi(1)(if1)(in1));
            }
         }   

      if (nbi == 8) // cas classique
        // on utilise des valeurs prédéfinies via deux methodes internes
        // fonctions d'interpolation globales
       { Phiphi();
        // derivees des fonctions d'interpolations
        DphiDphi();
        }
      else if ((nbi==1) || (nbi==27) || (nbi == 64))
        { // on utilise les méthodes internes pour calculer les fonctions
          // d'interpolation aux points d'intégrations
          for (int ptint=1;ptint<= Nbi(); ptint++)
            tabPhi(ptint) = Phi_point( ptInteg(ptint));
          for (int ptint=1;ptint<= Nbi(); ptint++)
            tabDPhi(ptint) = Dphi_point( ptInteg(ptint));
         }   
     
 // vérification suivant OK a priori
      // essai de calcul directe des fonctions d'interpolation
      for (int ptint=1;ptint<= Nbi(); ptint++)
       tabPhi(ptint) = Phi_point( ptInteg(ptint));
      for (int ptint=1;ptint<= Nbi(); ptint++)
       tabDPhi(ptint) = Dphi_point( ptInteg(ptint)); //*/
     // vérification des fonctions d'interpolation analytique et numériques
       for (int ptint=1;ptint<= Nbi(); ptint++)
       {
        Vecteur a = tabPhi(ptint);Vecteur b = Phi_point( ptInteg(ptint));
        for (int ne=1;ne<= NBNE;ne++)
         if (Dabs(a(ne) - b(ne)) >= 1.E-14)
           { 
             cout << (a(ne)) << "  " << (b(ne));
             cout << " erreur dans les points d'intégrations ";
             Sortie(1);
           }
        }   
      // vérification des dérivées des fonctions d'interpolation analytique 
      // et numériques
      for (int ptint=1;ptint<= Nbi(); ptint++)
       {
        Mat_pleine a = tabDPhi(ptint);Mat_pleine b = Dphi_point( ptInteg(ptint));
        for (int ne=1;ne<= NBNE;ne++)
         for (int ia =1; ia<= 3; ia++)
          if (Dabs(a(ia,ne) - b(ia,ne)) >= 1.E-14)
           { 
             cout << (a(ia,ne)) << "  " << (b(ia,ne));
             cout << " erreur dans les points d'intégrations ";
             Sortie(1);
           }
        }  // */
    // ---- constitution du tableau Extrapol ----- 
    Calcul_extrapol(nbi);      

};


// destructeur
GeomHexaQuad::~GeomHexaQuad() 
  { delete seg(1);
    delete face(1);
   };
// constructeur de copie
GeomHexaQuad::GeomHexaQuad(const GeomHexaQuad& a) :
  GeomHexaCom(a),phi_M(a.phi_M),dphi_M(a.dphi_M)
   { // la copie des parties pointées est à la charge de la classe spécifique
     // definition des  faces
     face(1) = new GeomQuadrangle(*((GeomQuadrangle*)(a.face(1))));
     // def des segments
     seg(1) = new GeomSeg(*((GeomSeg*)(a.seg(1)))) ;
     for (int il=2;il<= NBSE; il++)
           seg(il) = seg(1);
   };    
    
// création d'élément identiques : cette fonction est analogue à la fonction new
// elle y fait d'ailleurs appel. l'implantation est spécifique dans chaque classe
// dérivée
// pt est le pointeur qui est affecté par la fonction 
ElemGeomC0 * GeomHexaQuad::newElemGeomC0(ElemGeomC0 * pt) 
  { pt = new GeomHexaQuad(*this);
    return pt; 
   };

//--------- cas de coordonnees locales quelconques ---------------- 

    // retourne les fonctions d'interpolation au point M (en coordonnees locales)
const Vecteur& GeomHexaQuad::Phi_point(const Coordonnee& M)
  {
    #ifdef MISE_AU_POINT	 	   
     // verification de la dimension des coordonnees locales
     if (M.Dimension() != 3)
       { cout << "\n erreur la dimension des coordonnees locales :" << M.Dimension()
              <<"n\'est pas egale a 3 "
             << "\nGeomHexaQuad::Phi(Coordonnee& M)";
        Sortie(1);
      }
    #endif
//    Vecteur phi(NBNE); // tableau des fonctions d'interpolation
    // dimentionnement éventuelle du tableau des fonctions d'interpolation
    phi_M.Change_taille(NBNE); // si la taille est identique -> aucune action
    //------------------------------------------------------
    //        cas d'un Hexaedre trilquadratique incomplet
    //------------------------------------------------------
    if (NBNE == 20)
     {   // fonctions  pour les sommets
       double unsurhuit = 1./8.;
       for (int r=1; r<= 8; r++)
        { double x1 = M(1) * ptelem(r)(1);double  y1 = M(2) * ptelem(r)(2); 
          double z1 = M(3) * ptelem(r)(3);
          phi_M(r) = unsurhuit *(1. + x1)*(1. + y1)*(1. + z1)
                      *(-2.+ x1+y1+z1);
         }
       // pour les noeuds du plan xi=0
       double unsurquatre = 1./4.;
       Tableau<int> ind(4);
       ind(1) = 9; ind(2) = 11; ind(3) = 19; ind(4) = 17;
       for (int r=1;r<= 4; r++)
        {  double y1 = M(2) * ptelem(ind(r))(2);double  z1 = M(3) * ptelem(ind(r))(3);
          phi_M(ind(r)) = unsurquatre *(1. - M(1)*M(1))*(1. + y1)*(1. + z1);
         }
       // pour les noeuds du plan eta=0
       ind(1) = 10; ind(2) = 18; ind(3) = 20; ind(4) = 12;
       for (int r=1;r<= 4; r++)
        { double x1 = M(1) * ptelem(ind(r))(1);double  z1 = M(3) * ptelem(ind(r))(3);
          phi_M(ind(r)) = unsurquatre *(1. + x1)*(1. - M(2)*M(2))*(1. + z1);
         }
       // pour les noeuds du plan zeta=0
       ind(1) = 13; ind(2) = 14; ind(3) = 15; ind(4) = 16;
       for (int r=1;r<= 4; r++)
        { double  x1 = M(1) * ptelem(ind(r))(1);double  y1 = M(2) * ptelem(ind(r))(2);
          phi_M(ind(r)) = unsurquatre *(1. + x1)*(1. + y1)*(1. - M(3)*M(3));
         }
       }  
      else
       {cout << "\n erreur l'hexaedre de nombre de noeud NBNE = " << NBNE 
             << "\n n\'est pas implante !! ";
        cout << "\nGeomHexaQuad::Phi(Coordonnee& M) " << endl;     
        Sortie(1);
        }
    // retour de phi_M
    return phi_M;
   };
    // retourne les derivees des fonctions d'interpolation au point M (en coordonnees locales)
const Mat_pleine& GeomHexaQuad::Dphi_point(const Coordonnee& M)
  {  
   #ifdef MISE_AU_POINT	 	   
     // verification de la dimension des coordonnees locales
     if (M.Dimension() != 3)
       { cout << "\n erreur la dimension des coordonnees locales :" << M.Dimension()
              <<"n\'est pas egale a 3 "
             << "\nGeomHexaQuad::Dphi(Coordonnee& M)";
        Sortie(1);
      }
    #endif
//   Mat_pleine dphi(3,NBNE); // le tableau des derivees
   // le tableau des derivees: redimentionnement si nécessaire
   if ((dphi_M.Nb_ligne() != 3)&&(dphi_M.Nb_colonne() != NBNE))
      dphi_M.Initialise (3,NBNE,0.);
    //------------------------------------------------------
    //        cas d'un Hexaedre triquadratique incomplet
    //------------------------------------------------------
    if (NBNE == 20)
     { // fonctions et dérivées pour les sommets
       double unsurhuit = 1./8.;
       for (int r=1; r<= 8; r++)
        {double  x1 = M(1) * ptelem(r)(1);double  y1 = M(2) * ptelem(r)(2);
         double  z1 = M(3) * ptelem(r)(3);
          dphi_M(1,r)= unsurhuit *ptelem(r)(1)*(1. + y1)*(1. + z1)
                      *(-1+2.*x1+y1+z1);
          dphi_M(2,r)= unsurhuit *(1. + x1)*ptelem(r)(2)*(1. + z1)
                      *(-1+x1+2.*y1+z1);
          dphi_M(3,r)= unsurhuit *(1. + x1)*(1. + y1)*ptelem(r)(3)
                      *(-1+x1+y1+2.*z1);
         }
       // pour les noeuds du plan xi=0
       double unsurquatre = 1./4.;double unsur2 = 1./2.;
       Tableau<int> ind(4);
       ind(1) = 9; ind(2) = 11; ind(3) = 19; ind(4) = 17;
       for (int r=1;r<= 4; r++)
        { double y1 = M(2) * ptelem(ind(r))(2);double  z1 = M(3) * ptelem(ind(r))(3);
          dphi_M(1,ind(r))= -unsur2 * M(1)  *(1. + y1)*(1. + z1);
          dphi_M(2,ind(r))= unsurquatre *ptelem(ind(r))(2)*(1. - M(1)*M(1))*(1. + z1);
          dphi_M(3,ind(r))= unsurquatre *ptelem(ind(r))(3)*(1. - M(1)*M(1))*(1. + y1);
         }
       // pour les noeuds du plan eta=0
       ind(1) = 10; ind(2) = 18; ind(3) = 20; ind(4) = 12;
       for (int r=1;r<= 4; r++)
        { double  x1 = M(1) * ptelem(ind(r))(1);double  z1 = M(3) * ptelem(ind(r))(3);
          dphi_M(1,ind(r))= unsurquatre *ptelem(ind(r))(1)*(1. - M(2)*M(2))*(1. + z1);
          dphi_M(2,ind(r))= -unsur2 * M(2)  *(1. + x1)*(1. + z1);
          dphi_M(3,ind(r))= unsurquatre *ptelem(ind(r))(3)*(1. + x1)*(1. - M(2)*M(2));
         }
       // pour les noeuds du plan zeta=0
       ind(1) = 13; ind(2) = 14; ind(3) = 15; ind(4) = 16;
       for (int r=1;r<= 4; r++)
        { double  x1 = M(1) * ptelem(ind(r))(1);double  y1 = M(2) * ptelem(ind(r))(2);
          dphi_M(1,ind(r))= unsurquatre *ptelem(ind(r))(1)*(1. + y1)*(1. - M(3)*M(3));
          dphi_M(2,ind(r))= unsurquatre *ptelem(ind(r))(2)*(1. + x1)*(1. - M(3)*M(3));
          dphi_M(3,ind(r))= -unsur2 * M(3)  *(1. + x1)*(1. + y1);
         }
       }         
      else
       { cout << "\n erreur le nombre de noeud demande :" << NBNE <<"n\'est pas implante "
             << "\nGeomHexaQuad::Dphi(Coordonnee& M)";
        Sortie(1);
      }
    // retour des derivees  
    return dphi_M;
   };

// constitution du tableau Extrapol
void GeomHexaQuad::Calcul_extrapol(int nbi)
{ // cas de l'extrapolation de grandeur des points d'intégrations aux noeuds
  // def du tableau de pondération tab(i)(j) qu'il faut appliquer
  // aux noeuds pour avoir la valeur aux noeuds 
  // val_au_noeud(i) = somme_(de j=indir(i)(1) à indir(i)(taille(indir(i)) )) {tab(i)(j) * val_pt_integ(j) }
  // cas = 1: la valeur au noeud = la valeur au pt d'integ le plus près ou une moyenne des
  //          pt les plus près (si le nb de pt d'integ < nb noeud)
  //--- changement: 26 oct 2020: on fait une extrapolation tri-linéaire via les pti les plus proches
  


  switch  (nbi)
   { case 1: 
       { // cas avec un point d'intégration, on reporte la valeur au pt d'integ, telle quelle au noeud
         Tableau<Tableau<int> > & indir = extrapol(1).indir;  // pour simplifier
         Tableau<Tableau<double > > & tab = extrapol(1).tab;    // pour simplifier
         tab.Change_taille(NBNE);
         for (int ne=1;ne<=NBNE;ne++)
         	{ tab(ne).Change_taille(nbi);
            tab(ne)(1)=1.;
         	  indir(ne).Change_taille(1); indir(ne)(1)=1;
         	};
         break;
       }
   	 case 8:
   	   { // cas avec 8 points d'intégration ,
         //on définit un hexaèdre linéaire qui va nous permettre d'extrapoler
         // pour les nbi >= 8
         GeomHexalin hexa(8);

         // on a deux méthodes
         extrapol.Change_taille(2);
         // tab est supposé être initialisé à 0.
         Tableau <Coordonnee>  gi_B,gi_H; // bases naturelle et duale
         
         // --- méthode 1 (par défaut), on utilise une extrapolation via un hexaèdre bi-linéaire
         {Tableau<Tableau<int> > & indir = extrapol(1).indir;  // pour simplifier
          Tableau<Tableau<double > > & tab = extrapol(1).tab;    // pour simplifier
          Tableau <int> indirect(8); // tableau de travail: on a 8 pondérations
          tab.Change_taille(NBNE);indir.Change_taille(NBNE);

          // on extrapole bi-linéairement vers les  noeuds en considérant
          // les 8 pti dans l'ordre actuel: il faut utiliser une indirection
          // car la numérotation des pti ne correspond pas à celle des noeuds d'un hexa bi-linéaire
          indirect(1) = 8; indirect(2) = 4; indirect(3) = 2; indirect(4) = 6;
          indirect(5) = 7; indirect(6) = 3; indirect(7) = 1; indirect(8) = 5;
          

          for (int ne=1; ne <= NBNE; ne++)
           { // il nous faut calculer les coordonnées locales.
             // sachant que les pti ici considérés, sont aux sommets d'un hexa linéaire orthogonal
             // on peut traiter séparément chaque composante
             Coordonnee theta(3); // les coordonnées que l'on cherche
             Vecteur phi_x(2); // le conteneur pour les fonctions d'interpolation
             Vecteur phi_y(2); // le conteneur pour les fonctions d'interpolation
             Vecteur phi_z(2); // le conteneur pour les fonctions d'interpolation
             {Tableau <int> indirect_local(2); // tableau de travail
              Coordonnee theta_loc(1); // le conteneur pour les coordonnées locales
              Tableau <Coordonnee>  gi_loc_B,gi_loc_H; // bases naturelle et duale
              Vecteur phi_loc(2); // le conteneur pour les fonctions d'interpolation
              // suivant x
              {indirect_local(1) = indirect(1);indirect_local(2) = indirect(2);
               Bases_naturel_duales(indirect_local,gi_loc_B,gi_loc_H);
               Coordonnee O(0.5*(ptInteg(indirect_local(1))+ptInteg(indirect_local(2)))); // def de l'origine
               ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_x,theta_loc);
               theta(1)=theta_loc(1); // on enregistre
              }
              // suivant y
              {indirect_local(1) = indirect(1);indirect_local(2) = indirect(4);
               Bases_naturel_duales(indirect_local,gi_loc_B,gi_loc_H);
               Coordonnee O(0.5*(ptInteg(indirect_local(1))+ptInteg(indirect_local(2)))); //
               ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_y,theta_loc);
               theta(2)=theta_loc(1); // on enregistre
              }
              // suivant z
              {indirect_local(1) = indirect(1);indirect_local(2) = indirect(8);
               Bases_naturel_duales(indirect_local,gi_loc_B,gi_loc_H);
               Coordonnee O(0.5*(ptInteg(indirect_local(1))+ptInteg(indirect_local(2)))); //
               ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_z,theta_loc);
               theta(3)=theta_loc(1); // on enregistre
              }
             };

             // maintenant on calcule les fct d'interpolation au noeud ne
             // via ses coordonnées locales theta
             const Vecteur& phiphi = hexa.Phi_point(theta);
             // et on enregistre
             indir(ne).Change_taille(8);
             tab(ne).Change_taille(nbi);
             // on boucle sur les pti de l'hexa linéaire d'interpolation
             for (int i=1;i<9;i++)
              {tab(ne)(indirect(i))=phiphi(i);
               indir(ne)(i)=indirect(i);
              };
           }
         }; // fin méthode 1
         
         // --- deuxième méthode en utilisant un tétraèdre linéaire
         // le pb est que le choix du tétraèdre est vraiment arbitraire

         {// on extrapole linéairement vers les  noeuds en considérant à chaque fois les 4 pt d'intégration les plus près qui ne doivent pas être coplanaires
          Tableau<Tableau<int> > & indir = extrapol(2).indir;  // pour simplifier
          Tableau<Tableau<double > > & tab = extrapol(2).tab;    // pour simplifier
          tab.Change_taille(NBNE);indir.Change_taille(NBNE);
          //
          // -- le premier noeud
          {int ne = 1; indir(ne).Change_taille(4);
           indir(ne)(1)=8;indir(ne)(2)=4;indir(ne)(3)=6;indir(ne)(4)=7;
          }
          // -- le deuxième noeud
          {
           int ne = 2; indir(ne).Change_taille(4);
           indir(ne)(1)=4;indir(ne)(2)=2;indir(ne)(3)=6;indir(ne)(4)=3;
          }
          // -- le troisième noeud
          {int ne = 3; indir(ne).Change_taille(4);
           indir(ne)(1)=2;indir(ne)(2)=6;indir(ne)(3)=4;indir(ne)(4)=1;
          }
          // -- le quatrième noeud
          {int ne = 4; indir(ne).Change_taille(4);
           indir(ne)(1)=6;indir(ne)(2)=8;indir(ne)(3)=2;indir(ne)(4)=5;
          }
          // -- le cinquième noeud
          {int ne = 5; indir(ne).Change_taille(4);
           indir(ne)(1)=7;indir(ne)(2)=5;indir(ne)(3)=3;indir(ne)(4)=8;
          }
          // -- le sixième noeud
          {int ne = 6; indir(ne).Change_taille(4);
           indir(ne)(1)=3;indir(ne)(2)=7;indir(ne)(3)=1;indir(ne)(4)=4;
          }
          // -- le septième noeud
          {int ne = 7; indir(ne).Change_taille(4);
           indir(ne)(1)=1;indir(ne)(2)=3;indir(ne)(3)=5;indir(ne)(4)=2;
          }
          // -- le huitième noeud
          {int ne = 8; indir(ne).Change_taille(4);
           indir(ne)(1)=5;indir(ne)(2)=1;indir(ne)(3)=7;indir(ne)(4)=6;
          }
          // -- le neuvième noeud
          {int ne = 9; indir(ne).Change_taille(4);
           indir(ne)(1)=8;indir(ne)(2)=4;indir(ne)(3)=6;indir(ne)(4)=7;
          }
          // -- le dixième noeud
          {int ne = 10; indir(ne).Change_taille(4);
           indir(ne)(1)=4;indir(ne)(2)=2;indir(ne)(3)=6;indir(ne)(4)=3;
          }
          // -- le onzième noeud
          {int ne = 11; indir(ne).Change_taille(4);
           indir(ne)(1)=2;indir(ne)(2)=6;indir(ne)(3)=4;indir(ne)(4)=1;
          }
          // -- le douzième noeud
          {int ne = 12; indir(ne).Change_taille(4);
           indir(ne)(1)=6;indir(ne)(2)=8;indir(ne)(3)=4;indir(ne)(4)=5;
          }
          // -- le  noeud 13 .. idem noeud 1
          {int ne = 13; indir(ne).Change_taille(4);
           indir(ne)(1)=8;indir(ne)(2)=4;indir(ne)(3)=6;indir(ne)(4)=7;
          }
          // -- le noeud 14 .. idem noeud 2
          {int ne = 14; indir(ne).Change_taille(4);
           indir(ne)(1)=4;indir(ne)(2)=2;indir(ne)(3)=6;indir(ne)(4)=3;
          }
          // -- le  noeud 15 .. idem noeud 3
          {int ne = 15; indir(ne).Change_taille(4);
           indir(ne)(1)=2;indir(ne)(2)=6;indir(ne)(3)=4;indir(ne)(4)=1;
          }
          // -- le  noeud 16 .. idem le noeud 4
          {int ne = 16; indir(ne).Change_taille(4);
           indir(ne)(1)=6;indir(ne)(2)=8;indir(ne)(3)=2;indir(ne)(4)=5;
          }
          // -- le  noeud 17 .. idem le noeud 5
          {int ne = 17; indir(ne).Change_taille(4);
           indir(ne)(1)=7;indir(ne)(2)=5;indir(ne)(3)=3;indir(ne)(4)=8;
          }
          // -- le noeud 18 .. idem noeud 6
          {int ne = 18; indir(ne).Change_taille(4);
           indir(ne)(1)=3;indir(ne)(2)=7;indir(ne)(3)=1;indir(ne)(4)=4;
          }
          // -- le  noeud 19 .. idem noeud 7
          {int ne = 19; indir(ne).Change_taille(4);
           indir(ne)(1)=1;indir(ne)(2)=3;indir(ne)(3)=5;indir(ne)(4)=2;
          }
          // -- le  noeud 20 .. idem noeud 8
          {int ne = 20; indir(ne).Change_taille(4);
           indir(ne)(1)=5;indir(ne)(2)=1;indir(ne)(3)=7;indir(ne)(4)=6;
          }
          
          // 2) ---  on calcule l'interpolation
          Coordonnee theta(3); // les coordonnées que l'on cherche
          Tableau <Coordonnee>  gi_B,gi_H; // bases naturelle et duale
          Vecteur phi_(4); // le conteneur pour les fonctions d'interpolation
          for (int ne = 1;ne<=NBNE; ne++)
           { tab(ne).Change_taille(nbi);// 8 pti
             Tableau <int>& indirect= indir(ne) ; // tableau de travail: on a 4 pondérations
             Bases_naturel_duales(indirect,gi_B,gi_H);
             Coordonnee O(ptInteg(indirect(1))); // def de l'origine
             ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
             tab(ne)(indirect(1)) = phi_(1);tab(ne)(indirect(2)) = phi_(2);
             tab(ne)(indirect(3)) = phi_(3);tab(ne)(indirect(4)) = phi_(4);
           };
                
         }; // fin méthode 2

       
{//  ancienne méthode sans extrapolation --> donne des résultats vraiment moins bons !!
//         //on exporte directement la valeur du
//   	     // pt d'integ le plus proche ou d'une moyenne
//          int ne = 1; tab(ne)(8) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=8;
//          ne = 2; tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4;
//          ne = 3; tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
//          ne = 4; tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
//          ne = 5; tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;
//          ne = 6; tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
//          ne = 7; tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
//          ne = 8; tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
//          ne = 9;
//          tab(ne)(4) = 0.5;tab(ne)(8) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=4;indir(ne)(2)=8;
//          ne = 10;
//          tab(ne)(2) = 0.5;tab(ne)(4) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=2;indir(ne)(2)=4;
//          ne = 11;
//          tab(ne)(2) = 0.5;tab(ne)(6) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=2;indir(ne)(2)=6;
//          ne = 12;
//          tab(ne)(6) = 0.5;tab(ne)(8) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=6;indir(ne)(2)=8;
//          ne = 13;
//          tab(ne)(7) = 0.5;tab(ne)(8) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=7;indir(ne)(2)=8;
//          ne = 14;
//          tab(ne)(3) = 0.5;tab(ne)(4) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=3;indir(ne)(2)=4;
//          ne = 15;
//          tab(ne)(1) = 0.5;tab(ne)(2) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=1;indir(ne)(2)=2;
//          ne = 16;
//          tab(ne)(5) = 0.5;tab(ne)(6) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=5;indir(ne)(2)=6;
//          ne = 17;
//          tab(ne)(3) = 0.5;tab(ne)(7) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=3;indir(ne)(2)=7;
//          ne = 18;
//          tab(ne)(1) = 0.5;tab(ne)(3) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=1;indir(ne)(2)=3;
//          ne = 19;
//          tab(ne)(1) = 0.5;tab(ne)(5) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=1;indir(ne)(2)=5;
//          ne = 20;
//          tab(ne)(5) = 0.5;tab(ne)(7) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=5;indir(ne)(2)=7;
}
          break;  
   	 	}  // fin du cas avec 8 pt d'intégration
   	 case 27:
   	   { // cas avec 27 points d'intégration  

         // on va utiliser un hexaèdre quadratique complet et retenir que les 20 premier noeuds
         GeomHexaQuadComp hexa_inter;
         const ConteneurExtrapolation extrapol_inter =  hexa_inter.ExtrapolationNoeud(1);
         const Tableau<Tableau<int> > & indir_inter = extrapol_inter.indir;  // pour simplifier
         const Tableau<Tableau<double > > & tab_inter = extrapol_inter.tab;    // pour simplifier

         Tableau<Tableau<int> > & indir = extrapol(1).indir;  // pour simplifier
         Tableau<Tableau<double > > & tab = extrapol(1).tab;    // pour simplifier
         tab.Change_taille(NBNE);indir.Change_taille(NBNE);
         
         
         for (int ne=1; ne <= NBNE; ne++) // on boucle de 1 à 20
          {indir(ne)=indir_inter(ne);
           tab(ne).Change_taille(nbi);
           tab(ne)=tab_inter(ne);
          };



// programmation directe qui doit marcher, mais c'est plus simple d'utiliser
// un quadratique complet
{//         //on définit un hexaèdre linéaire qui va nous permettre d'extrapoler
//         // pour les nbi >= 8
//         GeomHexalin hexa(8);
//
//         // les pti sont numérotés à l'aide de la numérotation sur un segment
//         // contrairement à la numérotation des noeuds qui suit une autre logique
//         // du coup pour avoir le pti qui est en face (à peu près !) du noeud
//         // on utilise une numérotation indirecte: ind(i) = le numéro du noeud
//         // ou du pti i
//         Tableau<int> ind;
//         ind.Change_taille(NBNE);
//         ind(1) = 1; ind(2) = 3; ind(3) = 9; ind(4) = 7;ind(5) = 19;
//         ind(6) = 21; ind(7) = 27; ind(8) = 25; ind(9) = 2;ind(10) = 6;
//         ind(11) = 8; ind(12) = 4; ind(13) = 10; ind(14) = 12;ind(15) = 18;
//         ind(16) = 16; ind(17) = 20; ind(18) = 24; ind(19) = 26;ind(20) = 22;
//         ind(21) = 5; ind(22) = 11; ind(23) = 15; ind(24) = 17;ind(25) = 13;
//         ind(26) = 23; ind(27) = 14;
//
//         // et l'inverse
//         Tableau<int> jnd;
//         jnd.Change_taille(NBNE);
//         jnd(1) = 1; jnd(3) = 2; jnd(9) = 3; jnd(7) = 4;jnd(19) = 5;
//         jnd(21) = 6; jnd(27) = 7; jnd(25) = 8; jnd(2) = 9;jnd(6) = 10;
//         jnd(8) = 11; jnd(4) = 12; jnd(10) = 13; jnd(12) = 14;jnd(18) = 15;
//         jnd(16) = 16; jnd(20) = 17; jnd(24) = 18; jnd(26) = 19;jnd(22) = 20;
//         jnd(5) = 21; jnd(11) = 22; jnd(15) = 23; jnd(17) = 24;jnd(13) = 25;
//         jnd(23) = 26; jnd(14) = 27;
//
//         // on a deux méthodes
//         extrapol.Change_taille(2);
//         Tableau <Coordonnee>  gi_B,gi_H; // bases naturelle et duale
//
//         // --- méthode 1 (par défaut), on utilise une extrapolation via un hexaèdre tri-linéaire
//
//         {Tableau<Tableau<int> > & indir = extrapol(1).indir;  // pour simplifier
//          Tableau<Tableau<double > > & tab = extrapol(1).tab;    // pour simplifier
//          Tableau <int> indirect(8); // tableau de travail: on a 8 pondérations
//          tab.Change_taille(NBNE);indir.Change_taille(NBNE);
//
//          // on essaie de faire une boucle mais il y a des exceptions ...
//          int ni = 0;
//          for (int niz = 1; niz < 4; niz++)
//           for (int niy = 1; niy < 4; niy++)
//             for (int nix = 1; nix < 4; nix++)
//           {
//            ni++; // -- le noeud ne le plus proche
//
//            // on supprime les pti aux centre des faces et celui à l'origine
//            if ((ni != 5) && (ni != 11)&& (ni != 13)&& (ni != 14)&& (ni != 15)
//                && (ni != 17)&& (ni != 23)
//               )
//            {int ne = jnd(ni);
//             if (nix < 3)
//              {if (niy < 3)
//                {if (niz < 3) // nix < 3, niy < 3, niz < 3  OK
//                  {indirect(1)=ni;indirect(2)=ni+1;indirect(3)=ni+4;indirect(4)=ni+3;
//                   indirect(5)=ni+9;indirect(6)=ni+10;indirect(7)=ni+13;indirect(8)=ni+12;
//                  }
//                 else // nix < 3, niy < 3, niz = 3 OK
//                  {indirect(1)=ni-9;indirect(2)=ni-8;indirect(3)=ni-5;indirect(4)=ni-6;
//                   indirect(5)=ni;indirect(6)=ni+1;indirect(7)=ni+4;indirect(8)=ni+3;
//                  }
//                }
//               else // cas niy = 3
//                {if (niz < 3) // nix < 3, niy = 3, niz < 3 OK
//                  {indirect(1)=ni;indirect(2)=ni-3;indirect(3)=ni-2;indirect(4)=ni+1;
//                   indirect(5)=ni+9;indirect(6)=ni+6;indirect(7)=ni+7;indirect(8)=ni+10;
//                  }
//                 else // // nix < 3, niy = 3, niz = 3 OK
//                  {indirect(1)=ni-9;indirect(2)=ni-12;indirect(3)=ni-11;indirect(4)=ni-8;
//                   indirect(5)=ni;indirect(6)=ni-3;indirect(7)=ni-2;indirect(8)=ni+1;
//                  }
//                }
//              }
//             else // (nix == 3)
//              {if (niy < 3)
//                {if (niz < 3) // nix = 3, niy < 3, niz < 3 OK
//                  {indirect(1)=ni;indirect(2)=ni+3;indirect(3)=ni+2;indirect(4)=ni-1;
//                   indirect(5)=ni+9;indirect(6)=ni+12;indirect(7)=ni+11;indirect(8)=ni+8;
//                  }
//                 else // nix = 3, niy < 3, niz = 3 OK
//                  {indirect(1)=ni-9;indirect(2)=ni-6;indirect(3)=ni-7;indirect(4)=ni-10;
//                   indirect(5)=ni;indirect(6)=ni+3;indirect(7)=ni+2;indirect(8)=ni-1;
//                  }
//                }
//               else // cas niy = 3
//                {if (niz < 3) // nix = 3, niy = 3, niz < 3 OK
//                  {indirect(1)=ni;indirect(2)=ni-1;indirect(3)=ni-4;indirect(4)=ni-3;
//                   indirect(5)=ni+9;indirect(6)=ni+8;indirect(7)=ni+5;indirect(8)=ni+6;
//                  }
//                 else // // nix = 3, niy = 3, niz = 3 OK
//                  {indirect(1)=ni-9;indirect(2)=ni-10;indirect(3)=ni-13;indirect(4)=ni-12;
//                   indirect(5)=ni;indirect(6)=ni-1;indirect(7)=ni-4;indirect(8)=ni-3;
//                  }
//                }
//              }
//
//             indir(ne).Change_taille(8);
//             // il nous faut calculer les coordonnées locales.
//             // sachant que les pti ici considérés, sont aux sommets d'un hexa linéaire orthogonal
//             // on peut traiter séparément chaque composante
//             Coordonnee theta(3); // les coordonnées que l'on cherche
//             Vecteur phi_x(2); // le conteneur pour les fonctions d'interpolation
//             Vecteur phi_y(2); // le conteneur pour les fonctions d'interpolation
//             Vecteur phi_z(2); // le conteneur pour les fonctions d'interpolation
//             {Tableau <int> indirect_local(2); // tableau de travail
//              Coordonnee theta_loc(1); // le conteneur pour les coordonnées locales
//              Tableau <Coordonnee>  gi_loc_B,gi_loc_H; // bases naturelle et duale
//              Vecteur phi_loc(2); // le conteneur pour les fonctions d'interpolation
//              // suivant x
//              {indirect_local(1) = indirect(1);indirect_local(2) = indirect(2);
//               Bases_naturel_duales(indirect_local,gi_loc_B,gi_loc_H);
//               Coordonnee O(0.5*(ptInteg(indirect_local(1))+ptInteg(indirect_local(2)))); // def de l'origine
//               ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_x,theta_loc);
//               theta(1)=theta_loc(1); // on enregistre
//              }
//              // suivant y
//              {indirect_local(1) = indirect(1);indirect_local(2) = indirect(4);
//               Bases_naturel_duales(indirect_local,gi_loc_B,gi_loc_H);
//               Coordonnee O(0.5*(ptInteg(indirect_local(1))+ptInteg(indirect_local(2)))); // def
//               ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_y,theta_loc);
//               theta(2)=theta_loc(1); // on enregistre
//              }
//              // suivant z
//              {indirect_local(1) = indirect(1);indirect_local(2) = indirect(8);
//               Bases_naturel_duales(indirect_local,gi_loc_B,gi_loc_H);
//               Coordonnee O(0.5*(ptInteg(indirect_local(1))+ptInteg(indirect_local(2)))); // def de             ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_z,theta_loc);
//               ElemGeomC0::Coor_phi(O,gi_loc_H,ptelem(ne),phi_z,theta_loc);
//               theta(3)=theta_loc(1); // on enregistre
//              }
//             };
//             // maintenant on calcule les fct d'interpolation au noeud ne
//             // via ses coordonnées locales theta
//             const Vecteur& phiphi = hexa.Phi(theta);
//             // et on enregistre
//             // on boucle sur les pti de l'hexa linéaire d'interpolation
//             for (int i=1;i<9;i++)
//              {tab(ne)(indirect(i))=phiphi(i);
//               indir(ne)(i)=indirect(i);
//              };
//           };// fin d'exclusion des noeuds au centre des faces
//
//          };
//         }; // fin de la méthode 1
//
//        // --- deuxième méthode en utilisant un tétraèdre linéaire
//        // le pb est que le choix du tétraèdre est vraiment arbitraire
//
//
//        {Tableau<Tableau<int> > & indir = extrapol(2).indir;  // pour simplifier
//         Tableau<Tableau<double > > & tab = extrapol(2).tab;    // pour simplifier
//         Tableau <int> indirect(4); // tableau de travail: on a 4 pondérations
//         Coordonnee theta(3); // les coordonnées que l'on cherche
//         Tableau <Coordonnee>  gi_B,gi_H; // bases naturelle et duale
//         Vecteur phi_(4); // le conteneur pour les fonctions d'interpolation
//         tab.Change_taille(NBNE);indir.Change_taille(NBNE);
//         // on essaie de faire une boucle mais il y a des exceptions ...
//         int ni = 0;
//         for (int niz = 1; niz < 4; niz++)
//          for (int niy = 1; niy < 4; niy++)
//            for (int nix = 1; nix < 4; nix++)
//          {// -- le noeud ne le plus proche
//           ni++;
//           // on supprime les pti aux centre des faces et celui à l'origine
//           if ((ni != 5) && (ni != 11)&& (ni != 13)&& (ni != 14)&& (ni != 15)
//               && (ni != 17)&& (ni != 23)
//              )
//           {
//             int ne = ind(ni);
//             if (nix < 3)
//              {if (niy < 3)
//                {if (niz < 3) // nix < 3, niy < 3, niz < 3  OK
//                  {indirect(1)=ni;indirect(2)=ni+9;indirect(3)=ni+1;indirect(4)=ni+3;}
//                 else // nix < 3, niy < 3, niz = 3 OK
//                  {indirect(1)=ni;indirect(2)=ni+1;indirect(3)=ni-9;indirect(4)=ni+3;}
//                }
//               else // cas niy = 3
//                {if (niz < 3) // nix < 3, niy = 3, niz < 3 OK
//                  {indirect(1)=ni;indirect(2)=ni+1;indirect(3)=ni+9;indirect(4)=ni-3;}
//                 else // // nix < 3, niy = 3, niz = 3 OK
//                  {indirect(1)=ni;indirect(2)=ni-9;indirect(3)=ni+1;indirect(4)=ni-3;}
//                }
//              }
//             else // (nix == 3)
//              {if (niy < 3)
//                {if (niz < 3) // nix = 3, niy < 3, niz < 3 OK
//                  {indirect(1)=ni;indirect(2)=ni+9;indirect(3)=ni+3;indirect(4)=ni-1;}
//                 else // nix = 3, niy < 3, niz = 3  OK
//                  {indirect(1)=ni;indirect(2)=ni+3;indirect(3)=ni-9;indirect(4)=ni-1;}
//                }
//               else // cas niy = 3
//                {if (niz < 3) // nix = 3, niy = 3, niz < 3 OK
//                  {indirect(1)=ni;indirect(2)=ni+9;indirect(3)=ni-1;indirect(4)=ni-3;}
//                 else // // nix = 3, niy = 3, niz = 3  OK
//                  {indirect(1)=ni;indirect(2)=ni-1;indirect(3)=ni-9;indirect(4)=ni-3;}
//                }
//              }
//
//             Bases_naturel_duales(indirect,gi_B,gi_H);
//             Coordonnee O(ptInteg(indirect(1))); // def de l'origine
//             ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
//             tab(ne).Change_taille(nbi);
//             tab(ne)(indirect(1)) = phi_(1);tab(ne)(indirect(2)) = phi_(2);
//             tab(ne)(indirect(3)) = phi_(3);tab(ne)(indirect(4)) = phi_(4);
//             indir(ne).Change_taille(4);
//             indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);
//             indir(ne)(3)=indirect(3);indir(ne)(4)=indirect(4);
//           };// fin d'exclusion des noeuds au centre des faces
//          };
//         }; // -- fin méthode 2
}
{//         //  ancienne méthode sans extrapolation --> donne des résultats vraiment moins bons !!
//         // on exporte directement la valeur du pt d'integ le plus proche
//          int ne = 1; tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
//          ne = 2; tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
//          ne = 3; tab(ne)(9) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=9;
//          ne = 4; tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;
//          ne = 5; tab(ne)(19) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=19;
//          ne = 6; tab(ne)(21) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=21;
//          ne = 7; tab(ne)(27) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=27;
//          ne = 8; tab(ne)(25) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=25;
//
//          ne = 9; tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
//          ne = 10; tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
//          ne = 11; tab(ne)(8) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=8;
//          ne = 12; tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4;
//          ne = 13; tab(ne)(10) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=10;
//          ne = 14; tab(ne)(12) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=12;
//          ne = 15; tab(ne)(18) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=18;
//          ne = 16; tab(ne)(16) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=16;
//          ne = 17; tab(ne)(20) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=20;
//          ne = 18; tab(ne)(24) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=24;
//          ne = 19; tab(ne)(26) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=26;
//          ne = 20; tab(ne)(22) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=22;
}
          break;  
   	 	}  // fin du cas avec 27 pt d'intégration
   	 case 64:
   	   { // cas avec 64 points d'intégration
       
         // on va utiliser un hexaèdre quadratique complet et retenir que les 20 premier noeuds
         GeomHexaQuadComp hexa_inter;
         const ConteneurExtrapolation extrapol_inter =  hexa_inter.ExtrapolationNoeud(1);
         const Tableau<Tableau<int> > & indir_inter = extrapol_inter.indir;  // pour simplifier
         const Tableau<Tableau<double > > & tab_inter = extrapol_inter.tab;    // pour simplifier

         Tableau<Tableau<int> > & indir = extrapol(1).indir;  // pour simplifier
         Tableau<Tableau<double > > & tab = extrapol(1).tab;    // pour simplifier
         tab.Change_taille(NBNE);indir.Change_taille(NBNE);
         
         for (int ne=1; ne <= NBNE; ne++) // on boucle de 1 à 20
          {indir(ne)=indir_inter(ne);
           tab(ne).Change_taille(nbi);
           tab(ne)=tab_inter(ne);
          };
cout << "\n *** methode en developpement !! "
     << "\n void GeomHexaQuad::Calcul_extrapol(int nbi)" << flush;
Sortie(1);

// ---- ancienne méthode ... vraiment pas assez précise
{//         // on exporte directement la valeur du pt d'integ le plus proche
//          int ne = 1; tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
//          ne = 2; tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4;
//          ne = 3; tab(ne)(16) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=16;
//          ne = 4; tab(ne)(13) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=13;
//          ne = 5; tab(ne)(49) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=49;
//          ne = 6; tab(ne)(52) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=52;
//          ne = 7; tab(ne)(64) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=64;
//          ne = 8; tab(ne)(61) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=61;
//
//          ne = 9;
//          tab(ne)(2) = 0.5;tab(ne)(3) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=2;indir(ne)(2)=3;
//          ne = 10;
//          tab(ne)(8) = 0.5;tab(ne)(12) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=8;indir(ne)(2)=12;
//          ne = 11;
//          tab(ne)(14) = 0.5;tab(ne)(15) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=14;indir(ne)(2)=15;
//          ne = 12;
//          tab(ne)(5) = 0.5;tab(ne)(9) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=5;indir(ne)(2)=9;
//          ne = 13;
//          tab(ne)(17) = 0.5;tab(ne)(33) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=17;indir(ne)(2)=33;
//          ne = 14;
//          tab(ne)(20) = 0.5;tab(ne)(36) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=20;indir(ne)(2)=36;
//          ne = 15;
//          tab(ne)(32) = 0.5;tab(ne)(48) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=32;indir(ne)(2)=48;
//          ne = 16;
//          tab(ne)(29) = 0.5;tab(ne)(45) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=29;indir(ne)(2)=45;
//          ne = 17;
//          tab(ne)(50) = 0.5;tab(ne)(51) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=50;indir(ne)(2)=51;
//          ne = 18;
//          tab(ne)(56) = 0.5;tab(ne)(60) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=56;indir(ne)(2)=60;
//          ne = 19;
//          tab(ne)(62) = 0.5;tab(ne)(63) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=62;indir(ne)(2)=63;
//          ne = 20;
//          tab(ne)(53) = 0.5;tab(ne)(57) = 0.5;
//          indir(ne).Change_taille(2); indir(ne)(1)=53;indir(ne)(2)=57;
}
          break;  
   	 	}  // fin du cas avec 64 pt d'intégration
          
     default:
       { cout << "\n erreur le nombre de point d'integration demande :" << nbi <<"n\'est pas implante "
              << "\nGeomTriangle::Calcul_extrapol(..";
         Sortie(1);
       };
   };

};