// FICHIER : MetAxisymetrique3D.cp
// CLASSE : MetAxisymetrique3D

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



# include <iostream>
using namespace std;  //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"


#include "MetAxisymetrique3D.h"

// -Constructeur par defaut
MetAxisymetrique3D::MetAxisymetrique3D ():
  Met_abstraite()
     {};
// constructeur permettant de dimensionner uniquement certaines variables
// ici la dimension est 3
// des variables a initialiser
MetAxisymetrique3D::MetAxisymetrique3D (const  DdlElement& tabddl,
	                    const  Tableau<Enum_variable_metrique>& tab,int nomb_noeud) :
	           Met_abstraite(3,3,tabddl,tab, nomb_noeud)
	           {};
// constructeur de copie
MetAxisymetrique3D::MetAxisymetrique3D (const  MetAxisymetrique3D& a) :
    Met_abstraite (a) 
    {};

// destructeur
MetAxisymetrique3D::~MetAxisymetrique3D () {};

 //==================== METHODES PROTEGEES===============================

// calcul de la base  naturel a t0  
void MetAxisymetrique3D::Calcul_giB_0
 (const  Tableau<Noeud *>& tab_noeud,const   Mat_pleine& dphi, int nombre_noeud,const Vecteur& phi)
  { // on effectue le même calcul que dans le cas de met_abstraite mais pour les deux premiers vecteurs
    #ifdef MISE_AU_POINT	 	 
	  if  (giB_0 == NULL)
			{ cout << "\nErreur : la base a t=0 n'est pas dimensionne !\n";
			  cout << "void MetAxisymetrique3D::Calcul_giB_0 \n";
			  Sortie(1);
			};
    #endif
    // tout d'abord on calcul les deux pemiers vecteurs
    for (int i=1;i<= 2;i++)
	    {for (int a=1;a<= dim_base;a++)
		    { double & gib0_i_a = giB_0->CoordoB(i)(a);
		      gib0_i_a = 0.;
		      for (int r=1;r<=nombre_noeud;r++)
			       { gib0_i_a += tab_noeud(r)->Coord0()(a) * dphi(i,r);
          };
      };
	    };
    // puis on définit le dernier vecteur 
    // auparavant il faut le X1 du point, qui est le rayon de rotation
    double rho_0=0;
    for (int r=1;r<=nombre_noeud;r++)
	      rho_0 += tab_noeud(r)->Coord0()(1) * phi(r);
	   // puis on  peut calculer le dernier vecteur
    // cependant il y a un problème lorsque rho_0=0 , dans ce dernier cas on aurait avec la formule
    // normale une norme du deuxième vecteur nulle. On va donc la mettre de manière arbitraire à 1
    // et on considère que l'on ne peut pas avoir extension de ce point c-a-d: que c'est le centre de la matière
    // on ne peut pas avoir un trou qui se crée à ce niveau: on va mettre un test au cas où ce ne serait pas le cas
    if (Abs(rho_0) > ConstMath::pasmalpetit)
	     {giB_0->CoordoB(3)(1)=0.;giB_0->CoordoB(3)(2)=0.;giB_0->CoordoB(3)(3)=rho_0*2.*ConstMath::Pi;
      }
    else // cas d'un point sur l'axe de symétrie
	     {giB_0->CoordoB(3)(1)=0.;giB_0->CoordoB(3)(2)=0.;giB_0->CoordoB(3)(3)=1.;
      };
  };
  
