// 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 "Util.h"
#include "MathUtil.h"
#include "MathUtil2.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
//  =========================================================================================
		  

// routine (4) generale de calcul de la courbure
// cas 4: calcul de la courbure à partir d'un polynôme passant par les 6 points
Tenseur_ns2BB Met_Sfe1::Courbure4
 ( const Tableau<Coordonnee3>& tab_coor, const BaseB & aiB , const BaseH & aiH 
                 ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)

   { // calcul de la normale centrale
     // on utilise des expressions strictement identique ici et dans dcourbure
     // de manière à éviter les erreurs d'arrondis
     //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;
     Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
     Vecteur sm(3); // le second membre pour le calcul de la courbure
     // 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
     Coordonnee3 U1,CB,V,OpA;
     Coordonnee3 Mm,OpMm;  
     
     Coordonnee3 OpAp; // le vecteur horizontal correspondant à OpA
     for (int ncot=1;ncot<=3;ncot++)
       { // on utilise les notations : triangle principal DBC et secondaire BAC
          // 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);
         // examen de condition limite éventuelle
         // ncot % 3 +1: = le coté en face du noeud ncot
         if (tabTypeCL(ncot % 3 +1) != TANGENTE_CL)
          { // cas sans conditions limites particulières
            if (D != A) // cas ou le noeud externe existe
             {OpA = A - Op; 
              double theta3 = OpA * N; // cote du point externe
              OpAp = OpA - theta3 * N; // projection de OpA sur la facette 
              double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2); // les coordonnées locales de A
              // construction de la ligne de la matrice
              mat(ncot,1) = theta1 * (theta1 - 1.);
              mat(ncot,2) = theta2 * (theta2 - 1.);
              mat(ncot,3) = theta1 * theta2;
              // et du second membre correspondant
              sm(ncot) = theta3;
             }
            else // cas ou il n'y a pas de noeud externe
             {// On va calculer la normale non normée au coté externe
              // ce vecteur est dans le plan de la facette
              CB = B - C;
              U1 = CB ;
              V =  Util::ProdVec_coor(N,U1);
              //  maintenant on calcul les coordonnées locales de V
              double theta1 = V * aiH.Coordo(1); double theta2 = V * aiH.Coordo(2);
              // construction de la ligne de la matrice
              mat(ncot,1) = theta1 * theta1;
              mat(ncot,2) = theta2 * theta2;
              mat(ncot,3) = theta1 * theta2;
              // et du second membre correspondant
              sm(ncot) = 0;
             };
          }
         else
          { // cas avec une  condition de tangente imposée
            // on commence par définir la direction de la tangente
            CB = B - C;
            V =  Util::ProdVec_coor(CB,vplan(ncot % 3 +1));
            // calcul des coordonnées locales du vecteur tangent
            double d1 = V*aiH.Coordo(1); double d2=V*aiH.Coordo(2); double d3 = V*N;
            // on calcul les coordonnées du point milieu de BC que l'on nome Mm
            Mm = 0.5*(B+C);
            OpMm = Mm - Op;
            // les coordonnées locales de Mm
            double theta1 = OpMm * aiH.Coordo(1); double theta2 = OpMm * aiH.Coordo(2); 
            // construction de la ligne de la matrice
            mat(ncot,1) = (2.*theta1 - 1.) * d1;
            mat(ncot,2) = (2.*theta2 - 1.) * d2;
            mat(ncot,3) = (theta1 * d2 + theta2 * d1);
            // et du second membre correspondant
            sm(ncot) = d3;
          };
          
       };
           
     // calcul du tenseur de courbure dans la base ai
     Mat_pleine mati = mat.Inverse();
     Vecteur curbi = mati * sm;
     // retour de la courbure
     Tenseur_ns2BB curb;
     curb.Coor(1,1)=2. * curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(3); curb.Coor(2,2)= 2. * curbi(2);
     return curb;	 
   };
       
