// FICHIER : Met_Sfe1s3.cc
// CLASSE : Met_Sfe1

// 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-2021 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 <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"


#include "Met_Sfe1.h"


//  =========================================================================================
//  vu la taille des executables le fichier est decompose en quatre
// le premier : Met_Sfe1s1.cp concerne les constructeurs, destructeur et 
//              la gestion des variables
// le second :  Met_Sfe1s2.cp concerne  le calcul des grandeurs publics
// le troisieme : Met_Sfe1s3.cp concerne  le calcul des grandeurs protegees, et des courbures pour SFE1 et SFE2 et QSFE1
// le quatrième : Met_Sfe1s4.cp concerne  le calcul des courbures pour les éléments SFE3 et QSFE3
//  =========================================================================================
		  
		  	        
//==================== METHODES PROTEGEES===============================
   
// calcul des normales a la facette		  
void Met_Sfe1::Calcul_N_0 ()
   { N_0 = (Util::ProdVec_coorB((*aiB_0)(1),(*aiB_0)(2))).Normer();};
void Met_Sfe1::Calcul_N_t ()
   { N_t = (Util::ProdVec_coorB((*aiB_t)(1),(*aiB_t)(2)));
     norme_N_t = N_t.Norme();
     N_t /= norme_N_t;
   };
void Met_Sfe1::Calcul_N_tdt ()
   { N_tdt = (Util::ProdVec_coorB((*aiB_tdt)(1),(*aiB_tdt)(2)));
     norme_N_tdt = N_tdt.Norme();
     N_tdt /= norme_N_tdt;
   };

//== calcul des variations de la normales
// il faut avant le calcul des ap_alpha_t et de épaisseur à t
void Met_Sfe1::Cal_d_N_t()
  { int nbddlx = tab_ddl.NbDdl_famille(X1);
	   CoordonneeB d_NN; // un vecteur de travail
    for (int iddl=1;iddl<= nbddlx;iddl++)
      { // tout d'abord on calcul la variation de la normale non normée
    	   d_NN =   (Util::ProdVec_coorB((*d_aiB_t)(iddl)(1),(*aiB_t)(2)))
    	          + (Util::ProdVec_coorB((*aiB_t)(1),(*d_aiB_t)(iddl)(2)));
       	// maintenant on calcul la variation du vecteur normal unitaire
        (*d_N_t)(iddl) =  Util::VarUnVect_coorB(N_t,d_NN,norme_N_t);
       };
	};

// il faut avant le calcul des ap_alpha_tdt et de épaisseur à tdt
void Met_Sfe1::Cal_d_N_tdt() 
  { int nbddlx = tab_ddl.NbDdl_famille(X1);
	   CoordonneeB d_NN; // un vecteur de travail
    for (int iddl=1;iddl<= nbddlx;iddl++)
      { // tout d'abord on calcul la variation de la normale non normée
    	   d_NN =   (Util::ProdVec_coorB((*d_aiB_tdt)(iddl)(1),(*aiB_tdt)(2)))
    	          + (Util::ProdVec_coorB((*aiB_tdt)(1),(*d_aiB_tdt)(iddl)(2)));
    	   // maintenant on calcul la variation du vecteur normal unitaire
        (*d_N_tdt)(iddl) =  Util::VarUnVect_coorB(N_tdt,d_NN,norme_N_tdt);
       };
	};
	
// ---- calcul des points
	// calcul du point a t0
void Met_Sfe1::Calcul_M0( const Tableau<Noeud *>& , const  Vecteur& phiH, int )
      { // on considère qu'il y a deux points dans l'épaisseur, un a -h/2 et l'autre à h/2
          *M0 = *P0 +   (N_0 * (epais.epaisseur0 * 0.5 * (phiH(2)-phiH(1)))).Coor_const();};
void Met_Sfe1::Calcul_Mt( const Tableau<Noeud *>& , const  Vecteur& phiH, int )
      { // on considère qu'il y a deux points dans l'épaisseur, un a -h/2 et l'autre à h/2
          *Mt = *Pt +   (N_t * (epais.epaisseur_t * 0.5 * (phiH(2)-phiH(1)))).Coor_const();};
void Met_Sfe1::Calcul_Mtdt( const Tableau<Noeud *>& ,  const Vecteur& phiH, int )
      { // on considère qu'il y a deux points dans l'épaisseur, un a -h/2 et l'autre à h/2
          *Mtdt = *Ptdt +   (N_tdt * (epais.epaisseur_tdt * 0.5 * (phiH(2)-phiH(1)))).Coor_const();};

// cas de l'épaisseur interpolée à 0
void Met_Sfe1::Calcul_epaisseur_0(const Tableau<Noeud *>& tab_noeud, const Vecteur& phiS)
 {  double & epaisseur = epais.epaisseur0;
    epaisseur = 0.; // init
	for (int r=1;r<=nbNoeudCentral;r++)
       epaisseur += tab_noeud(r)->Valeur_0(EPAIS) * phiS(r);
	};
	    
// cas de l'épaisseur interpolée à t
void Met_Sfe1::Calcul_epaisseur_t(const Tableau<Noeud *>& tab_noeud, const Vecteur& phiS)
 {  double & epaisseur = epais.epaisseur_t;
    epaisseur = 0.; // init
	for (int r=1;r<=nbNoeudCentral;r++)
       epaisseur += tab_noeud(r)->Valeur_t(EPAIS) * phiS(r);
	};    

// cas de l'épaisseur interpolée à tdt
void Met_Sfe1::Calcul_epaisseur_tdt(const Tableau<Noeud *>& tab_noeud, const Vecteur& phiS)
 {  double & epaisseur = epais.epaisseur_tdt;
    epaisseur = 0.; // init
	for (int r=1;r<=nbNoeudCentral;r++)
       epaisseur += tab_noeud(r)->Valeur_tdt(EPAIS) * phiS(r);
	};    

// calcul de la variation de l'épaisseur par rapport aux ddl         
void Met_Sfe1::D_epaisseur_t(const  Vecteur& phiS)
	{ // l'épaisseur ne dépend que des ddl d'épaisseur
	  for (int r=1;r<=nbNoeudCentral;r++)
	    (*d_epais)(r).epaisseur_t = phiS(r);		
	};

// calcul de la variation de l'épaisseur par rapport aux ddl         
void Met_Sfe1::D_epaisseur_tdt(const  Vecteur& phiS)
	{ // l'épaisseur ne dépend que des ddl d'épaisseur
	  for (int r=1;r<=nbNoeudCentral;r++)
	    (*d_epais)(r).epaisseur_tdt = phiS(r);		
	};

//==  calcul de la variation de N,alpha
void Met_Sfe1::Cal_d_N_alpha_t()
  { int nbddlx = tab_ddl.NbDdl_famille(X1);
	   for (int iddl=1;iddl<= nbddlx;iddl++)
     { // variation de la derivee du vecteur normal
	      (*d_N_alpha_t_1)(iddl).ConstructionAPartirDe_H(
	                    -(dcurb_t(iddl)(1,1) * (*aiH_t)(1) + curb_t(1,1) * (*d_aiH_t)(iddl)(1) 
	                    +dcurb_t(iddl)(1,2) * (*aiH_t)(2) + curb_t(1,2) * (*d_aiH_t)(iddl)(2)));  
	      (*d_N_alpha_t_2)(iddl).ConstructionAPartirDe_H(
	                    -(dcurb_t(iddl)(2,1) * (*aiH_t)(1) + curb_t(2,1) * (*d_aiH_t)(iddl)(1) 
	                    + dcurb_t(iddl)(2,2) * (*aiH_t)(2) + curb_t(2,2) * (*d_aiH_t)(iddl)(2))); 
    };
 };

void Met_Sfe1::Cal_d_N_alpha_tdt()
  { int nbddlx = tab_ddl.NbDdl_famille(X1);
	   for (int iddl=1;iddl<= nbddlx;iddl++)
     { // variation de la derivee du vecteur normal
	      (*d_N_alpha_tdt_1)(iddl).ConstructionAPartirDe_H(
	                    -(dcurb_tdt(iddl)(1,1) * (*aiH_tdt)(1) + curb_tdt(1,1) * (*d_aiH_tdt)(iddl)(1) 
	                    +dcurb_tdt(iddl)(1,2) * (*aiH_tdt)(2) + curb_tdt(1,2) * (*d_aiH_tdt)(iddl)(2)));  
	      (*d_N_alpha_tdt_2)(iddl).ConstructionAPartirDe_H(
	                    -(dcurb_tdt(iddl)(2,1) * (*aiH_tdt)(1) + curb_tdt(2,1) * (*d_aiH_tdt)(iddl)(1) 
	                    + dcurb_tdt(iddl)(2,2) * (*aiH_tdt)(2) + curb_tdt(2,2) * (*d_aiH_tdt)(iddl)(2))); 
    };
 };

// ------ calcul des bases pour la surface médiane -------

// calcul de la base  naturel a t0  pour la surface médiane: partie a_alpha
void Met_Sfe1::Calcul_a_alpha_B_0
 (  const Tableau<Noeud *>& tab_noeud,  const Mat_pleine& dphi, int ,const Vecteur&)
  {
    #ifdef MISE_AU_POINT	 	 
	   if  (aiB_0 == NULL)
      { cout << "\nErreur : la base a t=0 n'est pas dimensionne !\n";
        cout << "void Met_Sfe1::Calcul_aiB_0 \n";
        Sortie(1);
      };
    #endif
	   for (int i=1;i<= 2;i++)
	    {for (int a=1;a<= dim_base;a++)
		     { double & aib0_i_a = aiB_0->CoordoB(i)(a);
		       aib0_i_a = 0.;
		       for (int r=1;r<=nbNoeudCentral;r++)
			        { aib0_i_a += tab_noeud(r)->Coord0()(a) * dphi(i,r);
			        };
		     };
	    };
 };