// calcul de la base  naturel a t  -- calcul classique
void MetAxisymetrique3D::Calcul_giB_t
 (const   Tableau<Noeud *>& tab_noeud,const   Mat_pleine& dphi, int nombre_noeud,const Vecteur& phi)
  { // on effectue le même calcul que dans le cas de met_abstraite mais pour les deux premiers vecteurs
    #ifdef MISE_AU_POINT	 	 
	   if  (giB_t == NULL)
     { cout << "\nErreur : la base a t n'est pas dimensionne !\n";
       cout << "void MetAxisymetrique3D::Calcul_giB_t \n";
       Sortie(1);
     };
    #endif
    for (int i=1;i<= 2;i++)
      {for (int a=1;a<= dim_base;a++)	  
        { double & gibt_i_a = giB_t->CoordoB(i)(a);
          gibt_i_a = 0.;
          for (int r=1;r<=nombre_noeud;r++)
           { gibt_i_a += tab_noeud(r)->Coord1()(a) * dphi(i,r);
			        };
		      };
	     };
    // puis on définit le dernier vecteur
    // auparavant il faut le X1 du point, qui est le rayon de rotation
    double rho_t=0;
    for (int r=1;r<=nombre_noeud;r++)
     rho_t += tab_noeud(r)->Coord1()(1) * phi(r);
    #ifdef MISE_AU_POINT	 	 
    double rho_0=0;
    for (int r=1;r<=nombre_noeud;r++)
      rho_0 += tab_noeud(r)->Coord0()(1) * phi(r);
    if ((Abs(rho_0)<=ConstMath::pasmalpetit)&&(Abs(rho_t)>ConstMath::pasmalpetit))
       { cout << "\n *** erreur: dans le calcul de la metrique axisymetrique en un point d'un element ligne "
              << " le rayon initiale du point est nul (="<<rho_0<<") tandis qu'il est non nul a l'instant t (="<<rho_t<<")"
              << " ce cas n'est pas pris en compte dans herezh actuellement !!";
         if (ParaGlob::NiveauImpression()>4)
          cout << "\n MetAxisymetrique3D::Calcul_giB_t(...";
         Sortie(1);
       };
     #endif
    
	   // puis on  peut calculer le dernier vecteur
    // cependant il y a un problème lorsque rho_0=0 , dans ce dernier cas on aurait avec la formule
    // normale une norme du deuxième vecteur nulle. On va donc la mettre de manière arbitraire à 1
    // et on considère que l'on ne peut pas avoir extension de ce point c-a-d: que c'est le centre de la matière
    // on ne peut pas avoir un trou qui se crée à ce niveau: on va mettre un test au cas où ce ne serait pas le cas
    if (Abs(rho_t) > ConstMath::pasmalpetit)
	 {giB_t->CoordoB(3)(1)=0.;giB_t->CoordoB(3)(2)=0.;giB_t->CoordoB(3)(3)=rho_t*2.*ConstMath::Pi;
     }
    else // cas d'un point sur l'axe de symétrie
	 {giB_t->CoordoB(3)(1)=0.; giB_t->CoordoB(3)(2)=0.; giB_t->CoordoB(3)(3)=1.;
     };
  };