// routine (4) 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 4: calcul de la courbure à partir d'un polynôme passant par les 6 points
void Met_Sfe1::Dcourbure4
(  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 & tabTypeCL,Tableau <Coordonnee3> const & vplan)
   { int nbddl = (nbNoeudCentral+3)*3; //18;   // nombre de ddl total
     Coordonnee3 Nul; // le vecteur nul de dim 3, pour les initialisations

     Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
     Tableau <Mat_pleine> dmat(nbddl,mat);  // sa variation
     Vecteur sm(3); // le second membre pour le calcul de la courbure
     Tableau <Vecteur> dsm(nbddl,sm);  // sa variation
     
     // 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 : 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);
       	};

     // 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 autres dérivées sont nulles
          
     // 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 CB,U1,V,dV;
     Coordonnee3 OpA,OpAp;
     Coordonnee3 dA,dOpA,dOpAp,dB,dC;
     Coordonnee3 Mm,dMm,OpMm,dOpMm; 
     
     for (int ncot=1;ncot<=3;ncot++)
      { // on utilise les notations : triangle principal DBC et secondaire BAC
           // 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);
       // examen de condition limite éventuelle
       int nbcote = ncot % 3 +1; // le numéro du coté en face du noeud ncot
       if (tabTypeCL(nbcote) != TANGENTE_CL)
        { // cas sans conditions limites particulières
          if (D != A) // cas ou le noeud externe existe
           {OpA = A - Op; 
            double theta3 = OpA * N; // cote du point externe
            OpAp = OpA - theta3 * N; // projection de OpA sur la facette 
            // les coordonnées locales de A
            double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2); 
            // construction de la ligne de la matrice
            mat(ncot,1) = theta1 * (theta1 - 1.);
            mat(ncot,2) = theta2 * (theta2 - 1.);
            mat(ncot,3) = theta1 * theta2;
            // et du second membre correspondant
            sm(ncot) = theta3;
            // calcul des variations
            for (int ia=1;ia<= 3;ia++)
             {{// indice de A: jnA, noeud extérieur
               int jnA = 9+(ncot-1)*3+ia;
               dA = Ia(ia);
               dOpA = dA; //  dOp(jnA)=0 // OpA = A - Op;
               // la coordonnée en 3
               double dtheta3 = dOpA * N; // dN(jnA)=0 car noeud externe // theta3 = OpA * N 
               dOpAp = dOpA - dtheta3 * N ; //dN(jnA)=0; // OpAp = OpA - theta3 * N 
               // les coordonnées locales
               double dtheta1 = dOpAp * aiH.Coordo(1); // DaiH(jnA).Coordo(1)=0; // theta1 = OpAp * aiH.Coordo(1)
               double dtheta2 = dOpAp * aiH.Coordo(2);// DaiH(jnA).Coordo(2)=0; // theta2 = OpAp * aiH.Coordo(2)
               // la matrice
               dmat(jnA)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
               dmat(jnA)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
               dmat(jnA)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
               // le second membre
               dsm(jnA)(ncot) = dtheta3; // sm(ncot) = theta3
              }
              {// indice de B: jnB, noeud de la facette centrale
               int  jnB = (indi(ncot+1)-1)*3+ia;
               // dA = 0 ;
               dOpA =  - dOp(jnB); // OpA = A - Op;
               // la coordonnée en 3
               double dtheta3 = dOpA * N + OpA * dN(jnB); // theta3 = OpA * N 
               dOpAp = dOpA - dtheta3 * N - theta3 * dN(jnB); // OpAp = OpA - theta3 * N 
               // les coordonnées locales
     // 1) à noter ici que OpAp est un vecteur de la facette, donc normale à N, et que DaiH ne contient que les
     // variations plane de aiH (il manque les variations suivant la normale) mais cette partie suffit car la partie
     // qui manque est suivant N et donc normale à OpAp (donc elle disparait au niveau du produit scalaire)          
               double dtheta1 = dOpAp * aiH.Coordo(1) + OpAp * DaiH(jnB).Coordo(1); // theta1 = OpAp * aiH.Coordo(1)
               double dtheta2 = dOpAp * aiH.Coordo(2) + OpAp * DaiH(jnB).Coordo(2); // theta2 = OpAp * aiH.Coordo(2)
               // la matrice
               dmat(jnB)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
               dmat(jnB)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
               dmat(jnB)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
               // le second membre
               dsm(jnB)(ncot) = dtheta3; // sm(ncot) = theta3
              }
              {// indice de C: jnC
               int jnC = (indi(ncot+2)-1)*3+ia;
               // dA = 0 ;
               dOpA =  - dOp(jnC); // OpA = A - Op;
               // la coordonnée en 3
               double dtheta3 = dOpA * N + OpA * dN(jnC); // theta3 = OpA * N 
               dOpAp = dOpA - dtheta3 * N - theta3 * dN(jnC); // OpAp = OpA - theta3 * N 
               // les coordonnées locales
     // mêmes remarques 1) que plus haut          
               double dtheta1 = dOpAp * aiH.Coordo(1) + OpAp * DaiH(jnC).Coordo(1); // theta1 = OpAp * aiH.Coordo(1)
               double dtheta2 = dOpAp * aiH.Coordo(2) + OpAp * DaiH(jnC).Coordo(2); // theta2 = OpAp * aiH.Coordo(2)
               // la matrice
               dmat(jnC)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
               dmat(jnC)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
               dmat(jnC)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
               // le second membre
               dsm(jnC)(ncot) = dtheta3; // sm(ncot) = theta3
              }
              {// indice de D: jnD
               int jnD = (indi(ncot)-1)*3+ia;
               // dA = 0 ;
               dOpA =  - dOp(jnD); // OpA = A - Op;
               // la coordonnée en 3
               double dtheta3 = dOpA * N + OpA * dN(jnD); // theta3 = OpA * N 
               dOpAp = dOpA - dtheta3 * N - theta3 * dN(jnD); // OpAp = OpA - theta3 * N 
               // les coordonnées locales
     // mêmes remarques 1) que plus haut          
               double dtheta1 = dOpAp * aiH.Coordo(1) + OpAp * DaiH(jnD).Coordo(1); // theta1 = OpAp * aiH.Coordo(1)
               double dtheta2 = dOpAp * aiH.Coordo(2) + OpAp * DaiH(jnD).Coordo(2); // theta2 = OpAp * aiH.Coordo(2)
               // la matrice
               dmat(jnD)(ncot,1)= dtheta1 * (theta1 - 1.) + theta1 * dtheta1; // mat(ncot,1) = theta1 * (theta1 - 1.)
               dmat(jnD)(ncot,2)= dtheta2 * (theta2 - 1.) + theta2 * dtheta2; // mat(ncot,2) = theta2 * (theta2 - 1.)
               dmat(jnD)(ncot,3)= dtheta1 * theta2 + theta1 * dtheta2; // mat(ncot,3) = theta1 * theta2
               // le second membre
               dsm(jnD)(ncot) = dtheta3; // sm(ncot) = theta3
              }
             }; // fin du for sur ia
           }
          else // cas ou il n'y a pas de noeud externe
          { // On va calculer la normale non normée au coté externe
            // ce vecteur est dans le plan de la facette
            CB = B - C;
            U1 = CB ;
            V =  Util::ProdVec_coor(N,U1);
            //  maintenant on calcul les coordonnées locales de V
            double theta1 = V * aiH.Coordo(1); double theta2 = V * aiH.Coordo(2);
            // construction de la ligne de la matrice
            mat(ncot,1) = theta1 * theta1;
            mat(ncot,2) = theta2 * theta2;
            mat(ncot,3) = theta1 * theta2;
            // et du second membre correspondant
            sm(ncot) = 0;
            // calcul des variations
            for (int ia=1;ia<= 3;ia++)
             {
//              {// indice de A: jnA, noeud extérieur
// comme il n'y a pas de noeud extérieur et que N ne dépend pas des ddl extérieur, 
// normalement jnA n'intervient pas !!              
//               int jnA = 9+(ncot-1)*3+ia;
//               // dCB = 0.; // CB = B - C
//               // dU1 = 0.; // U1 = CB 
//               dV = Util::ProdVec_coor(dN(jnA),U1);  // V = Util::ProdVec_coor(N,U1)
//               double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnA).Coordo(1); // theta1 = V * aiH.Coordo(1);
//               double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnA).Coordo(2); // theta2 = V * aiH.Coordo(2);
//               dmat(jnA)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
//               dmat(jnA)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
//               dmat(jnA)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
//               dsm(jnA)(ncot) = 0.; // sm(ncot) = 0
//              }
              {// indice de B: jnB, noeud de la facette centrale
               int  jnB = (indi(ncot+1)-1)*3+ia;
               dB = Ia(ia); 
               // d_CB = dB // CB = B - C
               // dU1 = dB; // U1 = CB 
               dV = Util::ProdVec_coor(dN(jnB),U1) + Util::ProdVec_coor(N,dB);  // V = Util::ProdVec_coor(N,U1)
     // 2) à noter ici que V est un vecteur de la facette, donc normale à N, et que DaiH ne contient que les
     // variations plane de aiH (il manque les variations suivant la normale) mais cette partie suffit car la partie
     // qui manque est suivant N et donc normale à V (donc elle disparait au niveau du produit scalaire)          
               double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnB).Coordo(1); // theta1 = V * aiH.Coordo(1);
               double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnB).Coordo(2); // theta2 = V * aiH.Coordo(2);
               dmat(jnB)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
               dmat(jnB)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
               dmat(jnB)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
               dsm(jnB)(ncot) = 0.; // sm(ncot) = 0
            	}
              {// indice de C: jnC
               int jnC = (indi(ncot+2)-1)*3+ia;
               dC = Ia(ia);
               //dCB = -dC; // CB = B - C
               // dU1 = -dC; // U1 = CB 
               dV = Util::ProdVec_coor(dN(jnC),U1) + Util::ProdVec_coor(N,-dC);  // V = Util::ProdVec_coor(N,U1)
      // mêmes remarques 2) que plus haut
               double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnC).Coordo(1); // theta1 = V * aiH.Coordo(1);
               double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnC).Coordo(2); // theta2 = V * aiH.Coordo(2);
               dmat(jnC)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
               dmat(jnC)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
               dmat(jnC)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
               dsm(jnC)(ncot) = 0.; // sm(ncot) = 0
            	}
              {// indice de D: jnD
               int jnD = (indi(ncot)-1)*3+ia;
               // dCB = 0.; // CB = B - C
               // dU1 = 0.; // U1 = CB 
               dV = Util::ProdVec_coor(dN(jnD),U1);  // V = Util::ProdVec_coor(N,U1)
      // mêmes remarques 2) que plus haut
               double dtheta1 = dV * aiH.Coordo(1) + V * DaiH(jnD).Coordo(1); // theta1 = V * aiH.Coordo(1);
               double dtheta2 = dV * aiH.Coordo(2) + V * DaiH(jnD).Coordo(2); // theta2 = V * aiH.Coordo(2);
               dmat(jnD)(ncot,1) = 2.* dtheta1 * theta1; // mat(ncot,1) = theta1 * theta1
               dmat(jnD)(ncot,2) = 2.* dtheta2 * theta2; // mat(ncot,2) = theta2 * theta2
               dmat(jnD)(ncot,3) = (dtheta1 * theta2 + theta1 * dtheta2); // mat(ncot,3) = 2. * theta1 * theta2
               dsm(jnD)(ncot) = 0.; // sm(ncot) = 0 
             }
            }; // fin du for sur ia 
          }; // fin du cas où il n'y a pas de noeud externe
        }
       else
        { // cas avec une  condition de tangente imposée
          // on commence par définir la direction de la tangente
          CB = B - C;
          V =  Util::ProdVec_coor(CB,vplan(nbcote));
          // on calcul les coordonnées du point milieu de BC que l'on nome Mm
          Mm = 0.5*(B+C);
          OpMm = Mm - Op;
          // calcul des coordonnées locales du vecteur tangent
          double d1 = V*aiH.Coordo(1); double d2=V*aiH.Coordo(2); double d3 = V*N;
          // les coordonnées locales de Mm
          double theta1 = OpMm * aiH.Coordo(1); double theta2 = OpMm * aiH.Coordo(2);
// les coordonnée sont ici des coordonnées matérielles qui dans le cas linéaires sont constantes (0.5,0) ou 
// (0.5,0.5) ou (0,0.5) leur variations est donc nulles 
// ***** donc on peut simplifier les expressions qui suivent en annulant les variations de theta_alpha          
          
           
          // construction de la ligne de la matrice
          mat(ncot,1) = (2.*theta1 - 1.) * d1;
          mat(ncot,2) = (2.*theta2 - 1.) * d2;
          mat(ncot,3) = (theta1 * d2 + theta2 * d1);
          // et du second membre correspondant
          sm(ncot) = d3;
          // calcul des variations
          for (int ia=1;ia<= 3;ia++)
          { {// indice de B: jnB, noeud de la facette centrale
             int  jnB = (indi(ncot+1)-1)*3+ia;
             dB = Ia(ia); 
             // d_CB = dB // CB = B - C
             dV = Util::ProdVec_coor(dB,vplan(nbcote));  // V = Util::ProdVec_coor(CB,vplan(nbcote))
      // mêmes remarques 2) que plus haut
             double dd1=dV*aiH.Coordo(1)+V*DaiH(jnB).Coordo(1); // d1 = V*aiH.Coordo(1)
             double dd2=dV*aiH.Coordo(2)+V*DaiH(jnB).Coordo(2); // d2=V*aiH.Coordo(2)
             double dd3=dV*N+V*dN(jnB);                // d3 = V*N
             // dMm = 0.5 * dB ;  // Mm = 0.5*(B+C);
             dOpMm = 0.5 * dB - dOp(jnB); // OpMm = Mm - Op;
             // les coordonnées locales
      // ici également, OpMm est dans le plan de la facette d'où mêmes remarques 2) que plus haut
             double dtheta1 = dOpMm * aiH.Coordo(1) + OpMm * DaiH(jnB).Coordo(1); // theta1 = OpMm * aiH.Coordo(1)
             double dtheta2 = dOpMm * aiH.Coordo(2) + OpMm * DaiH(jnB).Coordo(2); // theta2 = OpMm * aiH.Coordo(2)
             dmat(jnB)(ncot,1) = (2.*dtheta1-1.)*d1+(2.*theta1-1.)*dd1; // mat(ncot,1) = (2.*theta1 - 1.) * d1
             dmat(jnB)(ncot,2) = (2.*dtheta2-1.)*d2+(2.*theta2-1.)*dd2; // mat(ncot,1) = (2.*theta1 - 1.) * d1
             dmat(jnB)(ncot,3) = (dtheta1*d2+dtheta2*d1)+(theta1*dd2+theta2*dd1); // mat(ncot,3) = (theta1 * d2 + theta2 * d1);
             dsm(jnB)(ncot) = dd3; // sm(ncot) = d3
          	}
            {// indice de C: jnC
             int jnC = (indi(ncot+2)-1)*3+ia;
             dC = Ia(ia);
             //dCB = -dC; // CB = B - C
             dV = Util::ProdVec_coor(-dC,vplan(nbcote));  // V = Util::ProdVec_coor(CB,vplan(nbcote))
             double dd1=dV*aiH.Coordo(1)+V*DaiH(jnC).Coordo(1); // d1 = V*aiH.Coordo(1)
             double dd2=dV*aiH.Coordo(2)+V*DaiH(jnC).Coordo(2); // d2=V*aiH.Coordo(2)
             double dd3=dV*N+V*dN(jnC);                // d3 = V*N
             // dMm = 0.5 * dC ;  // Mm = 0.5*(B+C);
             dOpMm = 0.5 * dC - dOp(jnC); // OpMm = Mm - Op;
             // les coordonnées locales
      // ici également, OpMm est dans le plan de la facette d'où mêmes remarques 2) que plus haut
             double dtheta1 = dOpMm * aiH.Coordo(1) + OpMm * DaiH(jnC).Coordo(1); // theta1 = OpMm * aiH.Coordo(1)
             double dtheta2 = dOpMm * aiH.Coordo(2) + OpMm * DaiH(jnC).Coordo(2); // theta2 = OpMm * aiH.Coordo(2)
             dmat(jnC)(ncot,1) = (2.*dtheta1-1.)*d1+(2.*theta1-1.)*dd1; // mat(ncot,1) = (2.*theta1 - 1.) * d1
             dmat(jnC)(ncot,2) = (2.*dtheta2-1.)*d2+(2.*theta2-1.)*dd2; // mat(ncot,1) = (2.*theta1 - 1.) * d1
             dmat(jnC)(ncot,3) = (dtheta1*d2+dtheta2*d1)+(theta1*dd2+theta2*dd1); // mat(ncot,3) = (theta1 * d2 + theta2 * d1);
             dsm(jnC)(ncot) = dd3; // sm(ncot) = d3
          	}
            {// indice de D: jnD
             int jnD = (indi(ncot)-1)*3+ia;
             // dCB = 0.; // CB = B - C
             dV = 0. ;  // V = Util::ProdVec_coor(CB,vplan(nbcote))
             double dd1=V*DaiH(jnD).Coordo(1); // d1 = V*aiH.Coordo(1)
             double dd2=V*DaiH(jnD).Coordo(2); // d2=V*aiH.Coordo(2)
             double dd3=V*dN(jnD);                // d3 = V*N
             // dMm = 0.5 * dB ;  // Mm = 0.5*(B+C);
             dOpMm =  - dOp(jnD); // OpMm = Mm - Op;
             // les coordonnées locales
      // ici également, OpMm est dans le plan de la facette d'où mêmes remarques 2) que plus haut
             double dtheta1 = dOpMm * aiH.Coordo(1) + OpMm * DaiH(jnD).Coordo(1); // theta1 = OpMm * aiH.Coordo(1)
             double dtheta2 = dOpMm * aiH.Coordo(2) + OpMm * DaiH(jnD).Coordo(2); // theta2 = OpMm * aiH.Coordo(2)
             dmat(jnD)(ncot,1) = (2.*dtheta1-1.)*d1+(2.*theta1-1.)*dd1; // mat(ncot,1) = (2.*theta1 - 1.) * d1
             dmat(jnD)(ncot,2) = (2.*dtheta2-1.)*d2+(2.*theta2-1.)*dd2; // mat(ncot,1) = (2.*theta1 - 1.) * d1
             dmat(jnD)(ncot,3) = (dtheta1*d2+dtheta2*d1)+(theta1*dd2+theta2*dd1); // mat(ncot,3) = (theta1 * d2 + theta2 * d1);
             dsm(jnD)(ncot) = dd3; // sm(ncot) = d3
           }
          }; // fin du for sur ia
        }; // fin du test sur les CL: condition de tangence
      }; // fin de la boucle sur les 3 noeuds externes
      
     // calcul du tenseur de courbure dans la base ai
     Mat_pleine mati = mat.Inverse();
    
