// FICHIER : Met_Sfe1.cp
// CLASSE : Met_Sfe1


#include <iostream>
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "Util.h"
#include "MathUtil.h"


#include "Met_Sfe1.h"


//  =========================================================================================
//  vu la taille des executables le fichier est decompose en deux
// le premier : Met_Sfe1s1.cp concerne les constructeurs, destructeur et 
//              la gestion des variables
// le second :  Met_Sfe1s2.cp concerne  le calcul
//  =========================================================================================

// METHODES PUBLIQUES :
//  ------------------------ calculs ----------------------------------------
	    // cas explicite, toutes les grandeurs sont a 0 ou t	
/*  Met_Sfe1::Expli  Met_Sfe1::Cal_explicit
                ( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi,int nombre_noeud)
	{ Met_abstraite::Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
	  Calcul_giB_t(tab_noeud,dphi,nombre_noeud); // calcul de la base a t
	  Calcul_gijBB_0 ();  // metrique base naturelle
	  Calcul_gijBB_t ();  //      "
	  Calcul_gijHH_t ();        // composantes controvariantes
	  Jacobien_t();         // calcul du jacobien a t
	  
	  D_giB_t(dphi);   //
	  D_gijBB_t();    //variation de la  metrique en BB
	  // retour des infos
      Expli ex;
      ex.gijBB_0 = gijBB_0;
      ex.gijBB_t = gijBB_t;
      ex.gijHH_t = gijHH_t;
      ex.d_gijBB_t = &d_gijBB_t;
      ex.jacobien = jacobien_t;
      // liberation des tenseurs intermediaires
      LibereTenseur();   
      return ex;	  
	};   */
	
// cas implicite, toutes les grandeurs sont a 0 ou t+dt	
/*Met_abstraite::Impli Met_Sfe1::Cal_implicit
              ( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
	{ // tout d'abord on calcul les elements relatifs a la facette plane
	 Met_abstraite::Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
	 *aiB_0 = *giB_0;
	 Met_abstraite::Calcul_giB_t(tab_noeud,dphi,nombre_noeud); // calcul de la base a t
	 *aiB_t = *giB_t;
	 Met_abstraite::Calcul_giB_tdt(tab_noeud,dphi,nombre_noeud); // calcul de la base a tdt
	 *aiB_tdt = *giB_tdt;
	 Met_abstraite::Calcul_gijBB_0 ();    // metrique base naturelle
	 Met_abstraite::Calcul_gijBB_t ();    //      "
	 Met_abstraite::Calcul_gijBB_tdt ();  //      "
	 Met_abstraite::Calcul_gijHH_tdt ();        // composantes controvariantes
	 *aijBB_0 = *gijBB_0;
	 *aijBB_t = *gijBB_t;
	 *aijBB_tdt = *gijBB_tdt;
	 *aijHH_tdt = *gijHH_tdt;
	 Met_abstraite::Jacobien_tdt();         // calcul du jacobien a t
	 ajacobien_tdt = jacobien_tdt;
	 */
	  
/*	 D_giB_tdt(dphi);   //
	 D_gijBB_tdt();    //variation de la  metrique en BB
	 Calcul_giH_tdt();    // base duale
	 D_giH_tdt();         // variation de la base duale
	 D_gijHH_tdt();    //variation de la  metrique en HH
	 Djacobien_tdt();    // variation du jacobien */
	  // retour des infos
/*	 Impli ex;
	 ex.gijBB_0 = gijBB_0;
	 ex.gijBB_t = gijBB_t; ex.gijBB_tdt = gijBB_tdt;
	 ex.gijHH_tdt = gijHH_tdt;
	 ex.d_gijBB_tdt = &d_gijBB_tdt; ex.d_gijHH_tdt = &d_gijHH_tdt;
	 ex.jacobien = jacobien_tdt; ex.d_jacobien_tdt = &d_jacobien_tdt;
     // liberation des tenseurs intermediaires
     LibereTenseur();   
	 return ex;

	};*/
