2058 lines
96 KiB
C++
2058 lines
96 KiB
C++
|
|
// 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 "Debug.h"
|
|
# include "Deformation.h"
|
|
#include "ConstMath.h"
|
|
#include "ParaGlob.h"
|
|
#include "MathUtil.h"
|
|
#include "Tenseur3.h"
|
|
#include "Tenseur2.h"
|
|
#include "Tenseur1.h"
|
|
# include "TypeConsTens.h"
|
|
#include "Util.h"
|
|
#include "MathUtil2.h"
|
|
|
|
|
|
//================== initialisation des variables static ======================
|
|
// change 2 mai 2020: int Deformation::numInteg = 0;
|
|
// indicateur utilisé par VerifCal_deflog
|
|
int Deformation::indic_VerifCal_implicit = 0;
|
|
// def de deux tableaux d'indices pour avoir la suite: 12 21 13 31 23 32 à partir de 6 nombres
|
|
const Deformation::UtilIndexDeformation Deformation::ccdex;
|
|
// partie statique ou sont stocké tous les cas différents des variables de travail
|
|
list <Tableau <TenseurBH * > > Deformation::list_var_tens_BH;
|
|
list <Tableau <Tableau <TenseurBH* > > > Deformation::list_var_var_tens_BH;
|
|
list <Tableau <Coordonnee > > Deformation::list_tab_coor;
|
|
int Deformation::nombre_d_instance_deformation=0;
|
|
|
|
//================== fin d'initialisation des variables static ================
|
|
|
|
// def du constructeur par défaut de UtilIndexDeformation
|
|
Deformation::UtilIndexDeformation::UtilIndexDeformation() :
|
|
iii(6),jjj(6)
|
|
{ iii(6); iii(1)=1; iii(2)=2;iii(3)=1;iii(4)=3;iii(5)=2;iii(6)=3;
|
|
jjj(6); jjj(1)=2; jjj(2)=1;jjj(3)=3;jjj(4)=1;jjj(5)=3;jjj(6)=2;
|
|
};
|
|
|
|
// ---------- constructeur----------------------------------------
|
|
Deformation::Deformation () // constructeur ne devant pas etre utilise
|
|
//#ifndef ENLINUX
|
|
: tabnoeud((new Tableau <Noeud *> )) // dans le cas de linux on ne déclare pas le tableau
|
|
,saveDefResul(NULL),numInteg(0),sauve_numInteg(0)
|
|
//#endif
|
|
|
|
// : tabnoeud()
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
{ cout << "\nErreur : le constructeur par defaut ne doit pa etre utilise !\n";
|
|
cout << "Deformation::Deformation () \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
};
|
|
|
|
// constructeur normal dans le cas d'un ou de plusieurs pt d'integration
|
|
Deformation::Deformation (Met_abstraite & a,Tableau<Noeud *>& b,
|
|
Tableau <Mat_pleine> const & tabdphi,Tableau <Vecteur> const & tabphi
|
|
,Enum_type_deformation type_de_deformation) :
|
|
tabnoeud(&b),type_deformation(type_de_deformation)
|
|
,type_gradient_vitesse(GRADVITESSE_VCONST) // pour l'instant ********* en test
|
|
,saveDefResul(NULL)
|
|
// variables de travail pour log
|
|
,B_BH_tr(NULL),d_B_BH_tr(NULL),ki(),d_ki(NULL),Palpha_BH_tr(NULL),d_Palpha_BH_tr(NULL)
|
|
,dimensionnement_tableaux(0),numInteg(0)
|
|
,sauve_numInteg(0)
|
|
{ metrique = &a;
|
|
tabDphi = &tabdphi;
|
|
tabPhi = &tabphi;
|
|
nbNoeud = (*tabPhi)(1).Taille(); // par defaut
|
|
numInteg = 1; // 1 pt d'integ par defaut
|
|
nombre_d_instance_deformation++;
|
|
};
|
|
|
|
// constructeur de copie
|
|
Deformation::Deformation (const Deformation& a) :
|
|
tabnoeud(a.tabnoeud), type_deformation(a.type_deformation)
|
|
,type_gradient_vitesse(a.type_gradient_vitesse),saveDefResul(NULL)
|
|
// variables de travail pour log
|
|
,B_BH_tr(a.B_BH_tr),d_B_BH_tr(a.d_B_BH_tr),ki(a.ki),d_ki(a.d_ki)
|
|
,Palpha_BH_tr(a.Palpha_BH_tr),d_Palpha_BH_tr(a.d_Palpha_BH_tr)
|
|
,dimensionnement_tableaux(a.dimensionnement_tableaux)
|
|
,numInteg(a.numInteg),sauve_numInteg(a.sauve_numInteg)
|
|
{ metrique = a.metrique;
|
|
tabDphi = a.tabDphi;
|
|
tabPhi = a.tabPhi;
|
|
nbNoeud = a.nbNoeud;
|
|
nombre_d_instance_deformation++;
|
|
};
|
|
|
|
Deformation::~Deformation ()
|
|
{ LibereTenseur() ; nombre_d_instance_deformation--;
|
|
if (nombre_d_instance_deformation==0)
|
|
{// cas de la dernière instance
|
|
// ---- cas de list_var_tens_BH ----
|
|
list <Tableau <TenseurBH * > >::iterator ilv,ilvend= list_var_tens_BH.end();
|
|
for (ilv=list_var_tens_BH.begin();ilv!=ilvend;ilv++)
|
|
{ int tailletab=(*ilv).Taille();
|
|
for (int i=1;i<=tailletab;i++) delete (*ilv)(i);
|
|
};
|
|
list_var_tens_BH.erase(list_var_tens_BH.begin(),ilvend);
|
|
// ---- cas de list_var_var_tens_BH ----
|
|
list <Tableau <Tableau <TenseurBH* > > > ::iterator ilv1,ilvend1= list_var_var_tens_BH.end();
|
|
for (ilv1=list_var_var_tens_BH.begin();ilv1!=ilvend1;ilv1++)
|
|
{ int tailletab=(*ilv1).Taille();
|
|
for (int i=1;i<=tailletab;i++)
|
|
{ int tailletab2=(*ilv1)(i).Taille();
|
|
for (int j=1;j<=tailletab2;j++) delete (*ilv1)(i)(j);
|
|
};
|
|
};
|
|
list_var_var_tens_BH.erase(list_var_var_tens_BH.begin(),ilvend1);
|
|
// ---- cas de list_tab_coor ----
|
|
list_tab_coor.erase(list_tab_coor.begin(),list_tab_coor.end());
|
|
};
|
|
};
|
|
|
|
// réafectation du pointeur de tableau de noeuds interne
|
|
// il faut que le nouveau tableau ait la même taille que l'ancien (sert pour créer de nouveau éléments
|
|
// de même type, sinon ce n'est pas la bonne solution,
|
|
// il vaut mieux creer une nouvelle def
|
|
void Deformation::PointeurTableauNoeud(Tableau <Noeud *> & tabN)
|
|
{ // on vérifie la cohérence
|
|
if (tabN.Taille() != tabnoeud->Taille())
|
|
{ cout << "\n erreur de re-affectation de tableau de noeud dans la deformation ! "
|
|
<< " ancien nombre de noeud: " << tabnoeud->Taille() << " nouveau nombre: " << tabN.Taille()
|
|
<< "\n Deformation::PointeurTableauNoeud(..";
|
|
Sortie(1);
|
|
};
|
|
// sinon c'est ok
|
|
tabnoeud = &tabN;
|
|
};
|
|
|
|
// surcharge de l'opérateur d'affectation
|
|
Deformation& Deformation::operator= (const Deformation& def)
|
|
{ metrique = def.metrique;
|
|
tabnoeud = def.tabnoeud;
|
|
tabDphi = def.tabDphi;
|
|
tabPhi = def.tabPhi;
|
|
nbNoeud = def.nbNoeud;type_deformation=def.type_deformation;
|
|
numInteg = def.numInteg;
|
|
saveDefResul = def.saveDefResul;
|
|
// variables de travail pour log
|
|
B_BH_tr=def.B_BH_tr;d_B_BH_tr=def.d_B_BH_tr;ki=def.ki;d_ki=def.d_ki;
|
|
Palpha_BH_tr=def.Palpha_BH_tr;d_Palpha_BH_tr=def.d_Palpha_BH_tr;
|
|
|
|
return (*this);
|
|
};
|
|
|
|
// change le type de déformation
|
|
void Deformation::Change_type_de_deformation(Enum_type_deformation type_de_def)
|
|
{ // le changement de déformation doit-être cohérent avec la déformation actuelle d'où les
|
|
// différents choix restreinds
|
|
if ((type_de_def==DEFORMATION_STANDART) ||
|
|
(type_de_def==DEFORMATION_LOGARITHMIQUE) ||
|
|
(type_de_def==DEFORMATION_CUMU_LOGARITHMIQUE) ||
|
|
(type_de_def==DEF_CUMUL_CORROTATIONNEL) ||
|
|
(type_de_def==DEF_CUMUL_ROTATION_PROPRE) )
|
|
{ if ((type_deformation==DEFORMATION_STANDART) ||
|
|
(type_deformation==DEFORMATION_LOGARITHMIQUE) ||
|
|
(type_deformation==DEFORMATION_CUMU_LOGARITHMIQUE) ||
|
|
(type_deformation==DEF_CUMUL_CORROTATIONNEL) ||
|
|
(type_deformation==DEF_CUMUL_ROTATION_PROPRE) )
|
|
type_deformation = type_de_def; // ok cohérent
|
|
}
|
|
// sinon le type n'est pas pris en compte dans cette méthode donc erreur
|
|
else
|
|
{ cout << "\nErreur : valeur incorrecte du nouvrau type Enum_type_deformation !"
|
|
<< "\n type actuelle: " << Nom_type_deformation(type_deformation)
|
|
<< " type demande: " << Nom_type_deformation(type_de_def);
|
|
cout << "Deformation::Change_type_de_deformation(...) \n";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// affichage des informations
|
|
void Deformation::Affiche() const
|
|
{ cout << "\n -- deformation: --- "
|
|
<< "\n nb noeud concerne = "<< nbNoeud
|
|
<< " numInteg= "<< numInteg;
|
|
if (tabnoeud != NULL)
|
|
{cout << "\n -- information concernant les noeuds: \n";
|
|
for (int i = 1;i<= tabnoeud->Taille();i++)
|
|
(*tabnoeud)(i)->Affiche();
|
|
};
|
|
if (metrique != NULL)
|
|
{cout << "\n -- information concernant la metrique actuelle au pti --";
|
|
metrique->Affiche();
|
|
};
|
|
cout << flush;
|
|
};
|
|
|
|
|
|
// ============ METHODES PUBLIQUES : ==================
|
|
// ---------------- calcul des variables primaires pour la mécanique --------
|
|
// calcul explicit à t :tous les parametres sont des resultats
|
|
const Met_abstraite::Expli& Deformation::Cal_explicit_t
|
|
// ( bool gradV_instantane,TenseurBB & epsBB_t,Tableau <TenseurBB *> & d_epsBB
|
|
// ,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
|
|
// def_equi_t: est la def equi au debut du pas, et def_equi est la def finale
|
|
( const Tableau <double>& def_equi_t,TenseurBB & epsBB_t,Tableau <TenseurBB *> & d_epsBB
|
|
,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
|
|
{bool gradV_instantane = false; // ************ pour l'instant figé
|
|
// appel de la metrique: on le fait ici et non dans les méthodes Cal_explicit_Almansi_tdt etc
|
|
// car pour d'autre déformation, par exemple sfe ou plaque, la méthode Deformation::Cal_explicit_tdt
|
|
// est surchargée car différentes, par contre les méthodes Cal_explicit_Almansi_tdt etc sont les mêmes
|
|
const Met_abstraite::Expli& ex =
|
|
metrique->Cal_explicit_t(*tabnoeud,gradV_instantane,(*tabDphi)(numInteg),
|
|
nbNoeud,(*tabPhi)(numInteg),premier_calcul);
|
|
|
|
// ici il y a le choix entre les différents types de calcul de la déformation (Almansi, Green-Lagrange ..
|
|
switch (type_deformation)
|
|
{ case DEFORMATION_STANDART :
|
|
{Cal_explicit_Almansi
|
|
(gradV_instantane,epsBB_t,d_epsBB,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
case DEFORMATION_LOGARITHMIQUE:
|
|
{Cal_explicit_Logarithmique
|
|
(gradV_instantane,epsBB_t,d_epsBB,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
|
|
// case DEFORMATION_CUMU_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEF_CUMUL_ROTATION_PROPRE :
|
|
// {Cal_explicit_def_cumule
|
|
// (gradV_instantane,epsBB_t,d_epsBB,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
// break;
|
|
// }
|
|
default :
|
|
cout << "\nErreur : type de deformation ( " << Nom_type_deformation(type_deformation)
|
|
<< " ) non traite !\n";
|
|
cout << "Deformation::Cal_explicit(... \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
// --- on calcul la déformation cumulée
|
|
|
|
{TenseurBH * delta_epsBH = NevezTenseurBH(delta_epsBB.Dimension(), 0.);
|
|
*delta_epsBH = delta_epsBB * (*ex.gijHH_t);
|
|
double delta_eps_equi = sqrt(2./3. * ( ((*delta_epsBH) && (*delta_epsBH)) - Sqr(delta_epsBH->Trace()) /3. ));
|
|
def_equi(1) = def_equi_t(1) + delta_eps_equi;
|
|
def_equi(4) = delta_eps_equi;
|
|
delete delta_epsBH;
|
|
};
|
|
{TenseurBH * epsBH = NevezTenseurBH(epsBB_t.Dimension(), 0.);
|
|
*epsBH = epsBB_t * (*ex.gijHH_t);
|
|
def_equi(2) = sqrt(2./3. * ( ((*epsBH) && (*epsBH)) - Sqr(epsBH->Trace()) /3. ));
|
|
delete epsBH;
|
|
if (def_equi(2) > def_equi_t(3))
|
|
def_equi(3) = def_equi(2);
|
|
};
|
|
|
|
// sauvegarde des infos à 0 éventuellement et init par recopie des grandeurs sur t
|
|
if (premier_calcul) saveDefResul->MiseAJourGrandeurs_a_0(metrique);
|
|
// *** pour supprimer les warnings à la compilation *** mais on ne doit jamais passer ici
|
|
// const Met_abstraite::Expli& toto = * (new Met_abstraite::Expli()); return toto;
|
|
return ex;
|
|
};
|
|
|
|
// calcul explicit à tdt :tous les parametres sont des resultats
|
|
const Met_abstraite::Expli_t_tdt& Deformation::Cal_explicit_tdt
|
|
// ( bool gradV_instantane,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB,TenseurBB& DepsBB
|
|
// ,TenseurBB& delta_epsBB_tdt,bool premier_calcul)
|
|
( const Tableau <double>& def_equi_t,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB
|
|
,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& delta_epsBB_tdt,bool premier_calcul)
|
|
{ bool gradV_instantane = false; // ************ pour l'instant figé
|
|
// appel de la metrique: on le fait ici et non dans les méthodes Cal_explicit_Almansi_tdt etc
|
|
// car pour d'autre déformation, par exemple sfe ou plaque, la méthode Deformation::Cal_explicit_tdt
|
|
// est surchargée car différentes, par contre les méthodes Cal_explicit_Almansi_tdt etc sont les mêmes
|
|
const Met_abstraite::Expli_t_tdt& ex =
|
|
metrique->Cal_explicit_tdt(*tabnoeud,gradV_instantane,(*tabDphi)(numInteg),
|
|
nbNoeud,(*tabPhi)(numInteg),premier_calcul);
|
|
|
|
// ici il y a le choix entre les différents types de calcul de la déformation (Almansi, Green-Lagrange ..
|
|
switch (type_deformation)
|
|
{ case DEFORMATION_STANDART :
|
|
{Cal_explicit_Almansi_tdt
|
|
(gradV_instantane,epsBB_tdt,d_epsBB,DepsBB,delta_epsBB_tdt,premier_calcul,ex);
|
|
break;
|
|
}
|
|
case DEFORMATION_LOGARITHMIQUE:
|
|
{Cal_explicit_logarithmique_tdt
|
|
(gradV_instantane,epsBB_tdt,d_epsBB,DepsBB,delta_epsBB_tdt,premier_calcul,ex);
|
|
break;
|
|
}
|
|
// case DEFORMATION_CUMU_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEF_CUMUL_ROTATION_PROPRE :
|
|
// {Cal_explicit_def_cumule_tdt
|
|
// (gradV_instantane,epsBB_tdt,d_epsBB,DepsBB,delta_epsBB_tdt,premier_calcul,ex);
|
|
// break;
|
|
// }
|
|
default :
|
|
cout << "\nErreur : type de deformation ( " << Nom_type_deformation(type_deformation)
|
|
<< " ) non traite !\n";
|
|
cout << "Deformation::Cal_explicit_tdt(... \n";
|
|
Sortie(1);
|
|
};
|
|
// sauvegarde des infos à 0 éventuellement et init par recopie des grandeurs sur t
|
|
if (premier_calcul) saveDefResul->MiseAJourGrandeurs_a_0(metrique);
|
|
// sauvegarde des infos à t à chaque passage
|
|
saveDefResul->MiseAJourGrandeurs_a_tdt(metrique,DepsBB);
|
|
// --- on calcul la déformation cumulée
|
|
{TenseurBH * delta_epsBH = NevezTenseurBH(delta_epsBB_tdt.Dimension(), 0.);
|
|
*delta_epsBH = delta_epsBB_tdt * (*ex.gijHH_tdt);
|
|
double delta_eps_equi = sqrt(2./3. * ( ((*delta_epsBH) && (*delta_epsBH)) - Sqr(delta_epsBH->Trace()) /3. ));
|
|
def_equi(1) = def_equi_t(1) + delta_eps_equi;
|
|
def_equi(4) = delta_eps_equi;
|
|
delete delta_epsBH;
|
|
};
|
|
|
|
{TenseurBH * epsBH = NevezTenseurBH(epsBB_tdt.Dimension(), 0.);
|
|
*epsBH = epsBB_tdt * (*ex.gijHH_tdt);
|
|
def_equi(2) = sqrt(2./3. * ( ((*epsBH) && (*epsBH)) - Sqr(epsBH->Trace()) /3. ));
|
|
delete epsBH;
|
|
if (def_equi(2) > def_equi_t(3))
|
|
def_equi(3) = def_equi(2);
|
|
};
|
|
|
|
// *** pour supprimer les warnings à la compilation *** mais on ne doit jamais passer ici
|
|
// const Met_abstraite::Expli_t_tdt& toto = * (new Met_abstraite::Expli_t_tdt()); return toto;
|
|
return ex;
|
|
};
|
|
|
|
// calcul et récupération de la dérivée de la vitesse de déformation virtuelle à tdt
|
|
// dans le cas d'une discrétisation classique simple, ce calcul est indépendant des coordonnées
|
|
// donc peut être effectué indépendament du reste, il dépend cependant du point d'intégration
|
|
// dans le cas de discrétisation complexe, il faut s'assurer que les autres grandeurs utiles de la métrique
|
|
// sont déjà calculé (cf. D2_gijBB_tdt() )
|
|
// total = true : indique que le calcul est complet, à partir des données de base
|
|
// c-a-d calcul de la variation de vecteur de base puis calcul de la variation seconde
|
|
// = false: le calcul est simplifié, on considère que la variation des vecteurs de base vient
|
|
// juste d'être calculé, dans un calcul implicit, on ne calcul alors que la variation seconde
|
|
void Deformation::Cal_var_def_virtuelle(bool total,Tableau2 <TenseurBB *>& d2_epsBB_tdt)
|
|
{ const Tableau2 <TenseurBB * > & d2_gijBB_tdt
|
|
= metrique->CalVar2GijBBtdt(total, (*tabDphi)(numInteg),nbNoeud,(*tabPhi)(numInteg));
|
|
int d2_epsBB_tdtTaille1 = d2_epsBB_tdt.Taille1();
|
|
for (int i=1; i<= d2_epsBB_tdtTaille1; i++)
|
|
for (int j=1; j<= i; j++) // symétrie du tableau et du tenseur
|
|
{ *(d2_epsBB_tdt(i,j)) = 0.5 * (*((d2_gijBB_tdt)(i,j))) ;
|
|
*(d2_epsBB_tdt(j,i)) = *(d2_epsBB_tdt(i,j)) ;
|
|
};
|
|
};
|
|
|
|
// cas implicite : tous les parametres sont de resultats
|
|
//const Met_abstraite::Impli& Deformation::Cal_implicit( bool gradV_instantane
|
|
// ,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
|
|
// ,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
|
|
const Met_abstraite::Impli& Deformation::Cal_implicit(const Tableau <double>& def_equi_t
|
|
,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
|
|
,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
|
|
{ bool gradV_instantane = false; // ************ pour l'instant figé
|
|
// appel de la metrique: on le fait ici et non dans les méthodes Cal_explicit_Almansi_tdt etc
|
|
// car pour d'autre déformation, par exemple sfe ou plaque, la méthode Deformation::Cal_explicit_tdt
|
|
// est surchargée car différentes, par contre les méthodes Cal_explicit_Almansi_tdt etc sont les mêmes
|
|
// toutes les variables de passage a metrique apres l'appel
|
|
// pointeront sur des variables deja dimensionnees
|
|
// pour les Tableau <> il y a dim
|
|
//------debug
|
|
//premier_calcul=true;
|
|
//---- fin debug
|
|
const Met_abstraite::Impli& ex =metrique->Cal_implicit
|
|
( *tabnoeud,gradV_instantane,(*tabDphi)(numInteg),nbNoeud,(*tabPhi)(numInteg)
|
|
,premier_calcul);
|
|
// ici il y a le choix entre les différents types de calcul de la déformation (Almansi, Green-Lagrange ..
|
|
switch (type_deformation)
|
|
{ case DEFORMATION_STANDART :
|
|
{Cal_implicit_Almansi
|
|
(gradV_instantane,epsBB_tdt,d_epsBB_tdt,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
case DEFORMATION_LOGARITHMIQUE:
|
|
{Cal_implicit_Logarithmique
|
|
(gradV_instantane,epsBB_tdt,d_epsBB_tdt,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
case DEFORMATION_CUMU_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEF_CUMUL_ROTATION_PROPRE :
|
|
{Cal_implicit_def_cumule
|
|
(gradV_instantane,epsBB_tdt,d_epsBB_tdt,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
default :
|
|
cout << "\nErreur : type de deformation ( " << Nom_type_deformation(type_deformation)
|
|
<< " ) non traite !\n";
|
|
cout << "Deformation::Cal_implicit(... \n";
|
|
Sortie(1);
|
|
};
|
|
// sauvegarde des infos à 0 et à t
|
|
// *** non car dans le cas d'un restart ce n'est pas vrai !!!! // éventuellement et init par recopie des grandeurs sur t
|
|
if (premier_calcul) saveDefResul->MiseAJourGrandeurs_a_0(metrique);
|
|
// sauvegarde des infos à tdt à chaque passage
|
|
saveDefResul->MiseAJourGrandeurs_a_tdt(metrique,DepsBB);
|
|
// *** pour supprimer les warnings à la compilation *** mais on ne doit jamais passer ici
|
|
// const Met_abstraite::Impli& toto = * (new Met_abstraite::Impli()); return toto;
|
|
|
|
// --- on calcul la déformation cumulée
|
|
{TenseurBH * delta_epsBH = NevezTenseurBH(epsBB_tdt.Dimension(), 0.);
|
|
*delta_epsBH = delta_epsBB * (*ex.gijHH_tdt);
|
|
double x_inter = (2./3. * ( ((*delta_epsBH) && (*delta_epsBH)) - Sqr(delta_epsBH->Trace()) /3. ));
|
|
if (x_inter < 0)
|
|
{if (Dabs(x_inter) < ConstMath::trespetit)
|
|
{ //debug
|
|
//cout << "\n **************** x_inter= "<<x_inter<<endl;
|
|
//fin debug
|
|
x_inter = 0.;
|
|
}
|
|
else
|
|
{ if (ParaGlob::NiveauImpression() > 0)
|
|
cout << "\n erreur calcul def equi: delta_def_equi2 est negatif ("<<x_inter<<") la racine -> nan! "
|
|
<< " on met a 0, mais ce n'est pas normal !! ";
|
|
if (ParaGlob::NiveauImpression() > 5)
|
|
cout << "\n Deformation::Cal_implicit(.." << endl;
|
|
};
|
|
};
|
|
double delta_eps_equi = sqrt(x_inter);
|
|
def_equi(1) = def_equi_t(1) + delta_eps_equi;
|
|
def_equi(4) = delta_eps_equi;
|
|
delete delta_epsBH;
|
|
};
|
|
|
|
{TenseurBH * epsBH = NevezTenseurBH(epsBB_tdt.Dimension(), 0.);
|
|
*epsBH = epsBB_tdt * (*ex.gijHH_tdt);
|
|
|
|
double y_inter = 2./3. * ( ((*epsBH) && (*epsBH)) - Sqr(epsBH->Trace()) /3. );
|
|
if (y_inter < 0)
|
|
{if (Dabs(y_inter) < ConstMath::trespetit)
|
|
{ //debug
|
|
//cout << "\n **************** x_inter= "<<x_inter<<endl;
|
|
//fin debug
|
|
y_inter = 0.;
|
|
}
|
|
else
|
|
{ if (ParaGlob::NiveauImpression() > 0)
|
|
cout << "\n erreur calcul def equi: sqr(def_dual_mises est negatif) ("<<y_inter<<") la racine -> nan! "
|
|
<< " on met a 0, mais ce n'est pas normal !! ";
|
|
if (ParaGlob::NiveauImpression() > 5)
|
|
cout << "\n Deformation::Cal_implicit(.." << endl;
|
|
};
|
|
};
|
|
def_equi(2) = sqrt(y_inter);
|
|
delete epsBH;
|
|
if (def_equi(2) > def_equi_t(3))
|
|
def_equi(3) = def_equi(2);
|
|
};
|
|
|
|
return ex;
|
|
};
|
|
|
|
// ------------------------ flambage linéaire ------------------------------
|
|
const Met_abstraite::flambe_lin& Deformation::Cal_flambe_lin
|
|
// ( bool gradV_instantane,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt,Tableau2 <TenseurBB *>& d2_epsBB_tdt
|
|
// ,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
|
|
( const Tableau <double>& def_equi_t,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
|
|
,Tableau <double>& def_equi,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul)
|
|
{ bool gradV_instantane = false; // ************ pour l'instant figé
|
|
|
|
// --- le calcul est identique au cas implicite sauf que
|
|
// l'on doit absolument calculer la dérivée seconde de la déformation
|
|
|
|
// appel de la metrique: on le fait ici et non dans les méthodes Cal_explicit_Almansi_tdt etc
|
|
// car pour d'autre déformation, par exemple sfe ou plaque, la méthode Deformation::Cal_explicit_tdt
|
|
// est surchargée car différentes, par contre les méthodes Cal_explicit_Almansi_tdt etc sont les mêmes
|
|
// toutes les variables de passage a metrique apres l'appel
|
|
// pointeront sur des variables deja dimensionnees
|
|
// pour les Tableau <> il y a dimensionnement auto a l'affectation
|
|
const Met_abstraite::Impli& ex =metrique->Cal_implicit
|
|
( *tabnoeud,gradV_instantane,(*tabDphi)(numInteg),nbNoeud,(*tabPhi)(numInteg)
|
|
,premier_calcul);
|
|
// ici il y a le choix entre les différents types de calcul de la déformation (Almansi, Green-Lagrange ..
|
|
switch (type_deformation)
|
|
{ case DEFORMATION_STANDART :
|
|
{Cal_implicit_Almansi(gradV_instantane,epsBB_tdt,d_epsBB_tdt,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
case DEFORMATION_LOGARITHMIQUE:
|
|
{Cal_implicit_Logarithmique(gradV_instantane,epsBB_tdt,d_epsBB_tdt,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
break;
|
|
}
|
|
// case DEFORMATION_CUMU_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEF_CUMUL_ROTATION_PROPRE :
|
|
// {Cal_implicit_def_cumule
|
|
// (gradV_instantane,epsBB_tdt,d_epsBB_tdt,DepsBB,delta_epsBB,premier_calcul,ex);
|
|
// break;
|
|
// }
|
|
default :
|
|
cout << "\nErreur : type de deformation ( " << Nom_type_deformation(type_deformation)
|
|
<< " ) non traite !\n";
|
|
cout << "Deformation::Cal_flambe_lin(... \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
// --- on calcul la déformation cumulée
|
|
{TenseurBH * delta_epsBH = NevezTenseurBH(epsBB_tdt.Dimension(), 0.);
|
|
*delta_epsBH = delta_epsBB * (*ex.gijHH_tdt);
|
|
double delta_eps_equi = sqrt(2./3. * ( ((*delta_epsBH) && (*delta_epsBH)) - Sqr(delta_epsBH->Trace()) /3. ));
|
|
def_equi(1) = def_equi_t(1) + delta_eps_equi;
|
|
def_equi(4) = delta_eps_equi;
|
|
delete delta_epsBH;
|
|
};
|
|
|
|
{TenseurBH * epsBH = NevezTenseurBH(epsBB_tdt.Dimension(), 0.);
|
|
*epsBH = epsBB_tdt * (*ex.gijHH_tdt);
|
|
def_equi(2) = sqrt(2./3. * ( ((*epsBH) && (*epsBH)) - Sqr(epsBH->Trace()) /3. ));
|
|
delete epsBH;
|
|
if (def_equi(2) > def_equi_t(3))
|
|
def_equi(3) = def_equi(2);
|
|
};
|
|
|
|
// sauvegarde des infos à 0 éventuellement et init par recopie des grandeurs sur t
|
|
if (premier_calcul) saveDefResul->MiseAJourGrandeurs_a_0(metrique);
|
|
// sauvegarde des infos à tdt à chaque passage
|
|
saveDefResul->MiseAJourGrandeurs_a_tdt(metrique,DepsBB);
|
|
// *** pour supprimer les warnings à la compilation *** mais on ne doit jamais passer ici
|
|
// const Met_abstraite::flambe_lin& toto = * (new Met_abstraite::flambe_lin()); return toto;
|
|
return (Met_abstraite::flambe_lin&) ex;
|
|
};
|
|
|
|
// ------------------------- cas de la dilatation thermique ------------------------
|
|
// l'objectif est de calculer la déformation d'origine thermique et la déformation d'origine mécanique
|
|
// temperature_tdt,temperature_t : températures, temperature_0: température initiale
|
|
// atdt : booléen qui indique si les grandeurs à tdt sont disponible ou pas
|
|
// avec_repercution_sur_def_meca: si oui on met à jour la déformation méca, sinon, on ne calcul que les def thermiques
|
|
void Deformation::DeformationThermoMecanique(const double& temperature_0,const ThermoDonnee& dTP
|
|
,const TenseurBB & epsBB_totale,TenseurBB & epsBB_therm,const double& temperature_tdt,TenseurBB & epsBB_meca
|
|
,const TenseurBB & delta_epsBB_totale,TenseurBB & delta_epsBB_therm,TenseurBB & delta_epsBB_meca
|
|
,const double& temperature_t,TenseurBB & DepsBB_totale,TenseurBB & DepsBB_therm,TenseurBB & DepsBB_meca
|
|
,const Met_abstraite::Impli* ex_impli
|
|
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
|
|
,const Met_abstraite::Umat_cont* ex_umat
|
|
,bool atdt,bool avec_repercution_sur_def_meca)
|
|
{ // la dilatation est déterminée par l'élévation de température
|
|
double coef_dila_global = dTP.Dilatation() * (temperature_tdt-temperature_0);
|
|
|
|
// éléments de métrique
|
|
TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH;
|
|
BaseB* giB_0;BaseB* giB_t;
|
|
bool pas_de_metrique_dispo = false; // init
|
|
|
|
if (ex_impli != NULL)
|
|
{ gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt;
|
|
giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt;
|
|
giB_t = ex_impli->giB_t; giB_0 = ex_impli->giB_0;
|
|
}
|
|
else if (ex_expli_tdt != NULL)
|
|
{gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt;
|
|
giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt;
|
|
giB_t = ex_expli_tdt->giB_t; giB_0 = ex_expli_tdt->giB_0;
|
|
}
|
|
else if (ex_umat != NULL)
|
|
{gijHH = ex_umat->gijHH_t;gijBB = ex_umat->gijBB_t;
|
|
giB = giB_t = ex_umat->giB_t; giH_0 = ex_umat->giH_0;giH = ex_umat->giH_t;
|
|
giB_0 = ex_umat->giB_0;
|
|
}
|
|
else
|
|
{ pas_de_metrique_dispo = true;
|
|
};
|
|
|
|
// -- cas de la déformation
|
|
epsBB_therm = coef_dila_global * (*gijBB);
|
|
|
|
// -- cas de l'incrément de déformation
|
|
if (atdt)
|
|
{ // cas d'un calcul avec les grandeurs à tdt disponibles
|
|
double coef_dila_t_tdt = dTP.Dilatation() * (temperature_tdt-temperature_t);
|
|
delta_epsBB_therm = coef_dila_t_tdt * (*gijBB) ;
|
|
}
|
|
else
|
|
{ // cas où les valeurs à tdt ne sont pas dispo, on considère un coef moyen sur toute la durée
|
|
// et on considère que le delta = la valeur totale de la déformation
|
|
// donc c'est le même calcul que pour la déformation totale
|
|
delta_epsBB_therm = coef_dila_global * (*gijBB) ;
|
|
};
|
|
|
|
// -- cas de la vitesse de déformation
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{ unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
if (atdt)
|
|
{ // cas d'un calcul avec les grandeurs à tdt disponibles
|
|
double coef_dila_t_tdt = dTP.Dilatation() * (temperature_tdt-temperature_t);
|
|
DepsBB_therm = coef_dila_t_tdt * (*gijBB) * unSurDeltat;
|
|
}
|
|
else
|
|
{ // cas où les valeurs à tdt ne sont pas dispo, on considère un coef moyen sur toute la durée
|
|
// recup de l'incrément du temps global
|
|
double ttemps=ParaGlob::Variables_de_temps().TempsCourant();
|
|
double unSurttemps=0;
|
|
if (Abs(ttemps) >= ConstMath::trespetit)
|
|
{unSurttemps = 1./ttemps;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/ttemps par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurttemps < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurttemps = ConstMath::tresgrand;
|
|
};
|
|
DepsBB_therm = coef_dila_global * (*gijBB) * unSurttemps;
|
|
};
|
|
|
|
|
|
// switch (abs(epsBB_meca.Dimension()))
|
|
// { case 1: epsBB_therm = coef_dila * IdBB1;
|
|
// case 2: epsBB_therm = coef_dila * IdBB2;
|
|
// case 3: epsBB_therm = coef_dila * IdBB3;
|
|
// }
|
|
if (avec_repercution_sur_def_meca)
|
|
{ epsBB_meca = epsBB_totale - epsBB_therm; // pour la déformation
|
|
delta_epsBB_meca = delta_epsBB_totale - delta_epsBB_therm; // pour la déformation
|
|
DepsBB_meca = DepsBB_totale - DepsBB_therm; // pour la vitesse de déformation
|
|
};
|
|
|
|
|
|
//// -- debug
|
|
//if (!pas_de_metrique_dispo)
|
|
// {cout << "\n debug Deformation::DeformationThermoMecanique( ";
|
|
// cout << "\n coef_dila_global= " << coef_dila_global
|
|
// << " delta_temp total= "<< (temperature_tdt-temperature_0)
|
|
// << " increment temp= "<< (temperature_tdt-temperature_t);
|
|
// // on sort les tenseurs en absolue
|
|
// Tenseur3BB tiutiu;
|
|
// epsBB_therm.BaseAbsolue(tiutiu,*giH);
|
|
// cout << "\n epsBB_therm(1,1)=" << tiutiu(1,1);
|
|
// epsBB_meca.BaseAbsolue(tiutiu,*giH);
|
|
// cout << " epsBB_meca(1,1)=" << tiutiu(1,1);
|
|
// delta_epsBB_therm.BaseAbsolue(tiutiu,*giH);
|
|
// cout << " \n delta_epsBB_therm(1,1)= "<< tiutiu(1,1);
|
|
// delta_epsBB_meca.BaseAbsolue(tiutiu,*giH);
|
|
// cout << " delta_epsBB_meca(1,1)=" << tiutiu(1,1) << endl;
|
|
// };
|
|
//// -- fin debug
|
|
|
|
};
|
|
|
|
// ========== remontee aux informations =========================
|
|
|
|
// cas sortie d'un calcul implicit
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::InfoImp Deformation::RemontImp(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ Met_abstraite::InfoImp ex = metrique->Cal_InfoImp
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul implicit
|
|
// sans le calcul des matrices de passage
|
|
const Met_abstraite::InfoImp Deformation::RemontImp()
|
|
{ Met_abstraite::InfoImp ex = metrique->Cal_InfoImp
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_tdt);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul explicit à t
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::InfoExp_t Deformation::RemontExp_t(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ Met_abstraite::InfoExp_t ex = metrique->Cal_InfoExp_t
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_t);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_t);
|
|
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul explicit à t
|
|
// sans le calcul des matrices de passage
|
|
const Met_abstraite::InfoExp_t Deformation::RemontExp_t()
|
|
{ Met_abstraite::InfoExp_t ex = metrique->Cal_InfoExp_t
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_t);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_t);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul explicit à tdt
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::InfoExp_tdt Deformation::RemontExp_tdt(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ Met_abstraite::InfoExp_tdt ex = metrique->Cal_InfoExp_tdt
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul explicit à tdt
|
|
// sans le calcul des matrices de passage
|
|
const Met_abstraite::InfoExp_tdt Deformation::RemontExp_tdt()
|
|
{ Met_abstraite::InfoExp_tdt ex = metrique->Cal_InfoExp_tdt
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_tdt);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul à 0, t et tdt
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::Info0_t_tdt Deformation::Remont0_t_tdt(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aat,Mat_pleine& Aatdt)
|
|
{ Met_abstraite::Info0_t_tdt ex = metrique->Cal_Info0_t_tdt
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giBt = *(ex.giB_t);
|
|
BaseB & giBtdt = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giHt = *(ex.giH_t);
|
|
BaseH & giHtdt = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giBt,giBtdt,giH0,giHt,giHtdt,Aa0,Aat,Aatdt);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul à 0, t et tdt
|
|
// sans le calcul des matrices de passage
|
|
const Met_abstraite::Info0_t_tdt Deformation::Remont0_t_tdt()
|
|
{ Met_abstraite::Info0_t_tdt ex = metrique->Cal_Info0_t_tdt
|
|
( *tabnoeud,(*tabDphi)(numInteg),(*tabPhi)(numInteg),nbNoeud);
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giBt = *(ex.giB_t);
|
|
BaseB & giBtdt = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giHt = *(ex.giH_t);
|
|
BaseH & giHtdt = *(ex.giH_tdt);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul implicit
|
|
// mêmes fct que précédemment mais dans le cas où les éléments de métrique viennent justes
|
|
// d'être calculés, donc les seules grandeurs qui sont réellement calculées sont les matrices
|
|
// Aa0 et Aafin
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::InfoImp Deformation::RemontImpSansCalMet(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ Met_abstraite::InfoImp ex = metrique->Recup_InfoImp();
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul explicit à t
|
|
// mêmes fct que précédemment mais dans le cas où les éléments de métrique viennent justes
|
|
// d'être calculés, donc les seules grandeurs qui sont réellement calculées sont les matrices
|
|
// Aa0 et Aafin
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::InfoExp_t Deformation::RemontExp_tSansCalMet(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ Met_abstraite::InfoExp_t ex = metrique->Recup_InfoExp_t();
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_t);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_t);
|
|
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul explicit à tdt
|
|
// mêmes fct que précédemment mais dans le cas où les éléments de métrique viennent justes
|
|
// d'être calculés, donc les seules grandeurs qui sont réellement calculées sont les matrices
|
|
// Aa0 et Aafin
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::InfoExp_tdt Deformation::RemontExp_tdtSansCalMet(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ Met_abstraite::InfoExp_tdt ex = metrique->Recup_InfoExp_tdt();
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giB = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giH = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giB,giH0,giH,Aa0,Aafin);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul à 0, t et tdt
|
|
// mêmes fct que précédemment mais dans le cas où les éléments de métrique viennent justes
|
|
// d'être calculés, donc les seules grandeurs qui sont réellement calculées sont les matrices
|
|
// Aa0 et Aafin
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::Info0_t_tdt Deformation::Remont0_t_tdtSansCalMet(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aat,Mat_pleine& Aatdt)
|
|
{ Met_abstraite::Info0_t_tdt ex = metrique->Recup_Info0_t_tdt();
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giBt = *(ex.giB_t);
|
|
BaseB & giBtdt = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giHt = *(ex.giH_t);
|
|
BaseH & giHtdt = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giBt,giBtdt,giH0,giHt,giHtdt,Aa0,Aat,Aatdt);
|
|
return ex;
|
|
};
|
|
|
|
// cas sortie d'un calcul à 0, t et tdt avec tous les métriques associées
|
|
// mêmes fct que précédemment mais dans le cas où les éléments de métrique viennent justes
|
|
// d'être calculés, donc les seules grandeurs qui sont réellement calculées sont les matrices
|
|
// Aa0 et Aafin
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
const Met_abstraite::Info_et_metrique_0_t_tdt Deformation::Remont_et_metrique_0_t_tdtSansCalMet(bool absolue,Mat_pleine& Aa0,Mat_pleine& Aat,Mat_pleine& Aatdt)
|
|
{ Met_abstraite::Info_et_metrique_0_t_tdt ex = metrique->Recup_Info_et_metrique_0_t_tdt();
|
|
// determination des matrices de transformation de base
|
|
BaseB & giB0 = *(ex.giB_0);
|
|
BaseB & giBt = *(ex.giB_t);
|
|
BaseB & giBtdt = *(ex.giB_tdt);
|
|
BaseH & giH0 = *(ex.giH_0);
|
|
BaseH & giHt = *(ex.giH_t);
|
|
BaseH & giHtdt = *(ex.giH_tdt);
|
|
BasePassage(absolue,giB0,giBt,giBtdt,giH0,giHt,giHtdt,Aa0,Aat,Aatdt);
|
|
return ex;
|
|
};
|
|
|
|
// cas où l'on veut les matrices de passages à 0 et final
|
|
// voir Deformation.h pour plus d'info
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
// (voir Deformation.h pour un descriptif plus complet)
|
|
// si absolue = true: on fait un passage vers la base absolue
|
|
// si = false : on fait un passage vers une base ad hoc
|
|
void Deformation::BasePassage
|
|
(bool absolue,const BaseB & giB0,const BaseB & giB,const BaseH & giH0,const BaseH & giH,
|
|
Mat_pleine& Aa0,Mat_pleine& Aafin)
|
|
{ // determination des matrices de transformation de base
|
|
int dim = giB0.NbVecteur();
|
|
|
|
// dim_espace X nb_vecteur X taille_matrice
|
|
|
|
// cas d'une sortie en globale
|
|
// 1X1X1 ou 2X2X2 ou 3X3X3
|
|
if ((dim == ParaGlob::Dimension()) && (Aa0.Nb_ligne() == dim))
|
|
{ // --- ici il n'y a pas a priori de repère ad hoc, donc le calcul ne dépend pas de "absolue"
|
|
|
|
// petite vérif
|
|
#ifdef MISE_AU_POINT
|
|
if ((Aa0.Nb_ligne() != dim) && (Aa0.Nb_colonne() !=dim)
|
|
&& (Aafin.Nb_ligne() != dim) && (Aafin.Nb_colonne() !=dim))
|
|
{ cout << "\nErreur1 : les matrices Aa0 et /ou Aafin ne sont pas correctement"
|
|
<< " dimensionnees ! \n";
|
|
cout << "\n Aa0= "; Aa0.Affiche();
|
|
cout << "\n Aafin= "; Aafin.Affiche();
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int i=1;i<=dim;i++)
|
|
for (int j=1;j<=dim;j++)
|
|
{Aa0 (i,j) = giH0(i)(j);
|
|
Aafin(i,j) =giH(i)(j);
|
|
};
|
|
}
|
|
else if ((dim==2)&&(ParaGlob::Dimension()==3))
|
|
// cela comprend les cas: 3X2X2 -> un repère local ad hoc en 2D
|
|
// 3X2X3 -> un repère local ad hoc en 3D
|
|
// ==== cas d'élément de membrane en 3D ====
|
|
{// 1) si absolue = false: l'objectif est de determiner un repere pertinant dans le cas de membrane .
|
|
// choix : un repere qui appartiend a la facette, obtenu apres projection
|
|
// du repere global
|
|
//------ cas de la config initiale, on regarde si la projection de I1 n'est pas trop petite
|
|
// def de la normale a la facette
|
|
// 2) si absolue = true; on complète le repère des gi avec la normale et la matrice de passage
|
|
// correspond au passage des gi complètés -> I_a
|
|
|
|
Coordonnee N = (Util::ProdVec_coorBN(giB0(1),giB0(2))).Normer();
|
|
Coordonnee ve(0.,-N(3),N(2)); // = ProdVec(N,I1)
|
|
double norve = ve.Norme();
|
|
Tableau <Coordonnee> vi(2); // les vecteurs plans orthonormes de la facette
|
|
if (!absolue)
|
|
{Coordonnee ve(0.,-N(3),N(2)); // = ProdVec(N,I1)
|
|
double norve = ve.Norme();
|
|
Tableau <Coordonnee> vi(2); // les vecteurs plans orthonormes de la facette
|
|
if (norve >= 0.01)
|
|
{ vi(2) = ve.Normer();
|
|
vi(1) = Util::ProdVec_coor(vi(2),N);
|
|
}
|
|
else
|
|
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
|
|
vi(1) = ve.Normer();
|
|
vi(2) = Util::ProdVec_coor(N,vi(1));
|
|
};
|
|
for (int al=1 ;al<=2;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1) * vi(al) ;
|
|
Aa0(2,al) = giH0.Coordo(2) * vi(al) ;
|
|
};
|
|
// si on veut une matrice 3x3, il faut complèter la matrice avec un troisième vecteur
|
|
// comme les vi sont normaux à N, le troisième vecteur est le vecteur N
|
|
// qui ne peut pas s'exprimer dans le repère local des gi donc on considère qu'il s'exprime vis-à-vis
|
|
// du troisième vecteur qui n'est autre que N
|
|
if (Aa0.Nb_ligne()==3)
|
|
// ici on a Ip_3 == N et également g_3=g^3=N
|
|
{for (int al=1 ;al<=2;al++)
|
|
{Aa0(3,al) = N * vi(al) ;
|
|
Aa0(al,3) = giH0.Coordo(al) * N ;
|
|
};
|
|
Aa0(3,3)= 1.; // == N . N
|
|
};
|
|
}
|
|
else // cas on l'on veut le passage des gi et N vers I_a
|
|
{// g^i = Aa^i_{.a} * I^a -> donc Aa^i_{.a} c'est directement les coordonnées de g^i
|
|
for (int al=1 ;al < 4;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa0(2,al) = giH0.Coordo(2)(al) ;
|
|
};
|
|
// pour le troisième: N = Aa^3_{.a} * I^a -> c'est aussi directement les coordonnées de N
|
|
Aa0.RemplaceLigneSpe(3,N.Vect());
|
|
};
|
|
|
|
|
|
//------ cas de la config finale,
|
|
N = (Util::ProdVec_coorBN(giB(1),giB(2))).Normer();
|
|
if (!absolue)
|
|
{ve.Change_Coordonnee(3,0.,-N(3),N(2)); // = ProdVec(N,I1)
|
|
norve = ve.Norme();
|
|
Tableau <Coordonnee> Aa(2);
|
|
if (norve >= ConstMath::petit)
|
|
{ Aa(2) = ve.Normer();
|
|
Aa(1) = Util::ProdVec_coor(Aa(2),N);
|
|
}
|
|
else
|
|
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
|
|
Aa(1) = ve.Normer();
|
|
Aa(2) = Util::ProdVec_coor(N,Aa(1));
|
|
};
|
|
|
|
for (int be=1;be<=2;be++)
|
|
{ Aafin(1,be) = giH.Coordo(1) * Aa(be);
|
|
Aafin(2,be) = giH.Coordo(2) * Aa(be);
|
|
};
|
|
// si on veut une matrice 3x3, il faut complèter la matrice avec un troisième vecteur
|
|
// comme les vi sont normaux à N, le troisième vecteur est le vecteur N
|
|
// qui ne peut pas s'exprimer dans le repère local des gi donc on considère qu'il s'exprime vis-à-vis
|
|
// du troisième vecteur qui n'est autre que N
|
|
if (Aafin.Nb_ligne()==3)
|
|
// ici on a Ip_3 == N et également g_3=g^3=N
|
|
{for (int al=1 ;al<=2;al++)
|
|
{Aafin(3,al) = N * vi(al) ;
|
|
Aafin(al,3) = giH.Coordo(al) * N ;
|
|
};
|
|
Aafin(3,3)= 1.; // == N . N
|
|
};
|
|
}
|
|
else // cas on l'on veut le passage des gi et N vers I_a
|
|
{// g^i = Aa^i_{.a} * I^a -> donc Aa^i_{.a} c'est directement les coordonnées de g^i
|
|
for (int al=1 ;al < 4;al++)
|
|
{Aafin(1,al) = giH.Coordo(1)(al) ;
|
|
Aafin(2,al) = giH.Coordo(2)(al) ;
|
|
};
|
|
// pour le troisième: N = Aa^3_{.a} * I^a -> c'est aussi directement les coordonnées de N
|
|
Aafin.RemplaceLigneSpe(3,N.Vect());
|
|
};
|
|
}
|
|
else if (dim==1)
|
|
// ==== cas des éléments 1D en espace 1D, 2D ou 3D ====
|
|
{ if (!absolue)
|
|
{switch (ParaGlob::Dimension())
|
|
{ case 1: // normalement on ne devrait pas passer ici !! car ce cas est déjà traité au début
|
|
{cout << "\nErreur : bug!! on ne devrait pas passer ici \n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
case 3: // on est dans un espace 3D, ici on ne fait que le particulier à 3D et le reste est
|
|
// fait dans le case 2: qui suit
|
|
// 3X1X1 ou 3X1X2 ou 3X1X3, 2X1X1 , 2X1X2
|
|
if (Aa0.Nb_ligne() == 3 )
|
|
// 3X1X3
|
|
{ // tous les tenseurs vont utiliser que le premier vecteur, donc le troisième peut être quelconque
|
|
// et l'expression de giH(3) n'existant pas on dit que par défaut (et par simplicité)
|
|
// giH(3) = 1 * IP(3) d'où
|
|
Aa0 (3,3) = Aafin(3,3) = 1.;
|
|
// et on initialise les autres composantes au cas où
|
|
Aa0 (1,3) = Aa0 (3,1) = Aafin(1,3) = Aafin(3,1) = 0.;
|
|
Aa0 (2,3) = Aa0 (3,2) = Aafin(2,3) = Aafin(3,2) = 0.;
|
|
// pas de break, on continue sur le 2D
|
|
};
|
|
case 2: // on est dans un espace 2D, 3X1X1 ou 3X1X2 ou 3X1X3
|
|
{ // le premier vecteur intéressant c'est le vecteur collinéaire avec gH(1), mais normé
|
|
double norme0 = giH0(1).Norme();
|
|
double norme = giH(1).Norme();
|
|
Aa0 (1,1) = norme0;
|
|
Aafin(1,1) = norme;
|
|
// fin du cas 3X1X1 et 2X1X1
|
|
if (Aa0.Nb_ligne() > 1)
|
|
// 3X1X2 , 3X1X3 , 2X1X2
|
|
// tous les tenseurs vont utiliser que le premier vecteur, donc le second peut être quelconque
|
|
// et l'expression de giH(2) n'existant pas on dit que par défaut (et par simplicité)
|
|
// giH(2) = 1 * IP(2) d'où
|
|
{Aa0 (2,2) = Aafin(2,2) = 1.;
|
|
// et on initialise les autres composantes au cas où
|
|
Aa0 (1,2) = Aa0 (2,1) = Aafin(1,2) = Aafin(2,1) = 0.;
|
|
};
|
|
break;
|
|
}
|
|
default: cout << "\n erreur grossiere de dimension !! \n Deformation::BasePassage(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else // sinon on veut arriver en Ia
|
|
{ // On ne connait que le premier vecteur et normalement c'est seulement celui-là qui est
|
|
// utilisé, par contre il faut complèter la matrice pour qu'elle ne soit pas singulière
|
|
Aa0.Initialise(0.);Aafin.Initialise(0.);
|
|
switch (ParaGlob::Dimension())
|
|
{ case 1: // normalement on ne devrait pas passer ici !! car ce cas est déjà traité au début
|
|
{cout << "\nErreur : bug!! on ne devrait pas passer ici \n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
case 3: // on est dans un espace 3D,
|
|
{CoordonneeH V02(3); CoordonneeH V03(3);
|
|
// on calcule deux vecteurs arbitraires normaux, avec un résultat orthonormé V2,V3
|
|
// mais pas giH0(1) !
|
|
MathUtil2::Def_vecteurs_plan(giH0(1),V02,V03);
|
|
// idem pour le final
|
|
CoordonneeH V2(3); CoordonneeH V3(3);
|
|
MathUtil2::Def_vecteurs_plan(giH(1),V2,V3);
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aafin(1,al) = giH.Coordo(1)(al);
|
|
Aa0(2,al) = V02(al) ;
|
|
Aafin(2,al) = V2(al) ;
|
|
Aa0(3,al) = V03(al) ;
|
|
Aafin(3,al) = V3(al) ;
|
|
};
|
|
break;
|
|
}
|
|
case 2: // on est dans un espace 2D,
|
|
{CoordonneeH V02(2);
|
|
// on calcule le vecteur normal,
|
|
// mais pas giH0(1) !
|
|
MathUtil2::Def_vecteurs_plan(giH0(1),V02);
|
|
// idem pour le final
|
|
CoordonneeH V2(2);
|
|
MathUtil2::Def_vecteurs_plan(giH(1),V2);
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aafin(1,al) = giH.Coordo(1)(al);
|
|
Aa0(2,al) = V02(al) ;
|
|
Aafin(2,al) = V2(al) ;
|
|
};
|
|
break;
|
|
}
|
|
default: cout << "\n erreur2 grossiere de dimension !! \n Deformation::BasePassage(.. ";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
}
|
|
else
|
|
{// pour les autres cas on ramene la matrice identite
|
|
#ifdef MISE_AU_POINT
|
|
if ((Aa0.Nb_ligne() != dim) && (Aa0.Nb_colonne() !=dim) &&
|
|
(Aafin.Nb_ligne() != dim) && (Aafin.Nb_colonne() !=dim))
|
|
{ cout << "\nErreur : les matrices Aa0 et /ou Aafin ne sont pas correctement"
|
|
<< " dimensionnees !\n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int i=1;i<=dim;i++)
|
|
{ for (int j=1;j<=dim;j++)
|
|
{Aa0 (i,j) = 0.;Aafin(i,j) = 0.;}
|
|
Aa0 (i,i) = 1.; Aafin(i,i) = 1.;
|
|
};
|
|
};
|
|
// retour
|
|
return;
|
|
};
|
|
|
|
// cas où l'on veut les matrices de passages à 0 t et tdt
|
|
// voir Deformation.h pour plus d'info
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
// (voir Deformation.h pour un descriptif plus complet)
|
|
void Deformation::BasePassage(bool absolue,const BaseB & giB0,const BaseB & giB_t,const BaseB & giB_tdt
|
|
,const BaseH & giH0,const BaseH & giH_t,const BaseH & giH_tdt
|
|
,Mat_pleine& Aa0,Mat_pleine& Aa_t,Mat_pleine& Aa_tdt)
|
|
{
|
|
// determination des matrices de transformation de base
|
|
int dim = giB0.NbVecteur();
|
|
|
|
// dim_espace X nb_vecteur X taille_matrice
|
|
|
|
// cas d'une sortie en globale
|
|
// 1X1X1 ou 2X2X2 ou 3X3X3
|
|
if ((dim == ParaGlob::Dimension()) && (Aa0.Nb_ligne() == dim))
|
|
{ // --- ici il n'y a pas a priori de repère ad hoc, donc le calcul ne dépend pas de "absolue"
|
|
|
|
// petite vérif
|
|
#ifdef MISE_AU_POINT
|
|
if ((Aa0.Nb_ligne() != dim) && (Aa0.Nb_colonne() !=dim)
|
|
&& (Aa_t.Nb_ligne() != dim) && (Aa_t.Nb_colonne() !=dim)
|
|
&& (Aa_tdt.Nb_ligne() != dim) && (Aa_tdt.Nb_colonne() !=dim))
|
|
{ cout << "\nErreur1 : les matrices Aa0 et /ou Aa_t et / ou Aa_tdt ne sont pas correctement"
|
|
<< " dimensionnees ! \n";
|
|
cout << "\n Aa0= "; Aa0.Affiche();
|
|
cout << "\n Aa_t= "; Aa_t.Affiche();
|
|
cout << "\n Aa_tdt= "; Aa_tdt.Affiche();
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int i=1;i<=dim;i++)
|
|
for (int j=1;j<=dim;j++)
|
|
{Aa0 (i,j) = giH0(i)(j);
|
|
Aa_t(i,j) =giH_t(i)(j);
|
|
Aa_tdt(i,j) =giH_tdt(i)(j);
|
|
};
|
|
}
|
|
else if ((dim==2)&&(ParaGlob::Dimension()==3))
|
|
// cela comprend les cas: 3X2X2 -> un repère local ad hoc en 2D
|
|
// 3X2X3 -> un repère local ad hoc en 3D
|
|
// ==== cas d'élément de membrane en 3D ====
|
|
{// 1) si absolue = false: l'objectif est de determiner un repere pertinant dans le cas de membrane .
|
|
// choix : un repere qui appartiend a la facette, obtenu apres projection
|
|
// du repere global
|
|
//------ cas de la config initiale, on regarde si la projection de I1 n'est pas trop petite
|
|
// def de la normale a la facette
|
|
// 2) si absolue = true; on complète le repère des gi avec la normale et la matrice de passage
|
|
// correspond au passage des gi complètés -> I_a
|
|
|
|
Coordonnee N = (Util::ProdVec_coorBN(giB0(1),giB0(2))).Normer();
|
|
Coordonnee ve(0.,-N(3),N(2)); // = ProdVec(N,I1), donc normal à N
|
|
double norve = ve.Norme();
|
|
Tableau <Coordonnee> vi(2); // les vecteurs plans orthonormes de la facette
|
|
if (!absolue)
|
|
{if (norve >= 0.01)
|
|
{ vi(2) = ve.Normer();
|
|
vi(1) = Util::ProdVec_coor(vi(2),N);// donc normal à N
|
|
}
|
|
else
|
|
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N), donc normal à N
|
|
vi(1) = ve.Normer();
|
|
vi(2) = Util::ProdVec_coor(N,vi(1));// donc normal à N
|
|
};
|
|
// maintenant expression dans le repère local
|
|
// Aa^i_{.a} = g^i . Ip_a avec les Ip_a == vi(a)
|
|
for (int al=1 ;al<=2;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1) * vi(al) ;
|
|
Aa0(2,al) = giH0.Coordo(2) * vi(al) ;
|
|
};
|
|
// si on veut une matrice 3x3, il faut complèter la matrice avec un troisième vecteur
|
|
// comme les vi sont normaux à N, le troisième vecteur est le vecteur N
|
|
// qui ne peut pas s'exprimer dans le repère local des gi donc on considère qu'il s'exprime vis-à-vis
|
|
// du troisième vecteur qui n'est autre que N
|
|
if (Aa0.Nb_ligne()==3)
|
|
// ici on a Ip_3 == N et également g_3=g^3=N
|
|
{for (int al=1 ;al<=2;al++)
|
|
{Aa0(3,al) = N * vi(al) ;
|
|
Aa0(al,3) = giH0.Coordo(al) * N ;
|
|
};
|
|
Aa0(3,3)= 1.; // == N . N
|
|
};
|
|
}
|
|
else // cas on l'on veut le passage des gi et N vers I_a
|
|
{// g^i = Aa^i_{.a} * I^a -> donc Aa^i_{.a} c'est directement les coordonnées de g^i
|
|
for (int al=1 ;al < 4;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa0(2,al) = giH0.Coordo(2)(al) ;
|
|
};
|
|
// pour le troisième: N = Aa^3_{.a} * I^a -> c'est aussi directement les coordonnées de N
|
|
Aa0.RemplaceLigneSpe(3,N.Vect());
|
|
};
|
|
|
|
//------ cas de la config à t,
|
|
N = (Util::ProdVec_coorBN(giB_t(1),giB_t(2))).Normer();
|
|
Tableau <Coordonnee> Aa(2);
|
|
if (!absolue)
|
|
{ve.Change_Coordonnee(3,0.,-N(3),N(2)); // = ProdVec(N,I1)
|
|
norve = ve.Norme();
|
|
if (norve >= 0.01)
|
|
{ Aa(2) = ve.Normer();
|
|
Aa(1) = Util::ProdVec_coor(Aa(2),N);
|
|
}
|
|
else
|
|
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
|
|
Aa(1) = ve.Normer();
|
|
Aa(2) = Util::ProdVec_coor(N,Aa(1));
|
|
};
|
|
|
|
for (int be=1;be<=2;be++)
|
|
{ Aa_t(1,be) = giH_t.Coordo(1) * Aa(be);
|
|
Aa_t(2,be) = giH_t.Coordo(2) * Aa(be);
|
|
};
|
|
// si on veut une matrice 3x3, il faut complèter la matrice avec un troisième vecteur
|
|
// comme les vi sont normaux à N, le troisième vecteur est le vecteur N
|
|
// qui ne peut pas s'exprimer dans le repère local des gi donc on considère qu'il s'exprime vis-à-vis
|
|
// du troisième vecteur qui n'est autre que N
|
|
if (Aa_t.Nb_ligne()==3)
|
|
// ici on a Ip_3 == N et également g_3=g^3=N
|
|
{for (int al=1 ;al<=2;al++)
|
|
{Aa_t(3,al) = N * vi(al) ;
|
|
Aa_t(al,3) = giH_t.Coordo(al) * N ;
|
|
};
|
|
Aa_t(3,3)= 1.; // == N . N
|
|
};
|
|
}
|
|
else // cas on l'on veut le passage des gi et N vers I_a
|
|
{// g^i = Aa^i_{.a} * I^a -> donc Aa^i_{.a} c'est directement les coordonnées de g^i
|
|
for (int al=1 ;al < 4;al++)
|
|
{Aa_t(1,al) = giH_t.Coordo(1)(al) ;
|
|
Aa_t(2,al) = giH_t.Coordo(2)(al) ;
|
|
};
|
|
// pour le troisième: N = Aa^3_{.a} * I^a -> c'est aussi directement les coordonnées de N
|
|
Aa_t.RemplaceLigneSpe(3,N.Vect());
|
|
};
|
|
|
|
//// ---- vérification que le déterminant n'est pas nul
|
|
//{double det = Aa_t(1,1) * Aa_t(2,2) - Aa_t(1,2) * Aa_t(2,1);
|
|
// if (Abs(det) < 0.00001)
|
|
// cout << "\n erreur det0 null " << det << " dans Deformation::BasePassage( ";
|
|
//};
|
|
//// ---- fin vérification
|
|
//------ cas de la config finale à tdt,
|
|
N = (Util::ProdVec_coorBN(giB_tdt(1),giB_tdt(2))).Normer();
|
|
if (!absolue)
|
|
{ve.Change_Coordonnee(3,0.,-N(3),N(2)); // = ProdVec(N,I1)
|
|
norve = ve.Norme();
|
|
if (norve >= 0.01)
|
|
{ Aa(2) = ve.Normer();
|
|
Aa(1) = Util::ProdVec_coor(Aa(2),N);
|
|
}
|
|
else
|
|
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N)
|
|
Aa(1) = ve.Normer();
|
|
Aa(2) = Util::ProdVec_coor(N,Aa(1));
|
|
};
|
|
|
|
for (int be=1;be<=2;be++)
|
|
{ Aa_tdt(1,be) = giH_tdt.Coordo(1) * Aa(be);
|
|
Aa_tdt(2,be) = giH_tdt.Coordo(2) * Aa(be);
|
|
};
|
|
// si on veut une matrice 3x3, il faut complèter la matrice avec un troisième vecteur
|
|
// comme les vi sont normaux à N, le troisième vecteur est le vecteur N
|
|
// qui ne peut pas s'exprimer dans le repère local des gi donc on considère qu'il s'exprime vis-à-vis
|
|
// du troisième vecteur qui n'est autre que N
|
|
if (Aa_tdt.Nb_ligne()==3)
|
|
// ici on a Ip_3 == N et également g_3=g^3=N
|
|
{for (int al=1 ;al<=2;al++)
|
|
{Aa_tdt(3,al) = N * vi(al) ;
|
|
Aa_tdt(al,3) = giH_tdt.Coordo(al) * N ;
|
|
};
|
|
Aa_tdt(3,3)= 1.; // == N . N
|
|
};
|
|
}
|
|
else // cas on l'on veut le passage des gi et N vers I_a
|
|
{// g^i = Aa^i_{.a} * I^a -> donc Aa^i_{.a} c'est directement les coordonnées de g^i
|
|
for (int al=1 ;al < 4;al++)
|
|
{Aa_tdt(1,al) = giH_tdt.Coordo(1)(al) ;
|
|
Aa_tdt(2,al) = giH_tdt.Coordo(2)(al) ;
|
|
};
|
|
// pour le troisième: N = Aa^3_{.a} * I^a -> c'est aussi directement les coordonnées de N
|
|
Aa_tdt.RemplaceLigneSpe(3,N.Vect());
|
|
};
|
|
|
|
//// ---- vérification que le déterminant n'est pas nul
|
|
//{double det = Aa_tdt(1,1) * Aa_tdt(2,2) - Aa_tdt(1,2) * Aa_tdt(2,1);
|
|
// if (Abs(det) < 0.00001)
|
|
// cout << "\n erreur det0 null " << det << " dans Deformation::BasePassage( ";
|
|
//};
|
|
//// ---- fin vérification
|
|
}
|
|
else if (dim==1)
|
|
// ==== cas des éléments 1D en espace 1D, 2D ou 3D ====
|
|
{ if (!absolue)
|
|
{switch (ParaGlob::Dimension())
|
|
{ case 1: // normalement on ne devrait pas passer ici !! car ce cas est déjà traité au début
|
|
{cout << "\nErreur : bug!! on ne devrait pas passer ici \n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
case 3: // on est dans un espace 3D, ici on ne fait que le particulier à 3D et le reste est
|
|
// fait dans le case 2: qui suit
|
|
// 3X1X1 ou 3X1X2 ou 3X1X3, 2X1X1 , 2X1X2
|
|
if (Aa0.Nb_ligne() == 3 )
|
|
// 3X1X3
|
|
{ // tous les tenseurs vont utiliser que le premier vecteur, donc le troisième peut être quelconque
|
|
// et l'expression de giH(3) n'existant pas on dit que par défaut (et par simplicité)
|
|
// giH(3) = 1 * IP(3) d'où
|
|
Aa0 (3,3) = Aa_t(3,3)= Aa_tdt(3,3) = 1.;
|
|
// et on initialise les autres composantes au cas où
|
|
Aa0 (1,3) = Aa0 (3,1) = Aa_t(1,3) = Aa_t(3,1) = Aa_tdt(1,3) = Aa_tdt(3,1) = 0.;
|
|
Aa0 (2,3) = Aa0 (3,2) = Aa_t(2,3) = Aa_t(3,2) = Aa_tdt(2,3) = Aa_tdt(3,2) = 0.;
|
|
// pas de break, on continue sur le 2D
|
|
}
|
|
case 2: // on est dans un espace 2D, 3X1X1 ou 3X1X2 ou 3X1X3
|
|
{ // le premier vecteur intéressant c'est le vecteur collinéaire avec gH(1), mais normé
|
|
double norme0 = giH0(1).Norme();
|
|
double normet = giH_t(1).Norme();
|
|
double normetdt = giH_tdt(1).Norme();
|
|
Aa0 (1,1) = norme0;
|
|
Aa_t(1,1) = normet;
|
|
Aa_tdt(1,1) = normetdt;
|
|
// fin du cas 3X1X1 et 2X1X1
|
|
if (Aa0.Nb_ligne() > 1)
|
|
// 3X1X2 , 3X1X3 , 2X1X2
|
|
// tous les tenseurs vont utiliser que le premier vecteur, donc le second peut être quelconque
|
|
// et l'expression de giH(2) n'existant pas on dit que par défaut (et par simplicité)
|
|
// giH(2) = 1 * IP(2) d'où
|
|
{Aa0 (2,2) = Aa_t(2,2) = Aa_tdt(2,2)= 1.;
|
|
// et on initialise les autres composantes au cas où
|
|
Aa0 (1,2) = Aa0 (2,1) = Aa_t(1,2) = Aa_t(2,1)= Aa_tdt(1,2) = Aa_tdt(2,1) = 0.;
|
|
};
|
|
break;
|
|
}
|
|
default: cout << "\n erreur grossiere de dimension !! \n Deformation::BasePassage(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else // sinon on veut arriver en Ia
|
|
{ // On ne connait que le premier vecteur et normalement c'est seulement celui-là qui est
|
|
// utilisé, par contre il faut complèter la matrice pour qu'elle ne soit pas singulière
|
|
Aa0.Initialise(0.);Aa_t.Initialise(0.);Aa_tdt.Initialise(0.);
|
|
switch (ParaGlob::Dimension())
|
|
{ case 1: // normalement on ne devrait pas passer ici !! car ce cas est déjà traité au début
|
|
{cout << "\nErreur : bug!! on ne devrait pas passer ici \n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
case 3: // on est dans un espace 3D,
|
|
{CoordonneeH V02(3); CoordonneeH V03(3);
|
|
// on calcule deux vecteurs arbitraires normaux, avec un résultat orthonormé V2,V3
|
|
// mais pas giH0(1) !
|
|
MathUtil2::Def_vecteurs_plan(giH0(1),V02,V03);
|
|
// idem pour t et tdt
|
|
CoordonneeH V2_t(3); CoordonneeH V3_t(3);
|
|
MathUtil2::Def_vecteurs_plan(giH_t(1),V2_t,V3_t);
|
|
CoordonneeH V2_tdt(3); CoordonneeH V3_tdt(3);
|
|
MathUtil2::Def_vecteurs_plan(giH_tdt(1),V2_tdt,V3_tdt);
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa_t(1,al) = giH_t.Coordo(1)(al) ;
|
|
Aa_tdt(1,al) = giH_tdt.Coordo(1)(al) ;
|
|
Aa0(2,al) = V02(al) ;
|
|
Aa_t(2,al) = V2_t(al) ;
|
|
Aa_tdt(2,al) = V2_tdt(al) ;
|
|
Aa0(3,al) = V03(al) ;
|
|
Aa_t(3,al) = V3_t(al) ;
|
|
Aa_tdt(3,al) = V3_tdt(al) ;
|
|
};
|
|
break;
|
|
}
|
|
case 2: // on est dans un espace 2D,
|
|
{CoordonneeH V02(2);
|
|
// on calcule le vecteur normal,
|
|
// mais pas giH0(1) !
|
|
MathUtil2::Def_vecteurs_plan(giH0(1),V02);
|
|
// idem pour t et tdt
|
|
CoordonneeH V2_t(3);
|
|
MathUtil2::Def_vecteurs_plan(giH_t(1),V2_t);
|
|
CoordonneeH V2_tdt(3);
|
|
MathUtil2::Def_vecteurs_plan(giH_tdt(1),V2_tdt);
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa_t(1,al) = giH_t.Coordo(1)(al) ;
|
|
Aa_tdt(1,al) = giH_tdt.Coordo(1)(al) ;
|
|
Aa0(2,al) = V02(al) ;
|
|
Aa_t(2,al) = V2_t(al) ;
|
|
Aa_tdt(2,al) = V2_tdt(al) ;
|
|
};
|
|
break;
|
|
}
|
|
default: cout << "\n erreur2 grossiere de dimension !! \n Deformation::BasePassage(.. ";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
}
|
|
else
|
|
{// pour les autres cas on ramene la matrice identite
|
|
#ifdef MISE_AU_POINT
|
|
if ((Aa0.Nb_ligne() != dim) && (Aa0.Nb_colonne() !=dim) &&
|
|
(Aa_t.Nb_ligne() != dim) && (Aa_t.Nb_colonne() !=dim) &&
|
|
(Aa_tdt.Nb_ligne() != dim) && (Aa_tdt.Nb_colonne() !=dim))
|
|
{ cout << "\nErreur : les matrices Aa0 et /ou Aa_t et /ou Aa_tdt ne sont pas correctement"
|
|
<< " dimensionnees !\n";
|
|
cout << "void Deformation::RemontImp( \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int i=1;i<=dim;i++)
|
|
{ for (int j=1;j<=dim;j++)
|
|
{Aa0 (i,j) = 0.;Aa_t(i,j) = 0.;Aa_tdt(i,j) = 0.;}
|
|
Aa0 (i,i) = 1.; Aa_t(i,i) = 1.;Aa_tdt(i,i) = 1.;
|
|
};
|
|
};
|
|
// retour
|
|
return;
|
|
};
|
|
|
|
// cas où l'on veut la matrice de passage à 0 uniquement
|
|
// Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a
|
|
// tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer
|
|
// (voir Deformation.h pour un descriptif plus complet)
|
|
void Deformation::BasePassage(bool absolue,const BaseB & giB0,const BaseH & giH0,Mat_pleine& Aa0)
|
|
{
|
|
// determination des matrices de transformation de base
|
|
int dim = giB0.NbVecteur();
|
|
|
|
// dim_espace X nb_vecteur X taille_matrice
|
|
|
|
// cas d'une sortie en globale
|
|
// 1X1X1 ou 2X2X2 ou 3X3X3
|
|
if ((dim == ParaGlob::Dimension()) && (Aa0.Nb_ligne() == dim))
|
|
{ // petite vérif
|
|
#ifdef MISE_AU_POINT
|
|
if ((Aa0.Nb_ligne() != dim) && (Aa0.Nb_colonne() !=dim))
|
|
{ cout << "\nErreur1 : la matrice Aa0 n'est pas correctement"
|
|
<< " dimensionnee ! \n";
|
|
cout << "\n Aa0= "; Aa0.Affiche();
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int i=1;i<=dim;i++)
|
|
for (int j=1;j<=dim;j++)
|
|
{Aa0 (i,j) = giH0(i)(j);
|
|
};
|
|
}
|
|
else if ((dim==2)&&(ParaGlob::Dimension()==3))
|
|
// cela comprend les cas: 3X2X2 -> un repère local ad hoc en 2D
|
|
// 3X2X3 -> un repère local ad hoc en 3D
|
|
// ==== cas d'élément de membrane en 3D ====
|
|
{// 1) si absolue = false: l'objectif est de determiner un repere pertinant dans le cas de membrane .
|
|
// choix : un repere qui appartiend a la facette, obtenu apres projection
|
|
// du repere global
|
|
//------ cas de la config initiale, on regarde si la projection de I1 n'est pas trop petite
|
|
// def de la normale a la facette
|
|
// 2) si absolue = true; on complète le repère des gi avec la normale et la matrice de passage
|
|
// correspond au passage des gi complètés -> I_a
|
|
|
|
Coordonnee N ( (Util::ProdVec_coorBN(giB0(1),giB0(2))).Normer() );
|
|
Coordonnee ve(0.,-N(3),N(2)); // = ProdVec(N,I1), donc normal à N
|
|
double norve = ve.Norme();
|
|
Tableau <Coordonnee> vi(2); // les vecteurs plans orthonormes de la facette
|
|
if (!absolue)
|
|
{if (norve >= 0.01)
|
|
{ vi(2) = ve.Normer();
|
|
vi(1) = Util::ProdVec_coor(vi(2),N);// donc normal à N
|
|
}
|
|
else
|
|
{ ve.Change_Coordonnee(3,N(3),0.,-N(1)); // = ProdVec(I2,N), donc normal à N
|
|
vi(1) = ve.Normer();
|
|
vi(2) = Util::ProdVec_coor(N,vi(1));// donc normal à N
|
|
};
|
|
// maintenant expression dans le repère local
|
|
for (int al=1 ;al<=2;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1) * vi(al) ;
|
|
Aa0(2,al) = giH0.Coordo(2) * vi(al) ;
|
|
};
|
|
// si on veut une matrice 3x3, il faut complèter la matrice avec un troisième vecteur
|
|
// comme les vi sont normaux à N, le troisième vecteur est le vecteur N
|
|
// qui ne peut pas s'exprimer dans le repère local des gi donc on considère qu'il s'exprime vis-à-vis
|
|
// du troisième vecteur qui n'est autre que N
|
|
if (Aa0.Nb_ligne()==3)
|
|
// ici on a Ip_3 == N et également g_3=g^3=N
|
|
{for (int al=1 ;al<=2;al++)
|
|
{Aa0(3,al) = N * vi(al) ;
|
|
Aa0(al,3) = giH0.Coordo(al) * N ;
|
|
};
|
|
Aa0(3,3)= 1.; // == N . N
|
|
};
|
|
}
|
|
else // cas on l'on veut le passage des gi et N vers I_a
|
|
{// g^i = Aa^i_{.a} * I^a -> donc Aa^i_{.a} c'est directement les coordonnées de g^i
|
|
for (int al=1 ;al < 4;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa0(2,al) = giH0.Coordo(2)(al) ;
|
|
};
|
|
// pour le troisième: N = Aa^3_{.a} * I^a -> c'est aussi directement les coordonnées de N
|
|
Aa0.RemplaceLigneSpe(3,N.Vect());
|
|
};
|
|
|
|
//// ---- vérification que le déterminant n'est pas nul
|
|
//{double det = Aa0(1,1) * Aa0(2,2) - Aa0(1,2) * Aa0(2,1);
|
|
// if (Abs(det) < 0.00001)
|
|
// cout << "\n erreur det0 null " << det << " dans Deformation::BasePassage( ";
|
|
//};
|
|
//// ---- fin vérification
|
|
|
|
}
|
|
else if (dim==1)
|
|
// ==== cas des éléments 1D en espace 1D, 2D ou 3D ====
|
|
{ if (!absolue)
|
|
{switch (ParaGlob::Dimension())
|
|
{ case 1: // normalement on ne devrait pas passer ici !! car ce cas est déjà traité au début
|
|
{cout << "\nErreur : bug!! on ne devrait pas passer ici \n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
case 3: // on est dans un espace 3D, ici on ne fait que le particulier à 3D et le reste est
|
|
// fait dans le case 2: qui suit
|
|
// 3X1X1 ou 3X1X2 ou 3X1X3, 2X1X1 , 2X1X2
|
|
if (Aa0.Nb_ligne() == 3 )
|
|
// 3X1X3
|
|
{ // tous les tenseurs vont utiliser que le premier vecteur, donc le troisième peut être quelconque
|
|
// et l'expression de giH(3) n'existant pas on dit que par défaut (et par simplicité)
|
|
// giH(3) = 1 * IP(3) d'où
|
|
Aa0 (3,3) = 1.;
|
|
// et on initialise les autres composantes au cas où
|
|
Aa0 (1,3) = Aa0 (3,1) = 0.;
|
|
Aa0 (2,3) = Aa0 (3,2) = 0.;
|
|
// pas de break, on continue sur le 2D
|
|
}
|
|
case 2: // on est dans un espace 2D, 3X1X1 ou 3X1X2 ou 3X1X3
|
|
{ // le premier vecteur intéressant c'est le vecteur collinéaire avec gH(1), mais normé
|
|
double norme0 = giH0(1).Norme();
|
|
Aa0 (1,1) = norme0;
|
|
// fin du cas 3X1X1 et 2X1X1
|
|
if (Aa0.Nb_ligne() > 1)
|
|
// 3X1X2 , 3X1X3 , 2X1X2
|
|
// tous les tenseurs vont utiliser que le premier vecteur, donc le second peut être quelconque
|
|
// et l'expression de giH(2) n'existant pas on dit que par défaut (et par simplicité)
|
|
// giH(2) = 1 * IP(2) d'où
|
|
{Aa0 (2,2) = 1.;
|
|
// et on initialise les autres composantes au cas où
|
|
Aa0 (1,2) = Aa0 (2,1) = 0.;
|
|
};
|
|
break;
|
|
}
|
|
default: cout << "\n erreur grossiere de dimension !! \n Deformation::BasePassage(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else // sinon on veut arriver en Ia
|
|
{ // On ne connait que le premier vecteur et normalement c'est seulement celui-là qui est
|
|
// utilisé, par contre il faut complèter la matrice pour qu'elle ne soit pas singulière
|
|
Aa0.Initialise(0.);
|
|
switch (ParaGlob::Dimension())
|
|
{ case 1: // normalement on ne devrait pas passer ici !! car ce cas est déjà traité au début
|
|
{cout << "\nErreur : bug!! on ne devrait pas passer ici \n";
|
|
cout << "void Deformation::BasePassage( \n";
|
|
Sortie(1);
|
|
break;
|
|
}
|
|
case 3: // on est dans un espace 3D,
|
|
{CoordonneeH V02(3); CoordonneeH V03(3);
|
|
// on calcule deux vecteurs arbitraires normaux, avec un résultat orthonormé V2,V3
|
|
// mais pas giH0(1) !
|
|
MathUtil2::Def_vecteurs_plan(giH0(1),V02,V03);
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa0(2,al) = V02(al) ;
|
|
Aa0(3,al) = V03(al) ;
|
|
};
|
|
break;
|
|
}
|
|
case 2: // on est dans un espace 2D,
|
|
{CoordonneeH V02(2);
|
|
// on calcule le vecteur normal,
|
|
// mais pas giH0(1) !
|
|
MathUtil2::Def_vecteurs_plan(giH0(1),V02);
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
Aa0(2,al) = V02(al) ;
|
|
};
|
|
break;
|
|
}
|
|
default: cout << "\n erreur2 grossiere de dimension !! \n Deformation::BasePassage(.. ";
|
|
Sortie(1);
|
|
};
|
|
int bornesup = ParaGlob::Dimension() +1;
|
|
for (int al=1 ;al < bornesup;al++)
|
|
{Aa0(1,al) = giH0.Coordo(1)(al) ;
|
|
if (al!=1)
|
|
{Aa0 (al,al) = 1.;};
|
|
};
|
|
};
|
|
}
|
|
else
|
|
{// pour les autres cas on ramene la matrice identite
|
|
#ifdef MISE_AU_POINT
|
|
if ((Aa0.Nb_ligne() != dim) && (Aa0.Nb_colonne() !=dim) )
|
|
{ cout << "\nErreur : la matrices Aa0 n'est pas correctement"
|
|
<< " dimensionnee !\n";
|
|
cout << "void Deformation::RemontImp( \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
for (int i=1;i<=dim;i++)
|
|
{ for (int j=1;j<=dim;j++)
|
|
{Aa0 (i,j) = 0.;}
|
|
Aa0 (i,i) = 1.;
|
|
};
|
|
};
|
|
// retour
|
|
return;
|
|
};
|
|
|
|
// on suppose que gijHH et gijBB ont déjà été calculé et sont donc transmis
|
|
// on suppose que gijHH et gijBB ont déjà été calculé à t=0, et à temps
|
|
void Deformation::Cal_deformation (Enum_dure temps, TenseurBB & epsBB)
|
|
{ // ici il y a le choix entre les différents types de calcul de la déformation (Almansi, Green-Lagrange ..
|
|
switch (type_deformation)
|
|
{ case DEFORMATION_STANDART :
|
|
{return Cal_Almansi_auTemps(temps,epsBB);
|
|
break;
|
|
}
|
|
case DEFORMATION_LOGARITHMIQUE:
|
|
{return Cal_logarithmique_auTemps(temps,epsBB);
|
|
break;
|
|
}
|
|
case DEFORMATION_CUMU_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEF_CUMUL_ROTATION_PROPRE :
|
|
{return Cal_def_cumule_auTemps(temps,epsBB);
|
|
break;
|
|
}
|
|
default :
|
|
cout << "\nErreur : type de deformation ( " << Nom_type_deformation(type_deformation)
|
|
<< " ) non traite !\n";
|
|
cout << "Deformation::Cal_deformation (... \n";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// calcul de la deformation au temps donné dans le cas d'une déformation d'Almansi
|
|
// on suppose que les métriques en HH ou BB sont définis à 0 et au temps
|
|
void Deformation::Cal_def_cumule_auTemps (Enum_dure , TenseurBB & )
|
|
{ cout << "\n methode pas encore implante (dommage!) "
|
|
<< "\n Deformation::Cal_def_cumule_auTemps (..";
|
|
Sortie(1);
|
|
};
|
|
|
|
// gestion du parcours de tous les points d'integration
|
|
void Deformation::PremierPtInteg()
|
|
{ sauve_numInteg = numInteg; numInteg = 1; };
|
|
bool Deformation::DernierPtInteg()
|
|
{ sauve_numInteg = numInteg; // sauvegarde
|
|
if (numInteg <= tabPhi->Taille())
|
|
return true;
|
|
else
|
|
return false;
|
|
};
|
|
void Deformation::NevezPtInteg()
|
|
{ sauve_numInteg = numInteg; // sauvegarde
|
|
numInteg++;
|
|
};
|
|
|
|
// méthode pour revenir au pti précédant
|
|
void Deformation::Retour_pti_precedant()
|
|
{numInteg=sauve_numInteg;};
|
|
|
|
//==================================== méthodes protégées ===============================
|
|
|
|
|
|
|
|
// calcul implicite dans le cas d'une déformation tensorielle cumulée
|
|
const Met_abstraite::Impli& Deformation::Cal_implicit_def_cumule (bool
|
|
,TenseurBB & epsBB_tdt,Tableau <TenseurBB *> & d_epsBB_tdt
|
|
,TenseurBB& DepsBB,TenseurBB& delta_epsBB,bool premier_calcul
|
|
,const Met_abstraite::Impli& ex)
|
|
{
|
|
// dans le cas où ce n'est pas le premier calcul on récupère les datas qui ne sont pas recalculée
|
|
if (!premier_calcul)
|
|
{ Deformation::Stmet& a_0 = saveDefResul->meti_00; // pour simplifier
|
|
Deformation::Stmet& a_t = saveDefResul->meti_t; // pour simplifier
|
|
metrique->Recup_grandeur_0_t(*a_0.giB_,*a_0.giH_,*a_t.giB_,*a_t.giH_
|
|
,*a_0.gijBB_,*a_0.gijHH_,*a_t.gijBB_,*a_t.gijHH_
|
|
,*a_t.gradVmoyBB_,*a_t.jacobien_,*a_0.jacobien_);
|
|
};
|
|
// 1- calcul de la base duale et de la métrique à t=0, utiliser en 3-
|
|
// non déjà fait dans le calcul actuel de la métrique
|
|
// const Met_abstraite::gijHH_0_et_giH_0& ex1 = metrique->Cal_gijHH_0_et_giH_0_apres_im_expli();
|
|
// 2- calcul de l'opérateur gradient F(a,b)=d Xc^a/d X0^b, stockage sous forme d'un spseudo-tenseur
|
|
// alors que c'est un bi-tenseur !! permet ainsi une manipulation aisée
|
|
int dim = ParaGlob::Dimension();
|
|
TenseurHB * ftHB = NevezTenseurHB(dim);
|
|
for (int i=1;i<=dim;i++)
|
|
{ TenseurHB * interHB = Produit_tensorielHB((*ex.giH_0)(i),(*ex.giB_tdt)(i));
|
|
*ftHB += *interHB; delete interHB;
|
|
};
|
|
// 3- calcul du gradient de vitesse
|
|
// plusieurs possibilités,
|
|
// 1) . soit on calcul le gradient à partir du champ de vitesse au temps finale
|
|
// 2) . soit on calcul le gradient à partir d'une hypothèse de vitesse contante sur le pas de temps
|
|
// avec la vitesse = delta X / delta t, le gradient étant déterminé à t+ (deltat/2)
|
|
// 3) . soit on calcul le gradient à partir d'une hypothèse de vitesse moyenne au temps t+ (deltat/2)
|
|
// c-a-d en utilisant la vitesse moyenne entre t et t+dt, ceci dans le cas de l'utilisation
|
|
// directe des ddl vitesse
|
|
|
|
/* switch (type_gradient_vitesse) // traitement en fonction du cas
|
|
{ case GRADVITESSE_V_TDT : // le gradient est déterminé à partir de ddl de vitesse à t+dt
|
|
#ifdef MISE_AU_POINT
|
|
// on vérifie que le gradient a été calculé dans l'appel de la métrique
|
|
if (pas_de_gradV)
|
|
{ cout << "\nErreur : le gradient de vitesse n'a pas ete calcule dans la metrique"
|
|
<< " ce qui est contraire au type de calcul de gradient demande";
|
|
cout << "\n Deformation::Cal_implicit_def_cumule(... \n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// donc ici le gradient est dans (*ex.gradVBB_tdt)
|
|
break;
|
|
case GRADVITESSE_VCONST : // calcul en considérant V constant =delta X/delta t
|
|
// dans la config à t+(delta t)/2
|
|
|
|
// // 2) calcul du gradient à partir du champ de vitesse au temps final
|
|
// if (!pas_de_gradV)
|
|
// {// cas où le gradient a été calculé
|
|
// DepsBB = 0.5 * ((*ex.gradVBB_tdt) + ex.gradVBB_tdt->Transpose());}
|
|
// else
|
|
// {
|
|
// }
|
|
break;
|
|
case GRADVITESSE_VT_VTDT : //calcul à partir des ddl de vitesse : moyenne de t et t+dt,
|
|
//dans la config à t+(delta t)/2
|
|
break;
|
|
};*/
|
|
|
|
|
|
|
|
/* // epsBB_tdt et delta_epsBB est dimensionne auparavant
|
|
epsBB_tdt = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_0));
|
|
delta_epsBB = 0.5 * (*(ex.gijBB_tdt) - *(ex.gijBB_t));
|
|
if (pas_de_gradV)
|
|
{ // 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 on garde la vitesse actuelle
|
|
double deltat=vartemps.IncreTempsCourant();
|
|
if (deltat >= ConstMath::trespetit)
|
|
DepsBB = delta_epsBB/deltat;
|
|
}
|
|
else // cas où le gradient a été calculé
|
|
DepsBB = 0.5 * ((*ex.gradVBB_tdt) + ex.gradVBB_tdt->Transpose());
|
|
// variation de la déformation / au ddl
|
|
int d_epsBB_tdtTaille = d_epsBB_tdt.Taille();
|
|
for (int i=1; i<= d_epsBB_tdtTaille; i++)
|
|
*(d_epsBB_tdt(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
|
|
// calcul éventuelle de la dérivée seconde de la déformation
|
|
if (pa.Var_D())
|
|
{// recup des variations secondes de la déformation
|
|
int d2_epsBB_tdtTaille1 = d2_epsBB_tdt.Taille1();
|
|
for (int i=1; i<= d2_epsBB_tdtTaille1; i++)
|
|
for (int j=1; j<= i; j++) // symétrie du tableau et du tenseur
|
|
{ *(d2_epsBB_tdt(i,j)) = 0.5 * (*((*(ex.d2_gijBB_tdt))(i,j))) ;
|
|
*(d2_epsBB_tdt(j,i)) = *(d2_epsBB_tdt(i,j)) ;
|
|
}
|
|
} */
|
|
// vérification éventuelle
|
|
// VerifCal_implicit( ex);
|
|
// suppression des tenseurs de travail
|
|
delete ftHB;
|
|
// retour
|
|
return ex;
|
|
};
|
|
|
|
// calcul d'une vitesse de rotation cohérente avec un référentiel objectif particulier
|
|
void Deformation::Cal_vite_rota_objectif(EnumTypeViteRotat type_rotation,bool variation
|
|
,const TenseurBB& gradVBB,const Tableau <TenseurBB * >& d_gradVBB
|
|
,const TenseurBH& gij_HH_0,const TenseurBB& gij_BB,Tableau <TenseurBB *> & d_gijBB
|
|
,const TenseurHH& gij_HH,Tableau <TenseurHH *> & d_gijHH
|
|
,const TenseurBB& D_BB, Tableau <TenseurBB *> & d_D_BB
|
|
,TenseurBB& OmegaBB,Tableau <TenseurBB* > & d_OmegaBB)
|
|
{ // on calcul le pseudo vecteur W représentant le tenseur anti-symétrique: rotation du gradient de vitesse
|
|
// TenseurBB W_BB = 0.5 * (gradVBB - gradVBB.Transpose());
|
|
int dima=Abs(gradVBB.Dimension());
|
|
/* CoordonneeB W_B(dima);
|
|
switch (dima) { case 3: W_B(1) = W_BB(3,2);W_B(2) = W_BB(1,3);W_B(3) = W_BB(2,1);break;
|
|
case 2: W_B(1) = W_BB(1,2);W_B(2) = 0.;break; // en attente
|
|
case 1: W_B(1) = 0.;break;}; // en attente*/
|
|
int dimddl=d_gradVBB.Taille(); // dimension pour les variations
|
|
switch (type_rotation)
|
|
{ case R_CORROTATIONNEL:
|
|
OmegaBB =0.5 * (gradVBB - gradVBB.Transpose());
|
|
// dans le cas où l'on veut également des variations
|
|
if (variation)
|
|
{ for (int i=1;i<=dimddl;i++) (*(d_OmegaBB(i))) = 0.5*( (*(d_gradVBB(i))) - (d_gradVBB(i))->Transpose()); }
|
|
break;
|
|
case R_REF_ROT_PROPRE: case R_ROT_LOGARITHMIQUE:
|
|
// calcul des projections propres
|
|
switch (dima)
|
|
{ case 3:
|
|
{// simplification et transformation en tenseur de dim 3 (optimisation des calculs)
|
|
const Tenseur3HH & gijHH_0 = (const Tenseur3HH &) gij_HH_0;
|
|
const Tenseur3BB & gijBB = (const Tenseur3BB &) gij_BB;
|
|
const Tenseur3HH & gijHH = (const Tenseur3HH &) gij_HH;
|
|
// def des grandeurs de travail
|
|
DimensionementVarLog(dima,variation,dimddl);
|
|
// def du tenseur B_BH (cauchy_green gauche)
|
|
(*B_BH_tr) = gijBB * gijHH_0; Tenseur3BH& B_BH= *((Tenseur3BH*) B_BH_tr);
|
|
// calcul des valeurs propres et projection propres du tenseur B_BH
|
|
Tableau <Tenseur3BH * >& d_B_BH=*((Tableau <Tenseur3BH * >*) d_B_BH_tr);
|
|
if (variation) // dans le cas de variation, calcul des variation de B_BH
|
|
{ for (int ic=1;ic<=dimddl;ic++) {(*(*d_B_BH_tr)(ic))= (*(d_gijBB(ic))) * gijHH_0;};};
|
|
// def des grandeurs de travail
|
|
Tableau <Tenseur3BH* >& Palpha_BH=*((Tableau <Tenseur3BH* >*) Palpha_BH_tr); // les projections propres
|
|
Tableau <Tableau <Tenseur3BH* > >& d_Palpha_BH
|
|
=*((Tableau <Tableau <Tenseur3BH* > >*) d_Palpha_BH_tr); // variation des projections propres
|
|
// appel de la méthode calculant les grandeurs de travail
|
|
int cas_ki; // indique le cas des valeurs propres traitées (cf TenseurBH)
|
|
Val_et_projection_prop_tenseur(*B_BH_tr,*d_B_BH_tr,ki,*Palpha_BH_tr,variation,*d_ki,*d_Palpha_BH_tr,cas_ki);
|
|
// -- calcul de la rotation
|
|
// calcul de la vitesse de rotation en mixte
|
|
Tenseur3BH D_BH(D_BB * gijHH);
|
|
// tout d'abord: calcul de tenseurs intermédiaires
|
|
Tableau <Tenseur3BH> P_i_D_BH(3); Tableau <Tenseur3BH> D_P_i_BH(3);
|
|
P_i_D_BH(1)=((*Palpha_BH(1)) * D_BH); Tenseur3BH P_1_D_P_2_BH(P_i_D_BH(1) * (*Palpha_BH(2)));
|
|
P_i_D_BH(2)=((*Palpha_BH(2)) * D_BH); Tenseur3BH P_2_D_P_1_BH(P_i_D_BH(2) * (*Palpha_BH(1)));
|
|
D_P_i_BH(1)=(D_BH * (*Palpha_BH(1))); Tenseur3BH P_1_D_P_3_BH(P_i_D_BH(1) * (*Palpha_BH(3)));
|
|
P_i_D_BH(3)=((*Palpha_BH(3)) * D_BH); Tenseur3BH P_3_D_P_1_BH(P_i_D_BH(3) * (*Palpha_BH(1)));
|
|
D_P_i_BH(2)=(D_BH * (*Palpha_BH(2))); Tenseur3BH P_2_D_P_3_BH(P_i_D_BH(2) * (*Palpha_BH(3)));
|
|
D_P_i_BH(3)=(D_BH * (*Palpha_BH(3))); Tenseur3BH P_3_D_P_2_BH(P_i_D_BH(3) * (*Palpha_BH(2)));
|
|
|
|
double F_rot_1sur2= F_pour_rotat_objectif(type_rotation,ki(1)/ki(2));
|
|
double F_rot_1sur3= F_pour_rotat_objectif(type_rotation,ki(1)/ki(3));
|
|
double F_rot_2sur1= F_pour_rotat_objectif(type_rotation,ki(2)/ki(1));
|
|
double F_rot_3sur1= F_pour_rotat_objectif(type_rotation,ki(3)/ki(1));
|
|
double F_rot_2sur3= F_pour_rotat_objectif(type_rotation,ki(2)/ki(3));
|
|
double F_rot_3sur2= F_pour_rotat_objectif(type_rotation,ki(3)/ki(2));
|
|
|
|
OmegaBB =0.5 * (gradVBB - gradVBB.Transpose());
|
|
Tenseur3BH sup_OmegaBH ( F_rot_1sur2 * P_1_D_P_2_BH + F_rot_2sur1 * P_2_D_P_1_BH
|
|
+ F_rot_1sur3 * P_1_D_P_3_BH + F_rot_3sur1 * P_3_D_P_1_BH
|
|
+ F_rot_2sur3 * P_2_D_P_3_BH + F_rot_3sur2 * P_3_D_P_2_BH );
|
|
OmegaBB += sup_OmegaBH* gijBB ;
|
|
|
|
// calcul éventuel des variations
|
|
if (variation)
|
|
{ // mise sous forme indicée suivant l'ordre 1->12 2->21 3->13 4->31 5->23 6->32 à partir de 6 nombres
|
|
Tableau <Tenseur3BH*> P_i_D_j_BH(6);
|
|
P_i_D_j_BH(1)=&P_1_D_P_2_BH; P_i_D_j_BH(2)=&P_2_D_P_1_BH; P_i_D_j_BH(3)=&P_1_D_P_3_BH;
|
|
P_i_D_j_BH(4)=&P_3_D_P_1_BH; P_i_D_j_BH(5)=&P_2_D_P_3_BH; P_i_D_j_BH(6)=&P_3_D_P_2_BH;
|
|
Tableau <double *> F_rot_isurj(6);
|
|
F_rot_isurj(1)=&F_rot_1sur2; F_rot_isurj(2)=&F_rot_2sur1; F_rot_isurj(3)=&F_rot_1sur3;
|
|
F_rot_isurj(4)=&F_rot_3sur1; F_rot_isurj(5)=&F_rot_2sur3; F_rot_isurj(6)=&F_rot_3sur2;
|
|
// boucle sur les ddl
|
|
for (int ici=1;ici<=dimddl;ici++)
|
|
{ const Tenseur3BB & dgijBB = *((Tenseur3BB*)(d_gijBB(ici))) ; // pour simplifier l'ecriture
|
|
const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH(ici))) ; // pour simplifier l'ecriture
|
|
const Tenseur3BB & dD_BB = *((Tenseur3BB*)(d_D_BB(ici))) ; // pour simplifier l'ecriture
|
|
|
|
// variation du tenseur vitesse de déformation en mixte
|
|
Tenseur3BH dD_BH(dD_BB * gijHH + D_BB * dgijHH);
|
|
// variation du tenseur de rotation
|
|
Tenseur3BB & dOmegaBB = *((Tenseur3BB* ) (d_OmegaBB(ici))); // pour simplifier l'ecriture
|
|
dOmegaBB = 0.5*( (*(d_gradVBB(ici))) - (d_gradVBB(ici))->Transpose());
|
|
for (int m=1;m<=6;m++)
|
|
{ int i1=ccdex.iii(m); int i2=ccdex.jjj(m);
|
|
dOmegaBB += (( Fprime_pour_rotat_objectif(type_rotation,ki(i1)/ki(i2))
|
|
/ (ki(i2)*ki(i2)) * (ki(i2) * ((*d_ki)(ici))(i1) - ki(i1) * ((*d_ki)(ici))(i2)) ) * (* P_i_D_j_BH(m))
|
|
+ (*F_rot_isurj(m)) * ( (*d_Palpha_BH(ici)(i2)) * D_P_i_BH(i2) + (*Palpha_BH(i1)) * dD_BH * (*Palpha_BH(i2))
|
|
+ P_i_D_BH(i1) * (*d_Palpha_BH(ici)(i2)) )) * gijBB
|
|
+ sup_OmegaBH * dgijBB;
|
|
};
|
|
};
|
|
}; // -- fin du grand "if variation "
|
|
// effacement des tenseurs de travail
|
|
for (int it=1;it<=dima;it++) delete Palpha_BH(it);
|
|
if (variation) for (int ic=1;ic<=dimddl;ic++) for (int it=1;it<=dima;it++)
|
|
delete d_Palpha_BH(ic)(it);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n dimension (premier) non encore pris en compte: " << dima
|
|
<< "\n Deformation::Cal_vite_rota_objectif(...";
|
|
Sortie(1);
|
|
}
|
|
break;
|
|
default:
|
|
cout << "\n type de rotation non encore pris en compte"
|
|
<< "\n Deformation::Cal_vite_rota_objectif(...";
|
|
Sortie(1);
|
|
};
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
};
|
|
|
|
// fonction f et f' utilisé dans la méthode Cal_vite_rota_objectif:
|
|
double Deformation::F_pour_rotat_objectif(EnumTypeViteRotat type_rotation,const double& x)
|
|
{switch (type_rotation)
|
|
{ case R_CORROTATIONNEL: return 0.; break;
|
|
case R_REF_ROT_PROPRE: return (1.-sqrt(x))/(1.+sqrt(x)); break;
|
|
case R_ROT_LOGARITHMIQUE: return (1.+x)/(1.-x) + 2./log(x); break;
|
|
default:
|
|
cout << "\n type de rotation non encore pris en compte: " << Nom_TypeViteRotat(type_rotation)
|
|
<< "\n Deformation::F_pour_rotat_objectif(...";
|
|
Sortie(1);
|
|
}
|
|
return 0.; // pour éviter le warning du retour
|
|
};
|
|
double Deformation::Fprime_pour_rotat_objectif(EnumTypeViteRotat type_rotation,const double& x)
|
|
{switch (type_rotation)
|
|
{ case R_CORROTATIONNEL: return 0.; break;
|
|
case R_REF_ROT_PROPRE: return 1./(2.*x+sqrt(x)*(x+1.)); break;
|
|
case R_ROT_LOGARITHMIQUE:
|
|
{double logx = log(x); double log2_x=logx*logx;
|
|
return 2.*(x*log2_x-x*x+2.*x-1.)/((x*x*x-2.*x*x+x)*log2_x); break;}
|
|
default:
|
|
cout << "\n type de rotation non encore pris en compte: " << Nom_TypeViteRotat(type_rotation)
|
|
<< "\n Deformation::Fprime_pour_rotat_objectif(...";
|
|
Sortie(1);
|
|
}
|
|
return 0.; // pour éviter le warning du retour
|
|
};
|
|
|
|
// ----------------- méthodes de vérifications------- ----
|
|
|
|
void Deformation::VerifCal_implicit(bool gradV_instantane,const Met_abstraite::Impli & ex,TenseurBB& )
|
|
{ // 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_VerifCal_implicit == 0)
|
|
{ cout << "\n ****vérification du calcul de la déformation et des éléments de la métrique associé****";
|
|
cout << "\n Deformation::VerifCal_implicit \n";
|
|
}
|
|
indic_VerifCal_implicit++;
|
|
|
|
// on cré une seconde métrique pour éviter de détruire la métrique originale
|
|
Met_abstraite metrique_bis(*metrique);
|
|
// 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 = tabnoeud->Taille();
|
|
// le deltat pour les différences finis
|
|
double delta = ConstMath::unpeupetit;
|
|
double mini_val = ConstMath::pasmalpetit;
|
|
int numddl = 1; // le compteur de ddl
|
|
bool erreur = false; // indicateur d'erreur
|
|
bool premier_calcul=true;
|
|
for (int inoeud=1;inoeud<=nbnoeud;inoeud++)
|
|
{// on récupère les coordonnées du noeud
|
|
Coordonnee coordtdt = (*tabnoeud)(inoeud)->Coord2();
|
|
for (int ix= 1;ix<=dim;ix++,numddl++)
|
|
{ Coordonnee X(dim); X(ix) += delta;
|
|
(*tabnoeud)(inoeud)->Ajout_coord2(X);
|
|
// appel de la métrique
|
|
Met_abstraite::Expli_t_tdt ex_n =metrique_bis.Cal_explicit_tdt
|
|
(*tabnoeud,premier_calcul,(*tabDphi)(numInteg),nbNoeud,(*tabPhi)(numInteg),gradV_instantane);
|
|
premier_calcul=false;
|
|
|
|
// calcul des dérivées numériques et vérification
|
|
for (int j=1;j<=dim;j++)
|
|
{ // variation des vecteurs giB_tdt
|
|
CoordonneeB dgiB = ((*ex_n.giB_tdt)(j) -(*ex.giB_tdt)(j))/delta;
|
|
for (int i=1;i<=dim;i++)
|
|
if (diffpourcent(dgiB(i),(*ex.d_giB_tdt)(numddl)(j)(i),MaX(Dabs(dgiB(i)),Dabs((*ex.d_giB_tdt)(numddl)(j)(i))),0.05))
|
|
{if (MiN(Dabs(dgiB(i)),Dabs((*ex.d_giB_tdt)(numddl)(j)(i))) <= mini_val)
|
|
{if ( MaX(Dabs(dgiB(i)),Dabs((*ex.d_giB_tdt)(numddl)(j)(i))) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
// variation des vecteurs giH_tdt
|
|
CoordonneeH dgiH = ((*ex_n.giH_tdt)(j) - (*ex.giH_tdt)(j))/delta;
|
|
for (int i=1;i<=dim;i++)
|
|
if (diffpourcent(dgiH(i),(*ex.d_giH_tdt)(numddl)(j)(i),MaX(Dabs(dgiH(i)),Dabs((*ex.d_giH_tdt)(numddl)(j)(i))),0.05))
|
|
{if (MiN(Dabs(dgiH(i)),Dabs((*ex.d_giH_tdt)(numddl)(j)(i))) <= mini_val)
|
|
{if ( MaX(Dabs(dgiH(i)),Dabs((*ex.d_giH_tdt)(numddl)(j)(i))) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
}
|
|
// variation du tenseur gijBB_tdt
|
|
TenseurBB * gijBB = NevezTenseurBB(*ex_n.gijBB_tdt);
|
|
*gijBB = (*ex_n.gijBB_tdt - *ex.gijBB_tdt) / delta;
|
|
for (int i1=1;i1<=dim;i1++)
|
|
for (int j1=1;j1<=dim;j1++)
|
|
if (diffpourcent((*gijBB)(i1,j1),(*(*ex. d_gijBB_tdt)(numddl))(i1,j1),
|
|
MaX(Dabs((*gijBB)(i1,j1)),Dabs((*(*ex. d_gijBB_tdt)(numddl))(i1,j1))),0.05))
|
|
{if (MiN(Dabs((*gijBB)(i1,j1)),Dabs((*(*ex. d_gijBB_tdt)(numddl))(i1,j1))) <= mini_val)
|
|
{if( MaX(Dabs((*gijBB)(i1,j1)),Dabs((*(*ex. d_gijBB_tdt)(numddl))(i1,j1))) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
// variation du tenseur gijHH_tdt
|
|
TenseurHH * gijHH = NevezTenseurHH(*ex_n.gijHH_tdt);
|
|
*gijHH = (*ex_n.gijHH_tdt - *ex.gijHH_tdt) / delta;
|
|
for (int i1=1;i1<=dim;i1++)
|
|
for (int j1=1;j1<=dim;j1++)
|
|
if (diffpourcent((*gijHH)(i1,j1),(*(*ex. d_gijHH_tdt)(numddl))(i1,j1),
|
|
MaX(Dabs((*gijHH)(i1,j1)),Dabs((*(*ex. d_gijHH_tdt)(numddl))(i1,j1))),0.05))
|
|
{if (MiN(Dabs((*gijHH)(i1,j1)),Dabs((*(*ex. d_gijHH_tdt)(numddl))(i1,j1))) <= mini_val)
|
|
{if( MaX(Dabs((*gijHH)(i1,j1)),Dabs((*(*ex. d_gijHH_tdt)(numddl))(i1,j1))) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
// variation du jacobien
|
|
double djaco = ((*ex_n.jacobien_tdt) - (*ex.jacobien_tdt))/delta;
|
|
if (diffpourcent(djaco,(*ex.d_jacobien_tdt)(numddl),MaX(Dabs(djaco),Dabs((*ex.d_jacobien_tdt)(numddl))),0.05))
|
|
if (MiN(Dabs(djaco),Dabs((*ex.d_jacobien_tdt)(numddl))) <= mini_val)
|
|
{if( MaX(Dabs(djaco),Dabs((*ex.d_jacobien_tdt)(numddl))) > 50.*delta) erreur = true;
|
|
else erreur = true;
|
|
};
|
|
// effacement des tenseurs intermédiaires
|
|
delete gijBB; delete gijHH;
|
|
// maintenant on remet les coordonnées du noeud à l'état initial
|
|
(*tabnoeud)(inoeud)->Change_coord2(coordtdt);
|
|
} // fin de boucle sur la dimension de coordonnée
|
|
} // fin de boucle sur les noeuds
|
|
// message d'erreur si besoin
|
|
if (erreur)
|
|
{ cout << "\n erreur dans le calcul analytique des dérivees de la metrique";
|
|
cout << "\n Deformation::VerifCal_implicit(.."
|
|
<< " , numero d'increment = " << indic_VerifCal_implicit;
|
|
Sortie(1);
|
|
}
|
|
};
|
|
|
|
|
|
|
|
|