// calcul de la base  naturel a tdt  -- calcul classique
void MetAxisymetrique3D::Calcul_giB_tdt
 (const  Tableau<Noeud *>& tab_noeud,const   Mat_pleine& dphi, int nombre_noeud,const Vecteur& phi)
  { // on effectue le même calcul que dans le cas de met_abstraite mais pour les deux premiers vecteurs
    #ifdef MISE_AU_POINT	 	 
    if  (giB_tdt == NULL)
      { cout << "\nErreur : la base a tdt n'est pas dimensionne !\n";
        cout << "void MetAxisymetrique3D::Calcul_giB_tdt \n";
        Sortie(1);
      };
    #endif
    for (int i=1;i<= 2;i++)
      {for (int a=1;a<= dim_base;a++)	  
        { double & gibtdt_i_a = giB_tdt->CoordoB(i)(a);
          gibtdt_i_a = 0.;
          for (int r=1;r<=nombre_noeud;r++)
            { gibtdt_i_a += tab_noeud(r)->Coord2()(a) * dphi(i,r);
            };
        };
      };
   // puis on définit le dernier vecteur
   // auparavant il faut le X1 du point, qui est le rayon de rotation
   double rho_tdt=0;
   for (int r=1;r<=nombre_noeud;r++)
	   rho_tdt += tab_noeud(r)->Coord2()(1) * phi(r);
   #ifdef MISE_AU_POINT	 	 
   double rho_0=0;
   for (int r=1;r<=nombre_noeud;r++)
	    rho_0 += tab_noeud(r)->Coord0()(1) * phi(r);
   if ((Abs(rho_0)<=ConstMath::pasmalpetit)&&(Abs(rho_tdt)>ConstMath::pasmalpetit))
     { cout << "\n *** erreur: dans le calcul de la metrique axisymetrique en un point d'un element ligne "
            << " le rayon initiale du point est nul (="<<rho_0<<") tandis qu'il est non nul a l'instant t (="<<rho_tdt<<")"
            << " ce cas n'est pas pris en compte dans herezh actuellement !!";
       if (ParaGlob::NiveauImpression()>4)
        cout << "\n MetAxisymetrique3D::Calcul_giB_tdt(...";
       Sortie(1);
     };
   #endif
    
   // maintenant on peut calculer le dernier vecteur
   giB_tdt->CoordoB(3)(1)=0.; giB_tdt->CoordoB(3)(2)=0.;giB_tdt->CoordoB(3)(3)=rho_tdt*2.*ConstMath::Pi;
	  // puis on  peut calculer le dernier vecteur
   // cependant il y a un problème lorsque rho_0=0 , dans ce dernier cas on aurait avec la formule
   // normale une norme du deuxième vecteur nulle. On va donc la mettre de manière arbitraire à 1
   // et on considère que l'on ne peut pas avoir extension de ce point c-a-d: que c'est le centre de la matière
   // on ne peut pas avoir un trou qui se crée à ce niveau: on va mettre un test au cas où ce ne serait pas le cas
   if (Abs(rho_tdt) > ConstMath::pasmalpetit)
	    {giB_tdt->CoordoB(3)(1)=0.;giB_tdt->CoordoB(3)(2)=0.;giB_tdt->CoordoB(3)(3)=rho_tdt*2.*ConstMath::Pi;
     }
   else // cas d'un point sur l'axe de symétrie
	    {giB_tdt->CoordoB(3)(1)=0.; giB_tdt->CoordoB(3)(2)=0.; giB_tdt->CoordoB(3)(3)=1.;
     };
  };

//==  calcul de la variation des bases
 void MetAxisymetrique3D::D_giB_t( const  Mat_pleine& dphi, int nbnoeu,const  Vecteur & phi)
  { int indice; 
    // derivees des  gBi par rapport a XHbr, b=1 et 2 c'est-à-dire 
    int nb_ddl_en_var =  dim_base-1;  // on ne prend en compte que les variations suivant x1 et x2
    for (int b = 1; b<=2; b++) // uniquement pour les 2 premières coordonnées
      for (int r = 1; r <= nbnoeu; r++)
         { indice = (r-1)*nb_ddl_en_var + b; // sont  rangé 2 par 2
           BaseB & d_giB_t_indice=(*d_giB_t)(indice);
           // variation des deux premiers vecteurs
           for (int i =1;i<=2;i++) 
              { d_giB_t_indice.CoordoB(i).Zero();
                d_giB_t_indice.CoordoB(i)(b)=dphi(i,r);
               }
           // le troisième vecteur à une dérivée non nulle suivant XH1r
           d_giB_t_indice.CoordoB(3).Zero(); // init
           // comme g_3=X^{1r} * phi_r * I_3
           if (b==1)
             d_giB_t_indice.CoordoB(3)(3) = phi(r)*2.*ConstMath::Pi;
          } 
    // la dérivée de tous les vecteurs par rapport à XH3r est nulle                 
//    for (int r = 1; r <= nbnoeu; r++)
//         { indice = (r-1)*nb_ddl_en_var + 3; // sont quand même rangé 3 par 3
//           BaseB & d_giB_t_indice=(*d_giB_t)(indice);
//           d_giB_t_indice(1).Zero();
//           d_giB_t_indice(2).Zero();
//           d_giB_t_indice(3).Zero();
//          } 
   };
   
 void MetAxisymetrique3D::D_giB_tdt(  const Mat_pleine& dphi, int nbnoeu,const  Vecteur & phi)
  { int indice; 
    // derivees des  gBi par rapport a XHbr, b=1 et 2 c'est-à-dire    
    int nb_ddl_en_var =  dim_base-1;  // on ne prend en compte que les variations suivant x1 et x2
    for (int b = 1; b<=2; b++) // uniquement pour les 2 premières coordonnées
      for (int r = 1; r <= nbnoeu; r++)
         { indice = (r-1)*nb_ddl_en_var + b; // sont  rangé 2 par 2
           BaseB & d_giB_tdt_indice=(*d_giB_tdt)(indice);
           // variation des deux premiers vecteurs
           for (int i =1;i<=2;i++) 
              { d_giB_tdt_indice.CoordoB(i).Zero();
                d_giB_tdt_indice.CoordoB(i)(b)=dphi(i,r);
               }
           // le troisième vecteur à une dérivée non nulle suivant XH1r
           d_giB_tdt_indice.CoordoB(3).Zero(); // init
           // comme g_3=X^{1r} * phi_r * I_3
           if (b==1)
             d_giB_tdt_indice.CoordoB(3)(3) = phi(r)*2.*ConstMath::Pi;
          } 
    // la dérivée de tous les vecteurs par rapport à XH3r est nulle                 
//    for (int r = 1; r <= nbnoeu; r++)
//         { indice = (r-1)*nb_ddl_en_var + 3; // sont quand même rangé 3 par 3
//           BaseB & d_giB_tdt_indice=(*d_giB_tdt)(indice);
//           d_giB_tdt_indice(1).Zero();
//           d_giB_tdt_indice(2).Zero();
//           d_giB_tdt_indice(3).Zero();
//          } 
   };


  // calcul du gradient de vitesse à t