//---- pour le debug des cas particuliers ----
#ifdef MISE_AU_POINT
     if (ParaGlob::NiveauImpression() > 3)
      { int i,j;
        if (mati.MaxiValAbs(i,j) > ConstMath::grand)
         { cout << "\n *** Met_Sfe1::Dcourbure4, mati.MaxiValAbs(i,j) est trop grand  :\n mat= ";
           mat.Affiche();
           cout << "\n inverse de mat: mati= ";
           mati.Affiche();
           // affichage des éléments du triangle
           cout << "\n pour chaque cote : on utilise les notations : triangle principal DBC et     secondaire BAC";
           cout << "\n Op= "<< Op << "\n N= "<< N ;
           for (int ncot=1;ncot<=3;ncot++)
             { // on utilise les notations : triangle principal DBC et secondaire BAC
                  // simplification de l'ecriture
               cout << "\n cote : "<< ncot ;
               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);
               cout << "\n  D= "<< D << "\n B= " << B << "\n C= "<< C << "\n A= "<< A;
               cout << "\n aiH(1)= "<<aiH.Coordo(1) << "\n aiH(2)= "<<aiH.Coordo(2);
              
               int nbcote = ncot % 3 +1; // le numéro du coté en face du noeud ncot
               if (tabTypeCL(nbcote) != TANGENTE_CL) // seul cas envisagé pour l'instant
                 { // cas sans conditions limites particulières
                   if (D != A) // cas ou le noeud externe existe
                    {OpA = A - Op; 
                     double theta3 = OpA * N; // cote du point externe
                     OpAp = OpA - theta3 * N; // projection de OpA sur la facette
                     cout << "\n  OpA = A - Op= "<< OpA
                          << "\n OpAp = OpA - theta3 * N= " << OpAp ;
                     // les coordonnées locales de A
                     double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2);
                     cout << "\n theta1= OpAp * aiH(1)= "<< theta1 << " theta2= OpAp * aiH(2)= " << theta2 << " ";
                     cout << "\n  mat(ncot,1) = theta1 * (theta1 - 1.) ";
                     cout << "\n  mat(ncot,2) = theta2 * (theta2 - 1.)";
                     cout << "\n  mat(ncot,3) = theta1 * theta2";
                     cout << "\n  sm(ncot) = theta3";
                    };
                 };
             };
           cout << endl;
         };
      };