/* // pour la remontee au infos duaux	
Met_abstraite::InfoImp Met_Sfe1::Cal_InfoImp
   ( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi,Vecteur& phi, int nombre_noeud)
  {   Calcul_M0(tab_noeud,phi,nombre_noeud);
      Calcul_Mtdt(tab_noeud,phi,nombre_noeud);
      Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
	  Calcul_giB_tdt(tab_noeud,dphi,nombre_noeud); // calcul de la base a tdt
	  Calcul_gijBB_0 ();  //      "
	  Calcul_gijHH_0 ();        // composantes controvariantes
	  Calcul_gijBB_tdt ();  //      "
	  Calcul_gijHH_tdt ();        // composantes controvariantes
	  // bases duales
      Calcul_giH_0(); 
	  Calcul_giH_tdt(); 
	  
	  // retour des infos
      InfoImp ex;
      ex.M0 = M0;
      ex.Mtdt = Mtdt;
      ex.giB_0 = giB_0;
      ex.giB_tdt = giB_tdt;
      ex.gijBB_tdt = gijBB_tdt;
      ex.gijHH_tdt = gijHH_tdt;
      ex.giH_0 = giH_0;
      ex.giH_tdt = giH_tdt;
      // liberation des tenseurs intermediaires
      LibereTenseur();   
      return ex;  
   };
   
Met_abstraite::InfoExp Met_Sfe1::Cal_InfoExp
     (Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi,Vecteur& phi, int nombre_noeud)
  {   Calcul_M0(tab_noeud,phi,nombre_noeud);
      Calcul_Mt(tab_noeud,phi,nombre_noeud);
      Calcul_giB_0(tab_noeud,dphi,nombre_noeud); // calcul de la base a t=0
	  Calcul_giB_t(tab_noeud,dphi,nombre_noeud); // calcul de la base a t
	  Calcul_gijBB_0 ();  //      "
	  Calcul_gijHH_0 ();        // composantes controvariantes
	  Calcul_gijBB_t ();  //      "
	  Calcul_gijHH_t ();        // composantes controvariantes
	  // bases duales
      Calcul_giH_0(); 
	  Calcul_giH_t(); 
	  // retour des infos
      InfoExp ex;
      ex.M0 = M0;
      ex.Mt = Mt;
      ex.giB_0 = giB_0;
      ex.giB_t = giB_t;
      ex.gijBB_t = gijBB_t;
      ex.gijHH_t = gijHH_t;
      ex.giH_0 = giH_0;
      ex.giH_t = giH_t;
      // liberation des tenseurs intermediaires
      LibereTenseur();   
      return ex;     
   };
      
// ========== utilise par le contact =========================

// calcul d'un point M0 en fonction de phi et des coordonnees a 0
Coordonnee & Met_Sfe1::PointM_0 // ( stockage a t=0)
		  (Tableau<Noeud *>& tab_noeud,Vecteur& Phi)
  { Calcul_M0(tab_noeud,Phi,tab_noeud.Taille());
    return *M0;
   };

// calcul d'un point M0 en fonction de phi et des coordonnees a t
Coordonnee & Met_Sfe1::PointM_t // ( stockage a t)
		  (Tableau<Noeud *>& tab_noeud,Vecteur& Phi)
  { Calcul_Mt(tab_noeud,Phi,tab_noeud.Taille());
    return *Mt;
   };

// calcul d'un point Mt en fonction de phi et des coordonnees a tdt
Coordonnee & Met_Sfe1::PointM_tdt // ( stockage a tdt)
		  (Tableau<Noeud *>& tab_noeud,Vecteur& Phi)
  { Calcul_Mtdt(tab_noeud,Phi,tab_noeud.Taille());
    return *Mtdt;
   };
		  
// calcul de la base naturel ( stockage a t=0) au point correspondant  au dphi
BaseB& Met_Sfe1::BaseNat_0   // en fonction des coordonnees a t=0
		  (Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi)
  { Calcul_giB_0 (tab_noeud,dphi,tab_noeud.Taille());
    return *giB_0;
   };
		  
// calcul de la base naturel ( stockage a t) au point correspondant  au dphi
BaseB& Met_Sfe1::BaseNat_t   // en fonction des coordonnees a t
		  (Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi)
  { Calcul_giB_t (tab_noeud,dphi,tab_noeud.Taille());
    return *giB_t;
   };
		  
// calcul de la base naturel ( stockage a t=tdt) au point correspondant  au dphi
BaseB& Met_Sfe1::BaseNat_tdt   // en fonction des coordonnees a tdt
		  (Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi)
  { Calcul_giB_tdt (tab_noeud,dphi,tab_noeud.Taille());
    return *giB_tdt;
   };
		  
// calcul de la base naturelle et duale ( stockage a t=0) en fonction des coord a 0 
void Met_Sfe1::BaseND_0(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi,BaseB& bB,BaseH& bH)
  { Calcul_giB_0 (tab_noeud,dphi,tab_noeud.Taille());
    Calcul_gijBB_0 ();  //      "
	Calcul_gijHH_0 ();        // composantes controvariantes
    Calcul_giH_0 ();
    bB = *giB_0;
    bH =  *giH_0;
   };
// calcul de la base naturelle et duale ( stockage a t) en fonction des coord a t 
void Met_Sfe1::BaseND_t(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi,BaseB& bB,BaseH& bH)
  { Calcul_giB_t (tab_noeud,dphi,tab_noeud.Taille());
    Calcul_gijBB_t ();  //      "
	Calcul_gijHH_t ();        // composantes controvariantes
    Calcul_giH_t ();
    bB = *giB_t;
    bH =  *giH_t;
   };
// calcul de la base naturelle et duale ( stockage a tdt) en fonction des coord a tdt 
void Met_Sfe1::BaseND_tdt(Tableau<Noeud *>& tab_noeud,Mat_pleine& dphi,BaseB& bB,BaseH& bH)
  { Calcul_giB_tdt (tab_noeud,dphi,tab_noeud.Taille());
    Calcul_gijBB_tdt ();  //      "
	Calcul_gijHH_tdt ();        // composantes controvariantes
    Calcul_giH_tdt ();
    bB = *giB_tdt;
    bH =  *giH_tdt;
   };  */
		  	        