void MetAxisymetrique3D::Calcul_gradVBB_t
		   (const  Tableau<Noeud *>& tab_noeud, const  Mat_pleine& dphi,int nombre_noeud)
 {  
  #ifdef MISE_AU_POINT	 	 
	if  ((giB_t == NULL) || (gradVBB_t == NULL))
			{ cout << "\nErreur : la base et/ou le gradient de vitesse a t n'est pas defini !\n";
			  cout << "\nMetAxisymetrique3D::Calcul_gradVBB_t \n";
			  Sortie(1);
			};
  #endif 
  // Vi|j = V^{ar} \phi_{r,j} \vec I_a . \vec g_i = V,j . g_i 
  // Pour les  éléments axisymétriques, il n'y a pas de vitesse suivant z qui est la direction
  // tangente à l'axe de rotation qui est y
  // de plus la vitesse est normale à g_3
  for (int i=1;i<= nbvec_base;i++)
   {CoordonneeB & giBt = giB_t->CoordoB(i);
    // pour Vi|j, avec j variant de 1 à 2
    for (int j=1;j<= 2;j++) // j ici varie de 1 à 2 qui est la dimension d'interpolation
     {double& gradVij=(*gradVBB_t).Coor(i,j);
      gradVij = 0.;
	  for (int r=1;r<=nombre_noeud;r++)
       { // on ne tiens pas compte de la valeur V3 qui est nulle 
	     gradVij += tab_noeud(r)->Valeur_t(V2) * dphi(j,r) * giBt(2);
		 gradVij += tab_noeud(r)->Valeur_t(V1) * dphi(j,r) * giBt(1);      
       } // -- fin du for sur r 
     } // -- fin du for sur j
    // pour Vi|j, avec j=3   -> 0 car V,3 = 0
    (*gradVBB_t).Coor(i,3) = 0.;
    }// -- fin du for sur i   
  };  
		                         
  // calcul gradient de vitesse à t+dt