#endif
//---- fin pour le debug des cas particuliers ----
    
     Vecteur curbi = mati * sm;
     // retour de la courbure
     curb.Coor(1,1)=2. * curbi(1); curb.Coor(1,2)=curb.Coor(2,1)=curbi(3); curb.Coor(2,2)= 2. * curbi(2);
     // idem pour les variations
     Coordonnee3 dcurbi;
//     Mat_pleine dmata(mat); // une matrice intermédiaire
     for (int ddl=1;ddl<=nbddl;ddl++)
         { dcurbi = mati * (dsm(ddl)  - dmat(ddl) * curbi );
           Tenseur_ns2BB& dcur = dcurb(ddl);
// ----- pour le débug
//           cout << "\n dcurbi(1),(" << ddl << ")= " << dcurbi(1) 
//                << " dcurbi(2),(" << ddl << ")= " << dcurbi(2)
//                << " dcurbi(3),(" << ddl << ")= " << dcurbi(3) << endl;
// ----- fin pour le débug           
           dcur.Coor(1,1)=2. * dcurbi(1); dcur.Coor(1,2)=dcur.Coor(2,1)=dcurbi(3); dcur.Coor(2,2)= 2. * dcurbi(2);
         };
// ----- pour le débug
//         Sortie(1);
// ----- fin pour le débug           
   };