// calcul de la base  naturel a t    pour la surface médiane: partie a_alpha
void Met_Sfe1::Calcul_a_alpha_B_t
 (  const Tableau<Noeud *>& tab_noeud, const  Mat_pleine& dphi, int ,const Vecteur&)
  {
    #ifdef MISE_AU_POINT	 	 
	   if  (aiB_t == NULL)
			  { cout << "\nErreur : la base a t n'est pas dimensionne !\n";
			    cout << "void Met_Sfe1::Calcul_aiB_t \n";
			    Sortie(1);
			  };
    #endif
	   for (int i=1;i<= 2;i++)
	     {for (int a=1;a<= dim_base;a++)
		      { double & aibt_i_a = aiB_t->CoordoB(i)(a);
		        aibt_i_a = 0.;
          for (int r=1;r<=nbNoeudCentral;r++)
			         { aibt_i_a += tab_noeud(r)->Coord1()(a) * dphi(i,r);
			         };
		      };
	    };
 };

// calcul de la base  naturel a tdt    pour la surface médiane: partie a_alpha
void Met_Sfe1::Calcul_a_alpha_B_tdt
 (  const Tableau<Noeud *>& tab_noeud, const  Mat_pleine& dphi, int ,const Vecteur&)
 {
   #ifdef MISE_AU_POINT
	  if  (aiB_tdt == NULL)
     { cout << "\nErreur : la base a t+dt n'est pas dimensionne !\n";
       cout << "void Met_Sfe1::Calcul_aiB_tdt \n";
       Sortie(1);
     };
   #endif
   for (int i=1;i<= 2;i++)
     {for (int a=1;a<= dim_base;a++)
        { double & aibtdt_i_a = aiB_tdt->CoordoB(i)(a);
          aibtdt_i_a = 0.;
          for (int r=1;r<=nbNoeudCentral;r++)
            { aibtdt_i_a += tab_noeud(r)->Coord2()(a) * dphi(i,r);
            };
		     };
	    };
};
   		 

//==  calcul de la variation des bases pour la surface médiane: partie a_alpha
// ici il ne s'agit uniquement que de la variation de a_alpha et pas a_3
 void Met_Sfe1::D_a_alpha_B_t( const  Mat_pleine& dphi, int ,const  Vecteur & )
  { int indice; 
    // derivees de aBi par rapport a XHbr      
    for (int b = 1; b<= dim_base; b++)
      for (int r = 1; r <= nbNoeudCentral; r++)
         { indice = (r-1)*dim_base + b;
           BaseB & d_aiB_t_indice=(*d_aiB_t)(indice);
           for (int i =1;i<=2;i++)  // alpha = 1 et 2 uniquement
              { d_aiB_t_indice.CoordoB(i).Zero();
                d_aiB_t_indice.CoordoB(i)(b)=dphi(i,r);
               }
          };              
  };
// ici il ne s'agit uniquement que de la variation de a_alpha et pas a_3
 void Met_Sfe1::D_a_alpha_B_tdt(  const Mat_pleine& dphi, int ,const  Vecteur &)
  { int indice;       
    for (int b = 1; b<= dim_base; b++)
      for (int r = 1; r <= nbNoeudCentral; r++)
         { indice = (r-1)*dim_base + b;
           BaseB & d_aiB_tdt_indice=(*d_aiB_tdt)(indice);
           for (int i =1;i<=2;i++)  // alpha = 1 et 2 uniquement
              { d_aiB_tdt_indice.CoordoB(i).Zero();
                d_aiB_tdt_indice.CoordoB(i)(b)=dphi(i,r);
               }
          }              
  };
   
//==  calcul de la variation de a_3_B
void Met_Sfe1::D_a_3_B_t()
  { int nbddlx = tab_ddl.NbDdl_famille(X1);
    // calcul de la variation de a3B par rapport au ddlxi
    for (int iddl=1;iddl<= nbddlx;iddl++)
    	 (*d_aiB_t)(iddl).CoordoB(3) =  (epais.epaisseur_t * 0.5) * (*d_N_t)(iddl);
    // calcul de la variation en fonction des ddl d'épaisseur 
    int nbddl_total = nbddlx + nbNoeudCentral;int ih=1;
    for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
	 	 (*d_aiB_t)(iddl).CoordoB(3) = ((*d_epais)(ih).epaisseur_t * 0.5 ) * N_t;
 };

void Met_Sfe1::D_a_3_B_tdt()
  { int nbddlx = tab_ddl.NbDdl_famille(X1);
    // calcul de la variation de a3B par rapport au ddlxi
    for (int iddl=1;iddl<= nbddlx;iddl++)
    	 (*d_aiB_tdt)(iddl).CoordoB(3) =  (epais.epaisseur_tdt * 0.5) * (*d_N_tdt)(iddl);
    // calcul de la variation en fonction des ddl d'épaisseur 
    int nbddl_total = nbddlx + nbNoeudCentral;int ih=1;
    for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
	 	 (*d_aiB_tdt)(iddl).CoordoB(3) = ((*d_epais)(ih).epaisseur_tdt * 0.5 ) * N_tdt;
 };
		  
// --------- calcul de la base  naturel a t0 ------ 
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de la facette
// également l'épaisseur interpolée si on est avec des éléments 3D
void Met_Sfe1::Calcul_giB_0
 (   const  Tableau<Noeud *>& ,  const Mat_pleine& , int , const  Vecteur& phiH)
   { 
    #ifdef MISE_AU_POINT	 	 
	   if  (giB_0 == NULL)
     { cout << "\nErreur : la base a t=0 n'est pas dimensionne !\n";
       cout << "void Met_Sfe1::Calcul_giB_0 \n";
       Sortie(1);
     };
    #endif
    // vecteur de base gi,
    // 	theta^3 = (phiH(2)-phiH(1));
    double z0 = epais.epaisseur0 * 0.5 * (phiH(2)-phiH(1));
    // ici il n'y a pas le terme de variation de h par rapport au theta^alpha
    giB_0->CoordoB(1) = (*aiB_0)(1) + N_alpha_0_1 * z0;
    giB_0->CoordoB(2) = (*aiB_0)(2) + N_alpha_0_2 * z0;
    if (tCalEpais) // cas 3D
          giB_0->CoordoB(3) = (*aiB_0)(3);
  };

// calcul de la base  naturel a t  
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de la facette
void Met_Sfe1::Calcul_giB_t
 (  const  Tableau<Noeud *>& , const  Mat_pleine& , int, const  Vecteur& phiH)
   {	
    #ifdef MISE_AU_POINT	 	 
    if  (giB_t == NULL)
     { cout << "\nErreur : la base a t n'est pas dimensionne !\n";
       cout << "void Met_Sfe1::Calcul_giB_t \n";
       Sortie(1);
     };
    #endif
    // vecteur de base gi
    double z_t = epais.epaisseur_t * 0.5 * (phiH(2)-phiH(1));
    // ici il n'y a pas le terme de variation de h par rapport au theta^alpha
    giB_t->CoordoB(1) = (*aiB_t)(1) + N_alpha_t_1 * z_t;
    giB_t->CoordoB(2) = (*aiB_t)(2) + N_alpha_t_2 * z_t;
    if (tCalEpais) // cas 3D
          giB_t->CoordoB(3) = (*aiB_t)(3);
  };
// calcul de la base  naturel a tdt  
// neccessite le calcul prealable de la courbure et de la base naturelle et duale de la facette
void Met_Sfe1::Calcul_giB_tdt
 (  const Tableau<Noeud *>& , const  Mat_pleine& , int, const  Vecteur& phiH)
   {
    #ifdef MISE_AU_POINT	 	 
    if  (giB_tdt == NULL)
     { cout << "\nErreur : la base a t+dt n'est pas dimensionne !\n";
       cout << "void Met_Sfe1::Calcul_giB_tdt \n";
       Sortie(1);
     };
    #endif
	
    // vecteur de base gi
    double z_tdt = epais.epaisseur_tdt * 0.5 * (phiH(2)-phiH(1));
    // ici il n'y a pas le terme de variation de h par rapport au theta^alpha
    giB_tdt->CoordoB(1) = (*aiB_tdt)(1) + N_alpha_tdt_1 * z_tdt;
    giB_tdt->CoordoB(2) = (*aiB_tdt)(2) + N_alpha_tdt_2 * z_tdt;
    if (tCalEpais) // cas 3D
          giB_tdt->CoordoB(3) = (*aiB_tdt)(3);
  };
   