void MetAxisymetrique3D::Calcul_gradVBB_tdt
		   (const Tableau<Noeud *>& tab_noeud, const  Mat_pleine& dphi,int nombre_noeud)
 {  
  #ifdef MISE_AU_POINT	 	 
	if  ((giB_tdt == NULL) || (gradVBB_tdt == NULL))
			{ cout << "\nErreur : la base et/ou le gradient de vitesse a t+dt n'est pas defini !\n";
			  cout << "\nMetAxisymetrique3D::Calcul_gradVBB_tdt \n";
			  Sortie(1);
			};
  #endif
  // Vi|j = V^{ar} \phi_{r,j} \vec I_a . \vec g_i = V,j . g_i
  // Pour les  éléments axisymétriques, il n'y a pas de vitesse suivant z qui est la direction
  // tangente à l'axe de rotation qui est y
  // de plus la vitesse est normale à g_3
  for (int i=1;i<= nbvec_base;i++)
   {CoordonneeB & giBtdt = giB_tdt->CoordoB(i);
    // pour Vi|j, avec j variant de 1 à 2
    for (int j=1;j<= 2;j++) // j ici varie de 1 à 2 qui est la dimension d'interpolation
     {double& gradVij=(*gradVBB_tdt).Coor(i,j);
      gradVij = 0.;
	  for (int r=1;r<=nombre_noeud;r++)
       { // on ne tiens pas compte de la valeur V3 qui est nulle 
	     gradVij += tab_noeud(r)->Valeur_tdt(V2) * dphi(j,r) * giBtdt(2);
		 gradVij += tab_noeud(r)->Valeur_tdt(V1) * dphi(j,r) * giBtdt(1);      
       } // -- fin du for sur r 
     } // -- fin du for sur j
    // pour Vi|j, avec j=3   -> 0 car V,3 = 0
    (*gradVBB_tdt).Coor(i,3) = 0.;
    }// -- fin du for sur i   
  };  

// calcul du gradient de vitesse moyen en utilisant delta x^ar/delta t
// dans le cas où les ddl à tdt n'existent pas -> utilisation de la  vitesse sécante entre 0 et t !!
void MetAxisymetrique3D::Calcul_gradVBB_moyen_t
		   (const  Tableau<Noeud *>& tab_noeud, const  Mat_pleine& dphi,int nombre_noeud) 
 {  
  #ifdef MISE_AU_POINT
	if  ((giB_t == NULL) || (gradVmoyBB_t == NULL))
			{ cout << "\nErreur : la base et/ou le gradient de vitesse a t+dt n'est pas defini !\n";
			  cout << "\nMetAxisymetrique3D::Calcul_gradVBB_moyen_t \n";
			  Sortie(1);
			};
  #endif
  // on pose V^{ar} = delta X^{ar} / delta t
  // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
  const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
  // dans le cas où l'incrément de temps est nul la vitesse est nulle
  double deltat=vartemps.IncreTempsCourant();double unsurdeltat;
  double temps =vartemps.TempsCourant(); 
  // on regarde le premier noeud pour savoir s'il y a des coordonnées à t+dt
  if (tab_noeud(1)->ExisteCoord2())
   {// dans ce cas on utilise l'incrément de temps
    if (deltat >= ConstMath::trespetit)
        {unsurdeltat = 1./deltat;}
    else
        {unsurdeltat = 0.;} // o car on veut une vitesse nulle 
    }
  else // sinon on utilise la variation totale de 0 à t, le delta t vaudra le temps actuel      
   { if (temps >= ConstMath::trespetit)
       {unsurdeltat = 1./temps;}
     else
       {unsurdeltat = 0.;} // o car on veut une vitesse nulle
     }   
      // Vi|j = V^{ar} \phi_{r,j} \vec I_a . \vec g_i 
  gradVmoyBB_t->Inita(0.);
  for (int iv=1;iv<=nbvec_base;iv++)
    dmatV_moy_t->CoordoB(iv).Zero();
  for (int r=1;r<=nombre_noeud;r++)
    { // on construit d'abord le vecteur vitesse
      Coordonnee & V = V_moy_t(r);
      if (tab_noeud(r)->ExisteCoord2())
         V = (tab_noeud(r)->Coord2() - tab_noeud(r)->Coord1())*unsurdeltat; 
      else   // sinon utilisation des coordonnées sécantes !!!  
         V = (tab_noeud(r)->Coord1() - tab_noeud(r)->Coord0())*unsurdeltat;  
      // Pour les  éléments axisymétriques, il n'y a pas de vitesse suivant z qui est la direction
      // tangente à l'axe de rotation qui est y
      for (int i=1;i<= nbvec_base;i++)
        {CoordonneeB & giBt = giB_t->CoordoB(i);
         // pour Vi|j, avec j variant de 1 à 2
         for (int j=1;j<= 2;j++) // j ici varie de 1 à 2 qui est la dimension d'interpolation
          {double& gradVij=(*gradVmoyBB_t).Coor(i,j);
           // on ne tiens pas compte de V(3) qui doit être nul
		   gradVij += V(2) * dphi(j,r) * giBt(2);
		   gradVij += V(1) * dphi(j,r) * giBt(1);      
           } // -- fin du for sur j
         // calcul des vecteurs V,i moyens
         CoordonneeB & dmatV_moy_t_i= dmatV_moy_t->CoordoB(i);
         dmatV_moy_t_i(2) += V(2) * dphi(i,r);
         dmatV_moy_t_i(1) += V(1) * dphi(i,r);
         // pour Vi|j, avec j=3   -> 0 car V,3 = 0
         // ok car déjà initialisé, idem pour dmatV_moy_t
         }// -- fin du for sur i   
      } // -- fin du for sur r 
  };  