// 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_anormale4(double inf_normale, const Tableau<Coordonnee3>& tab_coor
                 , const BaseB & aiB , const BaseH & aiH)
   { // calcul de la normale centrale
     // on utilise des expressions strictement identique ici et dans dcourbure
     // de manière à éviter les erreurs d'arrondis
     //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;
     Mat_pleine mat(3,3); // matrice pour calculer le tenseur de courbure
     Vecteur sm(3); // le second membre pour le calcul de la courbure
     // 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
     Coordonnee3 U1,CB,V,OpA;
     Coordonnee3 Mm,OpMm;  
     
     Coordonnee3 OpAp; // le vecteur horizontal correspondant à OpA
     for (int ncot=1;ncot<=3;ncot++)
       { // on utilise les notations : triangle principal DBC et secondaire BAC
          // 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);
         // pas d'examen de condition limite éventuelle
         // on ne considère que le cas sans conditions limites particulières
         // ncot % 3 +1: = le coté en face du noeud ncot
          { // cas sans conditions limites particulières
            if (D != A) // cas ou le noeud externe existe
             {OpA = A - Op; 
              double theta3 = OpA * N; // cote du point externe
              OpAp = OpA - theta3 * N; // projection de OpA sur la facette 
              double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2); // les coordonnées locales de A
              // construction de la ligne de la matrice
              mat(ncot,1) = theta1 * (theta1 - 1.);
              mat(ncot,2) = theta2 * (theta2 - 1.);
              mat(ncot,3) = theta1 * theta2;
              // et du second membre correspondant
              sm(ncot) = theta3;
             }
            else // cas ou il n'y a pas de noeud externe
             {// On va calculer la normale non normée au coté externe
              // ce vecteur est dans le plan de la facette
              CB = B - C;
              U1 = CB ;
              V =  Util::ProdVec_coor(N,U1);
              //  maintenant on calcul les coordonnées locales de V
              double theta1 = V * aiH.Coordo(1); double theta2 = V * aiH.Coordo(2);
              // construction de la ligne de la matrice
              mat(ncot,1) = theta1 * theta1;
              mat(ncot,2) = theta2 * theta2;
              mat(ncot,3) = theta1 * theta2;
              // et du second membre correspondant
              sm(ncot) = 0;
             };
          }
       };
           
     // le principe du test est de regarder les valeurs de mat qui doit-être une matrice
     // non singulière
     double det = mat.Determinant();