//------------// variation des vecteurs de base
// il faut au par avant  avoir calculé la variation de l'épaisseur si on est en 3D
 void Met_Sfe1::D_giB_t( const  Mat_pleine& , int ,  const Vecteur & phiH)
  { double z_t = epais.epaisseur_t * 0.5 * (phiH(2)-phiH(1));
    int nbddlx = tab_ddl.NbDdl_famille(X1);
	
    for (int iddl=1;iddl<= nbddlx;iddl++)
     { // vecteur de base gi
       (*d_giB_t)(iddl).CoordoB(1) = (*d_aiB_t)(iddl)(1) + (*d_N_alpha_t_1)(iddl) * z_t;
       (*d_giB_t)(iddl).CoordoB(2) = (*d_aiB_t)(iddl)(2) + (*d_N_alpha_t_2)(iddl) * z_t;
     };
    // cas 3D
    if (tCalEpais) 
      {// variation de giB3 = en fait celle de aiB3
       int nbddl_total = nbddlx + nbNoeudCentral;
       for (int iddl=1;iddl<= nbddl_total;iddl++)
           (*d_giB_t)(iddl).CoordoB(3) = (*d_aiB_t)(iddl)(3);
       // ensuite les vecteurs g1B et g2B dépendent aussi des ddl d'épaisseur
       // on calcul donc cette variation, et on en profite pour calculer également la variation de g3B
       // en fonction des ddl d'épaisseur 
       int ih=1;
       for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
        {(*d_giB_t)(iddl).CoordoB(1) = ((*d_epais)(ih).epaisseur_t * 0.5 * (phiH(2)-phiH(1))) * N_alpha_t_1;
         (*d_giB_t)(iddl).CoordoB(2) = ((*d_epais)(ih).epaisseur_t * 0.5 * (phiH(2)-phiH(1))) * N_alpha_t_2;
    	   };
      };  
  };

 void Met_Sfe1::D_giB_tdt( const  Mat_pleine& , int , const  Vecteur & phiH)
  { double z_tdt = epais.epaisseur_tdt * 0.5 * (phiH(2)-phiH(1));
    int nbddlx = tab_ddl.NbDdl_famille(X1);
	
    for (int iddl=1;iddl<= nbddlx;iddl++)
     {// vecteur de base gi
      (*d_giB_tdt)(iddl).CoordoB(1) = (*d_aiB_tdt)(iddl)(1) + (*d_N_alpha_tdt_1)(iddl) * z_tdt;
      (*d_giB_tdt)(iddl).CoordoB(2) = (*d_aiB_tdt)(iddl)(2) + (*d_N_alpha_tdt_2)(iddl) * z_tdt;

      #ifdef MISE_AU_POINT
       if (ParaGlob::NiveauImpression() > 5)
        {if  (((*d_giB_tdt)(iddl)(1).Norme() > 10000.) || ((*d_giB_tdt)(iddl)(1).Norme() < -10000.)
           || ((*d_giB_tdt)(iddl)(2).Norme() > 10000.) || ((*d_giB_tdt)(iddl)(2).Norme() < -10000.) )
          {cout << "\n attention erreur , les  d_giB_tdt sont trop grand "
                << "\n Met_Sfe1::D_giB_tdt(..."
                << " (*d_giB_tdt)(iddl)(1).Norme()=" << (*d_giB_tdt)(iddl)(1).Norme()
                << " (*d_giB_tdt)(iddl)(2).Norme()=" << (*d_giB_tdt)(iddl)(2).Norme()
                << " (*d_aiB_tdt)(iddl)(1).Norme()=" << (*d_aiB_tdt)(iddl)(1).Norme()
                << " (*d_aiB_tdt)(iddl)(2).Norme()=" << (*d_aiB_tdt)(iddl)(2).Norme()
                << " (*d_N_alpha_tdt_1)(iddl).Norme()=" << (*d_N_alpha_tdt_1)(iddl).Norme()
                << " (*d_N_alpha_tdt_2)(iddl).Norme()=" << (*d_N_alpha_tdt_2)(iddl).Norme()
                << endl ;
          };
        };
      #endif         
     };
    // cas 3D
    if (tCalEpais) 
      {// variation de giB3 = en fait celle de aiB3
       int nbddl_total = nbddlx + nbNoeudCentral;
       for (int iddl=1;iddl<= nbddl_total;iddl++)
           (*d_giB_tdt)(iddl).CoordoB(3) = (*d_aiB_tdt)(iddl)(3);
       // ensuite les vecteurs g1B et g2B dépendent aussi des ddl d'épaisseur
       // on calcul donc cette variation, et on en profite pour calculer également la variation de g3B
       // en fonction des ddl d'épaisseur
       int ih=1; 
       for (int iddl=nbddlx+1;iddl <= nbddl_total;iddl++,ih++)
        {(*d_giB_tdt)(iddl).CoordoB(1) = ((*d_epais)(ih).epaisseur_tdt * 0.5 * (phiH(2)-phiH(1))) * N_alpha_tdt_1;
         (*d_giB_tdt)(iddl).CoordoB(2) = ((*d_epais)(ih).epaisseur_tdt * 0.5 * (phiH(2)-phiH(1))) * N_alpha_tdt_2;
        };
      };
  };

//-----------// calcul du tenseur de courbure  dans la base naturelle
// on utilise un tenseur non symétrique, bien qu'en théorie la courbure est symétrique
// suivant les différentes méthodes de calcul, le tenseur en sortie est symétrique ou pas
// on peut toujours le symétriser, si l'on veut un tenseur symétrique
// on a : N,al = - b_{al be} . a^{be}
// plusieurs cas sont etudies suivant l'instant considere
  // a l'instant t = 0
