
// 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 "GeomTriangle.h"
#include "GeomSeg.h"
#include "MathUtil.h"

// constructeur
// dimension 2, nbi pt integ, nbne noeuds, 1 face, 3 aretes
// lorsqu'il y a plusieurs cas avec un même nombre de points d'intégration 
// on les sépare par 1000: exemple: 3 ou 1003

GeomTriangle::GeomTriangle( int nbi, int nbne) :
   ElemGeomC0(2,(nbi%1000),nbne,1,3,TRIANGLE,RIEN_INTERPOL)
   ,phi_M(),dphi_M()
  { // definition de la face
    face(1) = this;
    NONF(1).Change_taille(NBNE);
    for (int i=1;i<=NBNE;i++)
        NONF(1)(i) = i;// def du pointeur de ligne
        
    // définition de la numérotation locale de l'élément de direction inverse
    switch  (NBNE)
        { case 3: 
            { INVCONNEC(1) = 1;INVCONNEC(2) = 3;
              INVCONNEC(3) = 2;
				          // le tableau des tranches
				          IND.Change_taille(1); IND(1)=3;
              break;
             }  
          case 6: 
             { INVCONNEC(1) = 1;INVCONNEC(2) = 3; INVCONNEC(3) = 2;
               INVCONNEC(4) = 6;INVCONNEC(5) = 5;INVCONNEC(6) = 4; 
				           // le tableau des tranches
				           IND.Change_taille(2); IND(1)=3;IND(2)=3;
               break;
              }  
          case 10: // élément cubique complet (cf dhatt et touzot, pg 113 ed 1984)
             { // avec une numérotation différente (dabord les sommets puis les noeuds internes
               // puis le noeud du centre
               INVCONNEC(1) = 1;INVCONNEC(2) = 3; INVCONNEC(3) = 2;
               INVCONNEC(4) = 9;INVCONNEC(5) = 8; INVCONNEC(6) = 7;
               INVCONNEC(7) = 6;INVCONNEC(8) = 5; INVCONNEC(9) = 4; 
               INVCONNEC(10) = 10; // noeud du centre
				           // le tableau des tranches
				           IND.Change_taille(3); IND(1)=3;IND(2)=6;IND(3)=1;
               break;
              }  
          case 9: // élément cubique incomplet (cf dhatt et touzot, pg 114 ed 1984)
             { INVCONNEC(1) = 1;INVCONNEC(2) = 3; INVCONNEC(3) = 2;
               INVCONNEC(4) = 9;INVCONNEC(5) = 8; INVCONNEC(6) = 7;
               INVCONNEC(7) = 6;INVCONNEC(8) = 5; INVCONNEC(9) = 4; 
				           // le tableau des tranches
				           IND.Change_taille(3); IND(1)=3;IND(2)=6;IND(3)=1;
				           // le tableau des tranches
				           IND.Change_taille(2); IND(1)=3;IND(2)=6;
               break;
              }  
          default :
                  { cout << "\n erreur le nombre de noeud demande :" << NBNE <<" n\'est"
                         << " pas implante "
                         << "\nGeomTriangle::GeomTriangle( int nbi, int nbne)";
                    Sortie(1);
                   }
        }
    double lambda; // variable de calcul intermediaire

   //      ------------------------------
   //      choix des points d'integration
   //      ------------------------------
    // dans le choix on utilise nbi plutôt que Nbi() pour différencier les deux cas avec 3 pt d'integ
    if (nbi ==1) 
     {
      KSI(1)= 1./3.;
      ETA(1)= KSI(1);
      WI(1)= 0.5 ;
      }
    else  if (nbi   == 3)
     { // cas avec des points d'intégrations sur les arrêtes
      KSI(1)= 0.5; 
      ETA(1)= 0.5; 
      WI(1)= 1. /6.; 
      KSI(2)= 0.;
      ETA(2)= 0.5 ;
      WI(2)= 1. /6.; 
      KSI(3)= 0.5 ;
      ETA(3)= 0. ;
      WI(3)= 1. /6. ;
      }
    else  if (nbi == 1003) // il n'y a cependant que 3 pt d'integ ==> Nbi() = 3
     { // cas avec des points d'intégrations qui ne sont pas sur les arrêtes
      KSI(1)= 1./6.; 
      ETA(1)= 1./6.; 
      WI(1)= 1. /6.; 
      KSI(2)= 2./3.;
      ETA(2)= 1./6. ;
      WI(2)= 1./6.; 
      KSI(3)= 1./6. ;
      ETA(3)= 2./3. ;
      WI(3)= 1./6. ;
     }
    else  if (nbi==4)
     {
      KSI(1)= 1. /3.; 
      ETA(1)= 1. /3. ;
      WI(1)= -27. /96.; 
      KSI(2)= 1. /5.; 
      ETA(2)= 1. /5. ;
      WI(2)= 25. /96. ;
      KSI(3)= 3. /5.;  
      ETA(3)= 1. /5.; 
      WI(3)= 25. /96.; 
      KSI(4)= 1. /5.;  
      ETA(4)= 3. /5.; 
      WI(4)= 25. /96.;
     }
    else  if (nbi==6)
     {double a = 0.445948490915965;
      double b = 0.091576213509771;
      KSI(1)= a; 
      ETA(1)= a ;
      WI(1)= 0.111690794839005; 
      KSI(2)= 1.-2.* a ; 
      ETA(2)= a ;
      WI(2)= 0.111690794839005;
      KSI(3)= a;  
      ETA(3)= 1.-2.* a; 
      WI(3)= 0.111690794839005; 
      KSI(4)= b;  
      ETA(4)= b; 
      WI(4)= 0.054975871827661;
      KSI(5)= 1. - 2. * b;  
      ETA(5)= b; 
      WI(5)= 0.054975871827661;
      KSI(6)= b;  
      ETA(6)= 1. - 2. * b; 
      WI(6)= 0.054975871827661;
     }
    else  if (nbi==7)
     {double a = (6.+sqrt(15.))/21.;
      double b = 4./7. -a;
      KSI(1)= 1./3.; 
      ETA(1)= 1./3. ;
      WI(1)= 9./80.;
      double A = (155. + sqrt(15.))/2400.; 
      KSI(2)= a; 
      ETA(2)= a ;
      WI(2)= A; 
      KSI(3)= 1.-2.* a ; 
      ETA(3)= a ;
      WI(3)= A;
      KSI(4)= a;  
      ETA(4)= 1.-2.* a; 
      WI(4)= A;
      A = 31./240.-A; 
      KSI(5)= b;  
      ETA(5)= b; 
      WI(5)= A;
      KSI(6)= 1. - 2. * b;  
      ETA(6)= b; 
      WI(6)= A;
      KSI(7)= b;  
      ETA(7)= 1. - 2. * b; 
      WI(7)= A;
     }
    else  if (nbi==12)
     {double a = 0.063089014491502;
      double b = 0.249286745170910;
      double c = 0.310352451033785;
      double d = 0.053145049844816;
      KSI(1)= a; 
      ETA(1)= a ;
      WI(1)= 0.025422453185103; 
      KSI(2)= 1.-2.* a ; 
      ETA(2)= a ;
      WI(2)= 0.025422453185103;
      KSI(3)= a;  
      ETA(3)= 1.-2.* a; 
      WI(3)= 0.025422453185103; 
      KSI(4)= b;  
      ETA(4)= b; 
      WI(4)= 0.058393137863189;
      KSI(5)= 1. - 2. * b;  
      ETA(5)= b; 
      WI(5)= 0.058393137863189;
      KSI(6)= b;  
      ETA(6)= 1. - 2. * b; 
      WI(6)= 0.058393137863189;
      double A = 0.041425537809187;
      KSI(7)= c;  
      ETA(7)= d; 
      WI(7)= A;
      KSI(8)= d;  
      ETA(8)= c; 
      WI(8)= A;
      KSI(9)= 1.-(c+d);  
      ETA(9)= c; 
      WI(9)= A;
      KSI(10)= 1.-(c+d);  
      ETA(10)= d; 
      WI(10)= A;
      KSI(11)= c;  
      ETA(11)= 1.-(c+d); 
      WI(11)= A;
      KSI(12)= d;  
      ETA(12)= 1.-(c+d); 
      WI(12)= A;
     }
    else
     {
      cout << "\n erreur, le nombre de point d\'integration demande n\'est pas implante ";
      cout << " nbi = " << nbi;
      cout << "\n GeomTriangle::GeomTriangle( int nbi, int nbne) " << endl;
      Sortie (1);
     };
    //--------------------------------
    //def des arretes
    //--------------------------------
    //def des aretes
    int nbil;  // nb de pt d'integ par ligne
    if (Nbi() == 1) 
          nbil = 1;
    else
          nbil = 2; 
    int nbnel;   // nb de noeud du segment
    if (NBNE == 3) nbnel = 2; else nbnel = 3;        
    seg(1) = new GeomSeg(nbil,nbnel);
    for (int il=2;il<= NBSE; il++)  // ici NBSE = 3
         seg(il) = seg(1);
    // def des tableaux de connection des noeuds  des cotes
    if (NBNE == 3)
         { for (int i =1;i<=NBSE;i++) NONS(i).Change_taille(2);
           NONS(1)(1) = 1;NONS(1)(2) = 2;NONS(2)(1) = 2;NONS(2)(2) = 3;
           NONS(3)(1) = 3;NONS(3)(2) = 1;
         }
    else 
         { for (int i =1;i<=NBSE;i++) NONS(i).Change_taille(3);
           NONS(1)(1) = 1;NONS(1)(2) = 4;NONS(1)(3) = 2;
           NONS(2)(1) = 2;NONS(2)(2) = 5;NONS(2)(3) = 3;
           NONS(3)(1) = 3;NONS(3)(2) = 6;NONS(3)(3) = 1;
         };
               
    //----------------------------------------------------------------------------------
    // définition des coordonnées des noeuds, de la triangulation linéaire de la face
    //  et valeurs aux points d'integration des fonctions et dérivées d'interpolation
    //---------------------------------------------------------------------------------- 
    //-----------------------------------------------------
    //------------------- cas lineaire --------------------
    //-----------------------------------------------------
    if (NBNE== 3)
     {id_interpol=LINEAIRE;
      // coordonnées des noeuds
      ptelem(1) = Coordonnee(0.,0.);
      ptelem(2) = Coordonnee(1.,0.);
      ptelem(3) = Coordonnee(0.,1.);
      // triangulation de la face
      NONFt(1).Change_taille(1); // une seule face pour la face
         // trois noeuds pour la face identique à ceux de l'élément
      NONFt(1)(1) = NONF(1);
      // puis
      for (int NI = 1; NI<= Nbi(); NI++)
         {
          //------------------------------
          //des fonctions d'interpolations
          //------------------------------
          tabPhi(NI)(1) = 1.  - KSI(NI) - ETA(NI);
          tabPhi(NI)(2) = KSI(NI);
          tabPhi(NI)(3) = ETA(NI);
          //-----------------
          //de leurs derivees
          //-----------------
          tabDPhi(NI)(1,1) = -1.; 
          tabDPhi(NI)(2,1) = -1.; 
          tabDPhi(NI)(1,2) = 1.; 
          tabDPhi(NI)(2,2) = 0.; 
          tabDPhi(NI)(1,3) = 0.; 
          tabDPhi(NI)(2,3) = 1.;
          };
     }

    //---------------------------------------------------------
    //--------------------- cas quadratique -------------------
    //---------------------------------------------------------
    else if (NBNE== 6)
     {id_interpol=QUADRACOMPL;
      // coordonnées des noeuds
      ptelem(1) = Coordonnee(0.,0.);
      ptelem(2) = Coordonnee(1.,0.);
      ptelem(3) = Coordonnee(0.,1.);
      ptelem(4) = Coordonnee(0.5,0.);
      ptelem(5) = Coordonnee(0.5,0.5);
      ptelem(6) = Coordonnee(0.,0.5);
      // triangulation de la face
      NONFt(1).Change_taille(4); // 4 face linéaires pour la face quadratique
         // trois noeuds pour chaque triangle linéaire
      NONFt(1)(1).Change_taille(3);NONFt(1)(2).Change_taille(3);
      NONFt(1)(3).Change_taille(3);NONFt(1)(4).Change_taille(3);
         // remplissage de la connexion par rapport à celle de l'élément
      NONFt(1)(1)(1) = 6; NONFt(1)(1)(2) = 4; NONFt(1)(1)(3) = 5; 
      NONFt(1)(2)(1) = 1; NONFt(1)(2)(2) = 4; NONFt(1)(2)(3) = 6; 
      NONFt(1)(3)(1) = 4; NONFt(1)(3)(2) = 2; NONFt(1)(3)(3) = 5; 
      NONFt(1)(4)(1) = 6; NONFt(1)(4)(2) = 5; NONFt(1)(4)(3) = 3; 
      // puis
      for (int NI = 1; NI<= Nbi(); NI++)
         {
          lambda= 1.  - KSI(NI) - ETA(NI);
         //------------------------------
         //des fonctions d'interpolations
         //------------------------------
          tabPhi(NI)(1) = -lambda * (1. -2. * lambda);
          tabPhi(NI)(2) = -KSI(NI) * (1. -2. * KSI(NI));
          tabPhi(NI)(3) = -ETA(NI) * (1. -2. * ETA(NI));
          tabPhi(NI)(4) = 4. * KSI(NI) * lambda;
          tabPhi(NI)(5) = 4. * KSI(NI) * ETA(NI);
          tabPhi(NI)(6) = 4. * ETA(NI) * lambda;
         //-----------------
         //de leurs derivees
         //-----------------
          tabDPhi(NI)(1,1) = 1. -4. * lambda;
          tabDPhi(NI)(2,1) = 1. -4. * lambda;
          tabDPhi(NI)(1,2) = -1. +4. * KSI(NI);
          tabDPhi(NI)(2,2) = 0.; 
          tabDPhi(NI)(1,3) = 0.; 
          tabDPhi(NI)(2,3) = -1. +4. * ETA(NI);
          tabDPhi(NI)(1,4) = 4. * (lambda-KSI(NI));
          tabDPhi(NI)(2,4) = -4. * KSI(NI);
          tabDPhi(NI)(1,5) = 4. * ETA(NI);
          tabDPhi(NI)(2,5) = 4. * KSI(NI);
          tabDPhi(NI)(1,6) = -4. * ETA(NI);
          tabDPhi(NI)(2,6) = 4. * (lambda-ETA(NI));
         };
     }
    //---------------------------------------------------------
    //--------------------- cas cubique -------------------
    //---------------------------------------------------------
    else if (NBNE== 9)
     {id_interpol=CUBIQUE_INCOMPL;
        // coordonnées des noeuds
        ptelem(1) = Coordonnee(0.,0.);
        ptelem(2) = Coordonnee(1.,0.);
        ptelem(3) = Coordonnee(0.,1.);
        double untiers= 1./3.; double deuxtiers = 2./3.;
        ptelem(4) = Coordonnee(untiers,0.);
        ptelem(5) = Coordonnee(deuxtiers,0.);
        ptelem(6) = Coordonnee(deuxtiers,untiers);
        ptelem(7) = Coordonnee(untiers,deuxtiers);
        ptelem(8) = Coordonnee(deuxtiers,0);
        ptelem(9) = Coordonnee(untiers,0.);
        // triangulation de la face
        NONFt(1).Change_taille(7); // 7 face linéaires
           // trois noeuds pour chaque triangle linéaire
        for (int i=1;i<= 7; i++ )  
           NONFt(1)(i).Change_taille(3); 
           // remplissage de la connexion par rapport à celle de l'élément
        NONFt(1)(1)(1) = 1;  NONFt(1)(1)(2) = 4;  NONFt(1)(1)(3) = 9; 
        NONFt(1)(2)(1) = 4;  NONFt(1)(2)(2) = 5;  NONFt(1)(2)(3) = 6; 
        NONFt(1)(3)(1) = 5;  NONFt(1)(3)(2) = 2;  NONFt(1)(3)(3) = 6; 
        NONFt(1)(4)(1) = 4;  NONFt(1)(4)(2) = 6;  NONFt(1)(4)(3) = 7; 
        NONFt(1)(5)(1) = 4;  NONFt(1)(5)(2) = 7;  NONFt(1)(5)(3) = 9; 
        NONFt(1)(6)(1) = 9;  NONFt(1)(6)(2) = 7;  NONFt(1)(6)(3) = 8; 
        NONFt(1)(7)(1) = 8;  NONFt(1)(7)(2) = 7;  NONFt(1)(7)(3) = 3; 
        // puis
        for (int NI = 1; NI<= Nbi(); NI++)
         {double ksi = KSI(NI); double eta = ETA(NI);
          double lambda= 1.  - ksi - eta;
          //------------------------------
          //des fonctions d'interpolations
          //------------------------------
          double troislambda= 3.*lambda;
          double moinsUnEtTroisLambda = -1.+ troislambda;
          double moinsDeuxEtTroisLambda = -2.+ troislambda;
          double troisKsi= 3.* ksi;
          double moinsUnEtTroisKsi = -1. + troisKsi;
          double moinsDeuxEtTroisKsi = -2. + troisKsi;
          double troisEta= 3.* eta;
          double moinsUnEtTroisEta = -1. + troisEta;
          double moinsDeuxEtTroisEta = -2. + troisEta;
         
          double N_10 = 0.5  * lambda * eta * ksi;
			 tabPhi(NI)(1) = 0.5 * lambda * moinsUnEtTroisLambda * moinsDeuxEtTroisLambda - N_10/6.;
          tabPhi(NI)(2) = 0.5 * ksi * moinsUnEtTroisKsi * moinsDeuxEtTroisKsi - N_10/6.;
          tabPhi(NI)(3) = 0.5 * eta * moinsUnEtTroisEta * moinsDeuxEtTroisEta - N_10/6.;
          tabPhi(NI)(4) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisLambda + N_10/4.;
          tabPhi(NI)(5) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisKsi + N_10/4.;
          tabPhi(NI)(6) = 0.5 * 9. * ksi * eta * moinsUnEtTroisKsi + N_10/4.;
          tabPhi(NI)(7) = 0.5 * 9. * ksi * eta * moinsUnEtTroisEta + N_10/4.;
          tabPhi(NI)(8) = 0.5 * 9. * lambda * eta * moinsUnEtTroisEta + N_10/4.;
          tabPhi(NI)(9) = 0.5 * 9. * lambda * eta * moinsUnEtTroisLambda + N_10/4.;

         //-----------------
         //de leurs derivees
         //-----------------
          double lambda2=lambda * lambda;
          double ksi2 = ksi * ksi;
          double eta2 = eta * eta;

          double d_N_10_1 = 27. * eta * (lambda - ksi);
          double d_N_10_2 = 27. * ksi * (lambda - eta);

          tabDPhi(NI)(1,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2) - d_N_10_1/6.;
          tabDPhi(NI)(2,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2) - d_N_10_2/6.;

          tabDPhi(NI)(1,2) = 0.5 * (2. - 18. * ksi + 27. * ksi2) - d_N_10_1/6.;
          tabDPhi(NI)(2,2) =  - d_N_10_2/6.; 

          tabDPhi(NI)(1,3) =  - d_N_10_1/6.; 
          tabDPhi(NI)(2,3) = 0.5 * (2. - 18. * eta + 27. * eta2) - d_N_10_2/6.;

          tabDPhi(NI)(1,4) = 4.5 * lambda*(moinsUnEtTroisLambda-6.*ksi) + 4.5 * ksi + d_N_10_1/4.;
          tabDPhi(NI)(2,4) = -4.5 * ksi * (-1. + 6. * lambda) + d_N_10_2/4.;
          
          tabDPhi(NI)(1,5) = 4.5 * ksi*(1.+6.*lambda-3.*ksi) - 4.5 * lambda + d_N_10_1/4.;
          tabDPhi(NI)(2,5) = -4.5 * ksi * (-1. + 3. * ksi) + d_N_10_2/4.;
          
          tabDPhi(NI)(1,6) = 4.5 * eta * (-1. + 6. * ksi) + d_N_10_1/4.;
          tabDPhi(NI)(2,6) = 4.5 * ksi * (-1. + 3. * ksi) + d_N_10_2/4.;

          tabDPhi(NI)(1,7) = 4.5 * eta * (-1. + 3. * eta) + d_N_10_1/4.;
          tabDPhi(NI)(2,7) = 4.5 * ksi * (-1. + 6. * eta) + d_N_10_2/4.;

          tabDPhi(NI)(1,8) = -4.5 * eta * (-1. + 3. * eta) + d_N_10_1/4.;
          tabDPhi(NI)(2,8) = 4.5 * eta * (1. + 6. * lambda  - 3. * eta) - 4.5 * lambda + d_N_10_2/4.;

          tabDPhi(NI)(1,9) = -4.5 * eta * (-1. + 6. * lambda) + d_N_10_1/4.;
          tabDPhi(NI)(2,9) = 4.5 * lambda * (-1. + 3. * lambda - 6. * eta) + 4.5 * eta + d_N_10_2/4.;

         }
        }  
    else if (NBNE== 10)
     {id_interpol=CUBIQUE;
        // coordonnées des noeuds
        ptelem(1) = Coordonnee(0.,0.);
        ptelem(2) = Coordonnee(1.,0.);
        ptelem(3) = Coordonnee(0.,1.);
        double untiers= 1./3.; double deuxtiers = 2./3.;
        ptelem(4) = Coordonnee(untiers,0.);
        ptelem(5) = Coordonnee(deuxtiers,0.);
        ptelem(6) = Coordonnee(deuxtiers,untiers);
        ptelem(7) = Coordonnee(untiers,deuxtiers);
        ptelem(8) = Coordonnee(deuxtiers,0);
        ptelem(9) = Coordonnee(untiers,0.);
        ptelem(10) = Coordonnee(untiers,untiers);
        // triangulation de la face
        NONFt(1).Change_taille(9); // 4 face linéaires pour la face quadratique
           // trois noeuds pour chaque triangle linéaire
        for (int i=1;i<= 9; i++ )  
           NONFt(1)(i).Change_taille(3); 
           // remplissage de la connexion par rapport à celle de l'élément
        NONFt(1)(1)(1) = 1;  NONFt(1)(1)(2) = 4;  NONFt(1)(1)(3) = 9; 
        NONFt(1)(2)(1) = 10; NONFt(1)(2)(2) = 9;  NONFt(1)(2)(3) = 4; 
        NONFt(1)(3)(1) = 4;  NONFt(1)(3)(2) = 5;  NONFt(1)(3)(3) = 10; 
        NONFt(1)(4)(1) = 6;  NONFt(1)(4)(2) = 10; NONFt(1)(4)(3) = 5; 
        NONFt(1)(5)(1) = 5;  NONFt(1)(5)(2) = 2;  NONFt(1)(5)(3) = 6; 
        NONFt(1)(6)(1) = 9;  NONFt(1)(6)(2) = 10; NONFt(1)(6)(3) = 8; 
        NONFt(1)(7)(1) = 7;  NONFt(1)(7)(2) = 8;  NONFt(1)(7)(3) = 10; 
        NONFt(1)(8)(1) = 10; NONFt(1)(8)(2) = 6;  NONFt(1)(8)(3) = 7; 
        NONFt(1)(9)(1) = 8;  NONFt(1)(9)(2) = 7;  NONFt(1)(9)(3) = 3; 
        // puis
        for (int NI = 1; NI<= Nbi(); NI++)
         {double ksi = KSI(NI); double eta = ETA(NI);
          double lambda= 1.  - ksi - eta;
          //------------------------------
          //des fonctions d'interpolations
          //------------------------------
          double troislambda= 3.*lambda;
          double moinsUnEtTroisLambda = -1.+ troislambda;
          double moinsDeuxEtTroisLambda = -2.+ troislambda;
          double troisKsi= 3.* ksi;
          double moinsUnEtTroisKsi = -1. + troisKsi;
          double moinsDeuxEtTroisKsi = -2. + troisKsi;
          double troisEta= 3.* eta;
          double moinsUnEtTroisEta = -1. + troisEta;
          double moinsDeuxEtTroisEta = -2. + troisEta;
         
          tabPhi(NI)(1) = 0.5 * lambda * moinsUnEtTroisLambda * moinsDeuxEtTroisLambda;
          tabPhi(NI)(2) = 0.5 * ksi * moinsUnEtTroisKsi * moinsDeuxEtTroisKsi;
          tabPhi(NI)(3) = 0.5 * eta * moinsUnEtTroisEta * moinsDeuxEtTroisEta;
          tabPhi(NI)(4) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisLambda;
          tabPhi(NI)(5) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisKsi;
          tabPhi(NI)(6) = 0.5 * 9. * ksi * eta * moinsUnEtTroisKsi;
          tabPhi(NI)(7) = 0.5 * 9. * ksi * eta * moinsUnEtTroisEta;
          tabPhi(NI)(8) = 0.5 * 9. * lambda * eta * moinsUnEtTroisEta;
          tabPhi(NI)(9) = 0.5 * 9. * lambda * eta * moinsUnEtTroisLambda;
          tabPhi(NI)(10) = 0.5  * lambda * eta * ksi;

         //-----------------
         //de leurs derivees
         //-----------------
          double lambda2=lambda * lambda;
          double ksi2 = ksi * ksi;
          double eta2 = eta * eta;
          
          tabDPhi(NI)(1,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2);
          tabDPhi(NI)(2,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2);

          tabDPhi(NI)(1,2) = 0.5 * (2. - 18. * ksi + 27. * ksi2);
          tabDPhi(NI)(2,2) = 0.; 

          tabDPhi(NI)(1,3) = 0.; 
          tabDPhi(NI)(2,3) = 0.5 * (2. - 18. * eta + 27. * eta2);

          tabDPhi(NI)(1,4) = 4.5 * lambda*(moinsUnEtTroisLambda-6.*ksi) + 4.5 * ksi;
          tabDPhi(NI)(2,4) = -4.5 * ksi * (-1. + 6. * lambda);
          
          tabDPhi(NI)(1,5) = 4.5 * ksi*(1.+6.*lambda-3.*ksi) - 4.5 * lambda;
          tabDPhi(NI)(2,5) = -4.5 * ksi * (-1. + 3. * ksi);
          
          tabDPhi(NI)(1,6) = 4.5 * eta * (-1. + 6. * ksi);
          tabDPhi(NI)(2,6) = 4.5 * ksi * (-1. + 3. * ksi);

          tabDPhi(NI)(1,7) = 4.5 * eta * (-1. + 3. * eta);
          tabDPhi(NI)(2,7) = 4.5 * ksi * (-1. + 6. * eta);

          tabDPhi(NI)(1,8) = -4.5 * eta * (-1. + 3. * eta);
          tabDPhi(NI)(2,8) = 4.5 * eta * (1. + 6. * lambda  - 3. * eta) - 4.5 * lambda;

          tabDPhi(NI)(1,9) = -4.5 * eta * (-1. + 6. * lambda);
          tabDPhi(NI)(2,9) = 4.5 * lambda * (-1. + 3. * lambda - 6. * eta) + 4.5 * eta;

          tabDPhi(NI)(1,10) = 27. * eta * (lambda - ksi);
          tabDPhi(NI)(2,10) = 27. * ksi * (lambda - eta);

         }
        }  
    else 
      {
       cout << "\n erreur, le nombre de noeud demande n\'est pas implante ";
       cout << " nbne = " << NBNE;
       cout << "\n GeomTriangle::GeomTriangle( int nbi, int nbne) " << endl;
       Sortie (1);
      };

    // ---- constitution du tableau Extrapol ----- 
    Calcul_extrapol(nbi);      
  };