// calcul du gradient de vitesse moyen en utilisant delta x^ar/delta t
void MetAxisymetrique3D::Calcul_gradVBB_moyen_tdt
		   (const  Tableau<Noeud *>& tab_noeud, const  Mat_pleine& dphi,int nombre_noeud) 
 {  
  #ifdef MISE_AU_POINT
	if  ((giB_tdt == NULL) || (gradVmoyBB_tdt == NULL))
			{ cout << "\nErreur : la base et/ou le gradient de vitesse a t+dt n'est pas defini !\n";
			  cout << "\nMetAxisymetrique3D::Calcul_gradVBB_moyen_tdt \n";
			  Sortie(1);
			};
  #endif
  // on pose V^{ar} = delta X^{ar} / delta t
  // on choisit une vitesse de déformation pour l'instant identique à deltaeps/t
  const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
  // dans le cas où l'incrément de temps est nul la vitesse est nulle
  double deltat=vartemps.IncreTempsCourant();
  double unsurdeltat;
  if (deltat >= ConstMath::trespetit)
     {unsurdeltat = 1./deltat;}
  else
     {unsurdeltat = 0.;} // o car on veut une vitesse nulle 
      // Vi|j = V^{ar} \phi_{r,j} \vec I_a . \vec g_i 
  gradVmoyBB_tdt->Inita(0.);
  for (int iv=1;iv<=nbvec_base;iv++)
    dmatV_moy_tdt->CoordoB(iv).Zero();
  for (int r=1;r<=nombre_noeud;r++)
    { // on construit d'abord le vecteur vitesse
      Coordonnee & V = V_moy_t(r);
      V = (tab_noeud(r)->Coord2() - tab_noeud(r)->Coord1())*unsurdeltat;  
      // Pour les  éléments axisymétriques, il n'y a pas de vitesse suivant z qui est la direction
      // tangente à l'axe de rotation qui est y
      for (int i=1;i<= nbvec_base;i++)
        {CoordonneeB & giBtdt = giB_tdt->CoordoB(i);
         // pour Vi|j, avec j variant de 1 à 2
         for (int j=1;j<= 2;j++) // j ici varie de 1 à 2 qui est la dimension d'interpolation
          {double& gradVij=(*gradVmoyBB_tdt).Coor(i,j);
           // on ne tiens pas compte de V(3) qui doit être nul
           gradVij += V(2) * dphi(j,r) * giBtdt(2);
           gradVij += V(1) * dphi(j,r) * giBtdt(1);
          }; // -- fin du for sur j
         // calcul des vecteurs V,i moyens
         CoordonneeB & dmatV_moy_tdt_i= dmatV_moy_tdt->CoordoB(i);
         dmatV_moy_tdt_i(2) += V(2) * dphi(i,r);
         dmatV_moy_tdt_i(1) += V(1) * dphi(i,r);
         // pour Vi|j, avec j=3   -> 0 car V,3 = 0
         // ok car déjà initialisé, idem pour dmatV_moy_tdt
         }// -- fin du for sur i   
      } // -- fin du for sur r 
  };  
  

   //== calcul de la variation du gradient instantannée	
   // par rapport aux composantes  V^ar (et non les X^ar )	