Tenseur_ns2BB& Met_Sfe1::courbure_0 ( const Tableau<Noeud *>& tab_noeud
       ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
  { Tableau<Coordonnee3> tab_coor(6);
    // actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
    // du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
    for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
      tab_coor(i) = tab_noeud(i)->Coord0();
    for (int i=1;i< 4;i++) // les noeuds externes
     tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord0();
//    for (int i=1;i<=6;i++)
//      tab_coor(i) = tab_noeud(i)->Coord0();
    curb_0 =  (this->*PtCourbure)(tab_coor,(*aiB_0),(*aiH_0),tabTypeCL,vplan);
    return curb_0;
  };
// a l'instant t 
Tenseur_ns2BB& Met_Sfe1::courbure_t ( const Tableau<Noeud *>& tab_noeud
       ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
  { Tableau<Coordonnee3> tab_coor(6);
    // actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
    // du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
    for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
      tab_coor(i) = tab_noeud(i)->Coord1();
    for (int i=1;i< 4;i++) // les noeuds externes
     tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord1();
//    for (int i=1;i<=6;i++)
//      tab_coor(i) = tab_noeud(i)->Coord1();
    curb_t =  (this->*PtCourbure)(tab_coor,(*aiB_t),(*aiH_t),tabTypeCL,vplan);
    return curb_t;
  };
// a l'instant t+dt
Tenseur_ns2BB& Met_Sfe1::courbure_tdt ( const Tableau<Noeud *>& tab_noeud
       ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
  { Tableau<Coordonnee3> tab_coor(6);
    // actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
    // du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
    for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
      tab_coor(i) = tab_noeud(i)->Coord2();
    for (int i=1;i< 4;i++) // les noeuds externes
     tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord2();
//    for (int i=1;i<=6;i++)
//      tab_coor(i) = tab_noeud(i)->Coord2();
    curb_tdt = (this->*PtCourbure)(tab_coor,(*aiB_tdt),(*aiH_tdt),tabTypeCL,vplan);
    return curb_tdt;
  };
       
//--------------// 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 
 ( const Tableau<Noeud *>& tab_noeud,Tenseur_ns2BB& curb,TabOper<Tenseur_ns2BB>& dcurb
       ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
  { Tableau<Coordonnee3> tab_coor(6);
    // actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
    // du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
    for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
      tab_coor(i) = tab_noeud(i)->Coord1();
    for (int i=1;i< 4;i++) // les noeuds externes
     tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord1();
//    for (int i=1;i<=6;i++)
//      tab_coor(i) = tab_noeud(i)->Coord1();
    (this->*PtDcourbure) (tab_coor,(*aiB_t),*d_aiB_t,(*aiH_t),*d_aiH_t,curb,dcurb,tabTypeCL,vplan);
  };
// a l'instant t+dt
void Met_Sfe1::Dcourbure_tdt 
 ( const Tableau<Noeud *>& tab_noeud,Tenseur_ns2BB& curb,TabOper<Tenseur_ns2BB>& dcurb
       ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
  { Tableau<Coordonnee3> tab_coor(6);
    // actuellement les courbures s'obtiennent à l'aide des noeuds sommets du triangle
    // du milieu et des noeuds externes: on est pour l'instant dans un cadre de triangle
    for (int i=1;i< 4;i++) // les 3 noeuds sommet interne
      tab_coor(i) = tab_noeud(i)->Coord2();
    for (int i=1;i< 4;i++) // les noeuds externes
     tab_coor(i+3) = tab_noeud(i+nbNoeudCentral)->Coord2();
//    for (int i=1;i<=6;i++)
//      tab_coor(i) = tab_noeud(i)->Coord2();
    (this->*PtDcourbure) (tab_coor,(*aiB_tdt),*d_aiB_tdt,(*aiH_tdt),*d_aiH_tdt,curb,dcurb,tabTypeCL,vplan);
  };

// ----------- routines internes pour calculer les différentes formes de courbures 

// routine generale de calcul de la courbure
// cas 1: calcul des 3 courbures suivant les 3 cotés
Tenseur_ns2BB Met_Sfe1::Courbure1
  ( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH & aiH 
    ,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )

   { Vecteur curbp(3),curbi(3); // : les tenseurs 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 
     Coordonnee3 N(Util::ProdVec_coorBN(aiB(1),aiB(2)).Normer());
     // calcul du centre de gravite
     Coordonnee3 G((tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.);
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- on utilise indi un tableau d'indice pour ne pas etre embete par les modulos 3
        // 2- on boucle sur le cotes du triangle
     Coordonnee3 Ni(3); // normale du triangle exterieur
     // def des différentes grandeurs
     Coordonnee3 G2;
     Coordonnee3 H2G2,CB,U1,CG2,C_G,CA,dNdV;
     Coordonnee3 U2p;   
     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 (vérifiée)
        const  Coordonnee3& D = tab_coor(indi(ncot));
        const  Coordonnee3& B = tab_coor(indi(ncot+1));
        const  Coordonnee3& C = tab_coor(indi(ncot+2));
        const  Coordonnee3& A = tab_coor(indi(ncot)+3);
         if (D != A) // cas ou le noeud externe existe
          { // calcul de G2, centre de gravite du triangle exterieur           
           G2 = (B + C + A) / 3.;
           // calcul de U1                      
           CB = B - C;
           U1 = CB / CB.Norme(); 
           // vecteur CG2
           CG2 = G2 - C;
           // longueur h2G2 et g_h1 et calcul du vecteur H2G2 norme
           // tout d'abord H2G2 ce situe dans le plan de la facette externe
           H2G2 = CG2 - (CG2 * U1) * U1;
           double h2G2 = H2G2.Norme();
           double unSurh2G2= 1./h2G2;
           // maintenant on considere la direction H2G2, ramenée dans le plan de la facette centrale: U2p
           U2p = Util::ProdVec_coor(N,U1);                    
           C_G = G - C; // on met un _ car ce nom existe pour une autre routine
           double g_h1 = sqrt ( C_G * C_G - ( C_G * U1)* ( C_G * U1)) ; 
           // calcul de la normale du triangle exterieur
           CA = A - C;
           Ni =  (Util::ProdVec_coor(CB,CA)).Normer(); 
           // calcul de la derivee de la normale
           dNdV = ( Ni - N) / ( g_h1 + h2G2);
           // courbure dans la direction  U2p
           curbp(ncot) = - dNdV * U2p;
           }
         else // cas ou il n'y a pas de noeud externe 
          { // la courbure est nulle
            curbp(ncot) = 0.;
            // la direction de la courbure nulle est normale au cote
            // c-a-d le produit vectoriel de la normale avec le cote
            CB = B - C;
            double norCB = CB.Norme();
            U1 = CB / norCB;
            U2p = Util::ProdVec_coor(N,U1);
           }
         // calcul de la matrice de passage de Ualpha a aalpha
         double Beta21 = U2p * aiH.Coordo(1);
         double Beta22 = U2p * aiH.Coordo(2);
         // remplissage de la matrice mat
         mat(ncot,1) = Beta21*Beta21;
         mat(ncot,2) = 2.* Beta21*Beta22;
         mat(ncot,3) = Beta22*Beta22;
        };
	 
      // calcul du tenseur de courbure dans la base ai
      Mat_pleine mati = mat.Inverse();
      curbi = mati * curbp;
      // retour de la courbure
      Tenseur_ns2BB curb;
      curb.Coor(1,1)=curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(2); curb.Coor(2,2)=curbi(3);
      return curb;	 
   };
       
// routine generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure  b11,b12,b22
//            dcurb , la variation de courbure.
// On ne considere dans les variations que les termes relatif aux ddl qui font varier
// les normales
// cas 1: calcul des 3 courbures suivant les 3 cotés
void Met_Sfe1::Dcourbure1
 (  const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB
    , const Tableau <BaseB>& DaiB,const BaseH & aiH
    , const Tableau <BaseH>& DaiH,Tenseur_ns2BB& curb,TabOper <Tenseur_ns2BB>& dcurb
    ,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
   { int nbddl = (nbNoeudCentral+3)*3; //18;   // nombre de ddl total
     Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations
     Vecteur curbi(3); // la courbure sous forme de vecteur
     Vecteur curbp(3); // : les tenseurs de courbures  dans les axes
                       // locaux des arretes du triangle
     double untiers=1./3.;
     Vecteur vecNul(3); // un vecteur nul pour l'initialisation                  
     Tableau <Vecteur> dcurbp(nbddl,vecNul); // 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 
     Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
     double nor = NN.Norme(); // la norme
     Coordonnee3 N = NN / nor;
     Coordonnee3 DNN; // vecteur intermediaire
     // calcul du centre de gravite
     Coordonnee3 G ((tab_coor(1) + tab_coor(2) + tab_coor(3)) * untiers);
     // et le tableau de ses variations
     Tableau <Coordonnee3> dG(nbddl,Nul); 
     // calcul de la variation de la normale centrale et du centre de gravité: vérifiée
     Tableau <Coordonnee3> dN(nbddl,Nul);
     // on s'intéresse ici qu'aux ddl de la face centrale
     for (int inc=1;inc<=3;inc++) // indice du numéro de noeud de la face centrale
       for (int ib=1;ib<=3;ib++)
       	{ // cas de la normale
       	  int i1=(inc-1)*3+ib;
          DNN = Util::ProdVec_coorBN(DaiB(i1)(1),aiB(2)) + Util::ProdVec_coorBN(aiB(1),DaiB(i1)(2));
          dN(i1) = Util::VarUnVect_coor(NN,DNN,nor);
          // cas du centre de gravité
          dG(i1)=untiers*Ia(ib);
       	};
          
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- on utilise un tableau d'indice pour ne pas etre embete par les modulos 3
        // 2- on boucle sur le cotes du triangle
     Coordonnee3 NNi,Ni; // normale du triangle exterieur 
     // -------- declaration pour la boucle sur les 3 cotes ------------
     Coordonnee3 G2(3); // init a zero et a 3 composantes 
     Tableau <Coordonnee3> dB(nbddl,Nul),dC(nbddl,Nul),dA(nbddl,Nul),dG2(nbddl,Nul) ; // variation des points
     Vecteur d_g_h1(nbddl);
     Coordonnee3 CB,U1,CG2,H2G2,U2p;
     Tableau <Coordonnee3> dCB(nbddl,Nul),dU1(nbddl,Nul);
     Coordonnee3 C_G,CA,dNdV;
     Tableau <Coordonnee3> dCA(nbddl,Nul),dNi(nbddl,Nul);
     Coordonnee3 ddNdv; 
     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
       const  Coordonnee3& D = tab_coor(indi(ncot));
       const  Coordonnee3& B = tab_coor(indi(ncot+1));
       const  Coordonnee3& C = tab_coor(indi(ncot+2));
       const  Coordonnee3& A = tab_coor(indi(ncot)+3);
        if (D != A) // cas ou le noeud externe existe
        {// calcul de G2, centre de gravite du triangle exterieur           
         G2 = (B + C + A) * untiers;
         // calcul de U1                      
         CB = B - C;
         double norCB = CB.Norme();
         U1 = CB / norCB;
         // vecteur CG2
         CG2 = G2 - C;
         // longueur h2G2 et g_h1 et calcul du vecteur H2G2 norme
         // tout d'abord H2G2 ce citue dans le plan de la facette externe
         H2G2 = CG2 - (CG2 * U1) * U1;
         double h2G2 = H2G2.Norme();
         double unSurh2G2= 1./h2G2;
         // maintenant on considere la direction H2G2, ramenée dans le plan de la facette centrale: U2p
         U2p = Util::ProdVec_coor(N,U1);                    
         C_G = G - C;
         double g_h1 = sqrt ( C_G * C_G - ( C_G * U1) * ( C_G * U1) ) ;
         // calcul de la normale du triangle exterieur
         CA = A - C;
         NNi =  Util::ProdVec_coor(CB,CA);
         double norNi = NNi.Norme();
         Ni =  NNi/norNi;
         // calcul de la derivee de la normale
         dNdV = ( Ni - N) / ( g_h1 + h2G2);
         // courbure dans la direction  U2p
         curbp(ncot) = - dNdV * U2p;
         // calcul de la matrice de passage de Ualpha a aalpha
         double Beta21 = U2p * aiH.Coordo(1);
         double Beta22 = U2p * aiH.Coordo(2);
         // remplissage de la matrice mat
         mat(ncot,1) = Beta21*Beta21;
         mat(ncot,2) = 2.* Beta21*Beta22;
         mat(ncot,3) = Beta22*Beta22;
         
         // calcul des variations
           // 1 - initialisation des vecteurs inter
         ddNdv.Zero();
           // 2 - init en fct des ddl
         for (int idi=1; idi<= nbddl; idi++)
           { dCA(idi).Zero(); dCB(idi).Zero();dNi(idi).Zero();
             dB(idi).Zero();dC(idi).Zero();dA(idi).Zero(); 
             dU1(idi).Zero();d_g_h1.Zero();dG2(idi).Zero();
            }; 
           // 3 - variation non nulle des Ni pour certaine valeur d'indice   
         for (int ia=1;ia<= 3;ia++)
           { // on a 3 indices pour chaque ia, pour lesquels la variation est non nulle
             // dans cette boucle on n'utilise pas systematiquement les 4 indices
             int  jnB = (indi(ncot+1)-1)*3+ia; // vérifié
             int jnC = (indi(ncot+2)-1)*3+ia;
             int jnA = 9+(ncot-1)*3+ia; //vérifié
             // tout d'abord les variations des points
             dB(jnB) = Ia(ia);
             dC(jnC) = Ia(ia);dA(jnA) = Ia(ia);
             dG2(jnA) = untiers * Ia(ia);dG2(jnB) = untiers * Ia(ia);dG2(jnC) = untiers * Ia(ia);
             // vecteur CB
             dCB(jnB) = dB(jnB); dCB(jnC) = - dC(jnC);
             // Vecteur CA = A - C
             dCA(jnC) = -dC(jnC); dCA(jnA) = dA(jnA);
             // Vecteur Ni =  (Util::ProdVec_coor(CB,CA)).Normer()
             dNi(jnC) = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dCB(jnC),CA)+Util::ProdVec_coor(CB,dCA(jnC))),norNi);
             dNi(jnB) = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dCB(jnB),CA)+Util::ProdVec_coor(CB,dCA(jnB))),norNi);
             dNi(jnA) = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dCB(jnA),CA)+Util::ProdVec_coor(CB,dCA(jnA))),norNi);
             // Vecteur U1 = CB / norCB;
             dU1(jnC) = Util::VarUnVect_coor(CB,dCB(jnC),norCB);
             dU1(jnB) = Util::VarUnVect_coor(CB,dCB(jnB),norCB);
            };
         for (int ib=1;ib<= 3;ib++) // indice de coordonnées
           { // 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
             // jnb[0] -> D, jnb[1] -> B, jnb[2] -> C, jnb[3] -> A,
             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++)
               { // variation de g_h1
                 Coordonnee3 d_CG =dC(jnb[jb])-dG(jnb[jb]);
                 double d_CGsur2 = (d_CG * C_G);
                 double d_CG_U1 = d_CG * U1 + C_G * dU1(jnb[jb]);
                 double d_g_h1 = (d_CGsur2 - d_CG_U1)/g_h1;
                 // variation de h2G2
                 Coordonnee3 d_CG2 =dC(jnb[jb])-dG2(jnb[jb]);
                 Coordonnee3 d_H2G2 =d_CG2 - (d_CG2 * U1 + CG2 * dU1(jnb[jb])) * U1
                                           - (CG2 * U1) *dU1(jnb[jb]);
                 double d_h2G2 = (H2G2 * d_H2G2)/h2G2;
                 // vecteur dNdV = ( Ni - N) / ( g_h1 + h2G2)
                 ddNdv = ((dNi(jnb[jb]) - dN(jnb[jb])) 
                         - (Ni-N) * (d_g_h1 + d_h2G2) / ( g_h1 + h2G2)) / ( g_h1 + h2G2); 
                 // variation de U2p
                 Coordonnee3 d_U2p = Util::ProdVec_coor(dN(jnb[jb]),U1)+Util::ProdVec_coor(N,dU1(jnb[jb]));
                 // courbure  curbp(ncot) = dNdV * U2p;
                 dcurbp(jnb[jb])(ncot) = - (ddNdv * U2p + dNdV * d_U2p); 
                 // -- variation des composantes de mat
                 double d_Beta21 = d_U2p * aiH.Coordo(1) + U2p * DaiH(jnb[jb]).Coordo(1);
                 double d_Beta22 = d_U2p * aiH.Coordo(2) + U2p * DaiH(jnb[jb]).Coordo(2);
                 dmat(jnb[jb])(ncot,1)= 2. * Beta21 * d_Beta21;     
                 dmat(jnb[jb])(ncot,2)= 2. * (d_Beta21 * Beta22 + Beta21 * d_Beta22);     
                 dmat(jnb[jb])(ncot,3)= 2. * Beta22 * d_Beta22;     
                };             
            } 
         }
        else // cas ou il n'y a pas de noeud externe 
        { // la courbure est nulle
          curbp(ncot) = 0.;
          // la direction de la courbure nulle est normale au cote
          // c-a-d le produit vectoriel de la normale avec le cote
          CB = B - C;
          double norCB = CB.Norme();
          U1 = CB / norCB;
          H2G2 = Util::ProdVec_coor(N,U1);
          // calcul de la matrice de passage de Ualpha a aalpha
          double Beta21 = H2G2 * aiH.Coordo(1);
          double Beta22 = H2G2 * aiH.Coordo(2);
          // remplissage de la matrice mat
          mat(ncot,1) = Beta21*Beta21;
          mat(ncot,2) = 2.* Beta21*Beta22;
          mat(ncot,3) = Beta22*Beta22;
          
        }
      }; //-- fin de la boucle sur les trois cotés
      
     // calcul du tenseur de courbure dans la base ai
     Mat_pleine mati = mat.Inverse();
     curbi = mati * curbp;
     curb.Coor(1,1)=curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(2); curb.Coor(2,2)=curbi(3);
     // idem pour les variations
     Mat_pleine dmata(mat); // une matrice intermédiaire
     for (int ddl=1;ddl<=nbddl;ddl++)
         { // tout d'abord le calcul de la variation de l'inverse de mat
           dmata = -(mati * dmat(ddl) * mati);
           // puis la variation de la courbure
           curbi = mati * dcurbp(ddl) + dmata * curbp;
           Tenseur_ns2BB& dcur = dcurb(ddl);
           dcur.Coor(1,1)=curbi(1); dcur.Coor(1,2)=dcur.Coor(2,1)=curbi(2); dcur.Coor(2,2)=curbi(3);
         };
   };