//==================== METHODES PROTEGEES===============================

/* // calcul de la base  naturel a t0  // a changer 
void Met_Sfe1::Calcul_giB_0
   {}; */
//-----------// calcul du tenseur de courbure  dans la base naturelle
// plusieurs cas sont etudies suivant l'instant considere
  // a l'instant t = 0
Vecteur& Met_Sfe1::courbure_0 (Tableau<Noeud *>& tab_noeud)
  { Vecteur & aiB1 = (*aiB_0)(1); Vecteur & aiB2 = (*aiB_0)(2); // pour simplifier
    Vecteur & aiH1 = (*aiH_0)(1); Vecteur & aiH2 = (*aiH_0)(2); // pour simplifier
    Tableau<Coordonnee> tab_coor(6);
    for (int i=1;i<=6;i++)
      tab_coor(i) = tab_noeud(i)->Coord0();
    return courbure1(tab_coor,aiB1,aiB2,aiH1,aiH2,tab_noeud);
   };      
// a l'instant t 
Vecteur& Met_Sfe1::courbure_t (Tableau<Noeud *>& tab_noeud)
  { Vecteur & aiB1 = (*aiB_t)(1); Vecteur & aiB2 = (*aiB_t)(2); // pour simplifier
    Vecteur & aiH1 = (*aiH_t)(1); Vecteur & aiH2 = (*aiH_t)(2); // pour simplifier
    Tableau<Coordonnee> tab_coor(6);
    for (int i=1;i<=6;i++)
      tab_coor(i) = tab_noeud(i)->Coord1();
    return courbure1(tab_coor,aiB1,aiB2,aiH1,aiH2,tab_noeud);
   };      
// a l'instant t+dt
Vecteur& Met_Sfe1::courbure_tdt (Tableau<Noeud *>& tab_noeud)
  { Vecteur & aiB1 = (*aiB_tdt)(1); Vecteur & aiB2 = (*aiB_tdt)(2); // pour simplifier
    Vecteur & aiH1 = (*aiH_tdt)(1); Vecteur & aiH2 = (*aiH_tdt)(2); // pour simplifier
    Tableau<Coordonnee> tab_coor(6);
    for (int i=1;i<=6;i++)
      tab_coor(i) = tab_noeud(i)->Coord2();
    return courbure1(tab_coor,aiB1,aiB2,aiH1,aiH2,tab_noeud);
   };      
       