void MetAxisymetrique3D::DgradVBB_t(const  Mat_pleine& dphi)
  { int indice; 
    // derivees par rapport a VHbr, b=1 et 2 c'est-à-dire     
    // uniquement pour les 2 premières coordonnées
    int nb_ddl_en_var =  dim_base-1;  // on ne prend en compte que les variations suivant x1 et x2
    // Pour les  éléments axisymétriques, il n'y a pas de vitesse suivant z qui est la direction
    // tangente à l'axe de rotation qui est y
    for (int b = 1; b<= 2; b++) 
      for (int r = 1; r <= nomb_noeud; r++)
         { indice = (r-1)*nb_ddl_en_var + b; // sont rangé 2 par 2
           TenseurBB & d_gradVBB_t_indice = (*(d_gradVBB_t(indice)));
           // les dérivées non-nulles
           for (int j =1;j<=2;j++)
             for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_point * gi
               d_gradVBB_t_indice.Coor(i,j) = dphi(j,r) * (*giB_t)(i)(b);
           // les dérivées nulles
           for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_point * gi
               d_gradVBB_t_indice.Coor(i,3) = 0.;
          };
    // les dérivées par rapport aux V3 sont nulles
//    for (int r = 1; r <= nomb_noeud; r++)
//      { indice = (r-1)*dim_base + 3; 
//        (*(d_gradVBB_t(indice))).Inita(0.);
//      };
   }; 

   //== calcul de la variation du gradient	
   // par rapport aux composantes  V^ar (et non les X^ar )	
void MetAxisymetrique3D::DgradVBB_tdt(const  Mat_pleine& dphi)
  { int indice; 
    // derivees par rapport a VHbr, b=1 et 2 c'est-à-dire     
    // uniquement pour les 2 premières coordonnées
    int nb_ddl_en_var =  dim_base-1;  // on ne prend en compte que les variations suivant x1 et x2
    // Pour les  éléments axisymétriques, il n'y a pas de vitesse suivant z qui est la direction
    // tangente à l'axe de rotation qui est y
    for (int b = 1; b<= 2; b++) 
      for (int r = 1; r <= nomb_noeud; r++)
         { indice = (r-1)*nb_ddl_en_var + b; // sont rangé 2 par 2
           TenseurBB & d_gradVBB_tdt_indice = (*(d_gradVBB_tdt(indice)));
           // les dérivées non-nulles
           for (int j =1;j<=2;j++)
             for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_pointdt * gi
               d_gradVBB_tdt_indice.Coor(i,j) = dphi(j,r) * (*giB_tdt)(i)(b);
           // les dérivées nulles
           for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_pointdt * gi
               d_gradVBB_tdt_indice.Coor(i,3) = 0.;
          };
    // les dérivées par rapport aux V3 sont nulles
//    for (int r = 1; r <= nomb_noeud; r++)
//      { indice = (r-1)*dim_base + 3; 
//        (*(d_gradVBB_tdt(indice))).Inita(0.);
//      };
   }; 
    