// test si la courbure est anormalement trop grande
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1:  si tout est ok,
//        0:  une courbure trop grande a été détecté
int Met_Sfe1::Test_courbure_anormale1(double inf_normale,const Tableau<Coordonnee3>& tab_coor
                 , const BaseB & aiB , const BaseH & aiH)
 {// pour l'instant cette méthode n'est pas alimentée on le signale
  if (ParaGlob::NiveauImpression()> 8)
    cout << "\n *** attention, cette methode n'est pas implante pour la courbure 1 "
         << "\n Met_Sfe1::Test_courbure_anormale1(... "<<endl;
  return 1; // retour par défaut
 };

// routine generale de calcul de la courbure
// cas 2: moyenne des normales de part et autres des arrêtes
Tenseur_ns2BB Met_Sfe1::Courbure2
  ( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH &  
    ,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
   { Tenseur_ns2BB curb; // : les tenseurs de courbures d'abord dans les axes
     // au lieu d'utiliser des expressions condensées (en commentaire)
     // on utilise exactement la même procédure que pour dcurb
     // car sinon on a de très petites différences qui ensuite
     // conduisent à des différences non négligeables sur les def
     // calcul de la normale centrale
     // Coordonnee3 N((Util::ProdVec_coorBN(aiB(1),aiB(2))).Normer());
    
     // calcul de la normale centrale 
     Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
     double nor = NN.Norme(); // la norme
     Coordonnee3 N = NN / nor;
    
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- on utilise indi un tableau d'indice pour ne pas etre embete par les modulos 3
        // 2- on boucle sur le cotes du triangle
     Tableau <Coordonnee3> t_Ni(3); // les 3 normales du triangle exterieur
//cout << "\n  N= "<<N;
     Coordonnee3 NNi,Ni,NietN; // 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 (vérifiée)
        const  Coordonnee3& D = tab_coor(indi(ncot));
        const  Coordonnee3& B = tab_coor(indi(ncot+1));
        const  Coordonnee3& C = tab_coor(indi(ncot+2));
        const  Coordonnee3& A = tab_coor(indi(ncot)+3);
         if (D != A) // cas ou le noeud externe existe
          {// calcul de la normale du triangle exterieur
           Coordonnee3 CB = B - C;
           Coordonnee3 CA = A - C;
           //t_Ni(ncot) = ((Util::ProdVec_coor(CB,CA)).Normer() + N).Normer();
           
          NNi =  Util::ProdVec_coor(CB,CA);
          double norNi = NNi.Norme();
          Ni =  NNi/norNi;
          NietN=(Ni + N);
          double norNietN = NietN.Norme();
          t_Ni(ncot) = NietN / norNietN; 
           }
         else // cas ou il n'y a pas de noeud externe 
          { // on utilise la normale actuelle comme normale sur l'arrête
            t_Ni(ncot) = N;
           };
        };
      // calcul de la variation de la normale
      Coordonnee3 dN_1 = 2.*(t_Ni(1)-t_Ni(2));
      Coordonnee3 dN_2 = 2.*(t_Ni(1)-t_Ni(3));
      // courbure
      curb.Coor(1,1)=-dN_1 * aiB.Coordo(1);
      curb.Coor(1,2)=-dN_1 * aiB.Coordo(2);
      curb.Coor(2,1)=-dN_2 * aiB.Coordo(1);
      curb.Coor(2,2)=-dN_2 * aiB.Coordo(2);
      // retour de la courbure
      return curb;	 
   };
       
// routine generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure  b11,b12,b22
//            dcurb , la variation de courbure.
// On ne considere dans les variations que les termes relatif aux ddl qui font varier
// les normales
// cas 2: moyenne des normales de part et autres des arrêtes
void Met_Sfe1::Dcourbure2
 (  const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB
    , const Tableau <BaseB>& DaiB,const BaseH & 
    , const Tableau <BaseH>& ,Tenseur_ns2BB& curb,TabOper <Tenseur_ns2BB>& dcurb
    ,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
   { int nbddl = (nbNoeudCentral+3)*3; //18;   // nombre de ddl total
     Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations
     
     double untiers=1./3.;                  
     // calcul de la normale centrale 
     Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
     double nor = NN.Norme(); // la norme
     Coordonnee3 N = NN / nor;
     Coordonnee3 DNN; // vecteur intermediaire
     // calcul de la variation de la normale centrale et du centre de gravité: vérifiée
     Tableau <Coordonnee3> dN(nbddl,Nul);
     // on s'intéresse ici qu'aux ddl de la face centrale
     for (int inc=1;inc<=3;inc++) // indice du numéro de noeud de la face centrale
       for (int ib=1;ib<=3;ib++)
       	{ // cas de la normale
       	  int i1=(inc-1)*3+ib;
          DNN = Util::ProdVec_coorBN(DaiB(i1)(1),aiB(2)) 
               + Util::ProdVec_coorBN(aiB(1),DaiB(i1)(2));
          dN(i1) = Util::VarUnVect_coor(NN,DNN,nor);
       	};
          
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- on utilise un tableau d'indice pour ne pas etre embete par les modulos 3
        // 2- on boucle sur le cotes du triangle
     Coordonnee3 NNi,Ni,NietN; // normale du triangle exterieur 
     // -------- declaration pour la boucle sur les 3 cotes ------------
     Coordonnee3 CB,CA;
     Tableau <Coordonnee3> t_Ni(3); // les 3 normales du triangle exterieur
     Tableau <Tableau <Coordonnee3> > d_t_Ni(nbddl,t_Ni); // le tableau des variations
     Coordonnee3 dA,dB,dC,dNi; 
//cout << "\n NN= "<<NN <<" nor= "<< nor <<" N= "<<N;
     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
       const  Coordonnee3& D = tab_coor(indi(ncot));
       const  Coordonnee3& B = tab_coor(indi(ncot+1));
       const  Coordonnee3& C = tab_coor(indi(ncot+2));
       const  Coordonnee3& A = tab_coor(indi(ncot)+3);
        if (D != A) // cas ou le noeud externe existe
        {CB = B - C;
         double norCB = CB.Norme();
         // calcul de la normale du triangle exterieur
         CA = A - C;
         NNi =  Util::ProdVec_coor(CB,CA);
         double norNi = NNi.Norme();
         Ni =  NNi/norNi;
         NietN=(Ni + N);
         double norNietN = NietN.Norme();
         t_Ni(ncot) = NietN / norNietN; 
//cout << "\n CB= "<<CB <<" CA= "<< CA
//     << "t_Ni("<<ncot<<")= "<<t_Ni(ncot);
//cout << "\n norNi= "<< norNi<< " norNietN= "<< norNietN;
         
         // calcul des variations
         for (int ia=1;ia<= 3;ia++)
          {{// indice de A: jnA
//            int jnA = 9+(ncot-1)*3+ia;
//modif pour généraliser             int jnA = 9+(ncot-1)*3+ia; //vérifié
            int jnA = nbNoeudCentral*3+(ncot-1)*3+ia;
            dA = Ia(ia);
            //dCA = dA;
            // Vecteur Ni =  (Util::ProdVec_coor(CB,CA)).Normer()
            dNi = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(CB,dA)),norNi);
            d_t_Ni(jnA)(ncot) = Util::VarUnVect_coor(NietN,(dNi+dN(jnA)),norNietN);
           }
           {// indice de B: jnB
            int  jnB = (indi(ncot+1)-1)*3+ia;
            dB = Ia(ia);
            //dCB = dB;
            // Vecteur Ni =  (Util::ProdVec_coor(CB,CA)).Normer()
            dNi = Util::VarUnVect_coor(NNi,(Util::ProdVec_coor(dB,CA)),norNi);
            d_t_Ni(jnB)(ncot) = Util::VarUnVect_coor(NietN,(dNi+dN(jnB)),norNietN);
           }
           {// indice de C: jnC
            int jnC = (indi(ncot+2)-1)*3+ia;
            dC = Ia(ia);
            //dCB = -dC; dCA= -dC;
            // Vecteur Ni =  (Util::ProdVec_coor(CB,CA)).Normer()
            dNi = Util::VarUnVect_coor(NNi,-(Util::ProdVec_coor(dC,CA)+Util::ProdVec_coor(CB,dC)),norNi);
            d_t_Ni(jnC)(ncot) = Util::VarUnVect_coor(NietN,(dNi+dN(jnC)),norNietN);
           }
           {// indice de D: jnD
            int jnD = (indi(ncot)-1)*3+ia;
            // Vecteur Ni =  (Util::ProdVec_coor(CB,CA)).Normer()
            // dNi = vecteur nul;
            d_t_Ni(jnD)(ncot) = Util::VarUnVect_coor(NietN,(dN(jnD)),norNietN);
           }
          }; // fin du for sur ia
        }
        else // cas ou il n'y a pas de noeud externe 
        { // on utilise la normale actuelle comme normale sur l'arrête
          t_Ni(ncot) = N;
//cout << "\n t_Ni("<<ncot<<")= "<<t_Ni(ncot);
          // et pour les variations, celle de N
          for (int ia=1;ia<= 3;ia++)
           { // indice de B: jnB
             int  jnB = (indi(ncot+1)-1)*3+ia;
             d_t_Ni(jnB)(ncot) = dN(jnB);
             // indice de C: jnC
             int jnC = (indi(ncot+2)-1)*3+ia;
             d_t_Ni(jnC)(ncot) = dN(jnC);
             // indice de D: jnD
             int jnD = (indi(ncot)-1)*3+ia;
             d_t_Ni(jnD)(ncot) = dN(jnD);
           }; // fin du for sur ia
         };
        
      }; //-- fin de la boucle sur les trois cotés
      
     // calcul de la variation de la normale
     Coordonnee3 dN_1 = 2.*(t_Ni(1)-t_Ni(2));
     Coordonnee3 dN_2 = 2.*(t_Ni(1)-t_Ni(3));
     // courbure
     curb.Coor(1,1)=-dN_1 * aiB.Coordo(1);
     curb.Coor(1,2)=-dN_1 * aiB.Coordo(2);
     curb.Coor(2,1)=-dN_2 * aiB.Coordo(1);
     curb.Coor(2,2)=-dN_2 * aiB.Coordo(2);