// destructeur
GeomTriangle::~GeomTriangle() 
  {   delete seg(1);
  };

// constructeur de copie
GeomTriangle::GeomTriangle(const GeomTriangle& a) :
  ElemGeomC0(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 de la face
     face(1) = this;
     // def des segments
     seg(1) = new GeomSeg(*((GeomSeg*)(a.seg(1)))) ;
     for (int il=2;il<= 3; 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 * GeomTriangle::newElemGeomC0(ElemGeomC0 * pt) 
  { pt = new GeomTriangle(*this);
    return pt; 
  };

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

    // retourne les fonctions d'interpolation au point M (en coordonnees locales)
const Vecteur& GeomTriangle::Phi_point(const Coordonnee& M)
  {
    #ifdef MISE_AU_POINT	 	   
     // verification de la dimension des coordonnees locales
     if (M.Dimension() != 2)
       { cout << "\n erreur la dimension des coordonnees locales :" << M.Dimension()
              <<"n\'est pas egale a 2 "
             << "\nGeomTriangle::Phi(Coordonnee& M)";
        Sortie(1);
      }
    #endif
//    Vecteur phi(NBNE); // tableau des fonctions d'interpolation
    phi_M.Change_taille(NBNE); // si la taille est identique -> aucune action
    // fonction  des fonctions d'interpolation
    if (NBNE== 3) 
         {
          //------------------------------
          //des fonctions d'interpolations
          //------------------------------
          phi_M(1) = 1.  - M(1) - M(2);
          phi_M(2) = M(1);
          phi_M(3) = M(2);
         }  
    else if (NBNE== 6)
       {
          double lambda= 1.  - M(1) - M(2);
         //------------------------------
         //des fonctions d'interpolations
         //------------------------------
          phi_M(1) = -lambda * (1. -2. * lambda);
          phi_M(2) = -M(1) * (1. -2. * M(1));
          phi_M(3) = -M(2) * (1. -2. * M(2));
          phi_M(4) = 4. * M(1) * lambda;
          phi_M(5) = 4. * M(1) * M(2);
          phi_M(6) = 4. * M(2) * lambda;
         } 
      else if (NBNE== 9)
       {  double ksi = M(1); double eta = M(2);
          double lambda= 1.  - ksi - eta;
          //------------------------------
          //des fonctions d'interpolations
          //------------------------------
          double troislambda= 3.*lambda;
          double moinsUnEtTroisLambda = -1.+ troislambda;
          double moinsDeuxEtTroisLambda = -2.+ troislambda;
          double troisKsi= 3.* ksi;
          double moinsUnEtTroisKsi = -1. + troisKsi;
          double moinsDeuxEtTroisKsi = -2. + troisKsi;
          double troisEta= 3.* eta;
          double moinsUnEtTroisEta = -1. + troisEta;
          double moinsDeuxEtTroisEta = -2. + troisEta;
         
          double N_10 = 0.5  * lambda * eta * ksi;
          phi_M(1) = 0.5 * lambda * moinsUnEtTroisLambda * moinsDeuxEtTroisLambda - N_10/6.;
          phi_M(2) = 0.5 * ksi * moinsUnEtTroisKsi * moinsDeuxEtTroisKsi - N_10/6.;
          phi_M(3) = 0.5 * eta * moinsUnEtTroisEta * moinsDeuxEtTroisEta - N_10/6.;
          phi_M(4) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisLambda + N_10/4.;
          phi_M(5) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisKsi + N_10/4.;
          phi_M(6) = 0.5 * 9. * ksi * eta * moinsUnEtTroisKsi + N_10/4.;
          phi_M(7) = 0.5 * 9. * ksi * eta * moinsUnEtTroisEta + N_10/4.;
          phi_M(8) = 0.5 * 9. * lambda * eta * moinsUnEtTroisEta + N_10/4.;
          phi_M(9) = 0.5 * 9. * lambda * eta * moinsUnEtTroisLambda + N_10/4.;

        }  
      else if (NBNE== 10)
       {  double ksi = M(1); double eta = M(2);
          double lambda= 1.  - ksi - eta;
          //------------------------------
          //des fonctions d'interpolations
          //------------------------------
          double troislambda= 3.*lambda;
          double moinsUnEtTroisLambda = -1.+ troislambda;
          double moinsDeuxEtTroisLambda = -2.+ troislambda;
          double troisKsi= 3.* ksi;
          double moinsUnEtTroisKsi = -1. + troisKsi;
          double moinsDeuxEtTroisKsi = -2. + troisKsi;
          double troisEta= 3.* eta;
          double moinsUnEtTroisEta = -1. + troisEta;
          double moinsDeuxEtTroisEta = -2. + troisEta;
         
          phi_M(1) = 0.5 * lambda * moinsUnEtTroisLambda * moinsDeuxEtTroisLambda;
          phi_M(2) = 0.5 * ksi * moinsUnEtTroisKsi * moinsDeuxEtTroisKsi;
          phi_M(3) = 0.5 * eta * moinsUnEtTroisEta * moinsDeuxEtTroisEta;
          phi_M(4) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisLambda;
          phi_M(5) = 0.5 * 9. * lambda * ksi * moinsUnEtTroisKsi;
          phi_M(6) = 0.5 * 9. * ksi * eta * moinsUnEtTroisKsi;
          phi_M(7) = 0.5 * 9. * ksi * eta * moinsUnEtTroisEta;
          phi_M(8) = 0.5 * 9. * lambda * eta * moinsUnEtTroisEta;
          phi_M(9) = 0.5 * 9. * lambda * eta * moinsUnEtTroisLambda;
          phi_M(10) = 0.5  * lambda * eta * ksi;

        }  
    else 
       {
         cout << "\n erreur, le nombre de noeud demande :" << NBNE <<"n\'est pas implante ";
         cout << "\n GeomTriangle::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& GeomTriangle::Dphi_point(const Coordonnee& M)
  {  
   #ifdef MISE_AU_POINT	 	   
     // verification de la dimension des coordonnees locales
     if (M.Dimension() != 2)
       { cout << "\n erreur la dimension des coordonnees locales :" << M.Dimension()
              <<"n\'est pas egale a 2 "
             << "\nGeomTriangle::Dphi(Coordonnee& M)";
        Sortie(1);
      }
    #endif
//   Mat_pleine dphi(2,NBNE); // le tableau des derivees
   // le tableau des derivees: redimentionnement si nécessaire
   if ((dphi_M.Nb_ligne() != 2)&&(dphi_M.Nb_colonne() != NBNE))
      dphi_M.Initialise (2,NBNE,0.);
   // fonction et derivees des fonctions d'interpolation
    if (NBNE== 3) 
         {
          //-----------------
          //de leurs derivees
          //-----------------
          dphi_M(1,1) = -1.;
          dphi_M(2,1) = -1.;
          dphi_M(1,2) = 1.;
          dphi_M(2,2) = 0.;
          dphi_M(1,3) = 0.;
          dphi_M(2,3) = 1.;
         }  
    else if (NBNE== 6)
       {
          double lambda= 1.  - M(1) - M(2);
         //-----------------
         //de leurs derivees
         //-----------------
          dphi_M(1,1) = 1. -4. * lambda;
          dphi_M(2,1) = 1. -4. * lambda;
          dphi_M(1,2) = -1. +4. * M(1);
          dphi_M(2,2) = 0.;
          dphi_M(1,3) = 0.;
          dphi_M(2,3) = -1. +4. * M(2);
          dphi_M(1,4) = 4. * (lambda-M(1));
          dphi_M(2,4) = -4. * M(1);
          dphi_M(1,5) = 4. * M(2);
          dphi_M(2,5) = 4. * M(1);
          dphi_M(1,6) = -4. * M(2);
          dphi_M(2,6) = 4. * (lambda-M(2));
         } 

      else if (NBNE== 9)
       {  double ksi = M(1); double eta = M(2);
          double lambda= 1.  - ksi - eta;

         //-----------------
          double lambda2=lambda * lambda;
          double ksi2 = ksi * ksi;
          double eta2 = eta * eta;
          
          double d_N_10_1 = 27. * eta * (lambda - ksi);
          double d_N_10_2 = 27. * ksi * (lambda - eta);

          dphi_M(1,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2) - d_N_10_1/6.;
          dphi_M(2,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2) - d_N_10_2/6.;

          dphi_M(1,2) = 0.5 * (2. - 18. * ksi + 27. * ksi2) - d_N_10_1/6.;
          dphi_M(2,2) =  - d_N_10_2/6.;

          dphi_M(1,3) =  - d_N_10_1/6.;
          dphi_M(2,3) = 0.5 * (2. - 18. * eta + 27. * eta2) - d_N_10_2/6.;

          dphi_M(1,4) = 4.5 * lambda*(-1.+3.*lambda-6.*ksi) + 4.5 * ksi + d_N_10_1/4.;
          dphi_M(2,4) = -4.5 * ksi * (-1. + 6. * lambda) + d_N_10_2/4.;
          
          dphi_M(1,5) = 4.5 * ksi*(1.+6.*lambda-3.*ksi) - 4.5 * lambda + d_N_10_1/4.;
          dphi_M(2,5) = -4.5 * ksi * (-1. + 3. * ksi) + d_N_10_2/4.;
          
          dphi_M(1,6) = 4.5 * eta * (-1. + 6. * ksi) + d_N_10_1/4.;
          dphi_M(2,6) = 4.5 * ksi * (-1. + 3. * ksi) + d_N_10_2/4.;

          dphi_M(1,7) = 4.5 * eta * (-1. + 3. * eta) + d_N_10_1/4.;
          dphi_M(2,7) = 4.5 * ksi * (-1. + 6. * eta) + d_N_10_2/4.;

          dphi_M(1,8) = -4.5 * eta * (-1. + 3. * eta) + d_N_10_1/4.;
          dphi_M(2,8) = 4.5 * eta * (1. + 6. * lambda  - 3. * eta) - 4.5 * lambda + d_N_10_2/4.;

          dphi_M(1,9) = -4.5 * eta * (-1. + 6. * lambda) + d_N_10_1/4.;
          dphi_M(2,9) = 4.5 * lambda * (-1. + 3. * lambda - 6. * eta) + 4.5 * eta + d_N_10_2/4.;
        }  
      else if (NBNE== 10)
       {  double ksi = M(1); double eta = M(2);
          double lambda= 1.  - ksi - eta;

         //-----------------
          double lambda2=lambda * lambda;
          double ksi2 = ksi * ksi;
          double eta2 = eta * eta;
          
          dphi_M(1,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2);
          dphi_M(2,1) = 0.5 * (-2. + 18.*lambda - 27.*lambda2);

          dphi_M(1,2) = 0.5 * (2. - 18. * ksi + 27. * ksi2);
          dphi_M(2,2) = 0.;

          dphi_M(1,3) = 0.;
          dphi_M(2,3) = 0.5 * (2. - 18. * eta + 27. * eta2);

          dphi_M(1,4) = 4.5 * lambda*(-1.+3.*lambda-6.*ksi) + 4.5 * ksi;
          dphi_M(2,4) = -4.5 * ksi * (-1. + 6. * lambda);
          
          dphi_M(1,5) = 4.5 * ksi*(1.+6.*lambda-3.*ksi) - 4.5 * lambda;
          dphi_M(2,5) = -4.5 * ksi * (-1. + 3. * ksi);
          
          dphi_M(1,6) = 4.5 * eta * (-1. + 6. * ksi);
          dphi_M(2,6) = 4.5 * ksi * (-1. + 3. * ksi);

          dphi_M(1,7) = 4.5 * eta * (-1. + 3. * eta);
          dphi_M(2,7) = 4.5 * ksi * (-1. + 6. * eta);

          dphi_M(1,8) = -4.5 * eta * (-1. + 3. * eta);
          dphi_M(2,8) = 4.5 * eta * (1. + 6. * lambda  - 3. * eta) - 4.5 * lambda;

          dphi_M(1,9) = -4.5 * eta * (-1. + 6. * lambda);
          dphi_M(2,9) = 4.5 * lambda * (-1. + 3. * lambda - 6. * eta) + 4.5 * eta;

          dphi_M(1,10) = 27 * eta * (lambda - ksi);
          dphi_M(2,10) = 27 * ksi * (lambda - eta);
        }  
    else 
      { cout << "\n erreur le nombre de noeud demande :" << NBNE <<"n\'est pas implante "
             << "\nGeomTriangle::Dphi(Coordonnee& M)";
        Sortie(1);
      }
    return dphi_M;
   };
          
// en fonction de coordonnees locales, retourne true si le point est a l'interieur
// de l'element, false sinon
bool GeomTriangle::Interieur(const Coordonnee& M)
  { if ((M(1) >= 0.) && (M(1) <= 1.) )
      if ((M(2) >= 0.) && (M(2) <= 1.-M(1)) )                
        return true;
      else 
        return false;  
    else
        return false;    
   };
              
// en fonction de coordonnees locales, retourne le point local P, maximum intérieur à l'élément, donc sur la frontière
// dont les coordonnées sont sur la droite GM: c-a-d GP = alpha GM, avec apha maxi et P appartenant à la frontière
// de l'élément, G étant le centre de gravité, sauf si GM est nul, dans ce cas retour de M
Coordonnee GeomTriangle::Maxi_Coor_dans_directionGM(const Coordonnee& M)
  { // le principe est de faire un traitement en fonction du secteur dans lequel se trouve le point
    double x=M(1); double y=M(2);
    // on appelle G le centre de gravité, A(1,0), B(0,1), C(0,0)
    double untiers = 1./3.;
    Coordonnee G(untiers,untiers);
    if ((x+2.*y-1. >=0.) && (2.*x + y -1 >=0.))
     { // cas où l'on est dans le secteur limité par les demies droites ]A,G] et ]B,G]
	      //1) on calcul la normale au plan: GAz->proportionnelle à (1,2,0), on doit être à l'avant de ce plan:-> x+2y-1 >=0 
		   //2) on calcul la normale au plan: GBz->proportionnelle à (2,1,0), on doit être à l'avant de ce plan-> 2x+y-1 >=0
       // on fait une rotation de 45 degré et le maxi des nouvelles coordonnées X Y est xmax*racine(2)/2
       double xmax =  MaX(Dabs(x+y-2./3.),Dabs(y-x)); 
       // ensuite sachant que le point que l'on cherche doit se trouver sur AB, c'est à dire pour
       // X = racine(2)/6, on doit donc avoir un coeff multiplicateur de 1/3/xmax (que l'on trouve
       // en ramenant de -45 degrés)
       if (xmax <= ConstMath::petit) 
        //si xmax est trop petit cela signifie que M est près de G,		 
        { return M;} 
       else
        { return (Coordonnee((untiers / xmax)* (M-G) + G));}
      }
     else if (((y-x) <= 0)&&(x+2.*y-1. <=0.))
      {// cas où l'on est dans le secteur limité par les demies droites ]C,G] et ]A,G]
       // on calcul le y-1/3 maxi et on limite à 1/3
       double xmax = untiers - y;
       if (xmax <= ConstMath::petit) 
        { return M;}
       else
        { return (Coordonnee((untiers / xmax)* (M-G) + G));}
       }
     else if (((y-x) >= 0)&&(2.*x + y -1 <=0.))
      {// cas où l'on est dans le secteur limité par les demies droites ]C,G] et ]B,G]
       // on calcul le x-1/3 maxi et on limite à 1/3
       double xmax = untiers - x;
       if (xmax <= ConstMath::petit) 
        { return M;}
       else
        { return (Coordonnee((untiers / xmax)* (M-G) + G));}
       }

      // cas qui ne doit jamais arriver car on se situe forcément dans un des trois secteurs
      // précédemment étudiés
      cout << "\n erreur inconnue !!"
           << "\n GeomTriangle::Maxi_Coor_dans_directionGM(..";
      Sortie(1);     
      return M; // pour éviter un warning
   }; 
  
// constitution du tableau Extrapol
void GeomTriangle::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)
  //  --- pour l'instant seul le cas 1 est implanté  ---
  Tableau<Tableau<int> > &  indir = extrapol(1).indir;  // pour simplifier
  Tableau<Tableau<double > > & tab = extrapol(1).tab;    // pour simplifier
  Tableau <int > indirect(3); 
  Tableau <Coordonnee>  gi_B,gi_H; 
  Vecteur phi_(3); // le conteneur pour les fonctions d'interpolation
  Coordonnee theta(2); // le conteneur pour les coordonnées locales

  switch  (nbi)
   { case 1: 
       { // cas avec un point d'intégration, quelque soit le nombre de noeuds, 
         // on reporte la valeur au pt d'integ, telle quelle au noeud
         for (int ne=1;ne<=NBNE;ne++)
         	{tab(ne)(1)=1.;
         	 indir(ne).Change_taille(1); indir(ne)(1)=1;
         	};
         break;
       }
   	 case 3:
   	   { // cas avec 3 points d'intégration (au milieu des arrêtes) 
         // on extrapole linéairement avec les valeurs aux pt d'integ, ceci idem pour tous les noeuds
         // --  extrapole en 2 3 1 
         indirect(1)=2;indirect(2)=3;indirect(3)=1; 
         Bases_naturel_duales(indirect,gi_B,gi_H);
         Coordonnee O(ptInteg(indirect(1))); // def de l'origine
         for (int ne=1;ne<=NBNE;ne++) 
          { 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);  
            indir(ne).Change_taille(3); 
            indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
           };
          break;  
   	 	}  // fin du cas avec 3 pt d'intégration
   	 case 1003:
   	   { // cas avec 3 points d'intégration chacun près d'un noeud sommet 
         // on extrapole linéairement 
         // --  extrapole en  1 2 3 
         indirect(1)=1;indirect(2)=2;indirect(3)=3; 
         Bases_naturel_duales(indirect,gi_B,gi_H);
         Coordonnee O(ptInteg(indirect(1))); // def de l'origine
         for (int ne=1;ne<=NBNE;ne++) 
          { 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);  
            indir(ne).Change_taille(3); 
            indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
           };
          break;  
   	 	}  // fin du cas avec 3 pt d'intégration
   	 case 4:
   	   { // cas avec 4 points d'intégration  
         // on extrapole linéairement vers les  noeuds en considérant à chaque fois les 3 pt d'intégration
         // extrèmes c-a-d 2,3,4 (on ne tiens pas compte du premier point qui est au centre: c'est un choix)
         //  tout d'abord on calcul la base duale associée au triangle constitué par ces 3 pt d'integ
         indirect(1)=2;indirect(2)=3;indirect(3)=4;
         Bases_naturel_duales(indirect,gi_B,gi_H);
         Coordonnee O(ptInteg(2)); // def de l'origine du triangle des pt d'integ
         // on boucle sur les noeuds
         for (int ne = 1; ne<= NBNE; ne++)
         	{ // calcul des fonctions d'interpolations au noeud
         	  ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	  // def du tableau tab pour ne
         	  tab(ne)(2) = phi_(1);
              tab(ne)(3) = phi_(2);
              tab(ne)(4) = phi_(3);
         	  indir(ne).Change_taille(3); 
              indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
         	};
          break;  
   	 	}  // fin du cas avec 4 pt d'intégration
   	 case 6:
   	   { // cas avec 6 points d'intégration  
         // on extrapole linéairement vers les  noeuds en considérant à chaque fois les 3 pt d'intégration
         // les plus près du noeud
         // -- def des conteneurs qui serviront pour chaque noeud
         // on différencie en fonction du nombre de noeud
         switch (NBNE)
         	{ case 3: // trois noeuds, 
         	     {// vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ
         	      int ne = 1;  // le pt d'integ le plus près = nb 4
         	      tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4; 
         	      ne = 2;  // le pt d'integ le plus près = nb 5
         	      tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
         	      ne = 3;  // le pt d'integ le plus près = nb 6
         	      tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
         		  break;
         	     }
         	  case 6:// comme pour le cas de 3 noeuds
         	    { // vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ
         	      int ne = 1;  // le pt d'integ le plus près = nb 4
         	      tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4;
         	      ne = 2;  // le pt d'integ le plus près = nb 5
         	      tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
         	      ne = 3;  // le pt d'integ le plus près = nb 6
         	      tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
         	      ne = 4;  // le pt d'integ le plus près = nb 3
         	      tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
         	      ne = 5;  // le pt d'integ le plus près = nb 1
         	      tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
         	      ne = 6;  // le pt d'integ le plus près = nb 2
         	      tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
         		  break;
         		}	
         	  case 10: // comme pour le cas de 3 noeuds
         	    { // vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ pour les 3 premiers noeuds ensuite on extrapole
         	      int ne = 1;  // le pt d'integ le plus près = nb 4
         	      tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4;
         	      ne = 2;  // le pt d'integ le plus près = nb 5
         	      tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
         	      ne = 3;  // le pt d'integ le plus près = nb 6
         	      tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
                  // -- maintenant on s'occupe des noeuds internes
                  //  - tout d'abord les noeuds sur les arrêtes
                  //    noeuds 4 5
                  indirect(1)=3;indirect(2)=4;indirect(3)=5;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  Coordonnee O(ptInteg(3)); // def de l'origine du triangle des pt d'integ
                  for (int ne = 4; ne<=5;ne++)
                  {ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(3) = phi_(1);tab(ne)(4) = phi_(2);tab(ne)(5) = phi_(3);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
                  //    noeuds 6 7
                  indirect(1)=1;indirect(2)=5;indirect(3)=6;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  O=ptInteg(1); // def de l'origine du triangle des pt d'integ
                  for (int ne = 6; ne<=7;ne++)
                  {ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(1) = phi_(1); tab(ne)(5) = phi_(2);tab(ne)(6) = phi_(3);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
                  //    noeuds 8 9
                  indirect(1)=2;indirect(2)=6;indirect(3)=4;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  O=ptInteg(2); // def de l'origine du triangle des pt d'integ
                  for (int ne = 8; ne<=9;ne++)
                  {ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(2) = phi_(1);
                   tab(ne)(4) = phi_(3);
                   tab(ne)(6) = phi_(2);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
                  //    noeud 10 : le noeud central
                  indirect(1)=1;indirect(2)=2;indirect(3)=3;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  O=ptInteg(1); // def de l'origine du triangle des pt d'integ
                  {int ne = 10;
                   ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(1) = phi_(1);
                   tab(ne)(2) = phi_(2);
                   tab(ne)(3) = phi_(3);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
         		  break;
         		}
         	  default:
         	    { cout << "\n erreur le nombre de noeud demande : " << NBNE <<" n\'est pas implante "
                       << "\nGeomTriangle::Calcul_extrapol(..";
                  Sortie(1);
                };
         	};
          break;  
   	 	}  // fin du cas avec 6 pt d'intégration
   	 case 7:
   	   { // cas avec 7 points d'intégration  
         // on extrapole linéairement vers les  noeuds en considérant à chaque fois les 3 pt d'intégration
         // les plus près du noeud
         // on différencie en fonction du nombre de noeud
         switch (NBNE)
         	{ case 3: // trois noeuds, 
         	     {// vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ
         	      int ne = 1;  // le pt d'integ le plus près = nb 5
         	      tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
         	      ne = 2;  // le pt d'integ le plus près = nb 6
         	      tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
         	      ne = 3;  // le pt d'integ le plus près = nb 7
         	      tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;
         		  break;
         	     }
         	  case 6:// comme pour le cas de 3 noeuds
         	    { // vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ
         	      int ne = 1;  // le pt d'integ le plus près = nb 5
         	      tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
         	      ne = 2;  // le pt d'integ le plus près = nb 6
         	      tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
         	      ne = 3;  // le pt d'integ le plus près = nb 7
         	      tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;
         	      ne = 4;  // le pt d'integ le plus près = nb 4
         	      tab(ne)(4) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=4;
         	      ne = 5;  // le pt d'integ le plus près = nb 2
         	      tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
         	      ne = 6;  // le pt d'integ le plus près = nb 3
         	      tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
         		  break;
         		}	
         	  case 10: // comme pour le cas de 3 noeuds
         	    { // vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ pour les 3 premiers noeuds ensuite on extrapole
         	      int ne = 1;  // le pt d'integ le plus près = nb 5
         	      tab(ne)(5) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=5;
         	      ne = 2;  // le pt d'integ le plus près = nb 6
         	      tab(ne)(6) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=6;
         	      ne = 3;  // le pt d'integ le plus près = nb 7
         	      tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;
                  // -- maintenant on s'occupe des noeuds internes
                  //  - tout d'abord les noeuds sur les arrêtes
                  //    noeuds 4 5
                  indirect(1)=4;indirect(2)=6;indirect(3)=5;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  Coordonnee O(ptInteg(4)); // def de l'origine du triangle des pt d'integ
                  for (int ne = 4; ne<=5;ne++)
                  {ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(4) = phi_(1);
                   tab(ne)(6) = phi_(2);
                   tab(ne)(5) = phi_(3);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
                  //    noeuds 6 7
                  indirect(1)=2;indirect(2)=7;indirect(3)=6;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  O=ptInteg(2); // def de l'origine du triangle des pt d'integ
                  for (int ne = 6; ne<=7;ne++)
                  {ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(2) = phi_(1);tab(ne)(6) = phi_(3);tab(ne)(7) = phi_(2);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
                  //    noeuds 8 9
                  indirect(1)=3;indirect(2)=5;indirect(3)=7;
                  Bases_naturel_duales(indirect,gi_B,gi_H);
                  O=ptInteg(3); // def de l'origine du triangle des pt d'integ
                  for (int ne = 8; ne<=9;ne++)
                  {ElemGeomC0::Coor_phi(O,gi_H,ptelem(ne),phi_,theta);
         	       tab(ne)(3) = phi_(1);
                   tab(ne)(5) = phi_(2);
                   tab(ne)(7) = phi_(3);
         	       indir(ne).Change_taille(3); 
                   indir(ne)(1)=indirect(1);indir(ne)(2)=indirect(2);indir(ne)(3)=indirect(3);
                  };
                  //    noeud 10 : le noeud central idem pt 7
         	      ne = 10;  
         	      tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;;
         		  break;
         	    }
         	  default:
         	    { cout << "\n erreur le nombre de noeud demande :" << NBNE <<"n\'est pas implante "
                       << "\nGeomTriangle::Calcul_extrapol(..";
                  Sortie(1);
                };
         	};
          break;  
   	 	}  // fin du cas avec 7 pt d'intégration
   	 case 12:
   	   { // cas avec 12 points d'intégration  
         // on extrapole linéairement vers les  noeuds en considérant à chaque fois les 3 pt d'intégration
         // les plus près du noeud
         // on différencie en fonction du nombre de noeud
         switch (NBNE)
         	{ case 3: // trois noeuds, 
         	     {// vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ
         	      int ne = 1;  // le pt d'integ le plus près = nb 1
         	      tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
         	      ne = 2;  // le pt d'integ le plus près = nb 2
         	      tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
         	      ne = 3;  // le pt d'integ le plus près = nb 3
         	      tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
         		  break;
         	     }
         	  case 6:// comme pour le cas de 3 noeuds
         	    { // vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ
         	      int ne = 1;  // le pt d'integ le plus près = nb 1
         	      tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
         	      ne = 2;  // le pt d'integ le plus près = nb 2
         	      tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
         	      ne = 3;  // le pt d'integ le plus près = nb 3
         	      tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
         	      ne = 4;  // on fait la moyenne des pt 7 et 10
         	      tab(ne)(7) = 0.5;tab(ne)(10) = 0.5;
         	      indir(ne).Change_taille(2); indir(ne)(1)=7;indir(ne)(2)=10;
         	      ne = 5;  // on fait la moyenne des pt 9 et 11
         	      tab(ne)(9) = 0.5;tab(ne)(11) = 0.5;
         	      indir(ne).Change_taille(2); indir(ne)(1)=9;indir(ne)(2)=11;
         	      ne = 6;  // on fait la moyenne des pt 12 et 8
         	      tab(ne)(8) = 0.5;tab(ne)(12) = 0.5;
         	      indir(ne).Change_taille(2); indir(ne)(1)=12;indir(ne)(2)=8;
         		  break;
         		}	
         	  case 10: // comme pour le cas de 3 noeuds
         	    { // vue la proximité de certain points d'intégration avec les noeuds on prend directement 
         	      // les valeurs à ces  pt d'integ pour les 3 premiers noeuds ensuite on extrapole
         	      int ne = 1;  // le pt d'integ le plus près = nb 1
         	      tab(ne)(1) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=1;
         	      ne = 2;  // le pt d'integ le plus près = nb 2
         	      tab(ne)(2) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=2;
         	      ne = 3;  // le pt d'integ le plus près = nb 3
         	      tab(ne)(3) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=3;
         	      // pour les arrêtes on a également des pt d'integ tous prés
          	      ne = 4;  // le pt d'integ le plus près = nb 7
         	      tab(ne)(7) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=7;
         	      ne = 5;  // le pt d'integ le plus près = nb 10
         	      tab(ne)(10) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=10;
         	      ne = 6;  // le pt d'integ le plus près = nb 9
         	      tab(ne)(9) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=9;
         	      ne = 7;  // le pt d'integ le plus près = nb 11
         	      tab(ne)(11) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=11;
         	      ne = 8;  // le pt d'integ le plus près = nb 12
         	      tab(ne)(12) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=12;
         	      ne = 9;  // le pt d'integ le plus près = nb 8
         	      tab(ne)(8) = 1.;indir(ne).Change_taille(1); indir(ne)(1)=8;
         	      // pour le noeud central on utilise les trois noeud 4 5 6 que l'on moyennes
         	      double untiers = 1./3.;
         	      ne = 10;
         	      tab(ne)(4)=untiers;tab(ne)(5)=untiers;tab(ne)(6)=untiers;
         	      indir(ne).Change_taille(3); indir(ne)(1)=4;indir(ne)(2)=5;indir(ne)(3)=6;
         		  break;
         	    }
         	  default:
         	    { cout << "\n erreur le nombre de noeud demande :" << NBNE <<"n\'est pas implante "
                       << "\nGeomTriangle::Calcul_extrapol(..";
                  Sortie(1);
                };
         	};
          break;  
   	 	}  // fin du cas avec 12 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);
       };
   };

}; 
        
  
  
  
  
  
