// 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 "ElemMeca.h"
#include <iomanip>
#include "ConstMath.h"
#include "Util.h"
#include "Coordonnee2.h"
#include "Coordonnee3.h"
#include "TypeQuelconqueParticulier.h"
#include "TypeConsTens.h"
#include "Loi_iso_elas3D.h"
#include "Loi_iso_elas1D.h"
#include "Loi_iso_elas2D_C.h"
#include "FrontPointF.h"
#include "MathUtil2.h"
#include "Enum_StabMembrane.h"
// -------------------- stabilisation d'hourglass -------------------
// calcul d'élément de contrôle d'hourglass associée à un comportement
// str_precision : donne le type particulier d'élément à construire
void ElemMeca::Init_hourglass_comp(const ElemGeomC0& elgeHour, const string & str_precision
,LoiAbstraiteGeneral * loiHourglass, const BlocGen & bloc)
{
//!!!! en fait pour l'instant on n'utilise pas elgehour, on considère que le nombre de pti complet est celui de l'élément sans info annexe !!!
// sauf pour l'hexaèdre quadratique pour lequel on a définit un str_precision particulier
// on met à jour les indicateurs
if (Type_Enum_StabHourglass_existe(bloc.Nom(1)))
{ type_stabHourglass = Id_Nom_StabHourglass(bloc.Nom(1).c_str());
coefStabHourglass = bloc.Val(1);
}
else
{ cout << "\n erreur, le type " << bloc.Nom(2) << " de gestion d'hourglass n'existe pas !!";
if (ParaGlob::NiveauImpression() > 5)
cout << "\n ElemMeca::Init_hourglass_comp(...";
cout << endl;
Sortie(1);
};
switch (type_stabHourglass)
{case STABHOURGLASS_PAR_COMPORTEMENT :
{ // -- choix de l'element pour le calcul de la raideur -----
// on recherche un élément de même type,
// par défaut : numéro de maillage = le numéro de this, car c'est bien rattaché à ce maillage, mais l'élément est caché, donc non accessible par le maillage
// numéro d'élément = 1000000 + numéro de this : a priori c'est deux numéros n'ont pas d'importance, car ils ne sont pas
// rattachés à un maillage existant en tant que tel
int nUmMail = this->Num_elt_const();
int nUmElem = 1000000 + this->Num_elt();
tab_elHourglass.Change_taille(1);
tab_elHourglass(1) = (ElemMeca*) Element::Choix_element(nUmMail,nUmElem,Id_geometrie(),Id_interpolation(),Id_TypeProblem(),str_precision);
if (tab_elHourglass(1) == NULL)
{ // l'operation a echouee
cout << "\n Erreur dans le choix d'un element interne utilise pour le blocage d'hourglass ****** ";
cout << "\n l\'element : " << Nom_interpol(id_interpol)
<< " " << Nom_geom(id_geom) << " " << str_precision
<< " n\'est pas present dans le programme ! " << endl;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n ElemMeca::Init_hourglass_comp(..." << endl;
Sortie (1);
};
// --- définition des noeuds de l'élément -----
// on construit à partir des mêmes noeuds que ceux de l'élément père
// affectation des noeuds au nouvel élément
tab_elHourglass(1)->Tab_noeud() = this->Tab_noeud();
//--- def de la loi de comportement ---
// affectation de la loi
tab_elHourglass(1)->DefLoi(loiHourglass);
break;
}
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // ici on a besoin de deux éléments. Le premier sert à faire le calcul complet
// le second sert à faire un calcul réduit, comme l'élément réel
tab_elHourglass.Change_taille(2);
// -- choix pour le calcul de la raideur -----
// on recherche un élément de même type,
// par défaut : numéro de maillage = le numéro de this, car c'est bien rattaché à ce maillage, mais l'élément est caché, donc non accessible par le maillage
// numéro d'élément = 1000000 + numéro de this : a priori c'est deux numéros n'ont pas d'importance, car ils ne sont pas
// rattachés à un maillage existant en tant que tel
//-- le premier élément:
int nUmMail = this->Num_elt_const();
int nUmElem = 1000000 + this->Num_elt();
tab_elHourglass(1) = (ElemMeca*) Element::Choix_element(nUmMail,nUmElem,Id_geometrie(),Id_interpolation(),Id_TypeProblem(),str_precision);
if (tab_elHourglass(1) == NULL)
{ // l'operation a echouee
cout << "\n Erreur dans le choix d'un element interne utilise pour le blocage d'hourglass ****** ";
cout << "\n l\'element : " << Nom_interpol(id_interpol)
<< " " << Nom_geom(id_geom) << " " << str_precision
<< " n\'est pas present dans le programme ! " << endl;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n ElemMeca::Init_hourglass_comp(..." << endl;
Sortie (1);
};
//-- le second élément identique à this: on utilise le même numéro de maillage
// même infos annexes
nUmElem = 20000000 + this->Num_elt(); // on ajoute une unité pour le numéro
tab_elHourglass(2) = (ElemMeca*) Element::Choix_element(nUmMail,nUmElem,Id_geometrie(),Id_interpolation(),Id_TypeProblem(),this->infos_annexes);
if (tab_elHourglass(2) == NULL)
{ // l'operation a echouee
cout << "\n Erreur dans le choix d'un element interne utilise pour le blocage d'hourglass ****** ";
cout << "\n l\'element : " << Nom_interpol(id_interpol)
<< " " << Nom_geom(id_geom) << " " << this->infos_annexes
<< " n\'est pas present dans le programme ! " << endl;
if (ParaGlob::NiveauImpression() > 5)
cout << "\n ElemMeca::Init_hourglass_comp(..." << endl;
Sortie (1);
};
// --- définition des noeuds de l'élément -----
// on construit à partir des mêmes noeuds que ceux de l'élément père
// affectation des noeuds au nouvel élément
tab_elHourglass(1)->Tab_noeud() = this->Tab_noeud();
tab_elHourglass(2)->Tab_noeud() = this->Tab_noeud();
//--- def de la loi de comportement ---
// affectation de la loi
tab_elHourglass(1)->DefLoi(loiHourglass);
tab_elHourglass(2)->DefLoi(loiHourglass);
// def de trois vecteurs de travail s'ils ne sont pas encore définit
if (vec_hourglass.Taille() < 3)
vec_hourglass.Change_taille(3);
break;
}
default :
cout << "\n*** Erreur : cas de gestop, d'hourglass non defini !\n";
cout << "\n ElemMeca::Cal_implicit_hourglass() \n";
Sortie(1);
};
};
// stabilisation pour un calcul implicit
void ElemMeca::Cal_implicit_hourglass()
{E_Hourglass=0.; // init par défaut
switch (type_stabHourglass)
{case STABHOURGLASS_PAR_COMPORTEMENT :
{ // --- calcul de la raideur de l'élément, dans le cas implicite ---
Element::ResRaid resraid = tab_elHourglass(1)->Calcul_implicit (ParaGlob::param->ParaAlgoControleActifs());
// on tiend compte du facteur de modération
(*(resraid.raid)) *= coefStabHourglass;
(*(resraid.res)) *= coefStabHourglass;
// on met à jour la raideur et le résidu de l'élément principal
(*raideur) += (*(resraid.raid));
(*residu) += (*(resraid.res));
// on récupère les énergies stockées à l'élément
const EnergieMeca& energieTotale = tab_elHourglass(1)->EnergieTotaleElement();
E_Hourglass = coefStabHourglass * (energieTotale.EnergieElastique()
+ energieTotale.DissipationPlastique() + energieTotale.DissipationVisqueuse());
//---- debug
//cout << "\n raideur locale après stabilisation d'hourglass ";
//matRaideur.Affiche();
//---- fin debug
};
break;
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // --- calcul de la raideur de l'élément, dans le cas implicite ---
// ici on calcul la partie supposée supplémentaire à celle déjà présente dans le calcul réduit
// au premier passage ou bien si on demande une grande précision de calcul, on calcul la raideur et le second membre,
// ensuite on ne calcul que le second membre
// if ((raid_hourglass_transitoire == NULL)||(prepa_niveau_precision > 8))
// en l'expérience prouve que c'est la forme des éléments du début qui est la moins tordu donc vaut mieux
// ne pas recalculer la raideur sur des éléments tordus, c'est dans ce cas que ça marche le mieux
if (raid_hourglass_transitoire == NULL)
{// 1) calcul de la partie complète
Element::ResRaid resraid1 = tab_elHourglass(1)->Calcul_implicit (ParaGlob::param->ParaAlgoControleActifs());
// 2) calcul de la partie réduite
Element::ResRaid resraid2 = tab_elHourglass(2)->Calcul_implicit (ParaGlob::param->ParaAlgoControleActifs());
// on soustrait en gardant le résultat dans resraid1
// **** non ce n'est pas une bonne idée, lorsque l'on fait certain test de stabilité: lorsque le facteur est très grand
// si on ne fait pas la soustraction cela converge, avec la soustraction cela ne converge pas donc il ne faut pas soustraire
// (*(resraid1.raid)) -= (*(resraid2.raid));
// (*(resraid1.res)) -= (*(resraid2.res));
// ---nouvelle idée:
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud.Taille();
int dim = ParaGlob::Dimension();
if(ParaGlob::AxiSymetrie())
dim--; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur& D_X = vec_hourglass(1); // pour simplier
Vecteur& X = vec_hourglass(2); // pour simplier
D_X.Change_taille(NBN*dim); // changement de la taille au cas où
X.Change_taille(NBN*dim); // changement de la taille au cas où
for (int ine=1;ine<=NBN;ine++)
{ Noeud& noe = *(tab_noeud(ine)); // pour simplifier
Coordonnee delta_X = noe.Coord2() - noe.Coord1();
Coordonnee Xnoeud = noe.Coord2();
int ipos = (ine-1)*dim;
switch (dim)
{ case 3: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);D_X(ipos+3)=delta_X(3);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);X(ipos+3)=Xnoeud(3);break;}
case 2: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);break;}
case 1: {D_X(ipos+1)=delta_X(1);
X(ipos+1)=Xnoeud(1);break;}
};
};
// on multiplie la matrice par le vecteur delta_X
(*(resraid1.res)) = resraid1.raid->Prod_mat_vec(D_X);
// on sauvegarde la raideur initiale
if (raid_hourglass_transitoire == NULL)
raid_hourglass_transitoire = new Mat_pleine((*(resraid1.raid)));
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * (D_X * (*(resraid1.res)));
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass;
if (E_Hourglass < 0.) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = -coefStabHourglass;};
// grandeurs finales avec le bon signe
(*(resraid1.raid)) *= coefMultiplicatif;
(*(resraid1.res)) *= coefMultiplicatif;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * (X * (*(resraid1.res)));
//debug
//cout << "\n premier passage, ";
//res.Affiche(); (resraid1.raid)->Affiche();
//fin debug
// on met à jour la raideur et le résidu de l'élément principal
(*raideur) += (*(resraid1.raid));
(*residu) -= (*(resraid1.res));
}
else
{// sinon on utilise la précédente raideur sauvegardée et on ne calcul que la partie résidue
// dans la nouvelle idée le vecteur doit exister à la bonne dimension, mais CalculResidu_tdt n'a pas besoin
// d'être appeler **** a améliorer
// Vecteur* resHourglass1 = NULL;
// resHourglass1 = tab_elHourglass(1)->CalculResidu_tdt (ParaGlob::param->ParaAlgoControleActifs());
/* Vecteur* resHourglass2 = NULL;
resHourglass2 = tab_elHourglass(2)->CalculResidu_tdt (ParaGlob::param->ParaAlgoControleActifs());
// on soustrait les résidus en gardant le résultat dans resraid1
(*(resHourglass1)) -= (*(resHourglass2));*/
// ---nouvelle idée:
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud.Taille();
int dim = ParaGlob::Dimension();
if(ParaGlob::AxiSymetrie())
dim--; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur& D_X = vec_hourglass(1); // pour simplier
Vecteur& X = vec_hourglass(2); // pour simplier
D_X.Change_taille(NBN*dim); // changement de la taille au cas où
X.Change_taille(NBN*dim); // changement de la taille au cas où
for (int ine=1;ine<=NBN;ine++)
{ Noeud& noe = *(tab_noeud(ine)); // pour simplifier
Coordonnee delta_X = noe.Coord2() - noe.Coord1();
Coordonnee Xnoeud = noe.Coord2();
int ipos = (ine-1)*dim;
switch (dim)
{ case 3: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);D_X(ipos+3)=delta_X(3);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);X(ipos+3)=Xnoeud(3);break;}
case 2: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);break;}
case 1: {D_X(ipos+1)=delta_X(1);
X(ipos+1)=Xnoeud(1);break;}
};
};
//debug
//cout << " res ddl, ";
//res.Affiche();
//fin debug
Vecteur& res = vec_hourglass(2); // pour simplier
res.Change_taille(NBN*dim); // changement de la taille au cas où
// on multiplie le vecteur delta_X par la matrice
res = raid_hourglass_transitoire->Prod_mat_vec(D_X);
// on définit un facteur de modération qui est fondé sur le ratio des énergies
// --- essai d'adaptation du coef, mais ça ne marche pas du tout, il n'y a que quand c'est fixe que c'est ok
/*
double energieHourglassDeBase = 0.5 * (D_X * res);
const EnergieMeca& energieActuelle = this->EnergieTotaleElement(); // récup de l'énergie actuelle
double val_ref_energie = (Dabs(energieActuelle.EnergieElastique())
+ Dabs(energieActuelle.DissipationPlastique())
+ Dabs(energieActuelle.DissipationVisqueuse()));
double coef_a_appliquer = coefStabHourglass;
if (Dabs(energieHourglassDeBase*coefStabHourglass) > MaX(ConstMath::petit,val_ref_energie))
coef_a_appliquer *= val_ref_energie / energieHourglassDeBase;
// on tiend compte du facteur de modération
res *= coef_a_appliquer;
// on met à jour le résidu de l'élément principal
(*residu) -= res;
// pour la partie raideur: on met à jour la raideur de l'élément principal
(*raideur) += coef_a_appliquer * (*raid_hourglass_transitoire);
// calcul de l'énergie
E_Hourglass = energieHourglassDeBase * coef_a_appliquer; //0.5 * (D_X * res);
*/
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * (D_X * res);
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass;
if (E_Hourglass < 0.) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = -coefStabHourglass;};
// grandeurs finales avec le bon signe
res *= coefMultiplicatif;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * (X * (res));
// on met à jour le résidu de l'élément principal
(*residu) -= res;
// pour la partie raideur: on met à jour la raideur de l'élément principal
(*raideur) += coefMultiplicatif * (*raid_hourglass_transitoire);
//debug
//if (E_Hourglass < 0.)
// { D_X.Affiche();
// res.Affiche();
// raid_hourglass_transitoire->Affiche();
// cout << "\n E_Hourglass= "<< E_Hourglass << endl ;
//
// };
//res.Affiche(); (raid_hourglass_transitoire)->Affiche();
//fin debug
};
//---- debug
//cout << "\n raideur locale après stabilisation d'hourglass ";
//matRaideur.Affiche();
//---- fin debug
};
break;
case STABHOURGLASS_NON_DEFINIE :
// on ne fait rien
break;
default :
cout << "\n*** Erreur : cas de gestop, d'hourglass non defini !\n";
cout << "\n ElemMeca::Cal_implicit_hourglass() \n";
Sortie(1);
};
};
// stabilisation pour un calcul explicit
void ElemMeca::Cal_explicit_hourglass(bool atdt)
{ switch (type_stabHourglass)
{case STABHOURGLASS_PAR_COMPORTEMENT :
{ // --- calcul de la raideur de l'élément, dans le cas implicite ---
Vecteur* resHourglass = NULL;
if (atdt)
{resHourglass = tab_elHourglass(1)->CalculResidu_tdt (ParaGlob::param->ParaAlgoControleActifs());}
else
{resHourglass = tab_elHourglass(1)->CalculResidu_t (ParaGlob::param->ParaAlgoControleActifs());};
// on tiend compte du facteur de modération
(*resHourglass) *= coefStabHourglass;
// on met à jour le résidu de l'élément principal
(*residu) += (*resHourglass);
//---- debug
//cout << "\n raideur locale après stabilisation d'hourglass ";
//matRaideur.Affiche();
//---- fin debug
};
break;
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{ // --- on calcul la raideur de l'élément juste à la définition ---
if ((raid_hourglass_transitoire == NULL)||(prepa_niveau_precision > 18)) // pour le debug
{ // 1) calcul de la partie complète
Element::ResRaid resraid1 = tab_elHourglass(1)->Calcul_implicit (ParaGlob::param->ParaAlgoControleActifs());
// on sauvegarde la raideur initiale
if (raid_hourglass_transitoire == NULL)
raid_hourglass_transitoire = new Mat_pleine((*(resraid1.raid)));
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud.Taille();
int dim = ParaGlob::Dimension();
if(ParaGlob::AxiSymetrie())
dim--; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur& D_X = vec_hourglass(1); // pour simplier
Vecteur& X = vec_hourglass(2); // pour simplier
D_X.Change_taille(NBN*dim); // changement de la taille au cas où
X.Change_taille(NBN*dim); // changement de la taille au cas où
for (int ine=1;ine<=NBN;ine++)
{ Noeud& noe = *(tab_noeud(ine)); // pour simplifier
Coordonnee delta_X = noe.Coord2() - noe.Coord1();
Coordonnee Xnoeud = noe.Coord2();
int ipos = (ine-1)*dim;
switch (dim)
{ case 3: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);D_X(ipos+3)=delta_X(3);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);X(ipos+3)=Xnoeud(3);break;}
case 2: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);break;}
case 1: {D_X(ipos+1)=delta_X(1);
X(ipos+1)=Xnoeud(1);break;}
};
};
// on multiplie la matrice par le vecteur delta_X
(*(resraid1.res)) = resraid1.raid->Prod_mat_vec(D_X);
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * (D_X * (*(resraid1.res)));
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass;
if (E_Hourglass < 0.) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = -coefStabHourglass;};
// grandeurs finales avec le bon signe
(*(resraid1.raid)) *= coefMultiplicatif;
(*(resraid1.res)) *= coefMultiplicatif;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * (X * (*(resraid1.res)));
// on met à jour le résidu de l'élément principal
(*residu) -= (*(resraid1.res));
// calcul de l'énergie
E_Hourglass = 0.5 * (D_X * (*(resraid1.res)));
}
else // cas courant
{// on utilise la précédente raideur sauvegardée et on ne calcul qu'une partie résidue linéaire
// on commence par récupérer dans le vecteur résidu les incréments de ddl correspondant à l'élément
int NBN = tab_noeud.Taille();
int dim = ParaGlob::Dimension();
if(ParaGlob::AxiSymetrie())
dim--; // on passe à une dimension 2 ! et il y a uniquement les x et les y qui sont concernés
Vecteur& D_X = vec_hourglass(1); // pour simplier
Vecteur& X = vec_hourglass(2); // pour simplier
D_X.Change_taille(NBN*dim); // changement de la taille au cas où
X.Change_taille(NBN*dim); // changement de la taille au cas où
for (int ine=1;ine<=NBN;ine++)
{ Noeud& noe = *(tab_noeud(ine)); // pour simplifier
Coordonnee delta_X = noe.Coord2() - noe.Coord1();
Coordonnee Xnoeud = noe.Coord2();
int ipos = (ine-1)*dim;
switch (dim)
{ case 3: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);D_X(ipos+3)=delta_X(3);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);X(ipos+3)=Xnoeud(3);break;}
case 2: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);
X(ipos+1)=Xnoeud(1);X(ipos+2)=Xnoeud(2);break;}
case 1: {D_X(ipos+1)=delta_X(1);
X(ipos+1)=Xnoeud(1);break;}
};
};
Vecteur& res = vec_hourglass(2); // pour simplier
res.Change_taille(NBN*dim); // changement de la taille au cas où
// on multiplie le vecteur delta_X par la matrice
res = raid_hourglass_transitoire->Prod_mat_vec(D_X);
// premier calcul de l'accroissement de l'énergie
E_Hourglass = 0.5 * (D_X * res);
// on tiend compte du facteur de modération
double coefMultiplicatif = coefStabHourglass;
if (E_Hourglass < 0.) // on veut un accroissement d'énergie positive (1/2 k deltaX ^2)
{ coefMultiplicatif = -coefStabHourglass;};
// grandeurs finales avec le bon signe
res *= coefMultiplicatif;
// calcul définitif: correspond à une grandeur qui peut être comparée avec l'énergie élastique totale
E_Hourglass = 0.5 * (X * (res));
// on met à jour le résidu de l'élément principal
(*residu) -= res;
};
};
break;
case STABHOURGLASS_NON_DEFINIE :
// on ne fait rien
break;
default :
cout << "\n*** Erreur : cas de gestop, d'hourglass non defini !\n";
cout << "\n ElemMeca::Cal_implicit_hourglass() \n";
Sortie(1);
};
};
/*
// récupération de l'energie d'hourglass éventuelle
double ElemMeca::Energie_Hourglass()
{double enerHourglass = 0.;
switch (type_stabHourglass)
{case STABHOURGLASS_PAR_COMPORTEMENT :
{ // on récupère les énergies stockées à l'élément
const EnergieMeca& energieTotale = tab_elHourglass(1)->EnergieTotaleElement();
enerHourglass = coefStabHourglass * (energieTotale.EnergieElastique()
+ energieTotale.DissipationPlastique() + energieTotale.DissipationVisqueuse());
break;
}
case STABHOURGLASS_PAR_COMPORTEMENT_REDUIT :
{// --- l'énergie est purement une énergie élastique du type 0.5 <D_X> [K] (D_X) ---
if (raid_hourglass_transitoire == NULL)
{ enerHourglass = 0.;} // car rien n'a été calculé
else
{ int NBN = tab_noeud.Taille(); // récup du nombre de noeuds
int dim = ParaGlob::Dimension();
Vecteur& D_X = vec_hourglass(1); // pour simplier
D_X.Change_taille(NBN*dim); // changement de la taille au cas où
for (int ine=1;ine<=NBN;ine++)
{ Noeud& noe = *(tab_noeud(ine)); // pour simplifier
Coordonnee delta_X = noe.Coord2() - noe.Coord1();
int ipos = (ine-1)*dim;
switch (dim)
{ case 3: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);D_X(ipos+3)=delta_X(3);break;}
case 2: {D_X(ipos+1)=delta_X(1);D_X(ipos+2)=delta_X(2);break;}
case 1: {D_X(ipos+1)=delta_X(1);break;}
};
};
Vecteur& res = vec_hourglass(2); // pour simplier
res.Change_taille(NBN*dim); // changement de la taille au cas où
// on multiplie la matrice par le vecteur delta_X
res = raid_hourglass_transitoire->Prod_mat_vec(D_X);
enerHourglass = 0.5 * (D_X * res);
//debug
//cout << " , " << enerHourglass;
//fin debug
};
break;
case STABHOURGLASS_NON_DEFINIE :
// on ne fait rien
break;
}
};
return enerHourglass;
};*/
// -------------------- stabilisation transversale de membrane ou biel -------------------
// mise à jour de "a_calculer" en fonction du contexte
// NB: pour l'instant ne concerne que le calcul en implicite
void ElemMeca::Mise_a_jour_A_calculer_force_stab()
{if (pt_StabMembBiel != NULL)
{pt_StabMembBiel->Aa_calculer() = true; // on init par défaut
Enum_StabMembraneBiel enu_stab = pt_StabMembBiel->Type_stabMembrane();
switch (enu_stab)
{ case STABMEMBRANE_BIEL_PREMIER_ITER: case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER:
case STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD:
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD:
{// on récupère le pointeur correspondant aux iteration
const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
if (pointe == NULL)
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< "et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible"
<< " on ne peut pas continuer "
<< "\nElemMeca::Mise_a_jour_A_calculer_force_stab(..."<<endl;
Sortie(1);
};
// récup du conteneur du compteur
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
Grandeur_scalaire_entier& gr
= *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
int compteur_iter = *(gr.ConteneurEntier());
if (compteur_iter > 1) // cela veut dire que la force de stabilisation a déjà été calculé
{pt_StabMembBiel->Aa_calculer() = false;}
else // sinon il faut la calculer
{pt_StabMembBiel->Aa_calculer() = true;};
break;
}
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR: case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR:
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR_NORMALE_AU_NOEUD:
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR_NORMALE_AU_NOEUD:
{int compteur_iter = 0;// init
{// on récupère le pointeur correspondant aux iteration
const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
if (pointe == NULL)
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< "et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible"
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// récup du conteneur du compteur
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
Grandeur_scalaire_entier& gr
= *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
compteur_iter = *(gr.ConteneurEntier());
}
// on récupère le pointeur correspondant aux incréments
int compteur_incr = 0; // init
{const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL));
if (pointe == NULL)
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est STABMEMBRANE_BIEL_PREMIER_ITER_INCR "
<< "et la variable globale COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL n'est pas disponible"
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// récup du conteneur du compteur
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
Grandeur_scalaire_entier& gr
= *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
compteur_iter = *(gr.ConteneurEntier());
}
if ((compteur_incr > 1) && (compteur_iter > 1)) // cela veut dire que la force de stabilisation a déjà été calculé
{pt_StabMembBiel->Aa_calculer() = false;}
else // sinon il faut la calculer
{pt_StabMembBiel->Aa_calculer() = true;};
break;
}
default:
cout << "\n *** erreur le type de stabilisation de membrane et/ou biel n'est pas definie "
<< " on ne peut pas continuer "
<< "\nElemMeca::Mise_a_jour_A_calculer_force_stab(..."<<endl;
Sortie(1);
break;
};
};
};
// stabilisation pour un calcul implicit,
// iteg -> donne le numéro du dernier pti sur lequel on a travaillé, auquel met est associé
// iteg == 0 : un seul calcul global, et on est à la suite d'une boucle
// iteg == -1 : fin d'un calcul avec boucle sur tous les pti, et on est à la suite de la boucle
// correspond au calcul d'alpha: == stab_ref
// iteg entre 1 et nbint: on est dans une boucle de pti
// nbint : le nombre maxi de pti
// poids_volume : si iteg > 0 : le poids d'intégration
// si iteg <= 0 : l'intégrale de l'élément (ici la surface totale)
void ElemMeca::Cal_implicit_StabMembBiel(int iteg,const Met_abstraite::Impli& ex, const int& nbint
,const double& poids_volume
,Tableau <int> * noeud_a_prendre_en_compte)
{// l'idée est de bloquer la direction normale à l'élément pour une plaque et
// les deux directions transverses pour une biel
int niveau_affichage = ParaGlob::NiveauImpression();
if (pt_StabMembBiel->Affichage_stabilisation() > 0)
niveau_affichage = pt_StabMembBiel->Affichage_stabilisation();
//---- détermination de l'intensité de stabilisation
// celle-ci n'est calculé que dans le cas ou iteg == 0
// c'est à dire que l'on n'est pas dans une boucle sur les pti, dans la méthode appelant
double alpha = 0.; // init
bool stab_utilisable = false; // init
// bool alpha_vient_etre_calculee = false;
Enum_StabMembraneBiel enu_stab = pt_StabMembBiel->Type_stabMembrane();
switch (enu_stab)
{ case STABMEMBRANE_BIEL_PREMIER_ITER: case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER:
case STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD :
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD:
{// on récupère le pointeur correspondant aux iteration
const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
if (pointe == NULL)
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< "et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible"
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// récup du conteneur du compteur
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
Grandeur_scalaire_entier& gr
= *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
int compteur_iter = *(gr.ConteneurEntier());
if ((iteg > 0) && (compteur_iter > 1)) // cela veut dire que la la référence de stabilisation a déjà été calculée
{alpha = pt_StabMembBiel->Stab_ref(); // récup de l'intensité de stabilisation
// ce n'est pas un pointeur, alpha est une variable locale
// dans le cas d'une fonction multiplicatrice ,
// et il faut appeler la fonction qui peut varier pendant les itérations
if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// cas où il y a une fonction nD
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_gamma(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,iteg,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,iteg,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
alpha *= tab_ret(1) ;
pt_StabMembBiel->Valgamma() = alpha; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe
{alpha *= pt_StabMembBiel->Valgamma();}
stab_utilisable=true;
}
else if ((pt_StabMembBiel->Aa_calculer())
&& (iteg == 0)
) // cas ou il faut calculer la référence de stabilisation, mais on ne le fait qu'ici car
// on est sortie de la boucle sur des pti, de manière à utiliser la raideur totale
{// on récupère la valeur maxi de la raideur
int i,j; // def
double& maxival = pt_StabMembBiel->Stab_ref(); // init
// là il s'agit d'un pointeur donc on travaille directement sur la grandeur stockée
if ( (enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER)
||(enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD))
{maxival = (*raideur).MaxiValAbs(i,j) / poids_volume;
// alpha_vient_etre_calculee = true;
}
else if ( (enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER)
||(enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD))
{maxival = (*raideur).Maxi_ligne_ValAbs(i) / poids_volume;
// alpha_vient_etre_calculee = true;
}
else
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< " pas pris en compte pour l'instant .... "
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// arrivée ici Stab_ref() a une valeur
// maintenant on va calculer la valeur locale d'alpha
// on récupère la stabilisation pour lui donner une valeur
if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// cas où il y a une fonction nD
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_gamma(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// ici c'est la métrique correspondante au dernier pti
// comme iteg == 0 on est donc sortie de la boucle sur les pti, le dernier pti c'est nbint
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,nbint,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,nbint,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
alpha = tab_ret(1) * maxival;
pt_StabMembBiel->Valgamma() = alpha; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe ;
{alpha = pt_StabMembBiel->Valgamma() * maxival; // on prend une proportion du maxi
};
stab_utilisable = true;
}
else if (iteg == -1)// cas de la fin d'une boucle, on passe car il n'y a rien à faire
{if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
{alpha = pt_StabMembBiel->Valgamma();} // on récupère la dernière valeur d'alpha calculée
else // sinon il s'agit d'une valeur fixe, on recalcule alpha
{alpha = pt_StabMembBiel->Valgamma() * pt_StabMembBiel->Stab_ref();};
stab_utilisable = true;
}
else // ce n'est pas prévu
{ cout << "\n *** erreur, iteg = "<< iteg //<< " local= "<< local
<< " cas non prevu dans la stabilisation "
<< "\n ElemMeca::Cal_implicit_StabMembBiel(..";
Sortie(1);
};
break;
}
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR: case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR:
case STABMEMBRANE_BIEL_PREMIER_ITER_INCR_NORMALE_AU_NOEUD:
case STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR_NORMALE_AU_NOEUD:
{int compteur_iter = 0;// init
// on récupère le pointeur correspondant aux iteration
{const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
if (pointe == NULL)
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< "et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible"
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// récup du conteneur du compteur
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
Grandeur_scalaire_entier& gr
= *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
compteur_iter = *(gr.ConteneurEntier());
}
// on récupère le pointeur correspondant aux incréments
int compteur_incr = 0; // init
{const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL));
if (pointe == NULL)
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< "et la variable globale COMPTEUR_INCREMENT_CHARGE_ALGO_GLOBAL n'est pas disponible"
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// récup du conteneur du compteur
TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
Grandeur_scalaire_entier& gr
= *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
compteur_iter = *(gr.ConteneurEntier());
}
if ((iteg != 0) && (compteur_incr > 1) && (compteur_iter > 1)) // cela veut dire que la la référence de stabilisation a déjà été calculée
{alpha = pt_StabMembBiel->Stab_ref(); // init de l'intensité de stabilisation
// dans le cas d'une fonction, la partie liée à la raideur est stockée dans Valgamma()
// et il faut appeler la fonction qui peut varier pendant les itérations
if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// cas où il y a une fonction nD
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_gamma(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,iteg,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,iteg,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
alpha *= tab_ret(1);
pt_StabMembBiel->Valgamma() = alpha; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe
{alpha *= pt_StabMembBiel->Valgamma();}
stab_utilisable=true;
}
else if ((pt_StabMembBiel->Aa_calculer())
&& (iteg == 0)
) // cas ou il faut calculer la référence de stabilisation, mais on ne le fait qu'ici car
// on est sortie de la boucle sur des pti, de manière à utiliser la raideur totale
{// on récupère la valeur maxi de la raideur
int i,j; // def
double& maxival = pt_StabMembBiel->Stab_ref(); // init
if ( (enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER_INCR)
|| (enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER_INCR_NORMALE_AU_NOEUD))
{maxival = (*raideur).MaxiValAbs(i,j);
// alpha_vient_etre_calculee = true;
}
else if ( (enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR)
||(enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_INCR_NORMALE_AU_NOEUD))
{maxival = (*raideur).Maxi_ligne_ValAbs(i);
// alpha_vient_etre_calculee = true;
}
else
{ cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
<< " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
<< " pas pris en compte pour l'instant .... "
<< " on ne peut pas continuer "
<< "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
Sortie(1);
};
// on récupère la stabilisation
if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// cas où il y a une fonction nD
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_gamma(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// ici c'est la métrique correspondante au dernier pti
// comme iteg == 0 on est donc sortie de la boucle sur les pti, le dernier pti c'est nbint
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,nbint,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,nbint,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
alpha = tab_ret(1) * maxival;
pt_StabMembBiel->Valgamma() = alpha; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe ;
{alpha = pt_StabMembBiel->Valgamma() * maxival; // on prend une proportion du maxi
};
stab_utilisable = true;
}
else if (iteg == -1)// cas de la fin d'une boucle, on passe car il n'y a rien à faire
{if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
{alpha = pt_StabMembBiel->Valgamma();} // on récupère la dernière valeur d'alpha calculée
else // sinon il s'agit d'une valeur fixe, on recalcule alpha
{alpha = pt_StabMembBiel->Valgamma() * pt_StabMembBiel->Stab_ref();};
stab_utilisable = true;
}
else // ce n'est pas prévu
{ cout << "\n *** erreur, iteg = "<< iteg //<< " local= "<< local
<< " cas non prevu dans la stabilisation "
<< "\n ElemMeca::Cal_implicit_StabMembBiel(..";
Sortie(1);
};
break;
}
/*
// case STAB_M_B_ITER_1_VIA_F_EXT: case STAB_M_B_ITER_1_VIA_F_EXT_NORMALE_AU_NOEUD:
// {// on récupère le pointeur correspondant aux iteration
// const void* pointe = (ParaGlob::param->GrandeurGlobal(COMPTEUR_ITERATION_ALGO_GLOBAL));
// if (pointe == NULL)
// { cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
// << " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
// << "et la variable globale COMPTEUR_ITERATION_ALGO_GLOBAL n'est pas disponible"
// << " on ne peut pas continuer "
// << "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
// Sortie(1);
// };
// // récup du conteneur du compteur
// TypeQuelconque* gr_quelc = (TypeQuelconque*) (pointe);
// Grandeur_scalaire_entier& gr
// = *((Grandeur_scalaire_entier*) gr_quelc->Grandeur_pointee()); // pour simplifier
// int compteur_iter = *(gr.ConteneurEntier());
// if (compteur_iter > 1) // cela veut dire que la force de stabilisation a déjà été calculé
// {// dans le cas d'une fonction, la partie liée à la raiseur est stockée dans Valgamma()
// // et il faut appeler la fonction qui peut varier pendant les itérations
// if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// // cas où il y a une fonction nD
// { Tableau <double>& tab_val = pt_StabMembBiel->Pt_fct_gamma()->Valeur_pour_variables_globales();
// // seule la première valeur est ok
// alpha = tab_val(1) * pt_StabMembBiel->Valgamma();
// pt_StabMembBiel->FF_StabMembBiel() = alpha; // on sauvegarde
// }
// else // sinon cela veut dire que l'intensité de stabilisation est fixe
// {alpha = pt_StabMembBiel->FF_StabMembBiel();}
// alpha_utilisable=true;
// }
// else if (!local) // sinon il faut la calculer, mais on ne le fait que
// // si on est sortie de la boucle sur des pti,
// {// ***** à faire ***** on récupère la valeur maxi de la raideur
// { cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
// << " partie pas operationnelle pour l'instant !! "
// << " on ne peut pas continuer "
// << "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
// Sortie(1);
// };
// int i,j; // def
// double maxival = 0.; // init
// if ( (enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER)
// ||(enu_stab == STABMEMBRANE_BIEL_PREMIER_ITER_NORMALE_AU_NOEUD))
// {maxival = (*raideur).MaxiValAbs(i,j);
// alpha_vient_etre_calculee = true;
// }
// else if ( (enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER)
// ||(enu_stab == STABMEMBRANE_BIEL_Gerschgorin_PREMIER_ITER_NORMALE_AU_NOEUD))
// {maxival = (*raideur).Maxi_ligne_ValAbs(i);
// alpha_vient_etre_calculee = true;
// }
// else
// { cout << "\n *** pb dans la stabilisation de membrane et/ou biel "
// << " le type de stabilisation est " << Nom_StabMembraneBiel(enu_stab)
// << " pas pris en compte pour l'instant .... "
// << " on ne peut pas continuer "
// << "\nElemMeca::Cal_implicit_StabMembBiel(..."<<endl;
// Sortie(1);
// };
// // on récupère la stabilisation
// alpha=0.; // init
// if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// // cas où il y a une fonction nD
// { Tableau <double>& tab_val = pt_StabMembBiel->Pt_fct_gamma()->Valeur_pour_variables_globales();
// // seule la première valeur est ok
// alpha = tab_val(1) * maxival;
// pt_StabMembBiel->Valgamma() = maxival; // sauvegarde
// }
// else // sinon cela veut dire que l'intensité de stabilisation est fixe ;
// {alpha = pt_StabMembBiel->Valgamma() * maxival; // on prend une proportion du maxi
// };
// pt_StabMembBiel->FF_StabMembBiel() = alpha; // on sauvegarde
// alpha_utilisable = true;
// };
// break;
// }
*/
default:
cout << "\n *** erreur le type de stabilisation de membrane et/ou biel n'est pas definie "
<< " on ne peut pas continuer " << endl;
Sortie(1);
break;
};
// si la stabilisation est utilisable on s'en sert
if (stab_utilisable)
// if ( (alpha_utilisable)
// // et soit la stabilisation vient juste d'être calculée, donc on est en dehors de la boucle sur les pti
// // donc on calcule globalement l'impact sur la raideur
// // ou soit on est en local : et on calcul localement l'impact sur la raideur, mais à la sortie de la boucle
// // on ne recalculera pas en globa
// && (alpha_vient_etre_calculee || local)
// )
{// choix en fonction de la géométrie
Enum_type_geom enu_type_geom = Type_geom_generique(this->Id_geometrie());
switch (enu_type_geom)
{ case SURFACE:
//#ifdef MISE_AU_POINT
// if (num_elt==25) // -----debug à virer
// {cout << "\n ElemMeca::Cal_implicit_StabMembBiel(..";
// cout << "\n iteg= "<< iteg << flush;
// };
//#endif
if (!Contient_Normale_au_noeud(enu_stab))
// cas où on utilise la normale au pti
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
Coordonnee normale = Util::ProdVec_coor( (*ex.giB_tdt).Coordo(1), (*ex.giB_tdt).Coordo(2));
// calcul de la variation de la normale
// 1) variation du produit vectoriel
Tableau <Coordonnee > D_pasnormale =
Util::VarProdVect_coor( (*ex.giB_tdt).Coordo(1),(*ex.giB_tdt).Coordo(2),(*ex.d_giB_tdt));
// 2) de la normale
Tableau <Coordonnee> D_normale = Util::VarUnVect_coor(normale,D_pasnormale,normale.Norme());
// calcul de la normale normée et pondérée par la stabilisation
normale.Normer();
int nbddl = (*residu).Taille();
int dimf = 3; // ne fonctionne qu'en dim 3
if(ParaGlob::AxiSymetrie())
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf--;
int nbne = tab_noeud.Taille();
double ponder_F_t = 1.; // init
// si iteg == 0 ==> un seul calcul, global: correspond au calcul d'alpha: == stab_ref
// si teg == 1 ==> premier calcul correspondant au premier pti
if ((iteg == 0) || (iteg == 1))
{ matD->Initialise(0.);
resD->Inita(0.);
maxi_F_t = 0.;
// on initialise les forces et énergie
int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) = 0.;
pt_StabMembBiel->EE_StabMembBiel(i) = pt_StabMembBiel->EE_StabMembBiel_t(i);
}
};
// on récupère l'élément géométrique
ElemGeomC0& elemGeom = this->ElementGeometrique();
// récup des fonctions d'interpolation
const Tableau <Vecteur>& taphi = elemGeom.TaPhi();
int nb_noeud_pris_en_compte = 0; //pour faire la moyenne dans le cas iteg == 0
// on applique une force de stabilisation
// uniquement suivant la direction de la normale
int num_pti=iteg; // init
if (iteg <= 0) num_pti = nbint;
// la force est proportionelle au déplacement suivant la normale
if (iteg > -1) // c-a-d == 0 ou plus
{ int ix=1;
for (int ne =1; ne<= nbne; ne++)
{bool continuer = true;
if (noeud_a_prendre_en_compte != NULL)
if (!(*noeud_a_prendre_en_compte)(ne))
continuer = false;
if (continuer)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
nb_noeud_pris_en_compte++;
// on va bloquer le déplacement dans la direction de la normale calculée au pti
//dans un calcul implicit on utilise la position à l'itération 0
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu(XI_ITER_0);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel = noe.Grandeur_quelconque(enu);
const Grandeur_coordonnee& cofo = *((Grandeur_coordonnee*) typquel.Const_Grandeur_pointee());
Coordonnee delta_X = noe.Coord2() - cofo.ConteneurCoordonnee_const();
double normal_fois_delta_X = normale * delta_X; // le déplacement au noeud suivant la normale au pti
// dans le cas où on veut pondérer F_t, on calcule la fonction
if (pt_StabMembBiel->Pt_fct_ponder_Ft() != NULL)
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_ponder_Ft(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi
(absolue,TEMPS_tdt,li_enu_scal,num_pti,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,num_pti,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
ponder_F_t = tab_ret(1) ;
}
// un - car la force doit s'opposer au déplacement
double intensite = -alpha * normal_fois_delta_X
+ ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX(maxi_F_t, pt_StabMembBiel->FF_StabMembBiel_t(num_pti));
// d'où une variation réelle d'intensité
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
pt_StabMembBiel->EE_StabMembBiel(num_pti) += 0.5 * delta_intensite * normal_fois_delta_X;
#ifdef MISE_AU_POINT
if (niveau_affichage>5)
// if (num_elt==25) // -----debug à virer
{cout << "\n ElemMeca::Cal_implicit_StabMembBiel(..";
cout << "\n noeud ("<<ne<<"), alpha= "<<alpha << " intensite= "<< intensite
<< " delta_intensite= " << delta_intensite
<< " delta_X= "<< delta_X
<< "\n normale= "<<normale << " normal_fois_delta_X= " << normal_fois_delta_X << flush;
};
#endif
// //----- debug -----
// if (num_elt==25) // -----debug à virer
// {int nb_pti = pt_StabMembBiel->Taille();
// for (int i=1;i<=nb_pti;i++)
// {cout << "\n avant update F_stab_t("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel_t(i)<< " ";
// cout << " F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
// cout << " ponder_F_t= " << ponder_F_t;
// }
// };
// //--- fin debug ----
// calcul de la raideur: on traite les cas :
// iteg >0 et iteg == 0
// le cas iteg == -1 n'est pas à prendre en compte car il correspond
// au cas où on vient de balayer tous les pti
if (iteg > 0) // cas où on est dans une boucle d'intégration
{// un premier coef pour optimiser
double phi_jacob_poids = taphi(iteg)(ne) * (*ex.jacobien_tdt) * poids_volume ;
pt_StabMembBiel->FF_StabMembBiel(iteg) += intensite * taphi(iteg)(ne); // * phi_jacob_poids;
// taphi(iteg)(ne) * intensite * (*ex.jacobien_tdt) * poids_volume ;
// NB: 1) comme on multiplie par taphi(iteg)(ne), et que l'on boucle sur les noeuds (ne = 1 à nbne)
// au sortir de la boucle on a l'intensité totale qui est due au déplacement de tous les noeuds
// 2) si alpha * normal_fois_delta_X = 0 c-a-d pas d'accroissement d'intensité, on va obtenir
// somme_noeud(ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(num_pti) * taphi(iteg)(ne))
// et comme somme_noeud(taphi(iteg)(ne))=1, on devrait obtenir une intensité globale identique
// à celle à t
pt_StabMembBiel->EE_StabMembBiel(iteg) += 0.5 * intensite * normal_fois_delta_X * phi_jacob_poids;
// 0.5 * delta_intensite * normal_fois_delta_X * taphi(iteg)(ne) * (*ex.jacobien_tdt) * poids_volume ;
// un second coef pour optimiser
double phi_poids = taphi(iteg)(ne) * poids_volume ;
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for (int i=1;i<= dimf;i++,ix++)
{ // = (- normale(i) * alpha * normal_fois_delta_X + F_t)
// * (taphi(iteg)(ne) * poids * jacobien) ;
(*resD)(ix) += normale(i) * intensite * phi_jacob_poids;
// taphi(iteg)(ne)* normale(i) * intensite
// * (*ex.jacobien_tdt) * poids_volume ; // ok
// contribution à la raideur: de signe inverse à celui du résidu
int j=1;
int me = 1; // numéro de noeud dans la direction des ddl
for (int jx = 1; jx <= nbddl; jx++,j++)
{if (j>dimf) {j=1;me++;};
(*matD)(ix,jx) += phi_poids // == taphi(iteg)(ne) * poids_volume
*( (*ex.jacobien_tdt) *
( D_normale(jx)(i) * alpha * normal_fois_delta_X
+ normale(i) * alpha * (D_normale(jx) * delta_X)
)
+ (*ex.d_jacobien_tdt)(jx) * normale(i) * (-intensite)
);
if (ne == me)
{(*matD)(ix,jx) += phi_poids * normale(i) * alpha * normale(j) * (*ex.jacobien_tdt);
// taphi(iteg)(ne)* alpha * normale(i) * normale(j) * poids_volume;
};
};
};
}
else if (iteg == 0)// on est hors boucle poids_volume == toute la surface
{// pour le cas particulier d'iteg == 0, on n'a qu'un calcul global
// et on affecte la même intensité à tous les pti
// du coup on sauvegarde uniquement pour le dernier
// ensuite après on répartira sur l'ensemble de pti
// --> en fait on cumulle l'influence de tous les normales
// il faudra ensuite prendre la moyenne
pt_StabMembBiel->FF_StabMembBiel(num_pti) += intensite ; //* poids_volume;
// taphi(iteg)(ne) * intensite * poids_volume ;
pt_StabMembBiel->EE_StabMembBiel(num_pti) += 0.5 * intensite * normal_fois_delta_X * poids_volume;
// 0.5 * intensite * normal_fois_delta_X * taphi(iteg)(ne) * poids_volume ;
for (int i=1;i<= dimf;i++,ix++)
{ // - normale(i) * alpha * normal_fois_delta_X;
(*resD)(ix) += normale(i) * intensite * poids_volume ; // ok
// contribution à la raideur: de signe inverse à celui du résidu
int j=1;
int me = 1; // numéro de noeud dans la direction des ddl
for (int jx = 1; jx <= nbddl; jx++,j++)
{if (j>dimf) {j=1;me++;};
(*matD)(ix,jx) += poids_volume
*(( D_normale(jx)(i) * alpha * normal_fois_delta_X
+ normale(i) * alpha * (D_normale(jx) * delta_X)
)
);
if (ne == me)
{(*matD)(ix,jx) += poids_volume * normale(i) * alpha * normale(j);
};
};
};
};
}; // fin du test sur les noeuds à considérer
}; // fin de la boucle sur les noeuds
};
// dans le cas particulier de iteg == 0, c-a-d correspond au calcul d'alpha: == stab_ref
// on n'a pas bouclé sur les pti, c'est uniquement un calcul global qui a été effectué
if (iteg == 0)
{ // on commence par moyenner: se rappeller que l'on a bouclé aussi sur les noeuds !
int taille = pt_StabMembBiel->Taille();
pt_StabMembBiel->FF_StabMembBiel(num_pti) /= (nb_noeud_pris_en_compte);
pt_StabMembBiel->EE_StabMembBiel(num_pti) /= (nb_noeud_pris_en_compte);
// ensuite on répartie la même moyenne sur tous les pti
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) = pt_StabMembBiel->FF_StabMembBiel(num_pti);
pt_StabMembBiel->EE_StabMembBiel(i) = pt_StabMembBiel->EE_StabMembBiel(num_pti);
};
// également sur le résidu et la raideur, car due à la boucle sur les noeuds
if (nb_noeud_pris_en_compte > 0)
{ double inverse = 1./nb_noeud_pris_en_compte;
(*matD) *= inverse;
(*resD) *= inverse;
};
};
////----- debug -----
// if (num_elt==25) // -----debug à virer
// {int nb_pti = pt_StabMembBiel->Taille();
// for (int i=1;i<=nb_pti;i++)
// cout << "\n final F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
// };
////--- fin debug ----
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
if ((iteg == 0) || (iteg == -1)) // signifie que l'on est en dehors de la boucle
{
if (pt_StabMembBiel->Gestion_maxi_mini())
{
#ifdef MISE_AU_POINT
if (niveau_affichage>9)
// if (num_elt==1) // -----debug à virer
{cout << "\n residu de stabilisation avant limitation: "<< (*resD)
<< "\n residu initial (sans stabilisation): "<< (*residu)
<< "\n maxi_F_t= " << maxi_F_t;
};
if (niveau_affichage > 5)
// if (num_elt==25) // -----debug à virer
{int nb_pti = pt_StabMembBiel->Taille();
for (int i=1;i<=nb_pti;i++)
cout << "\n F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
};
#endif
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false;
double coef=1.; // coef de limitation éventuelle si limitation devient true
// récup des limitations éventuelles
double beta = pt_StabMembBiel->Beta();
double F_mini = pt_StabMembBiel->F_mini();
double d_maxi = pt_StabMembBiel->D_maxi();
// on fait une boucle sur les noeuds
int ix=1;
for (int ne =1; ne<= nbne; ne++)
{bool continuer = true;
if (noeud_a_prendre_en_compte != NULL)
if (!(*noeud_a_prendre_en_compte)(ne))
continuer = false;
if (continuer)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab(dimf); // init
for (int i=1;i<= dimf;i++,ix++)
force_stab(i) = (*resD)(ix);
double intensite_generalise = force_stab.Norme();
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext(FORCE_GENE_EXT_t);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel_ext = noe.Grandeur_quelconque(enu_ext);
const Grandeur_coordonnee& cofo_ext = *((Grandeur_coordonnee*) typquel_ext.Const_Grandeur_pointee());
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs(cofo_ext.ConteneurCoordonnee_const() * normale);
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext;
double max_intens_via_Fext = maxi_intensite; // on sauve pour le cas où on est dep max
// on applique l'algo théorique cf. doc théorique
// erreur, maxi_intensité est une limitation, mais ce que l'on cherche à voir c'est
// est-ce que les forces externes sont non-nulles, donc c'est intensite_Fext qu'il faut tester
// if (maxi_intensite >= F_mini)
if (intensite_Fext >= F_mini)
// on est dans le cas où les forces externes sont à prendre en compte
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intensite et intensite très petites
)
{limitation = true;
if (niveau_affichage>3)
{cout << "\n limitation stabilisation due aux forces externes, F_stab gene= "<< intensite_generalise
<< " ==> " << maxi_intensite;
if (niveau_affichage > 5)
cout << "\n F_mini= "<< F_mini << " intensite_Fext= "<< intensite_Fext;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
// sinon pas de pb donc on ne fait rien
}
else // cas où l'intensité dans la direction normale, des forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{maxi_intensite = Dabs(-alpha*d_maxi + maxi_F_t);
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intensite et intensite très petites
&& (Dabs(alpha) > ConstMath::pasmalpetit) // il faut qu'alpha ne soit pas nul
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due deplace maxi acceptable, F_stab gene= "<< intensite_generalise
<< " ==> " << maxi_intensite
<< " (max_intens_via_Fext= "<<max_intens_via_Fext<<")" ;
if ((niveau_affichage > 4))
cout << "\n maxi_F_t= "<< maxi_F_t << " alpha*d_maxi= "
<< (alpha*d_maxi)
<< " alpha= " << alpha
<< " d_maxi= " << d_maxi;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
};
};
}; // fin du test sur les noeuds à considérer
};// fin de la boucle sur les noeuds
if (limitation)
{ (*resD) *= coef; // mise à jour de la force de stabilisation
// idem au niveau des pti
int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) *=coef;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel->EE_StabMembBiel(i) - pt_StabMembBiel->EE_StabMembBiel_t(i);
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef;
//ce qui donne au final:
pt_StabMembBiel->EE_StabMembBiel(i) =
delta_energie + pt_StabMembBiel->EE_StabMembBiel_t(i);
};
};
};
// 2) ---- on met à jour la raideur et résidu global
// en tenant compte de la stabilisation
(*residu) += (*resD);
(*raideur) += (*matD);
//et du coup on ne boucle pas sur les pti et on peut limiter la force sauf qu'une fois le calcul des forces
//aux noeuds fait, il faut ramener par interpolation, la force au pti pour la sauvegarder... bordel
//non... on sauvegarde dans un tableau indicé par les num de noeud, et seulement pour le post-traitement
//on calcul la valeur aux pti
//
#ifdef MISE_AU_POINT
if ((niveau_affichage > 9))
// if (num_elt==1) // -----debug à virer
{cout << "\n residu de stabilisation: "<< (*resD);
cout << "\n raideur de stabilisation: "; matD->Affiche();
cout << "\n fin ElemMeca::Cal_implicit_StabMembBiel(..";
};
#endif
}; // fin du cas iteg == 0 ou -1
} // fin du cas où on utilise la normale aux pti
else
{if (iteg < 1)
// on ne calcul que si iteg < 1, c-a-d qu'on est sortie
// de la boucle sur les pti
// sinon cas où on utilise la normale aux noeuds
{int nbddl = (*residu).Taille();
int dimf = 3; // ne fonctionne qu'en dim 3
if(ParaGlob::AxiSymetrie())
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf--;
int nbne = tab_noeud.Taille();
matD->Initialise(0.);
resD->Inita(0.);
maxi_F_t = 0.;
#ifdef MISE_AU_POINT
// on vérifie le bon dimensionnement
int taille = pt_StabMembBiel->Taille();
if (taille != nbne)
{cout << "\n erreur** ElemMeca::Cal_implicit_StabMembBiel(.."
<< " le stockage est mal dimensionne: taille = " << taille
<< " alors qu'il devrait etre : nbe = " << nbne << flush ;
Sortie(1);
};
#endif
// on applique une force de stabilisation aux noeuds
// uniquement suivant la direction de la normale
// la force est proportionelle au déplacement suivant la normale au noeud
int ix=1;
for (int ne =1; ne<= nbne; ne++)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// on initialise les forces et énergie
pt_StabMembBiel->FF_StabMembBiel(ne) = 0.;
pt_StabMembBiel->EE_StabMembBiel(ne) = pt_StabMembBiel->EE_StabMembBiel_t(ne);
// pour chaque noeud on va bloquer le déplacement dans la direction de la normale au noeud à t
// on récupère la normale au noeud
const TypeQuelconque& tiq = noe.Grandeur_quelconque(NN_SURF_t);
const Grandeur_coordonnee& gr= *((const Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
const Coordonnee& normale = gr.ConteneurCoordonnee_const();
//dans un calcul implicit on utilise la position à l'itération 0
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu(XI_ITER_0);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel = noe.Grandeur_quelconque(enu);
const Grandeur_coordonnee& cofo = *((Grandeur_coordonnee*) typquel.Const_Grandeur_pointee());
Coordonnee delta_X = noe.Coord2() - cofo.ConteneurCoordonnee_const();
double normal_fois_delta_X = normale * delta_X;
double ponder_F_t = 1.; // init
// dans le cas où on veut pondérer F_t, on calcule la fonction
if (pt_StabMembBiel->Pt_fct_ponder_Ft() != NULL)
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_ponder_Ft(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// au dernier pti !
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi
(absolue,TEMPS_tdt,li_enu_scal,nbint,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,nbint,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
ponder_F_t = tab_ret(1) ;
}
// un - car la force doit s'opposer au déplacement
double intensite = -alpha * normal_fois_delta_X
+ ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(ne);
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX(maxi_F_t, pt_StabMembBiel->FF_StabMembBiel_t(ne));
// d'où une variation réelle d'intensité
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(ne);
pt_StabMembBiel->EE_StabMembBiel(ne) += 0.5 * delta_intensite * normal_fois_delta_X;
#ifdef MISE_AU_POINT
if ((niveau_affichage > 5))
{cout << "\n ElemMeca::Cal_implicit_StabMembBiel(..";
cout << "\n alpha= "<<alpha << " intensite= "<< intensite
<< " delta_intensite= " << delta_intensite
<< " delta_X= "<< delta_X
<< "\n normale= "<<normale << " normal_fois_delta_X= " << normal_fois_delta_X << flush;
};
#endif
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for (int i=1;i<= dimf;i++,ix++)
{ // - normale(i) * stabili * normal_fois_delta_X;
(*resD)(ix) += normale(i) * intensite ;
// contribution à la raideur: de signe inverse à celui du résidu
int j=1;
int me = 1; // numéro de noeud dans la direction des ddl
for (int jx = 1; jx <= nbddl; jx++,j++)
{if (j>dimf) {j=1;me++;};
if (ne == me)
{(*matD)(ix,jx) += alpha * normale(i) * normale(j);
};
};
};
}; // fin de la boucle sur les noeuds
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
if (pt_StabMembBiel->Gestion_maxi_mini())
{
#ifdef MISE_AU_POINT
if ((niveau_affichage > 9))
{cout << "\n residu de stabilisation avant limitation: "<< (*resD)
<< "\n residu initial (sans stabilisation): "<< (*residu)
<< "\n maxi_F_t= " << maxi_F_t;
};
if ((niveau_affichage > 5))
{int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
cout << "\n F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
};
#endif
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false;
double coef=1.; // coef de limitation éventuelle si limitation devient true
// récup des limitations éventuelles
double beta = pt_StabMembBiel->Beta();
double F_mini = pt_StabMembBiel->F_mini();
double d_maxi = pt_StabMembBiel->D_maxi();
// on refait une boucle sur les noeuds
int ix=1;
for (int ne =1; ne<= nbne; ne++)
{bool continuer = true;
if (noeud_a_prendre_en_compte != NULL)
if (!(*noeud_a_prendre_en_compte)(ne))
continuer = false;
if (continuer)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab(dimf); // init
for (int i=1;i<= dimf;i++,ix++)
force_stab(i) = (*resD)(ix);
double intensite_generalise = force_stab.Norme();
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext(FORCE_GENE_EXT_t);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel_ext = noe.Grandeur_quelconque(enu_ext);
const Grandeur_coordonnee& cofo_ext = *((Grandeur_coordonnee*) typquel_ext.Const_Grandeur_pointee());
// on récupère la normale au noeud
const TypeQuelconque& tiq = noe.Grandeur_quelconque(NN_SURF_t);
const Grandeur_coordonnee& gr= *((const Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
const Coordonnee& normale = gr.ConteneurCoordonnee_const();
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs(cofo_ext.ConteneurCoordonnee_const() * normale);
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext;
double max_intens_via_Fext = maxi_intensite; // on sauve pour le cas où on est dep max
// on applique l'algo théorique cf. doc théorique
if (intensite_Fext >= F_mini)
// on est dans le cas où les forces externes existent
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intintensite très petites
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due aux forces externes, F_stab gene= "<< intensite_generalise
<< " ==> " << maxi_intensite;
if ((niveau_affichage>4)
|| (pt_StabMembBiel->Affichage_stabilisation() > 4))
cout << "\n F_mini= "<< F_mini << " intensite_Fext= "<< intensite_Fext;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
// sinon pas de pb donc on ne fait rien
}
else // cas où les forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{maxi_intensite = Dabs(-alpha*d_maxi + maxi_F_t);
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intintensite très petites
&& (Dabs(alpha) > ConstMath::pasmalpetit) // il faut qu'alpha ne soit pas nul
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due deplace maxi acceptable, F_s"<< intensite_generalise
<< " ==> " << maxi_intensite
<< " (max_intens_via_Fext= "<<max_intens_via_Fext<<")" ;
if ((niveau_affichage > 4))
cout << "\n maxi_F_t= "<< maxi_F_t << " alpha*d_maxi= "
<< (alpha*d_maxi)
<< " alpha= " << alpha
<< " d_maxi= " << d_maxi;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
};
};
}; // fin du test sur les noeuds à considérer
};// fin de la boucle sur les noeuds
if (limitation)
{ (*resD) *= coef; // mise à jour de la force de stabilisation
// idem au niveau des noeuds
int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) *=coef;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel->EE_StabMembBiel(i) - pt_StabMembBiel->EE_StabMembBiel_t(i);
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef;
//ce qui donne au final:
pt_StabMembBiel->EE_StabMembBiel(i) =
delta_energie + pt_StabMembBiel->EE_StabMembBiel_t(i);
};
// 2) ---- on met à jour la raideur et résidu global
// en tenant compte de la stabilisation
(*residu) += (*resD);
(*raideur) += (*matD);
}; // fin du cas if (limitation)
}; // fin du cas ou on applique éventuellement une limitation
#ifdef MISE_AU_POINT
if ((niveau_affichage > 9))
{cout << "\n residu de stabilisation: "<< (*resD);
cout << "\n raideur de stabilisation: "; matD->Affiche();
cout << "\n fin debug ElemMeca::Cal_implicit_StabMembBiel(..";
};
#endif
} // fin du cas où on utilise la normale aux noeuds
};
break;
case LIGNE:
break;
default: // on ne fait rien pour l'instant pour les autres types
break;
};
};
};
// stabilisation pour un calcul explicit
void ElemMeca::Cal_explicit_StabMembBiel(int iteg,const Met_abstraite::Expli_t_tdt ex
, const int& nbint,const double& poids_volume
,Tableau <int> * noeud_a_prendre_en_compte)
{
// la stabilisation en explicite n'est possible que s'il existe soit:
// - une force de stabilisation déjà enregistrée
// - un coefficient alpha non nulle que l'on peut calculer avec les données déjà enregistrées
// l'idée est de bloquer la direction normale à l'élément pour une plaque et
// les deux directions transverses pour une biel
int niveau_affichage = ParaGlob::NiveauImpression();
if (pt_StabMembBiel->Affichage_stabilisation() > 0)
niveau_affichage = pt_StabMembBiel->Affichage_stabilisation();
//---- détermination de l'intensité de stabilisation
// ici on ne s'occupe pas de l'itération ni de l'incrément, dans tous les cas
// on utilise les grandeurs déjà stockées
// ceci quelque soit le type de stabilisation
double alpha = 0.; // init
bool stab_utilisable = false; // init
Enum_StabMembraneBiel enu_stab = pt_StabMembBiel->Type_stabMembrane();
if (pt_StabMembBiel != NULL)
if (pt_StabMembBiel->Activite_en_explicite())
{ //alpha = pt_StabMembBiel->Stab_ref(); // récup de l'intensité de stabilisation
// ce n'est pas un pointeur, alpha est une variable locale
// stab_ref, par défaut == 1, mais dans le cas d'un passage
// implicit -> explicite : vaut la dernière valeur implicite calculée
// cela permet de profiter du passage en implicite, en récupérant l'intensité
// de raideur réelle
alpha = pt_StabMembBiel->Stab_ref(); // init
// dans le cas d'une fonction multiplicatrice ,
// et il faut appeler la fonction qui peut varier
if (pt_StabMembBiel->Pt_fct_gamma() != NULL)
// cas où il y a une fonction nD
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_gamma(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
int num_iteg_a_considerer; // init
if (iteg > 0) // on est dans la boucle
num_iteg_a_considerer = iteg;
else // sinon on est à la sortie de la boucle
num_iteg_a_considerer = nbint;
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi
(absolue,TEMPS_tdt,li_enu_scal,num_iteg_a_considerer,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,num_iteg_a_considerer,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
alpha *= tab_ret(1) ;
pt_StabMembBiel->Valgamma() = alpha; // on sauvegarde la valeur d'alpha
}
else // sinon cela veut dire que l'intensité de stabilisation est fixe
// en explicite c'est directement la valeur lue
{alpha *= pt_StabMembBiel->Valgamma() ;}
stab_utilisable=true;
};
// si la stabilisation est utilisable on s'en sert
if (stab_utilisable)
{// choix en fonction de la géométrie
Enum_type_geom enu_type_geom = Type_geom_generique(this->Id_geometrie());
switch (enu_type_geom)
{ case SURFACE:
if (!Contient_Normale_au_noeud(enu_stab))
// --- cas où on utilise la normale au pti
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
Coordonnee normale = Util::ProdVec_coor( (*ex.giB_tdt).Coordo(1), (*ex.giB_tdt).Coordo(2));
// calcul de la normale normée et pondérée par la stabilisation
normale.Normer();
int nbddl = (*residu).Taille();
int dimf = 3; // ne fonctionne qu'en dim 3
if(ParaGlob::AxiSymetrie())
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf--;
int nbne = tab_noeud.Taille();
double ponder_F_t = 1.; // init
// si iteg == 1 ==> premier calcul correspondant au premier pti
if (iteg == 1)
{ resD->Inita(0.);
maxi_F_t = 0.;
// on initialise les forces et énergie
int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) = 0.;
pt_StabMembBiel->EE_StabMembBiel(i) = pt_StabMembBiel->EE_StabMembBiel_t(i);
}
};
// on récupère l'élément géométrique
ElemGeomC0& elemGeom = this->ElementGeometrique();
// récup des fonctions d'interpolation
const Tableau <Vecteur>& taphi = elemGeom.TaPhi();
// on applique une force de stabilisation
// uniquement suivant la direction de la normale
// ici iteg peut soit == -1 : on est sortie de la boucle sur les pti
// soit > 0 : on et dans la boucle des pti
int num_pti=iteg; // init
if (iteg == -1) num_pti = nbint;
// la force est proportionelle au déplacement suivant la normale
if (iteg > 0)
{ int ix=1;
for (int ne =1; ne<= nbne; ne++)
{bool continuer = true;
if (noeud_a_prendre_en_compte != NULL)
if (!(*noeud_a_prendre_en_compte)(ne))
continuer = false;
if (continuer)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// pour chaque noeud on va rectifier la normale par rapport à celle calculée au pti
// on va bloquer le déplacement dans la direction de la normale au noeud à t
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu(XI_ITER_0);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel = noe.Grandeur_quelconque(enu);
const Grandeur_coordonnee& cofo = *((Grandeur_coordonnee*) typquel.Const_Grandeur_pointee());
Coordonnee delta_X = noe.Coord2() - cofo.ConteneurCoordonnee_const();
double normal_fois_delta_X = normale * delta_X; // le déplacement suivant la normale
// dans le cas où on veut pondérer F_t, on calcule la fonction
if (pt_StabMembBiel->Pt_fct_ponder_Ft() != NULL)
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_ponder_Ft(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi
(absolue,TEMPS_tdt,li_enu_scal,num_pti,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,num_pti,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
ponder_F_t = tab_ret(1) ;
}
// un - car la force doit s'opposer au déplacement
double intensite = -alpha * normal_fois_delta_X
+ ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX(maxi_F_t, pt_StabMembBiel->FF_StabMembBiel_t(num_pti));
// d'où une variation réelle d'intensité
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
pt_StabMembBiel->EE_StabMembBiel(num_pti) += 0.5 * delta_intensite * normal_fois_delta_X;
#ifdef MISE_AU_POINT
if ((niveau_affichage > 5))
{cout << "\n ElemMeca::Cal_explicit_StabMembBiel(..";
cout << "\n alpha= "<<alpha << " intensite= "<< intensite
<< " delta_intensite= " << delta_intensite
<< " delta_X= "<< delta_X
<< "\n normale= "<<normale << " normal_fois_delta_X= " << normal_fois_delta_X << flush;
};
#endif
// calcul de la raideur: on traite les cas :
// iteg >0 et iteg == 0
// le cas iteg == -1 n'est pas à prendre en compte car il correspond
// au cas où on vient de balayer tous les pti
if (iteg > 0) // cas où on est dans une boucle d'intégration
{// un premier coef pour optimiser
double phi_jacob_poids = taphi(iteg)(ne) * (*ex.jacobien_tdt) * poids_volume ;
pt_StabMembBiel->FF_StabMembBiel(iteg) += intensite * taphi(iteg)(ne); // * phi_jacob_poids;
// taphi(iteg)(ne) * intensite * (*ex.jacobien_tdt) * poids_volume ;
pt_StabMembBiel->EE_StabMembBiel(iteg) += 0.5 * intensite * normal_fois_delta_X * phi_jacob_poids;
// 0.5 * delta_intensite * normal_fois_delta_X * taphi(iteg)(ne) * (*ex.jacobien_tdt) * poids_volume ;
// un second coef pour optimiser
// double phi_poids = taphi(iteg)(ne) * poids_volume ;
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for (int i=1;i<= dimf;i++,ix++)
{ // = (- normale(i) * alpha * normal_fois_delta_X + F_t)
// * (taphi(iteg)(ne) * poids * jacobien) ;
(*resD)(ix) += normale(i) * intensite * phi_jacob_poids;
// taphi(iteg)(ne)* normale(i) * intensite
// * (*ex.jacobien_tdt) * poids_volume ; // ok
};
}
}; // fin du test sur les noeuds à considérer
}; // fin de la boucle sur les noeuds
};
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
// ici seul le cas iteg == -1 est pris en compte contrairement à l'implicite
if (iteg == -1) // signifie que l'on est en dehors de la boucle
{
if (pt_StabMembBiel->Gestion_maxi_mini())
{
#ifdef MISE_AU_POINT
if ((niveau_affichage > 5))
{cout << "\n residu de stabilisation avant limitation: "<< (*resD)
<< "\n residu initial (sans stabilisation): "<< (*residu)
<< "\n maxi_F_t= " << maxi_F_t;
};
if ((niveau_affichage > 5))
{int nb_pti = pt_StabMembBiel->Taille();
for (int i=1;i<=nb_pti;i++)
cout << "\n F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
};
#endif
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false;
double coef=1.; // coef de limitation éventuelle si limitation devient true
// récup des limitations éventuelles
double beta = pt_StabMembBiel->Beta();
double F_mini = pt_StabMembBiel->F_mini();
double d_maxi = pt_StabMembBiel->D_maxi();
// on refait une boucle sur les pti
int ix=1;
for (int ne =1; ne<= nbne; ne++)
{bool continuer = true;
if (noeud_a_prendre_en_compte != NULL)
if (!(*noeud_a_prendre_en_compte)(ne))
continuer = false;
if (continuer)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab(dimf); // init
for (int i=1;i<= dimf;i++,ix++)
force_stab(i) = (*resD)(ix);
double intensite_generalise = force_stab.Norme();
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext(FORCE_GENE_EXT_t);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel_ext = noe.Grandeur_quelconque(enu_ext);
const Grandeur_coordonnee& cofo_ext = *((Grandeur_coordonnee*) typquel_ext.Const_Grandeur_pointee());
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs(cofo_ext.ConteneurCoordonnee_const() * normale);
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext;
double max_intens_via_Fext = maxi_intensite; // on sauve pour le cas où on est dep max
// on applique l'algo théorique cf. doc théorique
if (intensite_Fext >= F_mini)
// on est dans le cas où les forces externes existent
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intensite et intensite très petites
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due aux forces externes, F_stab gene= "<< intensite_generalise
<< " ==> " << maxi_intensite;
if ((niveau_affichage > 5))
cout << "\n F_mini= "<< F_mini << " intensite_Fext= "<< intensite_Fext;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
// sinon pas de pb donc on ne fait rien
}
else // cas où les forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{maxi_intensite = Dabs(-alpha*d_maxi + maxi_F_t);
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intensite et intensite très petites
&& (Dabs(alpha) > ConstMath::pasmalpetit) // il faut qu'alpha ne soit pas nul
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due deplace maxi acceptable, F_stab gene= "<< intensite_generalise
<< " ==> " << maxi_intensite
<< " (max_intens_via_Fext= "<<max_intens_via_Fext<<")" ;
if ((niveau_affichage > 4))
cout << "\n maxi_F_t= "<< maxi_F_t << " alpha*d_maxi= "
<< (alpha*d_maxi)
<< " alpha= " << alpha
<< " d_maxi= " << d_maxi;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
};
};
}; // fin du test sur les noeuds à considérer
};// fin de la boucle sur les noeuds
if (limitation)
{ (*resD) *= coef; // mise à jour de la force de stabilisation
// idem au niveau des pti
int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) *=coef;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel->EE_StabMembBiel(i) - pt_StabMembBiel->EE_StabMembBiel_t(i);
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef;
//ce qui donne au final:
pt_StabMembBiel->EE_StabMembBiel(i) =
delta_energie + pt_StabMembBiel->EE_StabMembBiel_t(i);
};
};
}; // fin mise en place d'une limitation éventuelle
// 2) ---- on met à jour le résidu global
// en tenant compte de la stabilisation
(*residu) += (*resD);
#ifdef MISE_AU_POINT
if ((niveau_affichage > 9))
{cout << "\n residu de stabilisation: "<< (*resD);
cout << "\n fin ElemMeca::Cal_explicit_StabMembBiel(..";
};
#endif
}; // fin du cas iteg == -1
} // fin du cas où on utilise la normale aux pti
else //--- cas où on utilise la normale au noeud
{if (iteg == -1)
// on ne calcul que si iteg == -1, c-a-d qu'on est sortie
// de la boucle sur les pti
// sinon cas où on utilise la normale aux noeuds
{int nbddl = (*residu).Taille();
int dimf = 3; // ne fonctionne qu'en dim 3
if(ParaGlob::AxiSymetrie())
// cas axisymétrique, dans ce cas on ne prend en compte que les
// dimf-1 coordonnées donc on décrémente
dimf--;
int nbne = tab_noeud.Taille();
resD->Inita(0.);
maxi_F_t = 0.;
#ifdef MISE_AU_POINT
// on vérifie le bon dimensionnement
int taille = pt_StabMembBiel->Taille();
if (taille != nbne)
{cout << "\n erreur** ElemMeca::Cal_explicit_StabMembBiel(.."
<< " le stockage est mal dimensionne: taille = " << taille
<< " alors qu'il devrait etre : nbe = " << nbne << flush ;
Sortie(1);
};
#endif
// on applique une force de stabilisation aux noeuds
// uniquement suivant la direction de la normale
// la force est proportionelle au déplacement suivant la normale au noeud
int ix=1;
for (int ne =1; ne<= nbne; ne++)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// on initialise les forces et énergie
pt_StabMembBiel->FF_StabMembBiel(ne) = 0.;
pt_StabMembBiel->EE_StabMembBiel(ne) = pt_StabMembBiel->EE_StabMembBiel_t(ne);
// pour chaque noeud on va bloquer le déplacement dans la direction de la normale au noeud à t
// on récupère la normale au noeud
const TypeQuelconque& tiq = noe.Grandeur_quelconque(NN_SURF_t);
const Grandeur_coordonnee& gr= *((const Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
const Coordonnee& normale = gr.ConteneurCoordonnee_const();
//dans un calcul implicit on utilise la position à l'itération 0
// celle-ci peut-être différente de celle à t=0, lorsque l'on a une succession d'algo sur le même incrément
// on commence par récupérer les coordonnées pour l'itération 0
TypeQuelconque_enum_etendu enu(XI_ITER_0);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel = noe.Grandeur_quelconque(enu);
const Grandeur_coordonnee& cofo = *((Grandeur_coordonnee*) typquel.Const_Grandeur_pointee());
Coordonnee delta_X = noe.Coord2() - cofo.ConteneurCoordonnee_const();
double normal_fois_delta_X = normale * delta_X;
double ponder_F_t = 1.; // init
// dans le cas où on veut pondérer F_t, on calcule la fonction
if (pt_StabMembBiel->Pt_fct_ponder_Ft() != NULL)
{Fonction_nD* fct = pt_StabMembBiel->Pt_fct_ponder_Ft(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = -1; // signifie que l'on utilise la métrique qui vient juste d'être calculée
// au dernier pti !
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi
(absolue,TEMPS_tdt,li_enu_scal,nbint,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,nbint,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// seule la première valeur est ok
ponder_F_t = tab_ret(1) ;
}
// un - car la force doit s'opposer au déplacement
double intensite = -alpha * normal_fois_delta_X
+ ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(ne);
// on enregistre le maxi de la force à t
maxi_F_t = DabsMaX(maxi_F_t, pt_StabMembBiel->FF_StabMembBiel_t(ne));
// d'où une variation réelle d'intensité
double delta_intensite = intensite - ponder_F_t * pt_StabMembBiel->FF_StabMembBiel_t(ne);
pt_StabMembBiel->EE_StabMembBiel(ne) += 0.5 * delta_intensite * normal_fois_delta_X;
#ifdef MISE_AU_POINT
if ((niveau_affichage > 5))
{cout << "\n ElemMeca::Cal_explicit_StabMembBiel(..";
cout << "\n alpha= "<<alpha << " intensite= "<< intensite
<< " delta_intensite= " << delta_intensite
<< " delta_X= "<< delta_X
<< "\n normale= "<<normale << " normal_fois_delta_X= " << normal_fois_delta_X << flush;
};
#endif
// pour mémoire: intensite = -alpha * normal_fois_delta_X
// + pt_StabMembBiel->FF_StabMembBiel_t(num_pti);
for (int i=1;i<= dimf;i++,ix++)
{ // - normale(i) * stabili * normal_fois_delta_X;
(*resD)(ix) += normale(i) * intensite ;
};
}; // fin de la boucle sur les noeuds
// on est en dehors de la boucle, on effectue 2 choses
// 1) application d'une limitation éventuelle des efforts
// 2) prise en compte dans la raideur et résidu de l'élément, de la stabilisation
if (pt_StabMembBiel->Gestion_maxi_mini())
{
#ifdef MISE_AU_POINT
if ((niveau_affichage > 4))
{cout << "\n residu de stabilisation avant limitation: "<< (*resD)
<< "\n residu initial (sans stabilisation): "<< (*residu)
<< "\n maxi_F_t= " << maxi_F_t;
};
if ((niveau_affichage > 4))
{int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
cout << "\n F_stab("<<i<<")= "<<pt_StabMembBiel->FF_StabMembBiel(i)<< " ";
};
#endif
// 1) ---- on regarde s'il faut limiter la valeur de la stabilisation
// et on prend en compte dans la raideur globale et le résidu global
// la stabilisation
bool limitation = false;
double coef=1.; // coef de limitation éventuelle si limitation devient true
// récup des limitations éventuelles
double beta = pt_StabMembBiel->Beta();
double F_mini = pt_StabMembBiel->F_mini();
double d_maxi = pt_StabMembBiel->D_maxi();
// on refait une boucle sur les noeuds
int ix=1;
for (int ne =1; ne<= nbne; ne++)
{bool continuer = true;
if (noeud_a_prendre_en_compte != NULL)
if (!(*noeud_a_prendre_en_compte)(ne))
continuer = false;
if (continuer)
{Noeud& noe = *(tab_noeud(ne)); // pour simplifier
// récup de la force de stabilisation généralisée au noeud
// c'est différent de celle au pti car elle inclue la surface de l'élément
Coordonnee force_stab(dimf); // init
for (int i=1;i<= dimf;i++,ix++)
force_stab(i) = (*resD)(ix);
double intensite_generalise = force_stab.Norme();
// on récupère la force externe au noeud à l'instant t
TypeQuelconque_enum_etendu enu_ext(FORCE_GENE_EXT_t);
// récupération du conteneur pour lecture uniquement
const TypeQuelconque& typquel_ext = noe.Grandeur_quelconque(enu_ext);
const Grandeur_coordonnee& cofo_ext = *((Grandeur_coordonnee*) typquel_ext.Const_Grandeur_pointee());
// on récupère la normale au noeud
const TypeQuelconque& tiq = noe.Grandeur_quelconque(NN_SURF_t);
const Grandeur_coordonnee& gr= *((const Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee()));
const Coordonnee& normale = gr.ConteneurCoordonnee_const();
// on calcule l'intensité de la force externe selon la normale
double intensite_Fext = Dabs(cofo_ext.ConteneurCoordonnee_const() * normale);
// si l'intensité de la stabilisation est supérieure à beta * intensite_Fext
// on limite
double maxi_intensite = beta * intensite_Fext;
double max_intens_via_Fext = maxi_intensite; // on sauve pour le cas où on est en dep max
// on applique l'algo théorique cf. doc théorique
if (intensite_Fext >= F_mini)
// on est dans le cas où les forces externes existent
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intintensite très petites
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due aux forces externes, F_stab gene= "<< intensite_generalise
<< " ==> " << maxi_intensite;
if ((niveau_affichage > 4))
cout << "\n F_mini= "<< F_mini << " intensite_Fext= "<< intensite_Fext;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
// sinon pas de pb donc on ne fait rien
}
else // cas où les forces externes n'existent pas
// la limitation va s'exercée via un déplacement maxi
{maxi_intensite = Dabs(-alpha*d_maxi + maxi_F_t);
{if ((intensite_generalise > maxi_intensite)
&& (intensite_generalise > ConstMath::petit) // au cas où maxi_intintensite très petites
&& (Dabs(alpha) > ConstMath::pasmalpetit) // il faut qu'alpha ne soit pas nul
)
{limitation = true;
if ((niveau_affichage > 3))
{cout << "\n limitation stabilisation due deplace maxi acceptable, F_s"<< intensite_generalise
<< " ==> " << maxi_intensite << " (max_intens_via_Fext= "<<max_intens_via_Fext<<")" ;
if ((niveau_affichage > 4))
cout << "\n maxi_F_t= "<< maxi_F_t << " alpha*d_maxi= "
<< (alpha*d_maxi)
<< " alpha= " << alpha
<< " d_maxi= " << d_maxi;
};
coef = DabsMiN(coef, maxi_intensite/intensite_generalise);
};
};
};
}; // fin du test sur les noeuds à considérer
};// fin de la boucle sur les noeuds
if (limitation)
{ (*resD) *= coef; // mise à jour de la force de stabilisation
// idem au niveau des noeuds
int taille = pt_StabMembBiel->Taille();
for (int i=1;i<=taille;i++)
{pt_StabMembBiel->FF_StabMembBiel(i) *=coef;
// cas des énergies
double delta_energie_initiale =
pt_StabMembBiel->EE_StabMembBiel(i) - pt_StabMembBiel->EE_StabMembBiel_t(i);
// on coefficiente le delta
double delta_energie = delta_energie_initiale * coef;
//ce qui donne au final:
pt_StabMembBiel->EE_StabMembBiel(i) =
delta_energie + pt_StabMembBiel->EE_StabMembBiel_t(i);
};
}; // fin du cas if (limitation)
}; // fin du cas ou on applique éventuellement une limitation
// 2) ---- on met à jour la raideur et résidu global
// en tenant compte de la stabilisation
(*residu) += (*resD);
#ifdef MISE_AU_POINT
if ((niveau_affichage > 9))
{cout << "\n residu de stabilisation: "<< (*resD);
cout << "\n fin debug ElemMeca::Cal_explicit_StabMembBiel(..";
};
#endif
} // fin du cas où on utilise la normale aux noeuds
};
break;
case LIGNE:
break;
default: // on ne fait rien pour l'instant pour les autres types
break;
};
};
};
// fonction a renseigner par les classes dérivées, concernant les répercutions
// éventuelles due à la suppression de tous les frontières
// nums_i : donnent les listes de frontières supprimées
void ElemMeca::Prise_en_compte_des_consequences_suppression_tous_frontieres()
{// il faut supprimer toutes les déformations liés aux frontières
int tail_S = defSurf.Taille();
for (int i=1;i<=tail_S;i++)
{ delete defSurf(i);
defSurf(i)=NULL;
};
int tail_A = defArete.Taille();
for (int i=1;i<=tail_A;i++)
{ delete defArete(i);
defArete(i)=NULL;
};
};
// idem pour une frontière (avant qu'elle soit supprimée)
void ElemMeca::Prise_en_compte_des_consequences_suppression_une_frontiere(ElFrontiere* elemFront)
{ int taille = tabb.Taille();
// on choisit en fonction du type de frontière
if (( elemFront->Type_geom_front() == LIGNE)&& (defArete.Taille()!=0))
{int taille_ligne = taille - posi_tab_front_lin; // a priori = cas sans les points
if (posi_tab_front_point != 0) // cas où il y a des points
taille_ligne = posi_tab_front_point-posi_tab_front_lin;
for (int i=1; i<=taille_ligne; i++)
{if ((tabb(i+posi_tab_front_lin) == elemFront)&& (defArete(i) != NULL))
{delete defArete(i);defArete(i)=NULL;break;}
}
}
else if (( elemFront->Type_geom_front() == SURFACE) && (defSurf.Taille()!=0))
{if (posi_tab_front_lin != 0) // si == 0 cela signifie qu'il n'y a pas de surface à supprimer !!
{ for (int i=1; i<=posi_tab_front_lin; i++) // posi_tab_front_lin == le nombre de surfaces, qui sont obligatoirement en début de tableau
if ((tabb(i) == elemFront)&& (defSurf(i) != NULL))
{delete defSurf(i);defSurf(i)=NULL;break;}
};
};
};
// Calcul des frontieres de l'element
// creation des elements frontieres et retour du tableau de ces elements
// la création n'a lieu qu'au premier appel
// ou lorsque l'on force le paramètre force a true
// dans ce dernier cas seul les frontière effacées sont recréée
// cas :
// = 0 -> on veut toutes les frontières
// = 1 -> on veut uniquement les surfaces
// = 2 -> on veut uniquement les lignes
// = 3 -> on veut uniquement les points
// = 4 -> on veut les surfaces + les lignes
// = 5 -> on veut les surfaces + les points
// = 6 -> on veut les lignes + les points
Tableau <ElFrontiere*> const & ElemMeca::Frontiere_elemeca(int cas, bool force)
{// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
int taille = tabb.Taille(); // la taille initiales des frontières
if (force) // dans ce cas on commence par tout effacer
{ // on efface les surfaces (s'il y en a)
for (int i=1;i<=posi_tab_front_lin;i++)
{if (tabb(i) != NULL)
{delete tabb(i); //on commence par supprimer
tabb(i)=NULL;
// on supprime également éventuellement la déformation associée
if (defSurf.Taille() != 0)
if (defSurf(i) != NULL)
{delete defSurf(i);defSurf(i)=NULL;};
ind_front_surf = 0; // on indique qu'il ne reste plus de frontière surface
};
};
// on efface les lignes (s'il y en a)
for (int i=1;i<=posi_tab_front_point - posi_tab_front_lin;i++)
{if (tabb(i+posi_tab_front_lin) != NULL)
{delete tabb(i+posi_tab_front_lin); //on commence par supprimer
tabb(i+posi_tab_front_lin)=NULL;
// on supprime également éventuellement la déformation associée
if (defArete.Taille() != 0)
if (defArete(i) != NULL)
{delete defArete(i);defArete(i)=NULL;};
ind_front_lin = 0; // on indique qu'il ne reste plus de frontière ligne
};
};
// on efface les points (s'il y en a)
for (int i=1;i<=taille - posi_tab_front_point;i++)
{if (tabb(i+posi_tab_front_point) != NULL)
{delete tabb(i+posi_tab_front_point); //on commence par supprimer
tabb(i+posi_tab_front_point)=NULL;
ind_front_point = 0; // on indique qu'il ne reste plus de frontière ligne
};
};
};
// -- maintenant on s'occupe de la construction conditionnelle
bool built_surf = false;bool built_ligne = false; bool built_point = false;
switch (cas)
{case 0: built_surf = built_ligne = built_point = true; break;
case 1: built_surf = true; break;
case 2: built_ligne = true; break;
case 3: built_point = true; break;
case 4: built_surf = built_ligne = true; break;
case 5: built_surf = built_point = true; break;
case 6: built_ligne = built_point = true; break;
};
if ( ((ind_front_surf == 0)&& (ind_front_lin == 0) && (ind_front_point == 0)) || force )
{
// récup de l'élément géométrique
ElemGeomC0& el = ElementGeometrique();
int tail_ar = el.NbSe(); // nombre potentiel d'arêtes
int tail_fa = el.NbFe(); // nombre potentiel de faces
int tail_po = el.Nbne(); // nombre potentiel de points
// récup du tableau de ddl actuel
const DdlElement & tdd = TableauDdl();
// --- on va construire en fonction des indicateurs des tableaux intermédiaires
int new_posi_tab_front_point = 0; //init par défaut
int new_posi_tab_front_lin = 0; //init par défaut
int new_ind_front_point = 0;
int new_ind_front_lin = 0;
int new_ind_front_surf = 0;
// -- pour les surfaces
Tableau <ElFrontiere*> tabb_surf;
if ((built_surf)&& ((ind_front_surf == 0)||force))
{tabb_surf.Change_taille(tail_fa,NULL);// init par défaut
for (int num=1;num<=tail_fa;num++)
{ int nbnoe = el.Nonf()(num).Taille(); // nb noeud de la surface
Tableau <Noeud *> tab(nbnoe); // les noeuds de la frontiere
DdlElement ddelem(nbnoe); // les ddlelements des noeuds frontieres
for (int i=1;i<= nbnoe;i++)
{ tab(i) = tab_noeud(el.Nonf()(num)(i));
ddelem.Change_un_ddlNoeudElement(i,tdd(el.Nonf()(num)(i)));
};
tabb_surf(num) = new_frontiere_surf(num,tab,ddelem);
};
// nouveau indicateur d'existence
new_ind_front_surf = 1;
// on positionne les nouvelles positions
new_posi_tab_front_point += tail_fa;
new_posi_tab_front_lin += tail_fa;
};
// -- pour les lignes
Tableau <ElFrontiere*> tabb_ligne;
if ((built_ligne)&& ((ind_front_lin == 0)||force))
{tabb_ligne.Change_taille(tail_ar,NULL);// init par défaut
for (int num=1;num<=tail_ar;num++)
{ int nbnoe = el.NonS()(num).Taille(); // nb noeud de l'arête
Tableau <Noeud *> tab(nbnoe); // les noeuds de l'arête frontiere
DdlElement ddelem(nbnoe); // les ddlelements des noeuds frontieres
for (int i=1;i<= nbnoe;i++)
{ tab(i) = tab_noeud(el.NonS()(num)(i));
ddelem.Change_un_ddlNoeudElement(i,tdd(el.NonS()(num)(i)));
};
tabb_ligne(num) = new_frontiere_lin(num,tab,ddelem);
};
// nouveau indicateur d'existence
new_ind_front_lin = 1;
// on positionne les nouvelles positions
new_posi_tab_front_point += tail_ar;
};
// -- pour les points
Tableau <ElFrontiere*> tabb_point;
if ((built_point) && ((ind_front_point == 0)||force))
{tabb_point.Change_taille(tail_po,NULL);// init par défaut
// maintenant création des frontière point éventuellement
Tableau <Noeud *> tab(1); // les noeuds de la frontiere (tab de travail)
DdlElement ddelem(1); // les ddlelements des points frontieres (tab de travail)
for (int i=1;i<=tail_po;i++)
if (tabb_point(i) == NULL)
{ tab(1) = tab_noeud(i);
ddelem.Change_un_ddlNoeudElement(1,tdd(i));
tabb_point(i) = new FrontPointF(tab,ddelem);
};
// nouveau indicateur d'existence
new_ind_front_point = 1;
};
// --- mise à jour du tableau globale et des indicateurs ad hoc
int taille_finale = tabb_surf.Taille()+tabb_ligne.Taille()+tabb_point.Taille();
tabb.Change_taille(taille_finale);
// cas des points
if (new_ind_front_point) // là is s'agit de nouveaux éléments
{for (int i=tail_po;i>0;i--) // transfert pour les noeuds
{ tabb(i+new_posi_tab_front_point) = tabb_point(i);}
}
else if (ind_front_point) // là il s'agit d'anciens éléments
{for (int i=tail_po;i>0;i--) // transfert pour les noeuds en descendant
{ tabb(i+new_posi_tab_front_point) = tabb(i+posi_tab_front_point);}
};
// cas des lignes
if (new_ind_front_lin) // là il s'agit de nouveaux éléments
{for (int i=1;i<=tail_ar;i++) // transfert
{ tabb(i+new_posi_tab_front_lin) = tabb_ligne(i);}
}
else if (ind_front_lin) // là il s'agit d'anciens éléments
{for (int i=tail_ar;i>0;i--) // transfert en descendant
{ tabb(i+new_posi_tab_front_lin) = tabb(i+posi_tab_front_lin);}
};
// cas des surfaces
if (new_ind_front_surf) // là is s'agit de nouveaux éléments
{for (int i=1;i<=tail_fa;i++) // transfert
{ tabb(i) = tabb_surf(i);}
};
// dans le cas où il y avait des anciens éléments, il n'y a rien n'a faire
// car le redimentionnement de tabb ne change pas les premiers éléments
// mis à jour des indicateurs
ind_front_surf = new_ind_front_surf;
posi_tab_front_lin = new_posi_tab_front_lin;
ind_front_lin = new_ind_front_lin;
posi_tab_front_point = new_posi_tab_front_point;
ind_front_point = new_ind_front_point;
};
// retour du tableau
return (Tableau <ElFrontiere*>&)tabb;
};
// ramène la frontière point
// éventuellement création des frontieres points de l'element et stockage dans l'element
// si c'est la première fois sinon il y a seulement retour de l'elements
// a moins que le paramètre force est mis a true
// dans ce dernier cas la frontière effacéee est recréée
// num indique le numéro du point à créer (numérotation EF)
ElFrontiere* const ElemMeca::Frontiere_points(int num,bool force)
{ // -- tout d'abord on évacue le cas où il n'y a pas de frontière point à calculer
// récup de l'élément géométrique
ElemGeomC0& el = ElementGeometrique();
int tail_po = el.Nbne(); // nombre potentiel de points
if (num > tail_po)
return NULL;
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
// on regarde si les frontières points existent sinon on les crée
if (ind_front_point == 1)
{return (ElFrontiere*)tabb(posi_tab_front_point+num);}
else if ( ind_front_point == 2)
// cas où certaines frontières existent
{if (tabb(posi_tab_front_point+num) != NULL)
return (ElFrontiere*)tabb(posi_tab_front_point+num);
};
// arrivée ici cela veut dire que la frontière point n'existe pas
// on l'a reconstruit éventuellement
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
if ((ind_front_point == 0) || force )
{// récup du tableau de ddl
const DdlElement & tdd = TableauDdl();
int taille = tabb.Taille(); // la taille initiales des frontières
int tail_fa = el.NbFe(); // nombre potentiel de faces
int tail_ar = el.NbSe(); // nombre potentiel d'arêtes
// dimensionnement du tableau de frontières ligne si nécessaire
if (ind_front_point == 0)
{if ((ind_front_lin > 0) && (ind_front_surf == 0))
// cas où les frontières lignes existent seules, on ajoute les points
{ int taille_finale = tail_ar + tail_po;
tabb.Change_taille(taille_finale);
for (int i=1;i<= tail_ar;i++) // transfert pour les lignes
tabb(i) = tabb(i+tail_ar);
posi_tab_front_point=tail_ar;
posi_tab_front_lin = 0; // car on n'a pas de surface
}
else if ((ind_front_lin > 0) && (ind_front_surf > 0))
// cas où les frontières lignes existent et surfaces et pas de points, donc on les rajoutes
{ int taille_finale = tail_fa + tail_po+tail_ar;
tabb.Change_taille(taille_finale);// les grandeurs pour les surfaces et les lignes sont
// conservées, donc pas de transferts à prévoir
posi_tab_front_point=tail_ar+tail_fa; // après les faces et les lignes
posi_tab_front_lin = tail_fa; // après les surfaces
}
else
{ // cas où il n'y a pas de frontières lignes
if (ind_front_surf == 0) // cas où il n'y a pas de surface
{tabb.Change_taille(tail_po,NULL); // on n'a pas de ligne, pas de point et pas de surface
posi_tab_front_point = posi_tab_front_lin = 0;}
else {tabb.Change_taille(tail_po+tail_fa); // cas où les surfaces existent
// le redimensionnement n'affecte pas les surfaces qui sont en début de tableau
posi_tab_front_lin = posi_tab_front_point = tail_fa; // après les surfaces
};
};
// et on met les pointeurs de points en NULL
for (int i=1;i<= tail_po;i++)
{ tabb(i+posi_tab_front_point) = NULL;}
};
// maintenant création de la frontière point
Tableau <Noeud *> tab(1); // les noeuds de la frontiere
tab(1) = tab_noeud(num);
DdlElement ddelem(1); // les ddlelements des points frontieres
ddelem.Change_un_ddlNoeudElement(1,tdd(num));
tabb(posi_tab_front_point+num) = new FrontPointF(tab,ddelem);
// on met à jour l'indicateur ind_front_point
ind_front_point = 1; // a priori
for (int npoint=1;npoint<=tail_po;npoint++)
if (tabb(posi_tab_front_point+npoint) == NULL)
{ ind_front_point = 2; break;};
};
// maintenant normalement la frontière est créé on la ramène
return (ElFrontiere*)tabb(posi_tab_front_point+num);
};
// ramène la frontière linéique
// éventuellement création des frontieres linéique de l'element et stockage dans l'element
// si c'est la première fois et en 3D sinon il y a seulement retour de l'elements
// a moins que le paramètre force est mis a true
// dans ce dernier cas la frontière effacéee est recréée
// num indique le numéro de l'arête à créer (numérotation EF)
// nbneA: nombre de noeuds des segments frontières
// el : l'élément
ElFrontiere* const ElemMeca::Frontiere_lineique(int num,bool force)
{ // -- tout d'abord on évacue le cas où il n'y a pas de frontière linéique à calculer
// récup de l'élément géométrique
ElemGeomC0& el = ElementGeometrique();
int tail_ar = el.NbSe(); // nombre potentiel d'arêtes
if (num > tail_ar)
return NULL;
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
// on regarde si les frontières linéiques existent sinon on les crée
if (ind_front_lin == 1)
{return (ElFrontiere*)tabb(posi_tab_front_lin+num);}
else if ( ind_front_lin == 2)
// cas où certaines frontières existent
{if (tabb(posi_tab_front_lin+num) != NULL)
return (ElFrontiere*)tabb(posi_tab_front_lin+num);
};
// arrivée ici cela veut dire que la frontière ligne n'existe pas
// on l'a reconstruit
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
if ((ind_front_lin == 0) || force )
{// récup du tableau de ddl
const DdlElement & tdd = TableauDdl();
int taille = tabb.Taille(); // la taille initiales des frontières
int tail_fa = el.NbFe(); // nombre potentiel de faces
int tail_po = el.Nbne(); // nombre potentiel de points
// dimensionnement du tableau de frontières ligne si nécessaire
if (ind_front_lin == 0)
{if ((ind_front_point > 0) && (ind_front_surf == 0))
// cas où les frontières points existent seules, on ajoute les lignes
{ int taille_finale = tail_ar + tail_po;
tabb.Change_taille(taille_finale);
for (int i=1;i<= tail_po;i++)
tabb(i+tail_ar) = tabb(i);
posi_tab_front_point=tail_ar;
posi_tab_front_lin = 0; // car on n'a pas de surface
}
else if ((ind_front_point > 0) && (ind_front_surf > 0))
// cas où les frontières points existent et surfaces et pas de ligne, donc on les rajoutes
{ int taille_finale = tail_fa + tail_po+tail_ar;
tabb.Change_taille(taille_finale);// les grandeurs pour les surfaces sont conservées
for (int i=1;i<= tail_po;i++) // transfert pour les noeuds
{ tabb(i+tail_ar+tail_fa) = tabb(i+tail_fa);};
posi_tab_front_point=tail_ar+tail_fa; // après les faces et les lignes
posi_tab_front_lin = tail_fa; // après les surfaces
}
else
{ // cas où il n'y a pas de frontières points
if (ind_front_surf == 0) // cas où il n'y a pas de surface
{tabb.Change_taille(tail_ar,NULL); // on n'a pas de ligne, pas de point et pas de surface
posi_tab_front_lin = 0; posi_tab_front_point = tail_ar;
}
else {tabb.Change_taille(tail_ar+tail_fa); // cas où les surfaces existent
// le redimensionnement n'affecte pas les surfaces qui sont en début de tableau
posi_tab_front_lin = tail_fa; // après les surfaces
posi_tab_front_point = tail_fa+tail_ar;
};
};
// et on met les pointeurs de lignes en NULL
for (int i=1;i<= tail_ar;i++) // transfert pour les noeuds
{ tabb(i+posi_tab_front_lin) = NULL;}
};
// maintenant création de la ligne
int nbnoe = el.NonS()(num).Taille(); // nb noeud de l'arête
Tableau <Noeud *> tab(nbnoe); // les noeuds de l'arête frontiere
DdlElement ddelem(nbnoe); // les ddlelements des noeuds frontieres
for (int i=1;i<= nbnoe;i++)
{ tab(i) = tab_noeud(el.NonS()(num)(i));
ddelem.Change_un_ddlNoeudElement(i,tdd(el.NonS()(num)(i)));
};
tabb(posi_tab_front_lin+num) = new_frontiere_lin(num,tab,ddelem);
ind_front_lin = 1; // a priori
for (int nlign=1;nlign<=tail_ar;nlign++)
if (tabb(posi_tab_front_lin+nlign) == NULL)
{ ind_front_lin = 2; break;};
};
// maintenant normalement la frontière est créé on la ramène
return (ElFrontiere*)tabb(posi_tab_front_lin+num);
};
// ramène la frontière surfacique
// éventuellement création des frontieres surfacique de l'element et stockage dans l'element
// si c'est la première fois sinon il y a seulement retour de l'elements
// a moins que le paramètre force est mis a true
// dans ce dernier cas la frontière effacéee est recréée
// num indique le numéro de la surface à créer (numérotation EF)
ElFrontiere* const ElemMeca::Frontiere_surfacique(int num,bool force)
{ // -- tout d'abord on évacue le cas où il n'y a pas de frontière surfacique à calculer
// récup de l'élément géométrique
ElemGeomC0& el = ElementGeometrique();
int tail_fa = 0; // init
// test --- retour rapide si l'existance n'est pas possible
if (ParaGlob::AxiSymetrie())
// dans le cas axiSymétrique, le test est un peu différent car les surfaces sont générées
// par les lignes,
{tail_fa = el.NbSe(); // nombre potentiel de segments
if (num > tail_fa)
return NULL;
}
else // dans le cas où on n'est pas en axi, test simple
{tail_fa = el.NbFe(); // nombre potentiel de faces
if (num > tail_fa)
return NULL;
};
// --- la suite concerne le cas où on peut faire des surfaces (non axi ou axi !)
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
// on regarde si les frontières surfacique existent sinon on les crée
if (ind_front_surf == 1)
{return (ElFrontiere*)tabb(num);}
else if ( ind_front_surf == 2)
// cas où certaines frontières existent
{if (tabb(num) != NULL)
return (ElFrontiere*)tabb(num);
};
// arrivée ici cela veut dire que la frontière surface n'existe pas
// on l'a reconstruit
// le calcul et la création ne sont effectués qu'au premier appel
// ou lorsque l'on veut forcer une recréation
if ((ind_front_surf == 0) || force )
{// récup du tableau de ddl
const DdlElement & tdd = TableauDdl();
int taille = tabb.Taille(); // la taille initiales des frontières
int tail_ar = el.NbSe(); // nombre potentiel d'arêtes
int tail_po = el.Nbne(); // nombre potentiel de points
// dimensionnement du tableau de frontières surfaces si nécessaire
if (ind_front_surf == 0)
{if ((ind_front_point > 0) && (ind_front_lin == 0))
// cas où les frontières points existent seules, on ajoute les surfaces
{ int taille_finale = tail_fa + tail_po;
tabb.Change_taille(taille_finale);
for (int i=1;i<= tail_po;i++)
tabb(i+tail_fa) = tabb(i);
posi_tab_front_lin=posi_tab_front_point=tail_fa;
}
else if ((ind_front_point > 0) && (ind_front_lin > 0))
// cas où les frontières points existent et lignes et pas de surfaces, donc on les rajoutes
{ int taille_finale = tail_fa + tail_po+tail_ar;
tabb.Change_taille(taille_finale);
// --transfert pour les noeuds et les lignes
for (int i=1;i<= tail_po;i++) // transfert pour les noeuds
{ tabb(i+tail_ar+tail_fa) = tabb(i+tail_ar);};
for (int i=1;i<= tail_ar;i++) // transfert pour les lignes
{ tabb(i+tail_fa) = tabb(i);};
// --def des indicateurs
posi_tab_front_point=tail_ar+tail_fa; // après les faces et les lignes
posi_tab_front_lin = tail_fa; // après les surfaces
}
else
{ // cas où il n'y a pas de frontières points
if (ind_front_lin == 0) // cas où il n'y a pas de lignes
{tabb.Change_taille(tail_fa,NULL); // on n'a pas de ligne, pas de point et pas de surface
posi_tab_front_lin = posi_tab_front_point = tail_fa;
} // on peut tout mettre à NULL
else {tabb.Change_taille(tail_ar+tail_fa); // cas où les lignes existent
for (int i=1;i<= tail_ar;i++) // transfert pour les lignes
tabb(i+tail_fa) = tabb(i);
posi_tab_front_lin = tail_fa; // après les surfaces
posi_tab_front_point = tail_fa + tail_ar;
};
};
// --et on met les pointeurs de surfaces en NULL
for (int i=1;i<=tail_fa;i++)
tabb(i)=NULL;
};
// maintenant création de la surface
int nbnoe = el.Nonf()(num).Taille(); // b noeud de la surface
Tableau <Noeud *> tab(nbnoe); // les noeuds de la frontiere
DdlElement ddelem(nbnoe); // les ddlelements des noeuds frontieres
for (int i=1;i<= nbnoe;i++)
{ tab(i) = tab_noeud(el.Nonf()(num)(i));
ddelem.Change_un_ddlNoeudElement(i,tdd(el.Nonf()(num)(i)));
};
tabb(num) = new_frontiere_surf(num,tab,ddelem);
ind_front_surf = 1; // a priori
for (int nsurf=1;nsurf<=tail_fa;nsurf++)
if (tabb(nsurf) == NULL)
{ ind_front_surf = 2; break;};
};
// maintenant normalement la frontière est créé on la ramène
return (ElFrontiere*)tabb(num);
};
// calcul des intégrales de volume et de volume + temps
// cas = 1 : pour un appel après un calcul implicit
// cas = 2 : pour un appel après un calcul explicit
void ElemMeca::Calcul_Integ_vol_et_temps(int cas,const double& poid_jacobien,int iteg)
{ // 1) intégration de volume uniquement
//List_io <TypeQuelconque> integ_vol_typeQuel,integ_vol_typeQuel_t;
if (integ_vol_typeQuel != NULL)
// on balaie les conteneurs
{ int taille = integ_vol_typeQuel->Taille();
Tableau <TypeQuelconque> & tab_integ = *integ_vol_typeQuel; // pour simplifier
for (int il =1;il <= taille ;il++)
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ((*index_Integ_vol_typeQuel)(il))
{// différent choix
switch (tab_integ(il).Const_Grandeur_pointee()->Type_enumGrandeurParticuliere())
{case PARTICULIER_DDL_ETENDU:
{ Grandeur_Ddl_etendu& gr
= *((Grandeur_Ddl_etendu*) tab_integ(il).Grandeur_pointee()); // pour simplifier
Ddl_etendu& ddl = *(gr.ConteneurDdl_etendu());// récup du ddl
// on va utiliser la méthode Valeur_multi ce qui n'est pas optimal en vitesse
// mais 1) c'est cohérent avec l'implantation existante
// 2) identique avec les ddl déjà accessibles
List_io<Ddl_enum_etendu> inter; // une liste intermédiaire
inter.push_back(ddl.DdlEnumEtendu());
bool absolue = true; // on se place systématiquement en absolu
// calcul de la valeur et retour dans tab_d
Tableau <double> tab_d(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,inter,iteg,-cas));
// on cumule l'intégrale
ddl.Valeur() += tab_d(1) * poid_jacobien;
////---- debug
//cout << "\n debug: ElemMeca::Calcul_Integ_vol_et_temps(";
//cout << " tab_d(1)= "<< tab_d(1) << " ";ddl.Affiche();
//tab_noeud(1)->Affiche();
////---- fin debug
break;
}
case PARTICULIER_VECTEUR_NOMMER:
{ Grandeur_Vecteur_Nommer& gr
= *((Grandeur_Vecteur_Nommer*) tab_integ(il).Grandeur_pointee()); // pour simplifier
Fonction_nD* fct = gr.Fct(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,iteg,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,iteg,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// on cumule l'intégrale
Vecteur& V = gr.ConteneurVecteur();
int taille = V.Taille();
for (int i=1;i<=taille;i++)
V(i) += tab_ret(i) * poid_jacobien;
////----- debug
//cout << "\n debug *** ElemMeca::Grandeur_particuliere( ";
//cout << "\n li_enu_scal= "<< li_enu_scal << ", val_ddl_enum= " << val_ddl_enum;
//cout << "\n li_quelc= "<< li_quelc;
//V.Affiche();cout << " tab_ret= " << tab_ret << endl;;
////----- fin debug
break;
}
break;
default:
cout << " \n *** pas d'integration prevu pour la grandeur " << tab_integ(il);
if (ParaGlob::NiveauImpression() > 0)
cout << "\n ElemMeca::Calcul_Integ_vol_et_temps(..." << endl;
};
};
};
// 2) intégration de volume et en temps
//List_io <TypeQuelconque> integ_vol_t_typeQuel,integ_vol_t_typeQuel_t;
if (integ_vol_t_typeQuel != NULL)
// on balaie les conteneurs
{ int taille = integ_vol_t_typeQuel->Taille();
Tableau <TypeQuelconque> & tab_integ = *integ_vol_t_typeQuel; // pour simplifier
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
for (int il =1;il <= taille ;il++)
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ((*index_Integ_vol_t_typeQuel)(il))
{// différent choix
switch (tab_integ(il).Const_Grandeur_pointee()->Type_enumGrandeurParticuliere())
{case PARTICULIER_DDL_ETENDU:
{ Grandeur_Ddl_etendu& gr
= *((Grandeur_Ddl_etendu*) tab_integ(il).Grandeur_pointee()); // pour simplifier
Ddl_etendu& ddl = *(gr.ConteneurDdl_etendu());// récup du ddl
// on va utiliser la méthode Valeur_multi ce qui n'est pas optimal en vitesse
// mais 1) c'est cohérent avec l'implantation existante
// 2) identique avec les ddl déjà accessibles
List_io<Ddl_enum_etendu> inter; // une liste intermédiaire
inter.push_back(ddl.DdlEnumEtendu());
bool absolue = true; // on se place systématiquement en absolu
// calcul de la valeur et retour dans tab_d
Tableau <double> tab_d(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,inter,iteg,-cas));
// on cumule l'intégrale
ddl.Valeur() += deltat * tab_d(1) * poid_jacobien;
break;
}
case PARTICULIER_VECTEUR_NOMMER:
{ Grandeur_Vecteur_Nommer& gr
= *((Grandeur_Vecteur_Nommer*) tab_integ(il).Grandeur_pointee()); // pour simplifier
Fonction_nD* fct = gr.Fct(); // récupération de la fonction nD
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,iteg,-cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,iteg,-cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// on cumule l'intégrale
Vecteur& V = gr.ConteneurVecteur();
int taille = V.Taille();
for (int i=1;i<=taille;i++)
V(i) += deltat * tab_ret(i) * poid_jacobien;
break;
}
break;
default:
cout << " \n *** pas d'integration prevu pour la grandeur " << tab_integ(il);
if (ParaGlob::NiveauImpression() > 0)
cout << "\n ElemMeca::Calcul_Integ_vol_et_temps(..." << endl;
};
};
};
};
// au premier pti init éventuel des intégrales de volume et temps
void ElemMeca::Init_Integ_vol_et_temps()
{ // 1) intégration de volume uniquement
if (integ_vol_typeQuel != NULL)
{int taille = integ_vol_typeQuel->Taille();
Tableau <TypeQuelconque> & tab_integ = *integ_vol_typeQuel; // pour simplifier
for (int il =1;il <= taille ;il++)
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ((*index_Integ_vol_typeQuel)(il))
{tab_integ(il).Grandeur_pointee()->InitParDefaut(); // on initialise le conteneur
};
};
// 2) intégration de volume et en temps
if (integ_vol_t_typeQuel != NULL)
{int taille = integ_vol_t_typeQuel->Taille();
Tableau <TypeQuelconque> & tab_integ = *integ_vol_t_typeQuel; // pour simplifier
Tableau <TypeQuelconque> & tab_integ_t = *integ_vol_t_typeQuel_t; // ""
for (int il =1;il <= taille ;il++)
// l'intégrale n'est effectuée que si l'index est positif, sinon, il s'agit
// d'une valeur figée
if ((*index_Integ_vol_t_typeQuel)(il))
{(*tab_integ(il).Grandeur_pointee()) = (*tab_integ_t(il).Grandeur_pointee());
};
};
};
// définition d'un repère d'anisotropie à un point d'intégration
BaseH ElemMeca::DefRepereAnisotropie
(int iteg,LesFonctions_nD* lesFonctionsnD,const BlocGen & bloc)
{ // la définition du repère dépend du type de définition
// (cf. la def du bloc dans LesMaillages::Completer .... dans LesMaillages.cc)
if (bloc.Nom(4) == "par_fonction_nD_")
{// on doit simplement appeler la fonction nD
// Nom(5) est le nom de la fonction nD
Fonction_nD * fct = lesFonctionsnD->Trouve(bloc.Nom(5));
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Expli* ex_expli = NULL;
// pour l'instant on ne peut s'occuper que des grandeurs à t au maximum: mais pas
// à tdt
bool premier_calcul = true;
// def->ChangeNumInteg(iteg);
// const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
// ex_expli = &ex;
// // on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// // pour les grandeurs strictement scalaire
// Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
// (absolue,TEMPS_t,li_enu_scal,*def,ex_impli,ex_expli_tdt,ex_expli)
// );
// // on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// // pour les Coordonnees et Tenseur
// Valeurs_Tensorielles_interpoler_ou_calculer
// (absolue,TEMPS_t,li_quelc,*def,ex_impli,ex_expli_tdt,ex_expli);
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
int cas = 1; // on se met dans un cas normal avec utilisation de certaine partie de la métrique implicite
// a priori ce serait idem avec l'explicite
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_t,li_enu_scal,iteg,cas));
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
Valeurs_Tensorielles(absolue, TEMPS_t,li_quelc,iteg,cas);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_ret = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
// on doit récupérer 6 composantes en 3D (pour le reste c'est actuellement en attente)
int dim = ParaGlob::Dimension();
if (dim != 3)
{ cout << "\n arret: le type de construction de repere: par_fonction_nD_ "
<< " n'est implante que pour la dimension 3 pour l'instant !!! "
<< " affaire a suivre ....(et se plaindre! ) ";
Sortie(1);
};
int nb_composante = tab_ret.Taille();
if (nb_composante != 6)
{ cout << "\n erreur dans la construction d'un repere d'orthotropie par_fonction_nD_ "
<< " la fonction nD ramene "<<nb_composante << " composante(s) "
<< " or il en faut 6 uniquement ! ";
Sortie(1);
};
// arrivée ici, cela veut dire que tout est ok
// --- construction de la base dans I_a -------
BaseH baseH(3,3);
// premier vecteur:
CoordonneeH& c1 = baseH.CoordoH(1);
for (int i=1;i<4;i++)
c1(i) = tab_ret(i);
// on norme si c'est demandé
if ( bloc.Nom(3) == "orthonorme_")
{c1.Normer();
}
else
{ cout << "\n arret: le type de repere: "<<bloc.Nom(3)
<< " n'est pas encore implante !!! affaire a suivre ....(et se plaindre! ) ";
Sortie(1);
};
// le deuxième vecteur:
CoordonneeH& c2 = baseH.CoordoH(2);
for (int i=1;i<4;i++)
c2(i) = tab_ret(i+3);
// on norme
c2.Normer();
// le troisième vecteur
CoordonneeH& c3 = baseH.CoordoH(3);
c3 = Util::ProdVec_coorH(c1, c2);
/* ***** la suite est supprimée: elle correspond au calcul des coordonnées locales alpha_a^{.i}
du repère d'anisotropie. Le pb est que le repère local g_i à t= 0 peut changer pendant le calcul
par exemple dans le cas de l'utilisation des contraintes planes on va avoir un repère de travail
qui change. Sans doute que ce sera un pb identique avec les Umat
donc on passe en paramètre le repère en coordonnée absolue et au moment de l'utilisation
les coordonnées locales seront calculées
si c'est ok par la suite, il faudra supprimer la partie commentée qui suit
// --- maintenant on s'occupe de calculer la base en local -----
// en fait il s'agit de définir les coordonnées locales par projection de la base
// globale en local
def->ChangeNumInteg(iteg);
const Met_abstraite::InfoExp_t & ex_infoExp_t = def->RemontExp_t();
int nb_vec = ex_infoExp_t.giB_0->NbVecteur(); // nb vecteur local
BaseH base_ret_H(3,3);
// les coordonnées locales : théta^te = V * g^te
// soit O_a la base d'anisotropie: O_a = beta_a^{.i} g_i
// ona beta_a^{.i} = O_a . g^i
// et on stocke dans base_ret_H de la manière suivante
// base_ret_H(a)(i) = beta_a^{.i}
for (int a = 1;a <= 3; a++)
{CoordonneeH& inter = base_ret_H.CoordoH(a);
for (int i=1;i<nb_vec;i++)
inter(i)= baseB.CoordoB(a) * (*ex_infoExp_t.giH_0)(i);
};
// dans le cas où nb_vec n'est pas égale à 3 il faut complèter la base
switch (nb_vec)
{ case 3: break; // cas où c'est déjà ok
case 2: // cas où l'élément est une surface, le troisième vecteur est la normale
{CoordonneeH N_H = Util::ProdVec_coorH((*ex_infoExp_t.giH_0)(1),(*ex_infoExp_t.giH_0)(2));
for (int a = 1;a <= 3; a++)
base_ret_H.CoordoH(a)(3)= baseB.CoordoB(a) * N_H;
break;
}
case 1: // cas où l'élément est un segment, il faut définir deux autres vecteurs
{// on génère
la suite est problèmatique car lorsque l'on changera de repère local, les coordonnées locales
devraient changer donc on arrête
CoordonneeH N_H = Util::ProdVec_coorH((*ex_infoExp_t.giH_0)(1),(*ex_infoExp_t.giH_0)(2));
for (int a = 1;a <= 3; a++)
base_ret_H.CoordoH(a)(3)= baseB.CoordoB(a) * N_H;
break;
}
default:
break;
}
for (int i=1;i<4;i++)
{CoordonneeH& inter = base_ret_H.CoordoH(i);
for (int te = 1;te <= nb_vec; te++)
inter(te)= baseB.CoordoB(i) * (*ex_infoExp_t.giH_0)(te);
};
return base_ret_H;
*/
//cout << "\n debug ElemMeca::DefRepereAnisotropie ";
//for (int a = 1;a < 4; a++)
// {cout << "\n baseH("<<a<<"): ";
// baseH(a).Affiche();
// };
return baseH;
}
else if (bloc.Nom(4) == "Par_projection_et_normale_au_plan_")
{ cout << "\n arret: le type de construction de repere: "<<bloc.Nom(3)
<< " n'est pas encore implante !!! affaire a suivre ....(et se plaindre! ) ";
Sortie(1);
}
else
{ cout << "\n **** erreur: type de construction de repere: "<<bloc.Nom(3)
<< " inconnu ?? ";
Sortie(1);
};
return BaseH(); // pour faire taire le compilo
};
// calcul éventuel de la normale à un noeud
// ce calcul existe pour les éléments 2D, 1D axi, et aussi pour les éléments 1D
// qui possède un repère d'orientation
// en retour coor = la normale si coor.Dimension() est = à la dimension de l'espace
// si le calcul n'existe pas --> coor.Dimension() = 0
// ramène un entier :
// == 1 : calcul normal
// == 0 : problème de calcul -> coor.Dimension() = 0
// == 2 : indique que le calcul n'est pas licite pour le noeud passé en paramètre
// c'est le cas par exemple des noeuds exterieurs pour les éléments SFE
// mais il n'y a pas d'erreur, c'est seulement que l'élément n'est pas ad hoc pour
// calculer la normale à ce noeud là
// temps: indique à quel moment on veut le calcul
// pour des éléments particulier (ex: SFE) la méthode est surchargée
int ElemMeca::CalculNormale_noeud(Enum_dure temps,const Noeud& noe,Coordonnee& coor)
{
int retour = 1; // init du retour : on part d'un bon a priori
Enum_type_geom enutygeom = Type_geom_generique(this->Id_geometrie());
int dima = ParaGlob::Dimension();
// if ((enutygeom == LIGNE) || (enutygeom == SURFACE))
if ((enutygeom != SURFACE) // dans une première étape on ne s'occupe que des surfaces en 3D
|| ((enutygeom == SURFACE) && (dima == 2)) // si surface en 2D, ce sera les éléments frontières qui détermineront les normales
)
{coor.Libere(); // pas de calcul possible
retour = 0;
}
else // sinon le calcul est possible
{ // on commence par repérer le noeud dans la numérotation locale
int nuoe=0;
int borne_nb_noeud=tab_noeud.Taille()+1;
for (int i=1;i< borne_nb_noeud;i++)
{Noeud& noeloc = *tab_noeud(i);
if ( (noe.Num_noeud() == noeloc.Num_noeud())
&& (noe.Num_Mail() == noeloc.Num_Mail())
)
{nuoe = i; break;
};
};
// on ne continue que si on a trouvé le noeud
if (nuoe != 0)
{ ElemGeomC0& elemgeom = ElementGeometrique(); // récup de la géométrie
// récup des coordonnées locales du noeuds
const Coordonnee& theta_noeud = elemgeom.PtelemRef()(nuoe);
// récup des phi et dphi au noeud
const Vecteur & phi = elemgeom.Phi_point(theta_noeud);
const Mat_pleine& dphi = elemgeom.Dphi_point(theta_noeud);
switch (temps)
{case TEMPS_0 :
{const BaseB& baseB = met->BaseNat_0(tab_noeud,dphi,phi);
coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
coor.Normer();
break;
}
case TEMPS_t :
{const BaseB& baseB = met->BaseNat_t(tab_noeud,dphi,phi);
coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
coor.Normer();
break;
}
case TEMPS_tdt :
{const BaseB& baseB = met->BaseNat_tdt(tab_noeud,dphi,phi);
coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
coor.Normer();
break;
}
default :
cout << "\nErreur : valeur incorrecte du temps demande = "
<< Nom_dure(temps) << " !\n";
cout << "\n ElemMeca::CalculNormale_noeud(Enum_dure temps... \n";
retour = 0;
Sortie(1);
};
}
else
{cout << "\n *** erreur le noeud demande num= "<<noe.Num_noeud()
<< " du maillage "<< noe.Num_Mail()
<< " ne fait pas parti de l'element "
<< num_elt << " du maillage "<< noe.Num_Mail()
<< " on ne peut pas calculer la normale au noeud !!"
<< "\n ElemMeca::CalculNormale_noeud(...";
retour = 0;
Sortie(1);
}
};
// retour
return retour;
};