//cout << "\n 2) curb= "; curb.Ecriture(cout);
//cout << endl;
    
     // calcul des variations de la courbure par rapport au ddl
     for (int ddl=1;ddl<=nbddl;ddl++)
     	{Tenseur_ns2BB& dcur = dcurb(ddl);
     	 dcur.Coor(1,1) = (-2.) * ((d_t_Ni(ddl)(1)-d_t_Ni(ddl)(2)) * aiB.Coordo(1))  - dN_1 * DaiB(ddl).Coordo(1);
       dcur.Coor(1,2)= -(2.* (d_t_Ni(ddl)(1)-d_t_Ni(ddl)(2)) * aiB.Coordo(2) +      (dN_1 * DaiB(ddl).Coordo(2)));
       dcur.Coor(2,1)= -(2.*(d_t_Ni(ddl)(1)-d_t_Ni(ddl)(3)) * aiB.Coordo(1) +      (dN_2 * DaiB(ddl).Coordo(1)));
     	 dcur.Coor(2,2) = (-2.) * ((d_t_Ni(ddl)(1)-d_t_Ni(ddl)(3)) * aiB.Coordo(2)) - dN_2 * DaiB(ddl).Coordo(2);
     	};
   };

// test si la courbure est anormalement trop grande
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1:  si tout est ok,
//        0:  une courbure trop grande a été détecté
int Met_Sfe1::Test_courbure_anormale2(double inf_normale,const Tableau<Coordonnee3>& tab_coor
                 , const BaseB & aiB , const BaseH & aiH)
 {// pour l'instant cette méthode n'est pas alimentée on le signale
  if (ParaGlob::NiveauImpression()> 8)
    cout << "\n *** attention, cette methode n'est pas implante pour la courbure 2 "
         << "\n Met_Sfe1::Test_courbure_anormale2(... "<<endl;
  return 1; // retour par défaut
 };

// routine (3) generale de calcul de la courbure
// cas 3: moyenne pondérée des normales de part et autres des arrêtes, et localisation du point sur l'arrête
Tenseur_ns2BB Met_Sfe1::Courbure3
  ( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH & aiH 
    ,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )

   { // calcul de la normale centrale 
     Coordonnee3 N((Util::ProdVec_coorBN(aiB(1),aiB(2))).Normer());
     // calcul du centre de gravite
     Coordonnee3 G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) / 3.;
     // les coordonnées locales non nulles des points origines des normales sur les arrêtes
     // on les initialise comme si les points étaient au milieu des arrêtes
     double theta1_H1=0.5; double theta2_H1=0.5; 
     double theta2_H2=0.5; double theta1_H3=0.5;
     // origine
     const Coordonnee3& Op=tab_coor(1);// origine du repère local
     
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- on utilise indi un tableau d'indice pour ne pas etre embete par les modulos 3
        // 2- on boucle sur le cotes du triangle
     Tableau <Coordonnee3> t_Ni(3); // normale du triangle exterieur
     // def des différentes grandeurs
     Coordonnee3 Ne; // normale courante du triangle extérieur
     Coordonnee3 CB,CA,CD; Coordonnee3 G2;
     Coordonnee3 U1,HaA,OpHi;  
     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 (vérifiée)
        const  Coordonnee3& D = tab_coor(indi(ncot));
        const  Coordonnee3& B = tab_coor(indi(ncot+1));
        const  Coordonnee3& C = tab_coor(indi(ncot+2));
        const  Coordonnee3& A = tab_coor(indi(ncot)+3);
         if (D != A) // cas ou le noeud externe existe
          {// calcul de G2, centre de gravite du triangle exterieur           
           G2 = (B + C + A) / 3.;
           // calcul de U1                      
           CB = B - C;
           CA = A - C;
           U1 = CB / CB.Norme(); 
           // longueur HaA et HdD et calcul du vecteur HaA norme
           // tout d'abord HaA ce situe dans le plan de la facette externe
           HaA = CA - (CA * U1) * U1;
           double haA = HaA.Norme();
           double unSurhaA= 1./haA;
           CD = D - C; 
           double dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ; 
           // calcul de la normale du triangle exterieur
           Ne =  (Util::ProdVec_coor(CB,CA)).Normer(); 
           // calcul de la normale
           t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
           // calcul des coordonnées de la position des normales sur les arrêtes
           double beta = 0.5 * (G + G2 - 2. * B) * U1;
           OpHi=B+beta*U1-Op;
           switch (ncot)
           	{case 1: theta1_H1=OpHi*aiH.Coordo(1);theta2_H1=OpHi*aiH.Coordo(2); break;
           	 case 2: theta2_H2=OpHi*aiH.Coordo(2); break;
           	 case 3: theta1_H3=OpHi*aiH.Coordo(1); break;
           	};
           }
         else // cas ou il n'y a pas de noeud externe 
          {// on utilise comme normale extérieure la normale du triangle central
           t_Ni(ncot) = N; 
           // on utilise le point milieu de l'arrête pour le positionnement de la normale
           // normalement les points sont déjà initialisés
          };
       };
	 
     // calcul des coefficients des fonctions d'interpolation
     double s= theta1_H1 * theta2_H2 - theta1_H3 * theta2_H2 + theta2_H1 * theta1_H3;
     double a_r1= theta2_H2 / s; double b_r1 = theta1_H3 / s; double c_r1 = -theta1_H3 * theta2_H2/s;
     double a_r2= -theta2_H1 / s; double b_r2 = 1./theta2_H2 - theta1_H3 * theta2_H1 /theta2_H2/ s; 
                                               double c_r2 = -theta1_H3 * theta2_H1/s;
     double a_r3= 1./theta1_H3 - theta1_H1 * theta2_H2 /theta1_H3/ s; double b_r3 = - theta1_H1 / s; 
                                               double c_r3 = theta2_H2 * theta1_H1 /s;
     // calcul de la variation de la normale
     Coordonnee3 d_N_1= t_Ni(1) * a_r1 + t_Ni(2) * a_r2 + t_Ni(3) * a_r3 ;
     Coordonnee3 d_N_2= t_Ni(1) * b_r1 + t_Ni(2) * b_r2 + t_Ni(3) * b_r3 ;
     // calcul de la courbure
     Tenseur_ns2BB curb;
     curb.Coor(1,1) = - d_N_1 * aiB.Coordo(1); curb.Coor(1,2) = - d_N_1 * aiB.Coordo(2);
     curb.Coor(2,1) = - d_N_2 * aiB.Coordo(1); curb.Coor(2,2) = - d_N_2 * aiB.Coordo(2);
     // retour de la courbure
     return curb;	 
   };

// test si la courbure est anormalement trop grande
// inf_normale : indique en entrée le det mini pour la courbure en locale
// retour 1:  si tout est ok,
//        0:  une courbure trop grande a été détecté
int Met_Sfe1::Test_courbure_anormale3(double inf_normale,const Tableau<Coordonnee3>& tab_coor
                 , const BaseB & aiB , const BaseH & aiH)
 {// pour l'instant cette méthode n'est pas alimentée on le signale
  if (ParaGlob::NiveauImpression()> 8)
    cout << "\n *** attention, cette methode n'est pas implante pour la courbure 3 "
         << "\n Met_Sfe1::Test_courbure_anormale1(... "<<endl;
  return 1; // retour par défaut
 };