//--- debug
#ifdef MISE_AU_POINT
//cout << "\n debug Met_Sfe1::Test_courbure_anormale4 "
//     << " determinant= "<< det;
//if (Dabs(det) < ConstMath::petit)
     if ((Dabs(det) < inf_normale)&& (ParaGlob::NiveauImpression() > 5))
      { cout << "\n *** Met_Sfe1::Test_courbure_anormale4, determinant trop petit = "<< det
             << " \n mat= ";
        mat.Affiche();
        // affichage des éléments du triangle
        cout << "\n pour chaque cote : on utilise les notations : triangle principal DBC et     secondaire BAC";
        cout << "\n Op= "<< Op << "\n N= "<< N ;
        for (int ncot=1;ncot<=3;ncot++)
             { // on utilise les notations : triangle principal DBC et secondaire BAC
                  // simplification de l'ecriture
               cout << "\n cote : "<< ncot ;
               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);
               cout << "\n  D= "<< D << "\n B= " << B << "\n C= "<< C << "\n A= "<< A;
               cout << "\n aiH(1)= "<<aiH.Coordo(1) << "\n aiH(2)= "<<aiH.Coordo(2);
              
               int nbcote = ncot % 3 +1; // le numéro du coté en face du noeud ncot
                 { // cas sans conditions limites particulières
                   if (D != A) // cas ou le noeud externe existe
                    {OpA = A - Op; 
                     double theta3 = OpA * N; // cote du point externe
                     OpAp = OpA - theta3 * N; // projection de OpA sur la facette
                     cout << "\n  OpA = A - Op= "<< OpA
                          << "\n OpAp = OpA - theta3 * N= " << OpAp ;
                     // les coordonnées locales de A
                     double theta1 = OpAp * aiH.Coordo(1); double theta2 = OpAp * aiH.Coordo(2);
                     cout << "\n theta1= OpAp * aiH(1)= "<< theta1 << " theta2= OpAp * aiH(2)= " << theta2 << " ";
                     cout << "\n  mat(ncot,1) = theta1 * (theta1 - 1.) ";
                     cout << "\n  mat(ncot,2) = theta2 * (theta2 - 1.)";
                     cout << "\n  mat(ncot,3) = theta1 * theta2";
                     cout << "\n  sm(ncot) = theta3";
                    };
                 };
             };
        cout << endl;
      };
#endif
//--- fin debug
    
     int retour = 1; // init par défaut
//     if (det < ConstMath::petit)
     if (Dabs(det) < inf_normale)
       retour = 0;
     // retour de l'indicateur
     return retour;
   };