// calcul de la variation du gradient moyen à t
// par rapport aux composantes  X^ar (et non les V^ar )	
void MetAxisymetrique3D::DgradVmoyBB_t(const  Mat_pleine& dphi,const  Tableau<Noeud *>& tab_noeud)
  { // récup de delta t
    const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
    // dans le cas où l'incrément de temps est nul la vitesse est nulle
    double deltat=vartemps.IncreTempsCourant();double unsurdeltat;
    double temps =vartemps.TempsCourant(); 
    // on regarde le premier noeud pour savoir s'il y a des coordonnées à t+dt
    if (tab_noeud(1)->ExisteCoord2())
     {// dans ce cas on utilise l'incrément de temps
      if (deltat >= ConstMath::trespetit)
          {unsurdeltat = 1./deltat;}
      else
          {unsurdeltat = 0.;} // o car on veut une vitesse nulle 
       }
    else // sinon on utilise la variation totale de 0 à t, le delta t vaudra le temps actuel      
     { if (temps >= ConstMath::trespetit)
         {unsurdeltat = 1./temps;}
       else
         {unsurdeltat = 0.;} // o car on veut une vitesse nulle
       }       
    int indice; 
    // derivees par rapport a XHbr, b=1 à 2      
    int nb_ddl_en_var =  dim_base-1;  // on ne prend en compte que les variations suivant x1 et x2
    for (int b = 1; b<= 2; b++)
      for (int r = 1; r <= nomb_noeud; r++)
         { indice = (r-1)*nb_ddl_en_var + b;
           TenseurBB & d_gradVmoyBB_t_indice = (*(d_gradVmoyBB_t(indice)));
           BaseB & d_giB_t_indice=(*d_giB_t)(indice);
           // tout d'abord les variation impliquant V,alpha , d'ou j=1 à 2
           for (int j =1;j<=2;j++)
             for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_point * gi
               d_gradVmoyBB_t_indice.Coor(i,j) = dphi(j,r) *( unsurdeltat * (*giB_t)(i)(b)
                                            + (*dmatV_moy_t)(j).ScalBB(d_giB_t_indice(i)) );
           // puis  les variation impliquant V,3 ce qui donne 0  
           for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_point * gi
               d_gradVmoyBB_t_indice.Coor(i,3) = 0.;
          };              
    // derivees par rapport a XHbr, b=3, -> 0 car les composantes z sont fixes      
//    for (int r = 1; r <= nomb_noeud; r++)
//         { indice = (r-1)*dim_base + 3;
//           (*(d_gradVmoyBB_t(indice))).Inita(0.) ;
//          };              
   }; 
   
// calcul de la variation du gradient moyen à t dt
// par rapport aux composantes  X^ar (et non les V^ar )	
void MetAxisymetrique3D::DgradVmoyBB_tdt(const  Mat_pleine& dphi)
  { // récup de delta t
    const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
    // dans le cas où l'incrément de temps est nul la vitesse est nulle
    double deltat=vartemps.IncreTempsCourant();
    double unsurdeltat;
    if (deltat >= ConstMath::trespetit)
       {unsurdeltat = 1./deltat;}
    else
       {unsurdeltat = 0.;} // o car on veut une vitesse nulle 
    
    int indice; 
    // derivees par rapport a XHbr, b=1 à 2      
    int nb_ddl_en_var =  dim_base-1;  // on ne prend en compte que les variations suivant x1 et x2
    for (int b = 1; b<= 2; b++)
      for (int r = 1; r <= nomb_noeud; r++)
         { indice = (r-1)*nb_ddl_en_var + b;
           TenseurBB & d_gradVmoyBB_tdt_indice = (*(d_gradVmoyBB_tdt(indice)));
           BaseB & d_giB_tdt_indice=(*d_giB_tdt)(indice);
           // tout d'abord les variation impliquant V,alpha , d'ou j=1 à 2
           for (int j =1;j<=2;j++)
             for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_point * gi
               d_gradVmoyBB_tdt_indice.Coor(i,j) = dphi(j,r) *( unsurdeltat * (*giB_tdt)(i)(b)
                                            + (*dmatV_moy_tdt)(j).ScalBB(d_giB_tdt_indice(i)) );
           // puis  les variation impliquant V,3 ce qui donne 0  
           for (int i =1;i<=nbvec_base;i++)
               // grad(i,j) = Vi|j = gj_point * gi
               d_gradVmoyBB_tdt_indice.Coor(i,3) = 0.;
          };              
    // derivees par rapport a XHbr, b=3, -> 0 car les composantes z sont fixes      
//    for (int r = 1; r <= nomb_noeud; r++)
//         { indice = (r-1)*dim_base + 3;
//           (*(d_gradVmoyBB_tdt(indice))).Inita(0.) ;
//          };              
   };