// routine (3) generale de calcul de la courbure et de sa variation
// en sortie : curb , la courbure  b11,b12,b22
//            dcurb , la variation de courbure.
// On ne considere dans les variations que les termes relatif aux ddl qui font varier
// les normales
// cas 3: moyenne pondérée des normales de part et autres des arrêtes, et localisation du point sur l'arrête
void Met_Sfe1::Dcourbure3
 (  const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB
    , const Tableau <BaseB>& DaiB,const BaseH & aiH
    , const Tableau <BaseH>& DaiH,Tenseur_ns2BB& curb,TabOper <Tenseur_ns2BB>& dcurb
    ,Tableau <EnuTypeCL> const & ,Tableau <Coordonnee3> const & )
   { int nbddl = (nbNoeudCentral+3)*3; //18;   // nombre de ddl total
     Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations

     // les coordonnées locales non nulles des points origines des normales sur les arrêtes
     // on les initialise comme si les points étaient au milieu des arrêtes
     double theta1_H1=0.5; double theta2_H1=0.5; 
     double theta2_H2=0.5; double theta1_H3=0.5;
     // origine
     const Coordonnee3& Op=tab_coor(1);// origine du repère local
     Tableau <Coordonnee3> dOp(nbddl,Nul);
     dOp(1)=Ia(1);dOp(2)=Ia(2);dOp(3)=Ia(3);
     // et les variations
     Vecteur d_theta1_H1(nbddl), d_theta2_H1(nbddl), d_theta2_H2(nbddl), d_theta1_H3(nbddl); 
     double untiers=1./3.;                  
     // calcul du centre de gravite
     Coordonnee3 G = (tab_coor(1) + tab_coor(2) + tab_coor(3)) * untiers;
     // et le tableau de ses variations
     Tableau <Coordonnee3> dG(nbddl,Nul); 

     // calcul de la normale centrale 
     Coordonnee3 NN(Util::ProdVec_coorBN(aiB(1),aiB(2))); // normale non normee
     double nor = NN.Norme(); // la norme
     Coordonnee3 N = NN / nor;
     Coordonnee3 DNN; // vecteur intermediaire
     // calcul de la variation de la normale centrale et du centre de gravité: vérifiée
     Tableau <Coordonnee3> dN(nbddl,Nul);
     // on s'intéresse ici qu'aux ddl de la face centrale
     for (int inc=1;inc<=3;inc++) // indice du numéro de noeud de la face centrale
       for (int ib=1;ib<=3;ib++)
       	{ // cas de la normale
       	  int i1=(inc-1)*3+ib;
          DNN = Util::ProdVec_coorBN(DaiB(i1)(1),aiB(2)) 
               + Util::ProdVec_coorBN(aiB(1),DaiB(i1)(2));
          dN(i1) = Util::VarUnVect_coor(NN,DNN,nor);
          // cas du centre de gravité
          dG(i1)=untiers*Ia(ib);
       	};
          
     // le meme calcul est effectue pour les trois cotes du triangle principale
        //1- on utilise un tableau d'indice pour ne pas etre embete par les modulos 3
        // 2- on boucle sur le cotes du triangle
     Tableau <Coordonnee3> t_Ni(3); // normale du triangle exterieur
     Tableau < Tableau <Coordonnee3> > d_t_Ni(nbddl,t_Ni); // les dérivées, initialisées à celle de dN
     Coordonnee3 NNi; // normale du triangle exterieur 
     // def des différentes grandeurs
     Coordonnee3 Ne,NNe; // normale courante du triangle extérieur et sa valeur non normée
     Coordonnee3 CB,CA,CD; Coordonnee3 G2,dG2;
     Coordonnee3 U1,HaA,OpHi; Coordonnee3  BgimoinsGB;
     // -------- declaration pour la boucle sur les 3 cotes ------------     
     Coordonnee3 dA,dB,dC,dD,d_U1,dNe,d_HaA,d_NNi,d_Op,d_OpHi;
     Coordonnee3 t_NNi; // normale courante non normée sur l'arrête
     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
       const  Coordonnee3& D = tab_coor(indi(ncot));
       const  Coordonnee3& B = tab_coor(indi(ncot+1));
       const  Coordonnee3& C = tab_coor(indi(ncot+2));
       const  Coordonnee3& A = tab_coor(indi(ncot)+3);
        if (D != A) // cas ou le noeud externe existe
        {// calcul de G2, centre de gravite du triangle exterieur           
         G2 = (B + C + A) * untiers;
         // calcul de U1                      
         CB = B - C;
         CA = A - C;
         double norCB = CB.Norme();
         U1 = CB / norCB;
         // longueur HaA et HdD et calcul du vecteur HaA norme
         // tout d'abord HaA ce situe dans le plan de la facette externe
         HaA = CA - (CA * U1) * U1;
         double haA = HaA.Norme();
         double unSurhaA= 1./haA;
         CD = D - C; 
         double dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
//         double distH = dHd + haA;
//         double inv_distH = 1./distH; double inv_distH2 =  inv_distH * inv_distH;
         // calcul de la normale du triangle exterieur
         NNe =  (Util::ProdVec_coor(CB,CA)); 
         double norNe = NNe.Norme();
         Ne =  NNe/norNe; 
         // calcul de la normale
         t_NNi = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA)));
         double nor_t_NNi = t_NNi.Norme();
         t_Ni(ncot) = t_NNi / nor_t_NNi;
         // calcul des coordonnées de la position des normales sur les arrêtes
         BgimoinsGB = (G + G2 - 2. * B);
         double beta = 0.5 * BgimoinsGB * U1;
         OpHi=B+beta*U1-Op;
         switch (ncot)
         	{case 1: theta1_H1=OpHi*aiH.Coordo(1);theta2_H1=OpHi*aiH.Coordo(2); break;
         	 case 2: theta2_H2=OpHi*aiH.Coordo(2); break;
         	 case 3: theta1_H3=OpHi*aiH.Coordo(1); break;
         	}; 
         
         // calcul des variations
         for (int ia=1;ia<= 3;ia++)
          {{// indice de A: jnA, noeud extérieur
// modif pour généraliser            int jnA = 9+(ncot-1)*3+ia;
            int jnA = nbNoeudCentral*3+(ncot-1)*3+ia;
            dA = Ia(ia);
            //dCA = dA;
            // Vecteur Ne =  (Util::ProdVec_coor(CB,CA)).Normer()
            dNe = Util::VarUnVect_coor(NNe,(Util::ProdVec_coor(CB,dA)),norNe);
            // ici d_CB=0 donc d_U1=0
            // ici d_CD=0 donc d_dHd=0
            d_HaA = dA - (dA * U1) * U1;
            double d_haA = HaA * d_HaA * unSurhaA ;
            // t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
            d_NNi = (dNe*(dHd/(dHd+haA)) + dN(jnA)*(haA/(dHd+haA)))
                     + (Ne*(-dHd*d_haA/(dHd+haA)/(dHd+haA)) + N*(d_haA/(dHd+haA)))
                     + ( N*(-haA*d_haA/(dHd+haA)/(dHd+haA)));
            d_t_Ni(jnA)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
            // maintenant la variation des theta
            // G2 = (B + C + A) * untiers;
            dG2 = dA * untiers;
            //  beta = 0.5 * (G + G2 - 2. * B) * U1;
            // ici dG(jnA) = 0
            double d_beta = 0.5 * (dG2 * U1);
            d_OpHi=d_beta*U1;  // OpHi=B+beta*U1-Op;
            switch (ncot)
         	 {case 1: {d_theta1_H1(jnA)=d_OpHi*aiH.Coordo(1);d_theta2_H1(jnA)=d_OpHi*aiH.Coordo(2); break;}
         	  case 2:  d_theta2_H2(jnA)=d_OpHi*aiH.Coordo(2); break;
         	  case 3:  d_theta1_H3(jnA)=d_OpHi*aiH.Coordo(1); break;
         	};
           }
           {// indice de B: jnB, noeud de la facette centrale
            int  jnB = (indi(ncot+1)-1)*3+ia;
            dB = Ia(ia); // d_CB = dB, d_CA=0
            d_U1= Util::VarUnVect_coor(CB,dB,norCB);
            // HaA = CA - (CA * U1) * U1;
            d_HaA =  - ((CA*d_U1) * U1 +(CA * U1) * d_U1);
            double d_haA = HaA * d_HaA * unSurhaA ;
            // d_CD = 0; dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ;
            double d_dHd = - ( CD * d_U1)* ( CD * U1) / dHd ; 
            // Vecteur Ne =  (Util::ProdVec_coor(CB,CA)).Normer()
            dNe = Util::VarUnVect_coor(NNe,(Util::ProdVec_coor(dB,CA)),norNe);
            // t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
            d_NNi = (dNe*(dHd/(dHd+haA)) + dN(jnB)*(haA/(dHd+haA)))
                     + (Ne*(d_dHd/(dHd+haA)) + N*(d_haA/(dHd+haA)))
                     + (Ne*(-dHd*d_haA/(dHd+haA)/(dHd+haA))+ N*(-haA*d_haA/(dHd+haA)/(dHd+haA)))
                     + (Ne*(-dHd*d_dHd/(dHd+haA)/(dHd+haA)) + N*(-haA*d_dHd/(dHd+haA)/(dHd+haA)));
            d_t_Ni(jnB)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
            // maintenant la variation des theta
            // G2 = (B + C + A) * untiers;
            dG2 = dB * untiers;
            //  beta = 0.5 * (G + G2 - 2. * B) * U1;
            double d_beta = 0.5 * ((dG(jnB) + dG2 - 2. * dB) * U1 + BgimoinsGB * d_U1);
 //           if (ncot == 3) {d_Op=dB;} else { d_Op= Nul;};
//            d_OpHi=dB + d_beta*U1+beta*d_U1-d_Op;
            d_OpHi=dB + d_beta*U1+beta*d_U1-dOp(jnB);
            
            switch (ncot)
         	 {case 1: {d_theta1_H1(jnB)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnB).Coordo(1);
         	           d_theta2_H1(jnB)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnB).Coordo(2); break; }
         	  case 2:  d_theta2_H2(jnB)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnB).Coordo(2); break;
         	  case 3:  d_theta1_H3(jnB)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnB).Coordo(1); break;
         	};
           }
           {// indice de C: jnC
            int jnC = (indi(ncot+2)-1)*3+ia;
            dC = Ia(ia);
            //dCB = -dC; dCA= -dC; U1 = CB / norCB;
            d_U1= Util::VarUnVect_coor(CB,-dC,norCB);
            // HaA = CA - (CA * U1) * U1;
            d_HaA = -dC - ((-dC * U1)+(CA*d_U1)) * U1 -(CA * U1) * d_U1;
            double d_haA = HaA * d_HaA * unSurhaA ;
            // d_CD = -dC; dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ; 
            double d_dHd = - (CD*dC + ( CD * d_U1- dC*U1) * ( CD * U1) ) / dHd ; 
            // Vecteur Ne =  (Util::ProdVec_coor(CB,CA)).Normer()
            dNe = Util::VarUnVect_coor(NNe,(-(Util::ProdVec_coor(dC,CA)+Util::ProdVec_coor(CB,dC))),norNe);
            // t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
            d_NNi = (dNe*(dHd/(dHd+haA)) + dN(jnC)*(haA/(dHd+haA)))
                     + (Ne*(d_dHd/(dHd+haA)) + N*(d_haA/(dHd+haA)))
                     + (Ne*(-dHd*d_haA/(dHd+haA)/(dHd+haA)) + N*(-haA*d_haA/(dHd+haA)/(dHd+haA)))
                     + (Ne*(-dHd*d_dHd/(dHd+haA)/(dHd+haA)) + N*(-haA*d_dHd/(dHd+haA)/(dHd+haA)));
            d_t_Ni(jnC)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
            // maintenant la variation des theta
            // G2 = (B + C + A) * untiers;
            dG2 = dC * untiers;
            //  beta = 0.5 * (G + G2 - 2. * B) * U1;
            double d_beta = 0.5 * ((dG(jnC) + dG2 ) * U1 + BgimoinsGB * d_U1);
//            if (ncot == 2) {d_Op=dC;} else { d_Op= Nul;};
//            d_OpHi=d_beta*U1+beta*d_U1-d_Op;
            d_OpHi=d_beta*U1+beta*d_U1-dOp(jnC);
            
            switch (ncot)
         	 {case 1: {d_theta1_H1(jnC)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnC).Coordo(1);
         	           d_theta2_H1(jnC)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnC).Coordo(2); break; }
         	  case 2: d_theta2_H2(jnC)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnC).Coordo(2); break;
         	  case 3: d_theta1_H3(jnC)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnC).Coordo(1); break;
         	};
           }
           {// indice de D: jnD
            int jnD = (indi(ncot)-1)*3+ia;
            // d_CB = 0; d_CA = 0 ==> d_U1=0 
            // d_HaA = 0, et d_haA = 0
            // d_CD = dD; 
            // dHd = sqrt ( CD * CD - ( CD * U1)* ( CD * U1)) ; 
            double d_dHd = (CD*dD - (dD * U1) * ( CD * U1) ) / dHd ;
            // Vecteur Ne =  (Util::ProdVec_coor(CB,CA)).Normer()
            //dNe = 0
            // t_Ni(ncot) = (Ne*(dHd/(dHd+haA)) + N*(haA/(dHd+haA))).Normer();
            d_NNi = (dN(jnD)*(haA/(dHd+haA)))
                     + (Ne*(d_dHd/(dHd+haA)))
                     + (Ne*(-dHd*d_dHd/(dHd+haA)/(dHd+haA)) + N*(-haA*d_dHd/(dHd+haA)/(dHd+haA)));
            d_t_Ni(jnD)(ncot) = Util::VarUnVect_coor(t_NNi,d_NNi,nor_t_NNi);
            // maintenant la variation des theta
            // G2 = (B + C + A) * untiers;
            //dG2 = 0;
            //  beta = 0.5 * (G + G2 - 2. * B) * U1;
            double d_beta = 0.5 * (dG(jnD) * U1);
//            if (ncot == 1) {d_Op=dD;} else { d_Op= Nul;};
            // dB=0
//            d_OpHi=d_beta*U1-d_Op;
            d_OpHi=d_beta*U1-dOp(jnD);
            
            switch (ncot)
         	 {case 1: {d_theta1_H1(jnD)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnD).Coordo(1);
         	           d_theta2_H1(jnD)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnD).Coordo(2); break; }
         	  case 2:  d_theta2_H2(jnD)=d_OpHi*aiH.Coordo(2)+OpHi*DaiH(jnD).Coordo(2); break;
         	  case 3:  d_theta1_H3(jnD)=d_OpHi*aiH.Coordo(1)+OpHi*DaiH(jnD).Coordo(1); break;
         	};
           }
          }; // fin du for sur ia
         
         }
        else // cas ou il n'y a pas de noeud externe 
        { // on utilise la normale actuelle comme normale sur l'arrête
          t_Ni(ncot) = N;
          // et pour les variations, celle de N
          for (int ia=1;ia<= 3;ia++)
           { // indice de B: jnB
             int  jnB = (indi(ncot+1)-1)*3+ia;
             d_t_Ni(jnB)(ncot) = dN(jnB);
             // indice de C: jnC
             int jnC = (indi(ncot+2)-1)*3+ia;
             d_t_Ni(jnC)(ncot) = dN(jnC);
             // indice de D: jnD
             int jnD = (indi(ncot)-1)*3+ia;
             d_t_Ni(jnD)(ncot) = dN(jnD);
           }; // fin du for sur ia
          // concernant les variations des coordonnées des positions des normales,
          // ces variations sont nulles, on n'y touche pas
         };
      }; //-- fin de la boucle sur les trois cotés
      
	 
     // calcul des coefficients des fonctions d'interpolation
     double s= theta1_H1 * theta2_H2 - theta1_H3 * theta2_H2 + theta2_H1 * theta1_H3;
     double a_r1= theta2_H2 / s; double b_r1 = theta1_H3 / s; double c_r1 = -theta1_H3 * theta2_H2/s;
     double a_r2= -theta2_H1 / s; double b_r2 = 1./theta2_H2 - theta1_H3 * theta2_H1 /theta2_H2/ s; 
                                               double c_r2 = -theta1_H3 * theta2_H1/s;
     double a_r3= 1./theta1_H3 - theta1_H1 * theta2_H2 /theta1_H3/ s; double b_r3 = - theta1_H1 / s; 
                                               double c_r3 = theta2_H2 * theta1_H1 /s;
     // calcul de la variation de la normale
     Coordonnee3 d_N_1= t_Ni(1) * a_r1 + t_Ni(2) * a_r2 + t_Ni(3) * a_r3 ;
     Coordonnee3 d_N_2= t_Ni(1) * b_r1 + t_Ni(2) * b_r2 + t_Ni(3) * b_r3 ;
     // calcul de la courbure (non symétrique)
     curb.Coor(1,1) = - d_N_1 * aiB.Coordo(1); curb.Coor(1,2) = - d_N_1 * aiB.Coordo(2);
     curb.Coor(2,1) = - d_N_2 * aiB.Coordo(1); curb.Coor(2,2) = - d_N_2 * aiB.Coordo(2);
     // retour de la courbure

     double unsurs = 1./s; double unsurs2=unsurs * unsurs ;
     for (int ddl=1;ddl<=nbddl;ddl++)
      { // tout d'abord le calcul de la variation des coeff des fonctions d'interpolation
        double d_s = d_theta1_H1(ddl) * theta2_H2 + theta1_H1 * d_theta2_H2(ddl)
                   - d_theta1_H3(ddl) * theta2_H2 - theta1_H3 * d_theta2_H2(ddl)
                   + d_theta2_H1(ddl) * theta1_H3 + theta2_H1 * d_theta1_H3(ddl);
        double d_a_r1 = d_theta2_H2(ddl) * unsurs - theta2_H2 * d_s * unsurs2; 
        double d_b_r1 = d_theta1_H3(ddl) * unsurs - theta1_H3 * d_s * unsurs2; 
        double d_a_r2 = -d_theta2_H1(ddl) * unsurs + theta2_H1 * d_s * unsurs2; 
        double d_b_r2 = - d_theta2_H2(ddl)/theta2_H2/theta2_H2 
                        - d_theta1_H3(ddl) * theta2_H1 /theta2_H2* unsurs 
                        - theta1_H3 * d_theta2_H1(ddl) /theta2_H2* unsurs
                        + theta1_H3 * theta2_H1 * d_theta2_H2(ddl)/theta2_H2/theta2_H2 * unsurs 
                        + theta1_H3 * theta2_H1 * d_s /theta2_H2 * unsurs2; 
        double d_a_r3 = -d_theta1_H3(ddl)/theta1_H3/theta1_H3 
                        - d_theta1_H1(ddl) * theta2_H2 /theta1_H3* unsurs 
                        - theta1_H1 * d_theta2_H2(ddl) /theta1_H3* unsurs
                        + theta1_H1 * theta2_H2 * d_theta1_H3(ddl)/theta1_H3/theta1_H3 * unsurs 
                        + theta1_H1 * theta2_H2 * d_s/theta1_H3 * unsurs2; 
        double d_b_r3 = - d_theta1_H1(ddl) * unsurs  + theta1_H1 * d_s * unsurs2; 
        // calcul de la variation de la normale
        Coordonnee3 d_d_N_1=  d_t_Ni(ddl)(1) * a_r1 + t_Ni(1) * d_a_r1
                            + d_t_Ni(ddl)(2) * a_r2 + t_Ni(2) * d_a_r2 
                            + d_t_Ni(ddl)(3) * a_r3 + t_Ni(3) * d_a_r3;
        Coordonnee3 d_d_N_2=  d_t_Ni(ddl)(1) * b_r1 + t_Ni(1) * d_b_r1 
                            + d_t_Ni(ddl)(2) * b_r2 + t_Ni(2) * d_b_r2
                            + d_t_Ni(ddl)(3) * b_r3 + t_Ni(3) * d_b_r3;
        // variation de la courbure
        Tenseur_ns2BB& dcur = dcurb(ddl);
        dcur.Coor(1,1) = - (d_d_N_1 * aiB.Coordo(1)+ d_N_1 * DaiB(ddl).Coordo(1));
        dcur.Coor(1,2) = - (d_d_N_1 * aiB.Coordo(2)- d_N_1 * DaiB(ddl).Coordo(2));
        dcur.Coor(2,1) = - (d_d_N_2 * aiB.Coordo(1)- d_N_2 * DaiB(ddl).Coordo(2));
        dcur.Coor(2,2) = - (d_d_N_2 * aiB.Coordo(2)- d_N_2 * DaiB(ddl).Coordo(1));
      }; 

   };