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