// routine generale de calcul de la courbure
Vecteur& Met_Sfe1::courbure
 (Tableau<Coordonnee>& tab_coor,Vecteur & aiB1,Vecteur & aiB2,
             Vecteur & aiH1,Vecteur & aiH2,Tableau<Noeud *>& tab_noeud)
   { Vecteur curbp(3),curb(3); // : lestenseurs de courbures d'abord dans les axes
     // locaux des arretes du triangle puis dans le repere naturel: b11,b12,b22
     Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure 
     // calcul de la normale centrale 
     Vecteur N = (ProdVec(aiB1,aiB2)).Normer();
     // calcul du centre de gravite
     Coordonnee G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.;
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- tout d'abord un tableau d'indice pour ne pas etre embete par les modulos 3
     Vecteur indi(5);indi(1) = 1; indi(2) = 2; indi(3) = 3; indi(4) = 1; indi(5) = 2;
        // 2- on boucle sur le cotes du triangle
     Vecteur Ni(3); // normale du triangle exterieur   
     for (int ncot=1;ncot<=3;ncot++)
       { // on utilise les notations : triangle principal DBC et secondaire BAC
         // on cherche a calculer la longueur H2G2
           // simplification de l'ecriture
         Coordonnee& D = tab_noeud(indi(ncot))->Coord0();
         Coordonnee& B = tab_noeud(indi(ncot+1))->Coord0();
         Coordonnee& C = tab_noeud(indi(ncot+2))->Coord0();
         Coordonnee& A = tab_noeud(ncot+1)->Coord0();          
           // calcul de G2, centre de gravite du triangle exterieur           
         Coordonnee G2 = (B + C + A) / 3.;
         // calcul de U1                      
         Vecteur CB = B - C;
         Vecteur U1 = CB / CB.Norme(); 
         // vecteur CG2
         Vecteur CG2 = G2 - C;
         // longueur H2G2 et GH1
         Vecteur H2G2 = CG2 - (CG2 * U1) * U1;
         double h2G2 = CG2.Norme();         
         Vecteur CG = G - C;
         double GH1 = sqrt ( CG * CG - SQR( CG * U1)) ; 
         // calcul de la normale du triangle exterieur
         Vecteur CA = A - C;
         Ni =  (ProdVec(CB,CA)).Normer(); 
         // calcul de la derivee de la normale
         Vecteur dNdV = ( Ni - N) / ( GH1 + h2G2);
         // courbure dans la direction u2 c-a-d H2G2
         curbp(ncot) = dNdV * H2G2;
         // calcul de la matrice de passage de Ualpha a aalpha
         double Beta21 = H2G2 * aiH1;
         double Beta22 = H2G2 * aiH2;
         // remplissage de la matrice mat
         mat(ncot,1) = Beta21*Beta21;
         mat(ncot,2) = Beta21*Beta22;
         mat(ncot,3) = Beta22*Beta22;
        };
	 
      // calcul du tenseur de courbure dans la base ai
      Mat_pleine mati = mat.Inverse();
      curb = mati * curbp;
      // retour de la courbure
      return curb;	 
   };

//--------------// calcul du tenseur de courbure et de sa variation 
// plusieurs cas sont etudies suivant l'instant considere
// a l'instant t 
void Met_Sfe1::Dcourbure_t 
 (Tableau<Noeud *>& tab_noeud,Vecteur& curb,TabOper<Vecteur>& dcurb)
  { Vecteur & aiB1 = (*aiB_t)(1); Vecteur & aiB2 = (*aiB_t)(2); // pour simplifier
    Vecteur & aiH1 = (*aiH_t)(1); Vecteur & aiH2 = (*aiH_t)(2); // pour simplifier
    Tableau<Coordonnee> tab_coor(6);
    for (int i=1;i<=6;i++)
      tab_coor(i) = tab_noeud(i)->Coord1();
    Dcourbure1 (tab_coor,aiB1,aiB2,aiH1,aiH2,*d_aiB_t,*d_aiH_t,tab_noeud,curb,dcurb);
   };      
// a l'instant t+dt
void Met_Sfe1::Dcourbure_tdt 
 (Tableau<Noeud *>& tab_noeud,Vecteur& curb,TabOper<Vecteur>& dcurb)
  { Vecteur & aiB1 = (*aiB_tdt)(1); Vecteur & aiB2 = (*aiB_tdt)(2); // pour simplifier
    Vecteur & aiH1 = (*aiH_tdt)(1); Vecteur & aiH2 = (*aiH_tdt)(2); // pour simplifier
    Tableau<Coordonnee> tab_coor(6);
    for (int i=1;i<=6;i++)
      tab_coor(i) = tab_noeud(i)->Coord2();
    Dcourbure1 (tab_coor,aiB1,aiB2,aiH1,aiH2,*d_aiB_t,*d_aiH_t,tab_noeud,curb,dcurb);
   };      
       
