rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/Elements/Thermique/ElemThermi5.cc

848 lines
38 KiB
C++
Raw Permalink Normal View History

intégration des éléments de thermique (actuellement uniquement biel 1D)
2021-09-28 17:51:09 +02:00
// 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.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
intégration des éléments de thermique (actuellement uniquement biel 1D)
2021-09-28 17:51:09 +02:00
// 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
mise à jour 7.010
2023-05-03 17:23:49 +02:00
const Vecteur & phi = elemgeom.Phi_point(theta_noeud);
const Mat_pleine& dphi = elemgeom.Dphi_point(theta_noeud);
intégration des éléments de thermique (actuellement uniquement biel 1D)
2021-09-28 17:51:09 +02:00
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;
};
Reference in a new issue Copy permalink