877 lines
37 KiB
C++
877 lines
37 KiB
C++
// FICHIER : MetAxisymetrique2D.cc
|
|
// CLASSE : MetAxisymetrique2D
|
|
|
|
// 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 "MetAxisymetrique2D.h"
|
|
|
|
// -Constructeur par defaut
|
|
MetAxisymetrique2D::MetAxisymetrique2D ():
|
|
Met_abstraite()
|
|
{};
|
|
// constructeur permettant de dimensionner uniquement certaines variables
|
|
// ici la dimension est 2 ou 3, le nombre de vecteur par contre est fixe: 2
|
|
// des variables a initialiser
|
|
MetAxisymetrique2D::MetAxisymetrique2D (int dim_base,const DdlElement& tabddl,
|
|
const Tableau<Enum_variable_metrique>& tab,int nomb_noeud) :
|
|
Met_abstraite(dim_base,2,tabddl,tab, nomb_noeud)
|
|
{};
|
|
|
|
// constructeur de copie
|
|
MetAxisymetrique2D::MetAxisymetrique2D (const MetAxisymetrique2D& a) :
|
|
Met_abstraite (a)
|
|
{};
|
|
|
|
// destructeur
|
|
MetAxisymetrique2D::~MetAxisymetrique2D () {};
|
|
|
|
//==================== METHODES PROTEGEES===============================
|
|
|
|
// calcul de la base naturel a t0
|
|
void MetAxisymetrique2D::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 MetAxisymetrique2D::Calcul_giB_0 \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// on défini le X1 du point, qui est le rayon de rotation
|
|
// autour de y ou z
|
|
rho_0=0;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
rho_0 += tab_noeud(r)->Coord0()(1) * phi(r);
|
|
// calcul du premier vecteur
|
|
for (int a=1;a<= dim_base;a++) // le premier vecteur
|
|
{ double & gib0_1_a = giB_0->CoordoB(1)(a);
|
|
gib0_1_a = 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
gib0_1_a += tab_noeud(r)->Coord0()(a) * dphi(1,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)
|
|
{for (int k=1;k<=dim_base;k++)
|
|
{if (k==dim_base)
|
|
giB_0->CoordoB(2)(k)=rho_0*2.*ConstMath::Pi;
|
|
else
|
|
giB_0->CoordoB(2)(k)=0.;
|
|
};
|
|
}
|
|
else // cas d'un point sur l'axe de symétrie
|
|
{for (int k=1;k<=dim_base;k++)
|
|
{if (k==dim_base)
|
|
giB_0->CoordoB(2)(k)=1.;
|
|
else
|
|
giB_0->CoordoB(2)(k)=0.;
|
|
};
|
|
};
|
|
////debug
|
|
//tab_noeud(1)->Coord0().Affiche();
|
|
//tab_noeud(2)->Coord0().Affiche();
|
|
//(*(giB_0))(1).Affiche();
|
|
//(*(giB_0))(2).Affiche();
|
|
//cout << endl;
|
|
////fin debug
|
|
|
|
|
|
};
|
|
|
|
// calcul de la base naturel a t -- calcul classique
|
|
void MetAxisymetrique2D::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 MetAxisymetrique2D::Calcul_giB_t \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// on défini le X1 du point, qui est le rayon de rotation
|
|
// autour de y ou z
|
|
rho_t=0;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
rho_t += tab_noeud(r)->Coord1()(1) * phi(r);
|
|
#ifdef MISE_AU_POINT
|
|
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 MetAxisymetrique2D::Calcul_giB_t(...";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
// calcul du premier vecteur
|
|
for (int a=1;a<= dim_base;a++) // le premier vecteur
|
|
{ double & gibt_1_a = giB_t->CoordoB(1)(a);
|
|
gibt_1_a = 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
gibt_1_a += tab_noeud(r)->Coord1()(a) * dphi(1,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_t) > ConstMath::pasmalpetit)
|
|
{for (int k=1;k<=dim_base;k++)
|
|
{if (k==dim_base)
|
|
giB_t->CoordoB(2)(k)=rho_t*2.*ConstMath::Pi;
|
|
else
|
|
giB_t->CoordoB(2)(k)=0.;
|
|
};
|
|
}
|
|
else // cas d'un point sur l'axe de symétrie
|
|
{for (int k=1;k<=dim_base;k++)
|
|
{if (k==dim_base)
|
|
giB_t->CoordoB(2)(k)=1.;
|
|
else
|
|
giB_t->CoordoB(2)(k)=0.;
|
|
};
|
|
};
|
|
};
|
|
|
|
// calcul de la base naturel a tdt -- calcul classique
|
|
void MetAxisymetrique2D::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 MetAxisymetrique2D::Calcul_giB_tdt \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// on défini le X1 du point, qui est le rayon de rotation
|
|
// autour de y ou z
|
|
rho_tdt=0;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
rho_tdt += tab_noeud(r)->Coord2()(1) * phi(r);
|
|
#ifdef MISE_AU_POINT
|
|
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))
|
|
// if ((Abs(rho_0)== 0.)&&(Abs(rho_tdt)>ConstMath::petit))
|
|
{ cout << "\n *** attention: 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 est pris en compte de maniere arbitraire dans herezh actuellement !!";
|
|
if (ParaGlob::NiveauImpression()>4)
|
|
cout << "\n MetAxisymetrique2D::Calcul_giB_tdt(...";
|
|
// Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
// calcul du premier vecteur
|
|
for (int a=1;a<= dim_base;a++) // le premier vecteur
|
|
{ double & gibtdt_1_a = giB_tdt->CoordoB(1)(a);
|
|
gibtdt_1_a = 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
gibtdt_1_a += tab_noeud(r)->Coord2()(a) * dphi(1,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_tdt) > ConstMath::pasmalpetit)
|
|
{for (int k=1;k<=dim_base;k++)
|
|
{if (k==dim_base)
|
|
giB_tdt->CoordoB(2)(k)=rho_tdt*2.*ConstMath::Pi;
|
|
else
|
|
giB_tdt->CoordoB(2)(k)=0.;
|
|
};
|
|
}
|
|
else // cas d'un point sur l'axe de symétrie
|
|
{for (int k=1;k<=dim_base;k++)
|
|
{if (k==dim_base)
|
|
giB_tdt->CoordoB(2)(k)=1.;
|
|
else
|
|
giB_tdt->CoordoB(2)(k)=0.;
|
|
};
|
|
};
|
|
};
|
|
|
|
// calcul des variations du point a tdt
|
|
void MetAxisymetrique2D::Calcul_d_Mtdt
|
|
(const Tableau<Noeud *>& , const Vecteur& phi, int nbnoeu)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Mtdt == NULL)
|
|
{ cout << "\nErreur : le de variation du point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_Mtdt \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
int indice;
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b < dim_base; b++) //on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int r = 1; r <= nbnoeu; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
Coordonnee & d_Mtdt_indice=(*d_Mtdt)(indice);
|
|
d_Mtdt_indice.Zero();
|
|
d_Mtdt_indice(b)=phi(r);
|
|
}
|
|
};
|
|
// calcul des variations de la vitesse du point a t en fonction des ddl existants de vitesse
|
|
void MetAxisymetrique2D::Calcul_d_Vt
|
|
(const Tableau<Noeud *>& , const Vecteur& phi, int nombre_noeud)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_Vt \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
int indice;
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b < dim_base; b++) //on ne prend en compte que les variations jusqu'a dim-1
|
|
{for (int r = 1; r <= nombre_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
Coordonnee & d_Vt_indice=(*d_Vt)(indice);
|
|
d_Vt_indice.Zero();
|
|
d_Vt_indice(b)=phi(r);
|
|
};
|
|
};
|
|
};
|
|
// idem mais au noeud passé en paramètre
|
|
void MetAxisymetrique2D::Calcul_d_Vt (const Noeud* )
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_Vt(const Noeud* \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int b = 1; b < dim_base; b++) //on ne prend en compte que les variations jusqu'a dim-1
|
|
{ Coordonnee & d_Vt_indice=(*d_Vt)(b);
|
|
d_Vt_indice.Zero();
|
|
d_Vt_indice(b)=1;
|
|
};
|
|
};
|
|
|
|
// calcul des variations de la vitesse du point a tdt en fonction des ddl existants de vitesse
|
|
void MetAxisymetrique2D::Calcul_d_Vtdt
|
|
(const Tableau<Noeud *>& , const Vecteur& phi, int nombre_noeud)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vtdt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_Vtdt \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
int indice;
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b < dim_base; b++) //on ne prend en compte que les variations jusqu'a dim-1
|
|
{for (int r = 1; r <= nombre_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
Coordonnee & d_Vtdt_indice=(*d_Vtdt)(indice);
|
|
d_Vtdt_indice.Zero();
|
|
d_Vtdt_indice(b)=phi(r);
|
|
};
|
|
};
|
|
};
|
|
// idem mais au noeud passé en paramètre
|
|
void MetAxisymetrique2D::Calcul_d_Vtdt (const Noeud* )
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vtdt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse au point a tdt n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_Vtdt(const Noeud* \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int b = 1; b < dim_base; b++) //on ne prend en compte que les variations jusqu'a dim-1
|
|
{ Coordonnee & d_Vtdt_indice=(*d_Vtdt)(b);
|
|
d_Vtdt_indice.Zero();
|
|
d_Vtdt_indice(b)=1;
|
|
};
|
|
};
|
|
|
|
// calcul des variations de la vitesse moyenne du point a t en fonction des ddl de position
|
|
void MetAxisymetrique2D::Calcul_d_V_moyt
|
|
(const Tableau<Noeud *>& , const Vecteur& phi, int nombre_noeud)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse moyenne au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_V_moyt \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// la vitesse est calculée selon (M(t)-M(0)) / t, ici on se réfère par rapport au ddl de M(t)
|
|
// dans le cas où le temps est nul on pose que la vitesse est nulle et non infinie !
|
|
double temps = ParaGlob::Variables_de_temps().TempsCourant(); double unsurt;
|
|
if (temps >= ConstMath::trespetit) {unsurt = 1./temps;}
|
|
else {unsurt = 0.;} // o car on veut une vitesse nulle
|
|
int indice;
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b< dim_base; b++)
|
|
{for (int r = 1; r <= nombre_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
Coordonnee & d_Vt_indice=(*d_Vt)(indice);
|
|
d_Vt_indice.Zero();
|
|
d_Vt_indice(b)=unsurt * phi(r);
|
|
};
|
|
};
|
|
};
|
|
// idem mais au noeud passé en paramètre
|
|
void MetAxisymetrique2D::Calcul_d_V_moyt(const Noeud* )
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse moyenne au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_V_moyt(const Noeud* \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// la vitesse est calculée selon (M(t)-M(0)) / t, ici on se réfère par rapport au ddl de M(t)
|
|
// dans le cas où le temps est nul on pose que la vitesse est nulle et non infinie !
|
|
double temps = ParaGlob::Variables_de_temps().TempsCourant(); double unsurt;
|
|
if (temps >= ConstMath::trespetit) {unsurt = 1./temps;}
|
|
else {unsurt = 0.;} // o car on veut une vitesse nulle
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b< dim_base; b++)
|
|
{Coordonnee & d_Vt_indice=(*d_Vt)(b);
|
|
d_Vt_indice.Zero();
|
|
d_Vt_indice(b)=unsurt;
|
|
};
|
|
};
|
|
|
|
// calcul des variations de la vitesse moyenne du point a tdt en fonction des ddl de position
|
|
void MetAxisymetrique2D::Calcul_d_V_moytdt
|
|
(const Tableau<Noeud *>& , const Vecteur& phi, int nombre_noeud)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vtdt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse moyenne au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_V_moytdt \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// la vitesse est calculée selon (M(tdt)-M(t)) / deltat, les ddl sont ceux de M(tdt)
|
|
// dans le cas où l'incrément de temps est nul on pose que la vitesse est nulle et non infinie !
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();double unsurdeltat;
|
|
if (deltat >= ConstMath::trespetit) {unsurdeltat = 1./deltat;}
|
|
else {unsurdeltat = 0.;} // o car on veut une vitesse nulle
|
|
int indice;
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b< dim_base; b++)
|
|
{for (int r = 1; r <= nombre_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
Coordonnee & d_Vtdt_indice=(*d_Vtdt)(indice);
|
|
d_Vtdt_indice.Zero();
|
|
d_Vtdt_indice(b)=phi(r) * unsurdeltat;
|
|
};
|
|
};
|
|
};
|
|
// idem mais au noeud passé en paramètre
|
|
void MetAxisymetrique2D::Calcul_d_V_moytdt(const Noeud* )
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (d_Vtdt == NULL)
|
|
{ cout << "\nErreur : la variation de vitesse moyenne au point a t n'est pas dimensionne !\n";
|
|
cout << "void MetAxisymetrique2D::Calcul_d_V_moytdt(const Noeud* \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// la vitesse est calculée selon (M(tdt)-M(t)) / deltat, les ddl sont ceux de M(tdt)
|
|
// dans le cas où l'incrément de temps est nul on pose que la vitesse est nulle et non infinie !
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();double unsurdeltat;
|
|
if (deltat >= ConstMath::trespetit) {unsurdeltat = 1./deltat;}
|
|
else {unsurdeltat = 0.;} // o car on veut une vitesse nulle
|
|
for (int b = 1; b< dim_base; b++)
|
|
{ Coordonnee & d_Vtdt_indice=(*d_Vtdt)(b);
|
|
d_Vtdt_indice.Zero();
|
|
d_Vtdt_indice(b)= unsurdeltat;
|
|
};
|
|
};
|
|
|
|
|
|
//== calcul de la variation des bases
|
|
void MetAxisymetrique2D::D_giB_t( const Mat_pleine& dphi, int nbnoeu,const Vecteur & phi)
|
|
{int indice;
|
|
// derivees des gBi par rapport a XHbr , b=1 à 1 ou 2 suivant la dimension
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b<= nb_ddl_en_var; b++)
|
|
for (int r = 1; r <= nbnoeu; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
BaseB & d_giB_t_indice=(*d_giB_t)(indice);
|
|
// variation du premier vecteur
|
|
d_giB_t_indice.CoordoB(1).Zero();
|
|
d_giB_t_indice.CoordoB(1)(b)=dphi(1,r);
|
|
// le deuxième vecteur à une dérivée non nulle suivant XH1r
|
|
d_giB_t_indice.CoordoB(2).Zero(); // init
|
|
// comme g_2=X^{1r} * phi_r * I_{3 ou 2} suivant que dim_base= 3 ou 2
|
|
if (b==1)
|
|
d_giB_t_indice.CoordoB(2)(dim_base) = phi(r)*2.*ConstMath::Pi;
|
|
};
|
|
// en dim 3 la dérivée de tous les vecteurs par rapport à XH3r est nulle
|
|
// if (dim_base==3)
|
|
// for (int r = 1; r <= nbnoeu; r++)
|
|
// { indice = r*dim_base; // sont quand même rangé dim_base par dim_base
|
|
// BaseB & d_giB_t_indice=(*d_giB_t)(indice);
|
|
// for (int k=1;k<=2;k++)
|
|
// d_giB_t_indice(k).Zero();
|
|
// };
|
|
};
|
|
|
|
void MetAxisymetrique2D::D_giB_tdt( const Mat_pleine& dphi, int nbnoeu,const Vecteur & phi)
|
|
{int indice;
|
|
// derivees des gBi par rapport a XHbr , b=1 à 1 ou 2 suivant la dimension
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
for (int b = 1; b<= nb_ddl_en_var; b++)
|
|
for (int r = 1; r <= nbnoeu; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + b;
|
|
BaseB & d_giB_tdt_indice=(*d_giB_tdt)(indice);
|
|
// variation du premier vecteur
|
|
d_giB_tdt_indice.CoordoB(1).Zero();
|
|
d_giB_tdt_indice.CoordoB(1)(b)=dphi(1,r);
|
|
// le deuxième vecteur à une dérivée non nulle suivant XH1r
|
|
d_giB_tdt_indice.CoordoB(2).Zero(); // init
|
|
// comme g_2=X^{1r} * phi_r * I_{3 ou 2} suivant que dim_base= 3 ou 2
|
|
if (b==1)
|
|
d_giB_tdt_indice.CoordoB(2)(dim_base) = phi(r)*2.*ConstMath::Pi;
|
|
};
|
|
// en dim 3 la dérivée de tous les vecteurs par rapport à XH3r est nulle
|
|
// if (dim_base==3)
|
|
// for (int r = 1; r <= nbnoeu; r++)
|
|
// { indice = r*dim_base; // sont quand même rangé dim_base par dim_base
|
|
// BaseB & d_giB_tdt_indice=(*d_giB_tdt)(indice);
|
|
// for (int k=1;k<=2;k++)
|
|
// d_giB_tdt_indice(k).Zero();
|
|
// };
|
|
};
|
|
|
|
// calcul du gradient de vitesse à t
|
|
void MetAxisymetrique2D::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 << "\nMetAxisymetrique2D::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 l'axe dim qui est la direction
|
|
// tangente à l'axe de rotation qui est y (en 3D) ou z (en 2D)
|
|
// d'autre part en dim 3D, l'interpolation géométrique 1D nous permet d'obtenir le gradient spatial
|
|
// selon l'axe de l'élément (mais pas dans les autres directions évidemment !)
|
|
// de plus la vitesse est normale à g_{dim_base}
|
|
// --> tout d'abord le premier vecteur par rapport à theta 1
|
|
{CoordonneeB & giBt_1 = giB_t->CoordoB(1);
|
|
// pour Vi|j, avec j = 1 (dimension d'interpolation)
|
|
double& gradVij=(*gradVBB_t).Coor(1,1);
|
|
gradVij = 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
{ switch(dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
gradVij += tab_noeud(r)->Valeur_t(V2) * dphi(1,r) * giBt_1(2);
|
|
case 2:
|
|
gradVij += tab_noeud(r)->Valeur_t(V1) * dphi(1,r) * giBt_1(1);
|
|
};
|
|
};// -- fin du for sur r
|
|
}; // fin du bloc de calcul de V1|1
|
|
// --> maintenant le second vecteur, c'est à dire selon la direction tangente à l'axe de rotation
|
|
// la vitesse est nulle
|
|
// idem pour le premier vecteur
|
|
for (int i=1;i<=2;i++)
|
|
// pour Vi|j, avec j=2 -> 0 car V,2 = 0 (en local)
|
|
(*gradVBB_t).Coor(i,2) = 0.;
|
|
// pour la vitesse suivant 1 , elle est également normale à g_2 (du à l'axisymétrie)
|
|
// donc donne un gradient nulle
|
|
// V2|1 = V,1 . g_2 = 0
|
|
(*gradVBB_t).Coor(2,1) = 0.;
|
|
};
|
|
|
|
// calcul gradient de vitesse à t+dt
|
|
void MetAxisymetrique2D::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 << "\nMetAxisymetrique2D::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 l'axe dim qui est la direction
|
|
// tangente à l'axe de rotation qui est y (en 3D) ou z (en 2D)
|
|
// d'autre part en dim 3D, l'interpolation géométrique 1D nous permet d'obtenir le gradient spatial
|
|
// selon l'axe de l'élément (mais pas dans les autres directions évidemment !)
|
|
// de plus la vitesse est normale à g_{dim_base}
|
|
// --> tout d'abord le premier vecteur par rapport à theta 1
|
|
{CoordonneeB & giBtdt_1 = giB_tdt->CoordoB(1);
|
|
// pour Vi|j, avec j = 1 (dimension d'interpolation)
|
|
double& gradVij=(*gradVBB_tdt).Coor(1,1);
|
|
gradVij = 0.;
|
|
for (int r=1;r<=nombre_noeud;r++)
|
|
{ switch(dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
gradVij += tab_noeud(r)->Valeur_tdt(V2) * dphi(1,r) * giBtdt_1(2);
|
|
case 2:
|
|
gradVij += tab_noeud(r)->Valeur_tdt(V1) * dphi(1,r) * giBtdt_1(1);
|
|
};
|
|
};// -- fin du for sur r
|
|
}; // fin du bloc de calcul de V1|1
|
|
// --> maintenant le second vecteur, c'est à dire selon la direction tangente à l'axe de rotation
|
|
// la vitesse est nulle
|
|
// idem pour le premier vecteur
|
|
for (int i=1;i<=2;i++)
|
|
// pour Vi|j, avec j=2 -> 0 car V,2 = 0 (en local)
|
|
(*gradVBB_tdt).Coor(i,2) = 0.;
|
|
// pour la vitesse suivant 1 , elle est également normale à g_2 (du à l'axisymétrie)
|
|
// donc donne un gradient nulle
|
|
// V2|1 = V,1 . g_2 = 0
|
|
(*gradVBB_tdt).Coor(2,1) = 0.;
|
|
};
|
|
|
|
// 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 MetAxisymetrique2D::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 << "\nMetAxisymetrique2D::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 (en 3D) ou y (en 2D)
|
|
// = la direction tangente à l'axe de rotation qui est y en 3D et z en 2D
|
|
// ---> premier vecteur
|
|
{CoordonneeB & giBt_1 = giB_t->CoordoB(1);
|
|
// pour Vi|j, avec j = 1
|
|
double& gradVij=(*gradVmoyBB_t).Coor(1,1);
|
|
switch(dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
gradVij += V(2) * dphi(1,r) * giBt_1(2);
|
|
case 2:
|
|
gradVij += V(1) * dphi(1,r) * giBt_1(1);
|
|
};
|
|
// calcul du vecteurs V,1 moyens, sachant que V,2 = 0
|
|
CoordonneeB & dmatV_moy_t_1= dmatV_moy_t->CoordoB(1);
|
|
switch (dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
dmatV_moy_t_1(2) += V(2) * dphi(1,r);
|
|
case 2:
|
|
dmatV_moy_t_1(1) += V(1) * dphi(1,r);
|
|
}//-- fin du switch
|
|
}; // fin du bloc de calcul de V1|1
|
|
}; // fin de la boucle sur r
|
|
// --> maintenant le second vecteur, c'est à dire selon la direction tangente à l'axe de rotation
|
|
// la vitesse est nulle
|
|
// idem pour le premier vecteur
|
|
// pour la vitesse suivant 1 , elle est également normale à g_2 (du à l'axisymétrie)
|
|
// donc donne un gradient nulle
|
|
// V2|1 = V,1 . g_2 = 0
|
|
// rien n'a faire car gradVmoyBB_t et dmatV_moy_t ont été initialisés à 0.
|
|
};
|
|
|
|
// calcul du gradient de vitesse moyen en utilisant delta x^ar/delta t
|
|
void MetAxisymetrique2D::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 << "\nMetAxisymetrique2D::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_tdt(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 (en 3D) ou y (en 2D)
|
|
// = la direction tangente à l'axe de rotation qui est y en 3D et z en 2D
|
|
// ---> premier vecteur
|
|
{CoordonneeB & giBtdt_1 = giB_tdt->CoordoB(1);
|
|
// pour Vi|j, avec j = 1
|
|
double& gradVij=(*gradVmoyBB_tdt).Coor(1,1);
|
|
switch(dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
gradVij += V(2) * dphi(1,r) * giBtdt_1(2);
|
|
case 2:
|
|
gradVij += V(1) * dphi(1,r) * giBtdt_1(1);
|
|
};
|
|
// calcul du vecteurs V,1 moyens, sachant que V,2 = 0
|
|
CoordonneeB & dmatV_moy_tdt_1= dmatV_moy_tdt->CoordoB(1);
|
|
switch (dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
dmatV_moy_tdt_1(2) += V(2) * dphi(1,r);
|
|
case 2:
|
|
dmatV_moy_tdt_1(1) += V(1) * dphi(1,r);
|
|
}//-- fin du switch
|
|
}; // fin du bloc de calcul de V1|1
|
|
}; // fin de la boucle sur r
|
|
// --> maintenant le second vecteur, c'est à dire selon la direction tangente à l'axe de rotation
|
|
// la vitesse est nulle
|
|
// idem pour le premier vecteur
|
|
// pour la vitesse suivant 1 , elle est également normale à g_2 (du à l'axisymétrie)
|
|
// donc donne un gradient nulle
|
|
// V2|1 = V,1 . g_2 = 0
|
|
// rien n'a faire car gradVmoyBB_tdt et dmatV_moy_tdt ont été initialisés à 0.
|
|
};
|
|
|
|
//== calcul de la variation du gradient instantannée
|
|
|
|
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
// dans les dérivées des gradients il manque les variations par rapport à V2, cf les gi
|
|
// absolument les mettres à jour !!!!!!!!!!!!
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
//****************************************************************************************************
|
|
|
|
|
|
|
|
// par rapport aux composantes V^ar (et non les X^ar )
|
|
void MetAxisymetrique2D::DgradVBB_t(const Mat_pleine& dphi)
|
|
{ // derivees par rapport a VH1r
|
|
// Pour les éléments axisymétriques, il n'y a pas de vitesse suivant l'axe dim qui est la direction
|
|
// tangente à l'axe de rotation qui est y (en 3D) ou z (en 2D)
|
|
// seule la composante 1,1 du gradient est non nulle
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
int nddl = d_gradVBB_t.Taille();
|
|
for (int i=1;i<=nddl;i++)
|
|
(d_gradVBB_t(i))->Inita(0.);
|
|
// les dérivées non-nulles
|
|
int indice=0;
|
|
for (int r = 1; r <= nomb_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + 1;
|
|
// grad(i,j) = Vi|j = gj_point * gi
|
|
{ switch(dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
(*(d_gradVBB_t(indice+1))).Coor(1,1) = dphi(1,r) * (*giB_t)(1)(2);
|
|
case 2:
|
|
(*(d_gradVBB_t(indice))).Coor(1,1) = dphi(1,r) * (*giB_t)(1)(1);
|
|
};
|
|
};
|
|
};
|
|
};
|
|
|
|
//== calcul de la variation du gradient
|
|
// par rapport aux composantes V^ar (et non les X^ar )
|
|
void MetAxisymetrique2D::DgradVBB_tdt(const Mat_pleine& dphi)
|
|
{ // derivees par rapport a VH1r
|
|
// Pour les éléments axisymétriques, il n'y a pas de vitesse suivant l'axe dim qui est la direction
|
|
// tangente à l'axe de rotation qui est y (en 3D) ou z (en 2D)
|
|
// seule la composante 1,1 du gradient est non nulle
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
int nddl = d_gradVBB_tdt.Taille();
|
|
for (int i=1;i<=nddl;i++)
|
|
(d_gradVBB_tdt(i))->Inita(0.);
|
|
// les dérivées non-nulles
|
|
int indice=0;
|
|
for (int r = 1; r <= nomb_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + 1;
|
|
// grad(i,j) = Vi|j = gj_point * gi
|
|
{ switch(dim_base)
|
|
{case 3:
|
|
// on ne tiens pas compte de la valeur V3 qui est nulle
|
|
(*(d_gradVBB_tdt(indice+1))).Coor(1,1) = dphi(1,r) * (*giB_tdt)(1)(2);
|
|
case 2:
|
|
(*(d_gradVBB_tdt(indice))).Coor(1,1) = dphi(1,r) * (*giB_tdt)(1)(1);
|
|
};
|
|
};
|
|
};
|
|
};
|
|
|
|
// calcul de la variation du gradient moyen à t
|
|
// par rapport aux composantes X^ar (et non les V^ar )
|
|
void MetAxisymetrique2D::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
|
|
};
|
|
// initialisation
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
int nddl = d_gradVmoyBB_t.Taille();
|
|
for (int i=1;i<=nddl;i++)
|
|
(d_gradVmoyBB_t(i))->Inita(0.);
|
|
// seul les gradients (1,1) sont non nuls
|
|
// les dérivées non-nulles
|
|
int indice=0;
|
|
for (int r = 1; r <= nomb_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + 1;
|
|
// grad(i,j) = Vi|j = gj_point * gi
|
|
{ switch(dim_base)
|
|
{case 3:
|
|
(*(d_gradVmoyBB_t(indice+1))).Coor(1,1)
|
|
= dphi(1,r) *( unsurdeltat * (*giB_t)(1)(2)
|
|
+ (*dmatV_moy_t)(1).ScalBB((*d_giB_t)(indice+1)(1)) );
|
|
case 2:
|
|
(*(d_gradVmoyBB_t(indice))).Coor(1,1)
|
|
= dphi(1,r) *( unsurdeltat * (*giB_t)(1)(1)
|
|
+ (*dmatV_moy_t)(1).ScalBB((*d_giB_t)(indice)(1)) );
|
|
};
|
|
};
|
|
};
|
|
};
|
|
|
|
// calcul de la variation du gradient moyen à t dt
|
|
// par rapport aux composantes X^ar (et non les V^ar )
|
|
void MetAxisymetrique2D::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
|
|
// initialisation
|
|
int nb_ddl_en_var = dim_base-1; // on ne prend en compte que les variations jusqu'a dim-1
|
|
int nddl = d_gradVmoyBB_tdt.Taille();
|
|
for (int i=1;i<=nddl;i++)
|
|
(d_gradVmoyBB_tdt(i))->Inita(0.);
|
|
// seul les gradients (1,1) sont non nuls
|
|
// les dérivées non-nulles
|
|
int indice=0;
|
|
for (int r = 1; r <= nomb_noeud; r++)
|
|
{ indice = (r-1)*nb_ddl_en_var + 1;
|
|
// grad(i,j) = Vi|j = gj_point * gi
|
|
{ switch(dim_base)
|
|
{case 3:
|
|
(*(d_gradVmoyBB_tdt(indice+1))).Coor(1,1)
|
|
= dphi(1,r) *( unsurdeltat * (*giB_tdt)(1)(2)
|
|
+ (*dmatV_moy_tdt)(1).ScalBB((*d_giB_tdt)(indice+1)(1)) );
|
|
case 2:
|
|
(*(d_gradVmoyBB_tdt(indice))).Coor(1,1)
|
|
= dphi(1,r) *( unsurdeltat * (*giB_tdt)(1)(1)
|
|
+ (*dmatV_moy_tdt)(1).ScalBB((*d_giB_tdt)(indice)(1)) );
|
|
};
|
|
};
|
|
};
|
|
};
|