// routine generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure  b11,b12,b22
//            dcurb , la variation de courbure.
void Met_Sfe1::Dcourbure
 (Tableau<Coordonnee>& tab_coor,Vecteur & aiB1,Vecteur & aiB2,
          Vecteur & aiH1,Vecteur & aiH2,Tableau <BaseB>& DaiB,Tableau <BaseH>& DaiH,
          Tableau<Noeud *>& tab_noeud,Vecteur& curb,TabOper <Vecteur>& dcurb)
   { int nbddl = 18;   // nombre de ddl total
     int nbddlint = 9; // nombre de ddl interne
     Vecteur Nul(3); // le vecteur nul de dim 3, pour les initialisations
     
     Vecteur curbp(3); // : les tenseurs de courbures  dans les axes
                       // locaux des arretes du triangle
     Tableau <Vecteur> dcurbp(nbddl); // variation des courbures
     Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
     Tableau <Mat_pleine> dmat(nbddl,mat); 
     // calcul de la normale centrale 
     Vecteur NN = ProdVec(aiB1,aiB2); // normale non normee
     double nor = NN.Norme(); // la norme
     Vecteur N = NN / nor;
     Vecteur DNN; // vecteur intermediair
     Tableau <Vecteur> dN(nbddl,Nul);
     for (int i1=1;i1<= nbddlint;i1++)
       { DNN = ProdVec(DaiB(i1)(1),aiB2) + ProdVec(aiB1,DaiB(i1)(2));
         dN(i1) = VarUnVect(NN,DNN,nor);
        }
     // calcul du centre de gravite
     Coordonnee G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.;
     Tableau <Coordonnee> dG(nbddl);
     for (int ia=1;ia<= 3;ia++)
       {dG(ia) = Ia(ia)/3. ;dG(ia+3) = Ia(ia)/3. ;dG(ia+6) = Ia(ia)/3. ;}
     
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- tout d'abord un tableau d'indice pour ne pas etre embete par les modulos 3
     Tableau <int> indi(5);indi(1) = 1; indi(2) = 2; indi(3) = 3; indi(4) = 1; indi(5) = 2;
        // 2- on boucle sur le cotes du triangle
     Vecteur Ni(3); // normale du triangle exterieur 
     // -------- declaration pour la boucle sur les 3 cotes ------------
     Coordonnee G2(3); // init a zero et a 3 composantes 
     Tableau <Coordonnee> dD(nbddl,Nul),dB(nbddl,Nul),dC(nbddl,Nul),dA(nbddl,Nul) ; // variation des points
     Tableau <Coordonnee> dG2(nbddl,Nul);
     Vecteur CB(3),U1(3),CG2(3),H2G2(3);
     Tableau <Vecteur> dCB(nbddl,Nul),dU1(nbddl,Nul),dH2G2(nbddl,Nul);
     Tableau <double> dh2G2(nbddl);         
     Vecteur CG(3),CA(3),dNdV(3);
     Tableau <Vecteur> dCA(nbddl,Nul),dNi(nbddl,Nul);
     double under3 = 1./3.;
     Vecteur dCG(3),dGH1(3),ddNdv(3),dCG2_jnB,dCG2_jnC,dCG2_jnA; 
     for (int ncot=1;ncot<=3;ncot++)
       { // on utilise les notations : triangle principal DBC et secondaire BAC
         // on cherche a calculer la longueur H2G2
           // simplification de l'ecriture
         Coordonnee& D = tab_noeud(indi(ncot))->Coord0();
         Coordonnee& B = tab_noeud(indi(ncot+1))->Coord0();
         Coordonnee& C = tab_noeud(indi(ncot+2))->Coord0();
         Coordonnee& A = tab_noeud(ncot+1)->Coord0();          
         // calcul de G2, centre de gravite du triangle exterieur           
         G2 = (B + C + A) * under3;
         // calcul de U1                      
         CB = B - C;
         double norCB = CB.Norme();
         U1 = CB / norCB;
         // vecteur CG2
         CG2 = G2 - C;
         // longueur H2G2 et GH1
         H2G2 = CG2 - (CG2 * U1) * U1;
         double h2G2 = CG2.Norme();
         double unSurh2G2= 1./h2G2;
         CG = G - C;
         double GH1 = sqrt ( CG * CG - SQR( CG * U1)) ;
         // calcul de la normale du triangle exterieur
         CA = A - C;
         Ni =  ProdVec(CB,CA);
         double norNi = Ni.Norme();
         Ni /=  norNi;
         // calcul de la derivee de la normale
         dNdV = ( Ni - N) / ( GH1 + h2G2);
         // courbure dans la direction u2 c-a-d H2G2
         curbp(ncot) = dNdV * H2G2;
         // calcul de la matrice de passage de Ualpha a aalpha
         double Beta21 = H2G2 * aiH1;
         double Beta22 = H2G2 * aiH2;
         // remplissage de la matrice mat
         mat(ncot,1) = Beta21*Beta21;
         mat(ncot,2) = Beta21*Beta22;
         mat(ncot,3) = Beta22*Beta22;
         
         // calcul des variations
         for (int ia=1;ia<= 3;ia++)
           { // on a 4 indices pour chaque ia, pour lesquels la variation est non nulle
             // dans cette boucle on n'utilise pas systematiquement les 4 indices
             int jnD = (indi(ncot)-1)*3+ia; int  jnB = (indi(ncot+1)-1)*3+ia;
             int jnC = (indi(ncot+2)-1)*3+ia; int jnA = 9+(ncot-1)*3+ia;
             // tout d'abord les variations des points         
             dD(jnD) = Ia(ia);dB(jnB) = Ia(ia);
             dC(jnC) = Ia(ia);dA(jnA) = Ia(ia);
             // centre de gravite
             dG2(jnB) = Ia(ia)*under3;
             dG2(jnC) = Ia(ia)*under3;
             dG2(jnA) = Ia(ia)*under3;
             // vecteur CB
             dCB(jnB) = dB(jnB); dCB(jnC) = - dC(jnC);
             // vecteur U1
             dU1(jnB) = VarUnVect(CB,dCB(jnB),norCB);  
             dU1(jnC) = VarUnVect(CB,dCB(jnC),norCB);
             // vecteur CG2
             dCG2_jnB = dG2(jnB) - dC(jnB);  
             dCG2_jnC = dG2(jnC) - dC(jnC);  
             dCG2_jnA = dG2(jnA);
             // vecteur H2G2= CG2 - (CG2 * U1) * U1;
             dH2G2(jnB) = dCG2_jnB-(dCG2_jnB*U1 + CG2*dU1(jnB))* U1- (CG2 * U1)* dU1(jnB);
             dH2G2(jnC) = dCG2_jnC-(dCG2_jnC*U1 + CG2*dU1(jnC))* U1- (CG2 * U1)* dU1(jnC);
             dH2G2(jnA) = dCG2_jnA-(dCG2_jnA*U1 + CG2*dU1(jnA))* U1- (CG2 * U1)* dU1(jnA);
             // scalaire  h2G2 = CG2.Norme();
             dh2G2(jnB) = dH2G2(jnB) * H2G2 /unSurh2G2;
             dh2G2(jnC) = dH2G2(jnC) * H2G2 /unSurh2G2;
             dh2G2(jnA) = dH2G2(jnA) * H2G2 /unSurh2G2;
             // Vecteur CA = A - C
             dCA(jnC) = -dC(jnC); dCA(jnA) = dA(jnA);
             // Vecteur Ni =  (ProdVec(CB,CA)).Normer()
             dNi(jnC) = VarUnVect(Ni,(ProdVec(dCB(jnC),CA)+ProdVec(CA,dCA(jnC))),norNi);
             dNi(jnB) = VarUnVect(Ni,(ProdVec(dCB(jnB),CA)+ProdVec(CA,dCA(jnB))),norNi);
             dNi(jnA) = VarUnVect(Ni,(ProdVec(dCB(jnA),CA)+ProdVec(CA,dCA(jnA))),norNi);
            }
          for (int ib=1;ib<= 3;ib++)
           { // on a 4 indices pour chaque ia, pour lesquels la variation est non nulle
             // mais ici on utilise systematiquement les 4 indices d'ou une seconde boucle
             // egalement moins de stockage intermediaire
             int jnb[4] = {(indi(ncot)-1)*3+ib,(indi(ncot+1)-1)*3+ib,
                          (indi(ncot+2)-1)*3+ib,9+(ncot-1)*3+ib};
             for (int jb=0;jb<=3;jb++)
               { // vecteur CG = G - C
                 dCG = dG(jnb[jb]) - dC(jnb[jb]);
                 // scalaire GH1 = sqrt ( CG * CG - SQR( CG * U1)) 
                 dGH1 = (dCG * CG - (CG * U1) * (dCG * U1 + CG * dU1(jnb[jb])))/ GH1 ;
                 // vecteur dNdV = ( Ni - N) / ( GH1 + h2G2)
                 ddNdv = (dNi(jnb[jb]) - dN(jnb[jb])) / ( GH1 + h2G2)
                       - ( dGH1 + dh2G2(jnb[jb])) * (Ni - N) / ( SQR ( GH1 + h2G2));             
                 // courbure  curbp(ncot) = dNdV * H2G2;
                 dcurbp(jnb[jb])(ncot) = ddNdv * H2G2 + dNdV * dH2G2(jnb[jb]);                            
                 // matrice de passage : Beta21 = H2G2 * aiH1 et  Beta22 = H2G2 * aiH2
                 double dBeta21 = dH2G2(jnb[jb]) * aiH1 + H2G2 * DaiH(1)(jnb[jb]);
                 double dBeta22 = dH2G2(jnb[jb]) * aiH2 + H2G2 * DaiH(2)(jnb[jb]);
                 //  matrice mat : mat(ncot,1) = Beta21*Beta21 ;  mat(ncot,2) = Beta21*Beta22;
                 // mat(ncot,3) = Beta22*Beta22;
                 dmat(jnb[jb])(ncot,1) = 2. * dBeta21 * Beta21;
                 dmat(jnb[jb])(ncot,2) = dBeta21 * Beta22 + Beta21 * dBeta22;
                 dmat(jnb[jb])(ncot,3) = dBeta22 * Beta22 + Beta22 * dBeta22;               
                }             
            } 
        };
	 
      // calcul du tenseur de courbure dans la base ai
      Mat_pleine mati = mat.Inverse();
      curb = mati * curbp;
      // variations
      Mat_pleine dmati(mati);
      for (int ddl=1;ddl<=nbddl;ddl++)
        { dmati = - mat * ( dmat(ddl) * mati);
          dcurb(ddl) = dmati * curbp + mati * dcurbp(ddl);
         }
   };



