848 lines
38 KiB
C++
Executable file
848 lines
38 KiB
C++
Executable file
// 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 "ElemThermi.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"
|
|
|
|
|
|
// -------------------- 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 ElemThermi::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 ElemThermi::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) = (ElemThermi*) 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 ElemThermi::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) = (ElemThermi*) 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 ElemThermi::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) = (ElemThermi*) 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 ElemThermi::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);
|
|
break;
|
|
}
|
|
|
|
default :
|
|
cout << "\n*** Erreur : cas de gestop, d'hourglass non defini !\n";
|
|
cout << "\n ElemThermi::Cal_implicit_hourglass() \n";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// stabilisation pour un calcul implicit
|
|
void ElemThermi::Cal_implicit_hourglass()
|
|
{ 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));
|
|
|
|
//---- 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))
|
|
{// 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
|
|
(*(resraid1.raid)) -= (*(resraid2.raid));
|
|
(*(resraid1.res)) -= (*(resraid2.res));
|
|
// on tiend compte du facteur de modération
|
|
(*(resraid1.raid)) *= coefStabHourglass;
|
|
(*(resraid1.res)) *= coefStabHourglass;
|
|
// on sauvegarde la raideur
|
|
raid_hourglass_transitoire = new Mat_pleine((*(resraid1.raid)));
|
|
// 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
|
|
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));
|
|
// on tiend compte du facteur de modération
|
|
(*(resHourglass1)) *= coefStabHourglass;
|
|
// on met à jour le résidu de l'élément principal
|
|
(*residu) += (*(resHourglass1));
|
|
// pour la partie raideur: on met à jour la raideur de l'élément principal
|
|
(*raideur) += (*raid_hourglass_transitoire);
|
|
};
|
|
|
|
//---- 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 ElemThermi::Cal_implicit_hourglass() \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
};
|
|
|
|
// stabilisation pour un calcul explicit
|
|
void ElemThermi::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_NON_DEFINIE :
|
|
// on ne fait rien
|
|
break;
|
|
default :
|
|
cout << "\n*** Erreur : cas de gestop, d'hourglass non defini !\n";
|
|
cout << "\n ElemThermi::Cal_implicit_hourglass() \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
};
|
|
|
|
|
|
// récupération de l'energie d'hourglass éventuelle
|
|
double ElemThermi::Energie_Hourglass()
|
|
{ double enerHourglass = 0.;
|
|
switch (type_stabHourglass)
|
|
{case STABHOURGLASS_PAR_COMPORTEMENT :
|
|
{ // on récupère les énergies stockées à l'élément
|
|
const EnergieThermi& energieTotale = tab_elHourglass(1)->EnergieTotaleElement();
|
|
enerHourglass = coefStabHourglass * (energieTotale.EnergieElastique()
|
|
+ energieTotale.DissipationPlastique() + energieTotale.DissipationVisqueuse());
|
|
};
|
|
break;
|
|
};
|
|
return enerHourglass;
|
|
};
|
|
|
|
// 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 ElemThermi::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 ElemThermi::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 & ElemThermi::Frontiere_elethermi(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;
|
|
};
|
|
|
|
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 ElemThermi::Frontiere_points(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_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 ElemThermi::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 = posi_tab_front_point = 0;}
|
|
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 = posi_tab_front_point = tail_fa; // après les surfaces
|
|
};
|
|
};
|
|
// 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 ElemThermi::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 = el.NbFe(); // nombre potentiel de faces
|
|
if (num > tail_fa)
|
|
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 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 = posi_tab_front_point = tail_fa; // après les surfaces
|
|
};
|
|
};
|
|
// --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 é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 ElemThermi::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 ElemThermi::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 ElemThermi::CalculNormale_noeud(...";
|
|
retour = 0;
|
|
Sortie(1);
|
|
}
|
|
};
|
|
// retour
|
|
return retour;
|
|
};
|
|
|