// programme de vérification par différences finies des dérivées
void Met_Sfe1::Verif_Dcourbure
( const Tableau<Noeud *>& tab_noeud,bool gradV_instantane
		          ,Mat_pleine const & dphiS,int nombre_noeud,Vecteur const & phiS
		          ,Mat_pleine const & tabdphiH,Vecteur const& phiH,const Epai* epas
		          ,Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan)
 { // l'idée est de faire une vérification des dérivées à l'aide d'une méthode de différence finie
   int dim = ParaGlob::Dimension();
   // dans le cas du premier passage on indique qu'il y a vérification
   if (indic_Verif_DcourbureSFE == 0)
     { cout << "\n ****verification du calcul de la variation de courbure ****";
       cout << "\n Met_Sfe1::Verif_Dcourbure(.. \n";
      } 
   indic_Verif_DcourbureSFE++;
  
   // on définit 2 seconde métriques, car il faudra modifier les vecteurs ai et les vecteurs gi, donc
   // il vaut mieux travailler sur une copie
   Met_Sfe1 metrique_bis_1(*this);Met_Sfe1 metrique_bis_2(*this);
   // ici on considère que l'on a même nombre de ddl par noeud = dim
   // on va modifier chaque ddl de chaque noeud systématiquement
   int nbnoeud = tab_noeud.Taille();
   // le deltat pour les différences finis
   double delta = ConstMath::unpeupetit;
   double mini_val = ConstMath::pasmalpetit;
   int numddl = 1; // le compteur de ddl
   int nberr = 10;
   Tableau<int>  erreur(nberr);
   double diff_admi = 0.1;
   bool premier_calcul=true;
   for (int inoeud=1;inoeud<=nbnoeud;inoeud++)
   // on ne fait la vérification que si le noeud externe existe
    if (( inoeud <= 3) || ( (tab_noeud(inoeud-3)->Num_noeud()) != (tab_noeud(inoeud)->Num_noeud())))
    {// on récupère les coordonnées du noeud
     Coordonnee coordtdt = tab_noeud(inoeud)->Coord2();
     for (int ix= 1;ix<=dim;ix++,numddl++)
     { Coordonnee X(dim); X(ix) += delta;
       tab_noeud(inoeud)->Ajout_coord2(X);              
       // appel de la métrique
       metrique_bis_1.CalSfe1_explicit_tdt( tab_noeud,gradV_instantane,dphiS,nombre_noeud,phiS,premier_calcul
		                                  ,tabdphiH,phiH,epas,tabTypeCL,vplan);
       // maintenant on remet les coordonnées du noeud à l'état initial         
       tab_noeud(inoeud)->Change_coord2(coordtdt); 
       // et on décale dégativement
       X(ix) -= 2*delta;       
       tab_noeud(inoeud)->Ajout_coord2(X);              
       // appel de la métrique
       metrique_bis_2.CalSfe1_explicit_tdt( tab_noeud,gradV_instantane,dphiS,nombre_noeud,phiS,premier_calcul
		                                  ,tabdphiH,phiH,epas,tabTypeCL,vplan);
       // maintenant on remet les coordonnées du noeud à l'état initial         
       tab_noeud(inoeud)->Change_coord2(coordtdt); 

       premier_calcul=false;
               
       // calcul des dérivées numériques et  vérification
       int nb_vecteur = 2;
       for (int j=1;j<=nb_vecteur;j++)
        { // variation des vecteurs aiB_tdt
          CoordonneeB daiB = ((*metrique_bis_1.aiB_tdt)(j) -(*metrique_bis_2.aiB_tdt)(j))/(2*delta);
          CoordonneeB toto = ((*metrique_bis_1.aiB_tdt)(j) -(*metrique_bis_2.aiB_tdt)(j));
          
          for (int i=1;i<=dim;i++)
             if (diffpourcent(daiB(i),(*d_aiB_tdt)(numddl)(j)(i),MaX(Dabs(daiB(i)),Dabs((*d_aiB_tdt)(numddl)(j)(i))),diff_admi))
               { if (DabsMiN(daiB(i),(*d_aiB_tdt)(numddl)(j)(i)) <= mini_val)
                   {if ( DabsMaX(daiB(i),(*d_aiB_tdt)(numddl)(j)(i)) > 50.*delta) 
                	    {erreur(1) += 1;
                      cout << "\ndaiB("<<i<<")="<<daiB(i)<<",(*d_aiB_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiB_tdt)(numddl)(j)(i);
                      cout << "\n aiB(xi+delta)="<<(*metrique_bis_1.aiB_tdt)(j)(i) << " aiB(xi-delta)="<<(*metrique_bis_2.aiB_tdt)(j)(i);
                      cout << "\n toto= " << toto << " " << (*aiB_tdt)(j) << " " << (*metrique_bis_1.aiB_tdt)(j);
                	    }
                   }
                else 
                	{erreur(2) += 1;
                     cout << "\ndaiB("<<i<<")="<<daiB(i)<<",(*d_aiB_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiB_tdt)(numddl)(j)(i);
                	};
               };
        }; 
       for (int j=1;j<=nb_vecteur;j++)
        { // variation des vecteurs aiH_tdt
          CoordonneeH d_aiH = ((*metrique_bis_1.aiH_tdt)(j) -(*metrique_bis_2.aiH_tdt)(j))/(2*delta);
          Coordonnee daiH = d_aiH.Coor_const(); // un coordonné à la même place
          // en fait la variation de aiH calculée est uniquement celle dans le plan de la facette, il faut donc
          // retirer la partie hors facette pour la comparaison
          Met_Sfe1::Calcul_N_tdt(); // calcul de la normale
          daiH -= ((d_aiH * N_tdt) * N_tdt).Coor_const(); // on retire la partie portée par la normale
          
          for (int i=1;i<=dim;i++)
             if (diffpourcent(daiH(i),(*d_aiH_tdt)(numddl)(j)(i),MaX(Dabs(daiH(i)),Dabs((*d_aiH_tdt)(numddl)(j)(i))),diff_admi))
               {if (DabsMiN(daiH(i),(*d_aiH_tdt)(numddl)(j)(i)) <= mini_val)
                   {if ( DabsMaX(daiH(i),(*d_aiH_tdt)(numddl)(j)(i)) > 50.*delta) 
                	    {erreur(3) += 1;
                      cout << "\ndaiH("<<i<<")="<<daiH(i)<<",(*d_aiH_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiH_tdt)(numddl)(j)(i);
                      cout << "\n aiH_(xi+delta)="<<(*metrique_bis_1.aiH_tdt)(j)(i) << " aiH_(xi-delta)="<<(*metrique_bis_2.aiH_tdt)(j)(i);
                	    }
                   }
                else 
                	{erreur(4) += 1;
                     cout << "\ndaiH("<<i<<")="<<daiH(i)<<",(*d_aiH_tdt)("<<numddl<<")("<<j<<")("<<i<<")="<< (*d_aiH_tdt)(numddl)(j)(i);
                	};
               };
        }; 

// on ne peut pas faire la vérif des métriques car elles ne sont pas calculées !!
/*
       // variation du tenseur aijBB_tdt
       TenseurBB * aijBB = NevezTenseurBB(*metrique_bis_1.aijBB_tdt);
       *aijBB = (*metrique_bis_1.aijBB_tdt - *metrique_bis_2.aijBB_tdt) / (2*delta);
       for (int i1=1;i1<=nb_vecteur;i1++)
           for (int j1=1;j1<=nb_vecteur;j1++)
             if (diffpourcent((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1),
                DabsMaX((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1)),diff_admi))
                if (DabsMiN((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1)) <= mini_val)
                   {if( DabsMaX((*aijBB)(i1,j1),(*(d_aijBB_tdt(numddl)))(i1,j1)) > 50.*delta) 
                   	{erreur(5) += 1;
                     cout << "\ndaijBB("<<i1<<","<<j1<<")="<<(*aijBB)(i1,j1)<<",(*d_aijBB_tdt)("<<numddl<<")("<<i1<<","<<j1<<")="<< (*(d_aijBB_tdt(numddl)))(i1,j1);
                     cout << "\n(*d_aiB_tdt)("<<numddl<<")(1)="<< (*d_aiB_tdt)(numddl)(1) << " aiB(1)=" << (*aiB_tdt)(1);
                   	}}
                else erreur(6) += 1;
       // variation du tenseur aijHH_tdt
       TenseurHH * aijHH = NevezTenseurHH(*metrique_bis_1.aijHH_tdt);
       *aijHH = (*metrique_bis_1.aijHH_tdt - *metrique_bis_2.aijHH_tdt) / (2*delta);
       for (int i1=1;i1<=nb_vecteur;i1++)
           for (int j1=1;j1<=nb_vecteur;j1++)
             if (diffpourcent((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1),
                DabsMaX((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1)),diff_admi))
                if (DabsMiN((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1)) <= mini_val)
                   {if( DabsMaX((*aijHH)(i1,j1),(*(d_aijHH_tdt(numddl)))(i1,j1)) > 50.*delta) erreur(7) += 1;}
                else erreur(8) += 1;  
*/                
       // variation du tenseur de courbure 
       Tenseur_ns2BB d_curb(curb_tdt);
       // calcul de la variation de la métrique dans le repère absolu              
       d_curb = (metrique_bis_1.curb_tdt - metrique_bis_2.curb_tdt) / (2*delta);
       // on calcul the same pour la variation calculée
//       Tenseur_ns2BB d_curb_glob;
//       dcurb_tdt(numddl).BaseAbsolue(d_curb_glob,*aiH_tdt);
       
       for (int i1=1;i1<=nb_vecteur;i1++)
           for (int j1=1;j1<=nb_vecteur;j1++)
             if (diffpourcent(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1),
                DabsMaX(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1)),diff_admi))
               {if (DabsMiN(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1)) <= mini_val)
                   {if( DabsMaX(d_curb(i1,j1),dcurb_tdt(numddl)(i1,j1)) > 50.*delta) 
                   	 {erreur(9) += 1;
                      cout << "\nd_curb("<<i1<<","<<j1<<")="<<d_curb(i1,j1)<<",dcurb_tdt("<<numddl<<")("<<i1<<","<<j1<<")="<< dcurb_tdt(numddl)(i1,j1);
                      cout << "\n curbe= " << curb_tdt << ",dcurb_tdt("<<numddl<<")("<< dcurb_tdt(numddl) ;
                   	 }
                   }
                else 
             	    {erreur(10) += 1;
                   cout << "\nd_curb("<<i1<<","<<j1<<")="<<d_curb(i1,j1)<<",dcurb_tdt("<<numddl<<")("<<i1<<","<<j1<<")="<< dcurb_tdt(numddl)(i1,j1);
             	    }
              };
          // effacement des tenseurs intermédiaires
 //         delete aijHH; delete  aijBB;        
     };// fin de boucle sur la dimension de coordonnée
    } // fin du if: s'il y a un noeud externe
    else 
    {numddl += 3;}; // on passe les ddl
    
   // fin de boucle sur les noeuds
   
   // message d'erreur si besoin
   bool une_erreur=false;
   for (int l=1;l<=nberr;l++) if (erreur(l) != 0) une_erreur=true;
   if (une_erreur)
     { cout << "\n erreur dans le calcul analytique des derivees de la courbure";
       cout << "\n Met_Sfe1::Verif_Dcourbure(.."
            << " , numero de passage = " << indic_Verif_DcourbureSFE << " erreur = " << erreur << endl;
//       Sortie(1);
      } 
 };