/*// 
void Met_Sfe1::courbure2
 (Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
   { // tout d'abord on calcul la base naturelle du triangle central
     if (giB_0 =! NULL)
         Calcul_giB_0 (tab_noeud,dphi,3); // 3 pour les 3 noeuds centrals 
     Vecteur & aiB1 = (*giB_0)(1); Vecteur & aiB2 = (*giB_0)(2); // pour simplicite
     Calcul_gijBB_0 ();  // puis la metrique      "
	 Calcul_gijHH_0 ();  // composantes controvariantes
	 // les bases duales
     Calcul_giH_0(); 
     Vecteur & aiH1 = (*giH_0)(1); Vecteur & aiH2 = (*giH_0)(2); // pour simplicite 
     
     Vecteur curbp(3),curb(3); // lestenseurs de courbures d'abord dans les axes
     // locaux des arretes du triangle puis dans le repere naturel: b11,b12,b22
     Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure 
     // calcul de la normale centrale 
     Vecteur N = (ProdVec(aiB1,aiB2)).Normer();
     // calcul du centre de gravite
     Coordonnee G = (tab_noeud(1)->Coord0() + tab_noeud(2)->Coord0()
                    + tab_noeud(3)->Coord0()) / 3.;
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- tout d'abord un tableau d'indice pour ne pas etre embete par les modulos 3
     Vecteur indi(5);indi(1) = 1; indi(2) = 2; indi(3) = 3; indi(4) = 1; indi(5) = 2;
        // 2- on boucle sur le cotes du triangle
     Vecteur Ni(3); // normale du triangle exterieur   
     for (int ncot=1;ncot<=3;ncot++)
       { // on utilise les notations : triangle principal DBC et secondaire BAC
         // on cherche l'image de par rotation du plan BAC dans le plan DBC
         // le point A devient alors Ap
          // calcul des vecteurs cote
         Vecteur CA = tab_noeud(ncot+3)->Coord0() - tab_noeud(indi(ncot)+2)->Coord0();
         Vecteur CB = tab_noeud(indi(ncot)+1)->Coord0() - tab_noeud(indi(ncot)+2)->Coord0();
         Vecteur BA = tab_noeud(ncot+3)->Coord0() - tab_noeud(indi(ncot)+1)->Coord0();
         // calcul du terme (i) appele xi
         double CA2 = CA * CA; 
         double xi = CA2 + CB * CB - BA * BA;
         // calcul des coefficients a et b
         double g1CB = aiB1 * CB;
         double a = -(aiB2 * CB) / g1CB ;
         double b = xi * 0.5 / g1CB;
         // calcul des coefficients de l'equation du second degre
         Vecteur ag1g2 = a * aiB1 + aiB2;
         double aa = ag1g2 * ag1g2;
         double bb = (2 * b) * (aiB1 * ag1g2);
         double cc = (b*b) * (aiB1 * aiB1) - CA2;
         // resolution de l'equation du second degre
         double x1,x2; int isol;
            // ResoEqua2(aa,bb,cc,x1,x2,isol);
         // calcul de la position du point externe sur le plan de la facette centrale
         double teta1 = x1 * a + b; // choix de x1 a priori pour teta2
         Coordonnee Ap = teta1 * aiB1 + x1 * aiB2;
         // verif du point
         Vecteur DC = tab_noeud(indi(ncot)+2)->Coord0() - tab_noeud(indi(ncot))->Coord0();
         Vecteur CAp = Ap - tab_noeud(indi(ncot)+2)->Coord0();
         if (DC * CAp < 0)
          { // cas ou c'etait x2 le teta2
            teta1 = x2 * a + b;
            Ap = teta1 * aiB1 + x2 * aiB2;
            CAp = Ap - tab_noeud(indi(ncot)+2)->Coord0();
           }
         // calcul de la normale du triangle exterieur
         Ni =  (ProdVec(CB,CAp)).Normer(); 
         // calcul de GH1 et GpH2
         Vecteur CG = G - tab_noeud(indi(ncot)+2)->Coord0();
         Vecteur U1 = CB / CB.Norme(); 
         double GH1 = sqrt ( CG * CG - ( CG * U1) * ( CG * U2)) ; // a modifier !!!!
         double H2Gp = sqrt ( CG * CGp - ( CGp * U2) * ( CGp * U2)) ; // a modifier !!!!
         // calcul de la derivee de la normale
         Vecteur dNdV = ( Ni - N) / ( GH1 + H2Gp);
         // calcul du vecteur normal a U1, c'a dire U2, mais qui n'est pas norme
         Vecteur U2 = CAp - (CAp * U1) * U1;
         // courbure dans la direction u22
         curbp(ncot) = dNdV * U2;
         // calcul de la matrice de passage de Ualpha a aalpha
         double Beta21 = U2 * giH(1);
         double Beta22 = U2 * giH(2);
         // remplissage de la matrice mat
         mat(ncot,1) = Beta21*Beta21;
         mat(ncot,2) = Beta21*Beta22;
         mat(ncot,3) = Beta22*Beta22;
        };
	 
      // calcul du tenseur de courbure dans la base gi
      Mat_pleine mati = mat.Inverse();
      curb = mati * curbp;	 
   };
   
   
// calcul de la base  naturel a t 
void Met_Sfe1::Calcul_giB_t
 ( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
   {	
	 for (int a=1;a<= dim_base;a++)
		 { (*(giB_t))(1)(a)= 0.;
		   for (int r=1;r<=nombre_noeud;r++)
			{ (*(giB_t))(1)(a)=(*(giB_t))(1)(a) + tab_noeud(r)->Coord1()(a) * dphi(1,r);
			};
		};
   };
// calcul de la base  naturel a tdt  
void Met_Sfe1::Calcul_giB_tdt
 ( Tableau<Noeud *>& tab_noeud, Mat_pleine& dphi, int nombre_noeud)
   {
	for (int a=1;a<= dim_base;a++)
		 { (*(giB_tdt))(1)(a)= 0.;
		   for (int r=1;r<=nombre_noeud;r++)
			{ (*(giB_tdt))(1)(a)=(*(giB_tdt))(1)(a) + tab_noeud(r)->Coord2()(a) * dphi(1,r);
			};
		};
   };   */