// 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 "CharUtil.h"
// actualisation des ddl et des grandeurs actives de t+dt vers t
// appelé par les classes dérivées
void ElemThermi::TdtversT_()
{// on met à jour l'indicateur de premier calcul
// s'il y a des sauvegardes de grandeur aux déformations
// on ne regarde que le premier élément de tableau, a priori
// il y a toujours un pt d'integ et l'organisation est identique pour tous les pt d'integ
if (tabSaveDefDon(1) != NULL)
premier_calcul_thermi_impli_expli=false;
// cas des énergies
int nbi= tab_energ.Taille();
for (int ni=1;ni<= nbi; ni++)
tab_energ_t(ni) = tab_energ(ni);
E_elem_bulk_t = E_elem_bulk_tdt; // énergie due au bulk viscosity
};
// actualisation des ddl et des grandeurs actives de t vers tdt
// appelé par les classes dérivées
void ElemThermi::TversTdt_()
{// on met à jour l'indicateur de premier calcul
// on considère que si l'on revient en arrière, il vaut mieux re-initialiser les
// grandeurs correspondantes au premier_calcul
// s'il y a des sauvegardes de grandeur aux déformations
if (tabSaveDefDon(1) != NULL)
premier_calcul_thermi_impli_expli=true;
// cas des énergies
int nbi= tab_energ.Taille();
for (int ni=1;ni<= nbi; ni++)
tab_energ(ni) = tab_energ_t(ni);
E_elem_bulk_tdt = E_elem_bulk_t; // énergie due au bulk viscosity
};
// calcul du résidu et de la matrice de raideur pour le calcul d'erreur
// cas d'une intégration suivant une seule liste
void ElemThermi::FluxAuNoeud_ResRaid(const int nbne,const Tableau <Vecteur>& taphi
,const Vecteur& poids,Tableau <Vecteur *>& resErr,Mat_pleine& raidErr
,const Tableau <Vecteur>& taphiEr,const Vecteur& poidsEr)
{ int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
// // dimension des tenseurs
// int dim = lesPtIntegThermiInterne->DimTens();
// // inialisation du second membre et de la raideur
// int nbSM = resErr.Taille(); // nombre de second membre
// for (int j =1; j<= nbSM; j++) // boucle sur les seconds membres
// (*resErr(j)).Zero();
// raidErr.Zero();
// // Il faut déterminer l'ordre dans lequel on parcours les contraintes qui doit
// // être compatible avec l'ordre des ddl
// Tableau2 <int> ordre = OrdreContrainte(nbSM);
// // création d'un tenseur au dimension absolu pour le calcul des contraintes
// // dans la base absolue
// TenseurHH& sigHH = *(NevezTenseurHH(dimAbsolue)) ;
// // ====== calcul du second membre =======
// int ni; // compteur globale de point d'integration
// bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
// // donc il faut calculer tous les éléments de la métrique
// for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
// {PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
// // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
// // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
// const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
// // passage dans le repère absolu du tenseur contrainte final
// sigHH = (*(ptIntegThermi.SigHH_t())).BaseAbsolue(sigHH,*ex.giB_t);
// for (int ne =1; ne<= nbne; ne++) // 1ere boucle sur les noeuds
// { // les résidus : mais il faut suivre l'ordre de l'enregistrement des ddl
// for (int itot = 1; itot<= nbSM; itot++)
// { int ix = (int) (ordre(itot,1)); int iy = (int) (ordre(itot,2));
// (*resErr(itot))(ne) += taphi(ni)(ne)*sigHH(ix,iy) * (poids(ni) * (*ex.jacobien_t));
// }
// }
// }
// // ====== calcul de la raideur c'est à dire du hessien ========
// // boucle sur les pt d'integ spécifiques à l'erreur
// for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
// {
// // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
// // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
// const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
// for (int ne =1; ne<= nbne; ne++) // 1ere boucle sur les noeuds
// for (int me =1; me<= nbne; me++) // 2ere boucle sur les noeuds
// raidErr(ne,me) += taphiEr(ni)(ne) * taphiEr(ni)(me) * poidsEr(ni) * (*ex.jacobien_t);
// }
// // liberation des tenseurs intermediaires
// TenseurHH * ptsigHH = &sigHH; delete ptsigHH;
// LibereTenseur();
// // calcul de l'erreur relative
};
// calcul de l'erreur sur l'élément. Ce calcul n'est disponible
// qu'une fois la remontée aux contraintes effectuées sinon aucune
// action. En retour la valeur de l'erreur sur l'élément
// = 1 : erreur = (int (delta sigma):(delta sigma) dv)/(int sigma:sigma dv)
// le numerateur et le denominateur sont tel que :
// errElemRelative = numerateur / denominateur , si denominateur different de 0
// sinon denominateur = numerateur si numerateur est different de 0, sinon
// tous sont nuls mais on n'effectue pas la division , les autres variables sont spécifiques
// a l'element.
void ElemThermi::Cal_ErrElem(int type,double& errElemRelative,double& numerateur
, double& denominateur,const int nbne,const Tableau <Vecteur>& taphi
,const Vecteur& poids,const Tableau <Vecteur>& taphiEr,const Vecteur& poidsEr)
{ int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
// // dimension des tenseurs
// int dim = lesPtIntegThermiInterne->DimTens();
// bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
// // donc il faut calculer tous les éléments de la métrique
// switch (type)
// {case 1 : // cas du calcul aux moindres carrés
//
// {// création d'un tenseur au dimension absolu pour le calcul des contraintes
// // dans la base absolue, on le choisit HB pour le double produit contracté
// // mais en absolu la variance n'a pas d'importance
// TenseurHB& sigHB = *(NevezTenseurHB(dimAbsolue)) ;
// // idem au point d'intégration et un tenseur nul pour l'initialisation
// TenseurHB& signiHB = *(NevezTenseurHB(dimAbsolue)) ;
// TenseurHB& sigiHB = *(NevezTenseurHB(dimAbsolue)) ; // tenseur de travail
// TenseurHB& nulHB = *(NevezTenseurHB(dimAbsolue)) ;
// // ====== calcul des termes de l'erreur =======
// // ---- tout d'abord on parcourt les points d'intégration de la mécanique
// numerateur = 0.; denominateur = 0.; // initialisation
// int ni; // compteur globale de point d'integration
// for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
// {PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(ni);
// // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
// // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
// const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
// // passage dans le repère absolu du tenseur contrainte initiale
// sigHB = (*(ptIntegThermi.SigHH_t())).BaseAbsolue(sigHB,*ex.giB_t);
// // calcul du denominateur
// denominateur += (sigHB && sigHB) * (poids(ni) * (*ex.jacobien_t));
// // calcul de la premiere partie du numerateur, celle qui dépend des points d'intégration
// // mécanique.
// // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds,
// signiHB = nulHB; // initialisation
// for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds
// signiHB += taphi(ni)(ne) * ((tab_noeud(ne))->Contrainte(sigiHB));
// // 2) intégrale de la partie dépendant de ni
// numerateur += denominateur - 2 * (signiHB && sigHB) * (poids(ni) * (*ex.jacobien_t)) ;
// }
// // ---- on parcourt maintenant les points d'intégration pour le calcul d'erreur
// // ---- ce qui permet de calculer la deuxième partie du numérateur
// // boucle sur les pt d'integ spécifiques à l'erreur
// for (defEr->PremierPtInteg(), ni = 1;defEr->DernierPtInteg();defEr->NevezPtInteg(),ni++)
// {
// // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
// // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
// const Met_abstraite::Expli& ex = defEr->Cal_explicit_t(premier_calcul);
// // 1) calcul au point d'intégration du tenseur des contraintes défini aux noeuds,
// signiHB = nulHB; // initialisation
// for (int ne =1; ne<= nbne; ne++) // boucle sur les noeuds
// signiHB += taphiEr(ni)(ne)*(tab_noeud(ne))->Contrainte(sigiHB);
// // 2) intégrale de la partie dépendant de ni
// numerateur += (signiHB && signiHB) * poidsEr(ni) * (*ex.jacobien_t);
// }
// // liberation des tenseurs intermediaires
// TenseurHB * ptsigHB = &sigHB; delete ptsigHB;
// ptsigHB = &signiHB; delete ptsigHB;ptsigHB = &sigiHB; delete ptsigHB;
// ptsigHB = &nulHB; delete ptsigHB;
// LibereTenseur();
// // enregistrement de l'erreur pour l'élément
// if (sigErreur == NULL) sigErreur = new double;
// *sigErreur = numerateur;
// break;
// }
// default :
// { cout << "\nErreur : valeur incorrecte du type de calcul d'erreur !\n";
// cout << "ElemThermi::ErreurElement(int type , type = \n" << type;
// Sortie(1);
// }
// };
// // --- calcul de l'erreur relative
// if (denominateur <= ConstMath::trespetit)
// // cas d'un champ de contraintes nulles initialement
// if (numerateur <= ConstMath::trespetit)
// // cas également du numérateur nul
// errElemRelative = 0.;
// else
// // on fait denominateur = numérateur -> erreur relative = 1
// errElemRelative = 1.;
// else
// // cas du dénominateur non nul
// errElemRelative = numerateur / denominateur;
};
// calcul de l'erreur aux noeuds. Contrairement au cas des contraintes
// seul le résidu est calculé. Cas d'une intégration suivant une liste
void ElemThermi::Cal_ErrAuxNoeuds(const int nbne, const Tableau <Vecteur>& taphi,
const Vecteur& poids,Tableau <Vecteur *>& resErr )
{ int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
// // inialisation du second membre, on utilise le premier vecteur uniquement
// (*resErr(1)).Zero();
// // ====== calcul du second membre =======
// int ni; // compteur globale de point d'integration
// double volume = 0.; // le volume de l'élément
// bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
// // donc il faut calculer tous les éléments de la métrique
// for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
// {
// // calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
// // que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
// const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
// for (int ne =1; ne<= nbne; ne++) // 1ere boucle sur les noeuds
// (*resErr(1))(ne) += taphi(ni)(ne)*(*sigErreur) * (poids(ni) * (*ex.jacobien_t));
// // calcul du volume
// volume += (poids(ni) * (*ex.jacobien_t));
// }
// // on relativise par rapport au volume, car initialement sigErreur représente l'erreur
// // totale sur tout l'élément. Donc on divise par le volume pour retrouver après résolution
// // une valeur au noeud qui représente une valeur ponctuelle et non une valeur qui
// // qui est relative à un volume
// *resErr(1) /= volume;
//// *resErr(1) = -sigErreur;
// LibereTenseur();
// // calcul de l'erreur relative
};
// ajout des ddl relatif aux contraintes pour les noeuds de l'élément
void ElemThermi::Ad_ddl_Flux(const DdlElement& tab_ddlErr)
{ int nbne = tab_noeud.Taille();
// for (int ne=1; ne<= nbne; ne++) // pour chaque noeud
// { // création du tableau de ddl
// int tab_ddlErr_Taille = tab_ddlErr(ne).tb.Taille(); // nb de ddl du noeud Thermi
// Tableau <Ddl> ta(tab_ddlErr_Taille);
// for (int i =1; i<= tab_ddlErr_Taille; i++)
// ta(i) = Ddl (tab_ddlErr(ne).tb(i),0.,LIBRE);
// // ajout du tableau dans le noeud
// tab_noeud(ne)->PlusTabDdl(ta);
// }
};
// inactive les ddl du problème primaire de mécanique
void ElemThermi::Inact_ddl_primaire(DdlElement& tab_ddl)
{ int nbne = tab_noeud.Taille();
for (int ne=1; ne<= nbne; ne++)
tab_noeud(ne)->Met_hors_service(tab_ddl(ne).tb);
};
// active les ddl du problème primaire de mécanique
void ElemThermi::Act_ddl_primaire(DdlElement& tab_ddl)
{ int nbne = tab_noeud.Taille();
for (int ne=1; ne<= nbne; ne++)
tab_noeud(ne)->Met_en_service(tab_ddl(ne).tb);
};
// inactive les ddl du problème de recherche d'erreur : les contraintes
void ElemThermi::Inact_ddl_Flux(DdlElement& tab_ddlErr)
{ int nbne = tab_noeud.Taille();
for (int ne=1; ne<= nbne; ne++)
tab_noeud(ne)->Met_hors_service(tab_ddlErr(ne).tb);
};
// active les ddl du problème de recherche d'erreur : les contraintes
void ElemThermi::Act_ddl_Flux(DdlElement& tab_ddlErr)
{ int nbne = tab_noeud.Taille();
for (int ne=1; ne<= nbne; ne++)
tab_noeud(ne)->Met_en_service(tab_ddlErr(ne).tb);
};
// active le premier ddl du problème de recherche d'erreur : FLUX1
void ElemThermi::Act_premier_ddl_Flux()
{ int nbne = tab_noeud.Taille();
for (int ne=1; ne<= nbne; ne++)
tab_noeud(ne)->Met_en_service(FLUXD1);
};
// retourne un tableau de ddl element, correspondant à la
// composante de densité de flux -> FLUXD1, pour chaque noeud qui contiend
// des ddl de contrainte
// -> utilisé pour l'assemblage de la raideur d'erreur
DdlElement& ElemThermi::Tableau_de_Flux1() const
{ cout << "\n erreur, fonction non defini pour cette element "
<< "\n ElemThermi::Tableau_de_Flux1()" << endl;
Sortie(1);
DdlElement * toto = new DdlElement();
return *toto;
};
// lecture des flux sur le flot d'entrée
void ElemThermi::LectureDesFlux(bool cas,UtilLecture * entreePrinc,Tableau <CoordonneeH *>& tabfluxH)
{ // dimensionnement de la metrique identique au cas d'un calcul explicite
if( cas)
{ Tableau<Enum_variable_metrique> tab(7);
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
met->PlusInitVariables(tab) ;
};
int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
// dimension des flux
int dim = (*(tabfluxH(1))).Dimension();
// création d'un Coordonnee au dimension absolu pour récupérer les flux
// dans la base absolue
Coordonnee flux(dimAbsolue) ;
int ni; // compteur globale de point d'integration
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
// donc il faut calculer tous les éléments de la métrique
for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
{ // calcul des éléments de la métrique, on utilise le même calcul
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
for (int i=1;i<= dimAbsolue;i++)
// récupération des coordonnées du flux en absolu
*(entreePrinc->entree) >> flux(i);
// passage dans le repère local du flux final
CoordonneeH & fluxH = *tabfluxH(ni); // pour simplifier
fluxH.Change_dim(dim);
for (int i = 1;i<=dim;i++)
fluxH(i) = ex.giH_t->Coordo(i) * flux;
}
};
// retour des flux en absolu retour true si ils existes sinon false
void ElemThermi::FluxEnAbsolues
(bool cas,Tableau <CoordonneeH *>& tabfluxH,Tableau <Vecteur>& tabflux)
{ // dimensionnement de la metrique identique au cas d'un calcul explicite
if( cas)
{ Tableau<Enum_variable_metrique> tab(7);
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
met->PlusInitVariables(tab) ;
};
int dimAbsolue = ParaGlob::Dimension(); // dimension absolue
// dimension des tenseurs
int dim = (*(tabfluxH(1))).Dimension();
// redimensionnement éventuel du tableau de sortie
int nbi = tabfluxH.Taille();
if (tabflux.Taille() != nbi)
tabflux.Change_taille(nbi);
for (int ni=1;ni<= nbi;ni++)
if ( tabflux(ni).Taille() != dimAbsolue)
tabflux(ni).Change_taille(dimAbsolue);
int ni; // compteur globale de point d'integration
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
// donc il faut calculer tous les éléments de la métrique
for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
{ // calcul des éléments de la métrique, on utilise le même calcul
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
const Met_abstraite::Expli& ex = def->Cal_explicit_t(premier_calcul);
// passage dans le repère global du flux local
Coordonnee flux; // init à 0
for (int i=1;i<= dim;i++)
flux += (*tabfluxH(ni))(i) * (ex.giB_t->Coordo(i));
};
};
// ---------------- lecture écriture dans base info ----------------
// programmes utilisés par les classes dérivées
// cas donne le niveau de la récupération
// = 1 : on récupère tout
// = 2 : on récupère uniquement les données variables (supposées comme telles)
void ElemThermi::Lecture_bas_inf
(istream& ent,const Tableau<Noeud *> * tabMaillageNoeud,const int cas)
{ // appel de la routine d'élément
Element::Lect_bas_inf_element(ent,tabMaillageNoeud,cas);switch (cas)
{ case 1 : // ------- on récupère tout -------------------------
{ string toto,nom;
// récup de la masse volumique
ent >> toto >> masse_volumique ;
// données thermique
ent >> toto >> dilatation;
// blocage éventuelle d'hourglass
ent >> toto >> nom;
type_stabHourglass=Id_Nom_StabHourglass(nom.c_str());
break;
}
case 2 : // ----------- lecture uniquement de se qui varie --------------------
{ break;
}
default :
{ cout << "\nErreur : valeur incorrecte du type de lecture !\n";
cout << "ElemThermi::Lecture_bas_inf (istream& ent,const int cas)"
<< " cas= " << cas << endl;
Sortie(1);
}
}
// ------ lecture dans tous les cas -------
// résultat d'erreur
string toto;
ent >> toto;
if (toto == "erreur_de_densite_flux")
{ if (fluxErreur == NULL)
fluxErreur = new double;
ent >> (*fluxErreur) ;
}
// données particulière pour les lois de comportement mécanique
int tabSaveDonTaille = tabSaveDon.Taille();
if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) ent >> toto;
int num ; // numéro du pt d'integ, n'est pas vraiment utilisé mais c'est mieux que de lire un string
for (int i=1; i<= tabSaveDonTaille; i++)
if (tabSaveDon(i) != NULL)
{ ent >> toto >> num;
tabSaveDon(i)->Lecture_base_info(ent,cas);}
// données particulière pour les lois de comportement thermo physique
int tabSaveTPTaille = tabSaveTP.Taille();
if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) ent >> toto;
for (int i=1; i<= tabSaveTPTaille; i++)
if (tabSaveTP(i) != NULL)
{ ent >> toto >> num;
tabSaveTP(i)->Lecture_base_info(ent,cas);}
// données particulière pour la déformation mécanique
int tabSaveDefDonTaille = tabSaveDefDon.Taille();
if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL)) ent >> toto;
for (int i=1; i<= tabSaveDefDonTaille; i++)
if (tabSaveDefDon(i) != NULL)
{ ent >> toto >> num;
tabSaveDefDon(i)->Lecture_base_info(ent,cas);}
// lecture des énergies
int energ_taille = tab_energ_t.Taille();
for (int i=1;i<= energ_taille;i++) ent >> tab_energ_t(i) ;
// // énergie et puissance éventuelle de la partie bulk viscosity
// ent >> toto >> E_elem_bulk_t >> toto >> P_elem_bulk;
// E_elem_bulk_tdt = E_elem_bulk_t;
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void ElemThermi::Ecriture_bas_inf(ostream& sort,const int cas)
{ // appel de la routine d'élément
Element::Ecri_bas_inf_element(sort,cas);
// en fait ici on sauvegarde la même chose dans tous les cas, par contre la sortie
// totale est documentée.
switch (cas)
{ case 1 : // ------- on sauvegarde tout -------------------------
{ // écriture de la masse volumique,
sort << "masse_volumique " << masse_volumique <<" " << "\n";
// données thermique
sort << "dilatation_thermique " << dilatation << " ";
// blocage éventuelle d'hourglass
sort << "\n hourglass: " << Nom_StabHourglass(type_stabHourglass) << " ";
break;
}
case 2 : // ----------- sauvegarde uniquement de se qui varie --------------------
{ break;
}
default :
{ cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
cout << "ElemThermi::Ecriture_bas_inf(ostream& sort,const int cas)"
<< " cas= " << cas << endl;
Sortie(1);
}
}
// --- informations à sortir dans tous les cas ----
// résultat d'erreur
if (fluxErreur != NULL)
sort << "erreur_de_densite_de_flux " << (*fluxErreur) <<" " << "\n";
else
sort << "pas_d'erreur_de_densite_de_flux \n";
// données particulière pour les lois de comportement mécaniques
int tabSaveDonTaille = tabSaveDon.Taille();
if ((tabSaveDonTaille != 0) && (tabSaveDon(1) != NULL)) sort << " data_spec_loi_comp_meca ";
for (int i=1; i<= tabSaveDonTaille; i++)
if (tabSaveDon(i) != NULL)
{ if (i==1) {sort << "pt_int= " <<i << " ";} else {sort << "\n pt_int= " <<i << " ";};
tabSaveDon(i)->Ecriture_base_info(sort,cas);}
// données particulière pour les lois de comportement thermo physiques
int tabSaveTPTaille = tabSaveTP.Taille();
if ((tabSaveTPTaille != 0) && (tabSaveTP(1) != NULL)) sort << " data_spec_loi_comp_ThemoPhysique ";
for (int i=1; i<= tabSaveTPTaille; i++)
if (tabSaveTP(i) != NULL)
{ if (i==1) {sort << "pt_int= " <<i << " ";} else {sort << "\n pt_int= " <<i << " ";};
tabSaveTP(i)->Ecriture_base_info(sort,cas);}
// données particulière pour la déformation mécanique
int tabSaveDefDonTaille = tabSaveDefDon.Taille();
if ((tabSaveDefDonTaille != 0) && (tabSaveDefDon(1) != NULL)) sort << " data_spec_def ";
for (int i=1; i<= tabSaveDefDonTaille; i++)
if (tabSaveDefDon(i) != NULL)
{ if (i==1) {sort << "pt_int= " <<i << " ";} else {sort << "\n pt_int= " <<i << " ";};
tabSaveDefDon(i)->Ecriture_base_info(sort,cas);}
// sortie des énergies
int energ_taille = tab_energ_t.Taille();
sort << "\n ";
for (int i=1;i<= energ_taille;i++) sort << tab_energ_t(i) << " " ;
// // énergie et puissance éventuelle de la partie bulk viscosity
// sort << "\n E_el_bulk= " << E_elem_bulk_t << " P_el_bulk= " << P_elem_bulk;
};
// calcul de la longueur d'arrête de l'élément minimal
// divisé par la célérité la plus rapide dans le matériau
// appelé par les classes dérivées
// nb_noeud : =0 indique que l'on utilise tous les noeuds du tableau de noeuds
// = un nombre > 0, indique le nombre de noeuds à utiliser au début du tableau
double ElemThermi::Interne_Long_arrete_mini_sur_c(Enum_dure temps,int nb_noeud)
{ int nbne = tab_noeud.Taille(); // récup du nombre de noeud
if (nb_noeud != 0)
nbne = nb_noeud;
// tout d'abord l'objectif est de déterminer la distance minimum
// entre les différents noeuds
// initialisation de la distance entre les deux noeuds
double dist = ConstMath::tresgrand;
for (int i=1;i<= nbne;i++)
// on itère sur les noeuds restants
switch (temps)
{ case TEMPS_0:
for (int j=i+1;j<= nbne;j++)
{ double dist_new = (tab_noeud(i)->Coord0() - tab_noeud(j)->Coord0()).Norme();
if (dist_new < dist) dist = dist_new;
}
break;
case TEMPS_t:
for (int j=i+1;j<= nbne;j++)
{ double dist_new = (tab_noeud(i)->Coord1() - tab_noeud(j)->Coord1()).Norme();
if (dist_new < dist) dist = dist_new;
}
break;
case TEMPS_tdt:
for (int j=i+1;j<= nbne;j++)
{ double dist_new = (tab_noeud(i)->Coord2() - tab_noeud(j)->Coord2()).Norme();
if (dist_new < dist) dist = dist_new;
}
break;
default :
cout << "\n cas du temps non implante temps= " << Nom_dure(temps)
<< "\n ElemThermi::Interne_Long_arrete_mini_sur_c(Enum_dure temps)";
Sortie(1);
}
// traitement d'une erreur éventuelle
if (dist <= ConstMath::petit)
{ cout << "\n **** ERREUR une longueur d'arrete de l'element est nulle"
<< "\n ElemThermi::Interne_Long_arrete_mini_sur_c(..."
<< "\n element: "; this->Affiche(1);
#ifdef MISE_AU_POINT
cout << "\n *** version mise au point: on continue neanmoins avec une longueur "
<< " arbitrairement tres petite (" <<ConstMath::petit <<") ";
#else
Sortie(1);
#endif
};
// on calcul la célérité
double cc = sqrt(ElemThermi::Calcul_maxi_E_qui(temps)/masse_volumique);
// calcul du résultat
double l_sur_c = dist/cc;
return l_sur_c;
};
// initialisation pour le calcul de la matrice masse dans le cas de l'algorithme
// de relaxation dynamique avec optimisation en continu de la matrice masse
// casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse
void ElemThermi::InitCalculMatriceMassePourRelaxationDynamique(int )
{ // on active le calcul des invariants de contraintes
int taille = lesPtIntegThermiInterne->NbPti();
// on définit le ddl étendu correspondant, s'il n'existe pas
// on définit le Ddl_enum_etendu correspondant à la masse au noeud pour la méthode
if (masse_relax_dyn.Nom_vide() )
masse_relax_dyn = (Ddl_enum_etendu("masse_relax_dyn"));
};
// phase de calcul de la matrice masse dans le cas de l'algo de relaxation dynamique
// mi= lambda * delta_t**2 / 2 * (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises)
// avec delta_t=1 par défaut a priori, donc n'intervient pas ici
// ep: epaisseur, K module de compressibilite, mu: module de cisaillement, Isig trace de sigma,
// Sig_mises la contrainte de mises
// casMass_relax: permet de choisir entre différentes méthodes de calcul de la masse
void ElemThermi::CalculMatriceMassePourRelaxationDynamique
(const double& alph, const double& beta, const double & lambda
,const double & gamma,const double & theta, int casMass_relax)
{ // 1)
// on calcule la moyenne des grandeurs qui servent pour le calcul de la pseudo-masse que l'on distribue
// également sur les noeuds
double trace_moyenne= 0.; double mises_moyen= 0.;
double compress_moy=0.; double cisaille_moy = 0.;
int taille = lesPtIntegThermiInterne->NbPti();
cout << "\n stop *** pour l'instant la methode : ElemThermi::CalculMatriceMassePourRelaxationDynamique"
<< " n'est pas encore operationnelle " << endl;
Sortie(1);
// // on boucle sur les points d'intégrations
// for (int i= 1; i<= taille; i++)
// {PtIntegThermiInterne & ptIntegThermi = (*lesPtIntegThermiInterne)(i);
// const Vecteur & invariant = ptIntegThermi.SigInvar();
// trace_moyenne += invariant(1) ; // le premier invariant c'est la trace
// if (invariant.Taille() > 1)
// mises_moyen += sqrt((3.*invariant(2)-invariant(1)*invariant(1))/2.);
// compress_moy += ptIntegThermi.ModuleCompressibilite();
// cisaille_moy += ptIntegThermi.ModuleCisaillement();
// };
// --- tout d'abord la partie contrainte
// double ensemble = (val_propre_sup + val_propre_inf + val_cisaillement) / taille;
// on divise par taille, ce qui conduit à la moyenne sur les points d'intégration
double ensemble = (gamma*trace_moyenne /3.+ 0.5*theta*mises_moyen ) / taille;
// --- maintenant la partie raideur
ensemble += (alph * compress_moy + beta * cisaille_moy) / taille;
// --- pour avoir l'intégration totale on utilise le volume
// ensemble *= volume;
// --- et maintenant on distribue sur tous les noeuds de manière identique
// ne fonctionne que si tous les noeuds sont actifs !!! (ce qui est un cas particulier courant)
// la masse élémentaire pour chaque noeud
// double masse_elementaire = ensemble * lambda * 0.5 / tab_noeud.Taille();
// pour des membrannes on a:
// mi= lambda * delta_t**2 / 2 * epaisseur/4 * (alpha*K+beta*mu+gamma*Isig/3+theta/2*Sig_mises)
// avec delta_t=1 par défaut a priori, donc n'intervient pas ici
// d'où le facteur 0.125, que l'on garde également pour le volume
double masse_elementaire = ensemble * lambda * 0.125 ;
// prise en compte du type d'élément:
double ParNoeud = 1.; // n'est utilisé que par la méthode de Barnes classique
int nbn=tab_noeud.Taille();
switch (Type_geom_generique(ElementGeometrique_const().TypeGeometrie()))
{ case VOLUME : // par rapport à la surface, on utilise la racine cubique du volume
{ double long_caracteristique = pow(Volume(),1./3.);
masse_elementaire *= long_caracteristique;
// ParNoeud = volume / nombre de noeuds, donc c'est la partie du volume
// attribuée à chaque noeud
ParNoeud = Volume() / nbn;
break;
}
case SURFACE : // l'algorithme historique
{ double epaisseur_moyenne = EpaisseurMoyenne(TEMPS_tdt);
masse_elementaire *= epaisseur_moyenne;
// ParNoeud = section / nombre de noeuds, donc c'est la partie de la section
// attribuée à chaque noeud
ParNoeud = Volume() / epaisseur_moyenne/ nbn;
break;
}
case LIGNE :
{ double section_moyenne = SectionMoyenne(TEMPS_tdt);
double long_caracteristique = Volume()/section_moyenne;
masse_elementaire *= long_caracteristique;
// ParNoeud = longueur / nombre de noeuds, donc c'est la partie de la longueur
// attribuée à chaque noeud
ParNoeud = Volume() / section_moyenne/ nbn;
break;
}
default :
cout << "\nErreur : pour l'instant les types autres que volume, surface, ligne ne sont pas pris en compte dans la relaxation "
<< " dynamique selon barnes !\n";
cout << "\n ElemThermi::CalculMatriceMassePourRelaxationDynamique(.. \n";
Sortie(1);
};
// on parcours les noeuds de l'élément et on remplace le cas échéent, la valeur de la masse au noeud.
switch (casMass_relax)
{ case 1: case 3: // cas 1 : on cumule aux noeuds, cas 3 on fait la moyenne (pas effectué ici, on ne fait
// que préparer le travail pour le cas 3)
{for (int inoe=1;inoe<=nbn;inoe++)
{Noeud& noe = *tab_noeud(inoe);
// on cumule
noe.ModifDdl_etendu(masse_relax_dyn).Valeur() += masse_elementaire;
};
break;
}
case 2: // cas 2 : on prend la masse maxi
{for (int inoe=1;inoe<=nbn;inoe++)
{Noeud& noe = *tab_noeud(inoe);
const Ddl_etendu& masse_actuel =noe.DdlEtendue(masse_relax_dyn);
// dans le cas où la masse actuelle est plus petite on la remplace
if (masse_actuel.ConstValeur() <= masse_elementaire)
noe.ModifDdl_etendu(masse_relax_dyn).Valeur() = masse_elementaire;
};
break;
}
case 4: case 5: // on cumule (puis on moyennera pour le cas 4, autre part) et on divise par la section comme dans la formule de barnes
{for (int inoe=1;inoe<=nbn;inoe++)
{Noeud& noe = *tab_noeud(inoe);
const Ddl_etendu& masse_actuel =noe.DdlEtendue(masse_relax_dyn);
// on cumule
if (masse_actuel.ConstValeur() <= masse_elementaire)
// !!!! je pense que ce cas est très bizarre il faudra revoir ce que cela signifie.... !!
noe.ModifDdl_etendu(masse_relax_dyn).Valeur() += (masse_elementaire/ParNoeud);
};
break;
}
}; //-- fin du switch
};
// METHODES VIRTUELLES:
// recuperation des coordonnées du point de numéro d'ordre = iteg pour
// la grandeur enu
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
// utilisé par les classes dérivées
Coordonnee ElemThermi::CoordPtInt(Enum_dure temps,Enum_ddl enu,int iteg,bool& err)
{ Coordonnee ptret;err = false;
// récupération de l'élément géométrique correspondant à Enu
ElemGeomC0& ele = this->ElementGeometrie(enu);
// récupération du tableau des fonctions d'interpolations
const Tableau <Vecteur> & tabphi = ele.TaPhi();
if ((iteg < 0) || (iteg > tabphi.Taille()))
{ err = true;}
else
{ switch (temps)
{ case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
}
};
return ptret; // retour
};
// recuperation des coordonnées du point d'intégration numéro = iteg pour
// la face : face
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
Coordonnee ElemThermi::CoordPtIntFace(int face, Enum_dure temps,int iteg,bool& err)
{ Coordonnee ptret;err = false;
if (SurfExiste(face))
{// récupération de l'élément géométrique correspondant à la face
const ElFrontiere* elf = this->Frontiere_surfacique(face,false);
if (elf != NULL) // cas où la frontière existe déjà
{ const ElemGeomC0 & elegeom = elf->ElementGeometrique();
// récupération du tableau des fonctions d'interpolations
const Tableau <Vecteur> & tabphi = elegeom.TaPhi();
if ((iteg < 0) || (iteg > tabphi.Taille()))
{ err = true;}
else
{ switch (temps)
{ case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
};
};
}
}
else // cas où la frontière n'existe pas, on renvoie une erreur
// et on laisse la valeur par défaut pour ptret
{err = true;
};
return ptret; // retour
};
// recuperation des coordonnées du point d'intégration numéro = iteg pour
// la face : face
// temps: dit si c'est à 0 ou t ou tdt
// si erreur retourne erreur à true
Coordonnee ElemThermi::CoordPtIntArete(int arete, Enum_dure temps,int iteg,bool& err)
{ Coordonnee ptret;err = false;
if (AreteExiste(arete))
{// récupération de l'élément géométrique correspondant à l'arête
const ElFrontiere* elf = this->Frontiere_lineique(arete,false);
if (elf != NULL) // cas où la frontière existe déjà
{ const ElemGeomC0 & elegeom = elf->ElementGeometrique();
// récupération du tableau des fonctions d'interpolations
const Tableau <Vecteur> & tabphi = elegeom.TaPhi();
if ((iteg < 0) || (iteg > tabphi.Taille()))
{ err = true;}
else
{ switch (temps)
{ case TEMPS_0 : ptret = met->PointM_0(tab_noeud,tabphi(iteg)); break;
case TEMPS_t : ptret = met->PointM_t(tab_noeud,tabphi(iteg)); break;
case TEMPS_tdt : ptret = met->PointM_tdt(tab_noeud,tabphi(iteg)); break;
};
};
}
}
else // cas où la frontière n'existe pas, on renvoie une erreur
// et on laisse la valeur par défaut pour ptret
{err = true;
};
return ptret; // retour
};
// retourne le numero du pt d'ing le plus près ou est exprimé la grandeur enum
// temps: dit si c'est à 0 ou t ou tdt
// utilisé par les classes dérivées
int ElemThermi::PtLePlusPres(Enum_dure temps,Enum_ddl enu, const Coordonnee& M)
{ int iret;
// récupération de l'élément géométrique correspondant à Enu
ElemGeomC0& ele = this->ElementGeometrie(enu);
// récupération du tableau des fonctions d'interpolations
const Tableau <Vecteur> & tabphi = ele.TaPhi();
// on boucle sur les pt d'integ pour avoir le point le plus près
int tabphitaille = tabphi.Taille();
Coordonnee P; iret=1; double dist= ConstMath::tresgrand;
for (int ipt = 1;ipt<=tabphitaille;ipt++)
{ switch (temps)
{ case TEMPS_0 : P = met->PointM_0(tab_noeud,tabphi(ipt)); break;
case TEMPS_t : P = met->PointM_t(tab_noeud,tabphi(ipt)); break;
case TEMPS_tdt : P = met->PointM_tdt(tab_noeud,tabphi(ipt)); break;
}
double di=(M-P).Norme();
if (di < dist) { dist = di; iret = ipt;};
}
return iret; // retour
};
// ==== >>>> methodes virtuelles redéfini éventuellement dans les classes dérivées ============
// ramene l'element geometrique correspondant au ddl passé en paramètre
ElemGeomC0& ElemThermi::ElementGeometrie(Enum_ddl ddl) const
{ Enum_ddl enu = PremierDdlFamille(ddl);
switch (enu)
{ case X1 : return ElementGeometrique(); break;
case TEMP : return ElementGeometrique(); break;
case FLUXD1 : return ElementGeometrique(); break;
case GRADT1 : return ElementGeometrique(); break;
case DGRADT1 : return ElementGeometrique(); break;
case ERREUR : return ElementGeometrique(); break;
default :
{cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl)
<< "\n ElemThermi::ElementGeometrie(Enum_ddl ddl) " ;
Sortie(1);}
}
return ElementGeometrique(); // pour taire le compilo
};
// ramène le nombre de grandeurs génératrices pour un pt d'integ, correspondant à un type enuméré
// peut-être surchargé pour des éléments particuliers
int ElemThermi::NbGrandeurGene(Enum_ddl ddl) const
{ Enum_ddl enu = PremierDdlFamille(ddl);
int nbGG = 0.; // par défaut
switch (enu)
{ case TEMP : case FLUXD1 :case GRADT1 :case DGRADT1 :
{ // cas d'un calcul de thermique classique
switch (Type_geom_generique(id_geom))
{ case LIGNE : case POINT_G : nbGG = 1; break; // FLUXD1
case SURFACE : nbGG = 2; break; // FLUXD1, FLUXD2
case VOLUME : nbGG = 3; break; // FLUXD1 FLUXD2 FLUXD3
default :
cout << "\nErreur : cas non traite, id_geom= :" << Nom_type_geom(Type_geom_generique(id_geom))
<< "ElemThermi::NbGrandeurGene(.. \n";
Sortie(1);
};
break;
}
default :
{cout << "\n cas non prevu ou non encore implante: ddl= " << Nom_ddl(ddl)
<< "\n ElemThermi::NbGrandeurGene(Enum_ddl ddl) " ;
Sortie(1);}
};
return nbGG;
};
// modification de l'orientation de l'élément en fonction de cas_orientation
// =0: inversion simple (sans condition) de l'orientation
// si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration
// 1. si tous les jacobiens sont négatifs on change d'orientation
// 2. si tous les jacobiens sont positifs on ne fait rien
// 3. si certains jacobiens sont positifs et d'autres négatifs message
// d'erreur et on ne fait rien
// ramène true: s'il y a eu changement effectif, sinon false
bool ElemThermi::Modif_orient_elem(int cas_orientation)
{ // retour:
bool retour=false; // par défaut pas de changement
if (cas_orientation == 0) // cas où on inverse l'orientation sans condition particulière
{ // on change l'orientation de l'élément
retour = true;
int nbnoe = tab_noeud.Taille();
Tableau<Noeud *> tab_inter(tab_noeud); // on crée un tableau intermédiaire
// on récupère la numérotation locale inverse
const Tableau<int> tabi = ElementGeometrique().InvConnec();
// on met à jour le tableau actuel
for (int n=1;n<=nbnoe;n++)
tab_noeud(n)=tab_inter(tabi(n));
}
else
{ // si cas_orientation est diff de 0: on calcul le jacobien aux différents points d'intégration
// 1. si tous les jacobiens sont négatifs on change d'orientation
// 2. si tous les jacobiens sont positifs on ne fait rien
// 3. si certains jacobiens sont positifs et d'autres négatifs message
// d'erreur et on ne fait rien
int cas=1; // a priori tout est ok
// boucle sur les pt d'integ
for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg())
{ double jacobien_0 = def->JacobienInitial();
if (jacobien_0 < 0)
{if (cas ==1) // on a trouvé un jacobien négatif
{cas =0; }
}// si c'était positif --> négatif
else // cas positif
{if (cas == 0) // on a déjà changé
{ cas = 2; break;} // on sort de la boucle
};
};
// gestion de pb
if (cas == 2)
{ cout << "\n **** Attention **** element nb= "<< this->Num_elt() << " peut-etre trop distordu ?"
<< " pt d'integ Thermi positif et negatif !! ";
if (ParaGlob::NiveauImpression() > 7)
{// on va essayer de sortir plus d'info
int pti=1;cout << "\n les jacobiens initiaux \n";
// on sort les valeurs des jacobiens
for (def->PremierPtInteg();def->DernierPtInteg();def->NevezPtInteg(),pti++)
{ cout << " pti= "<<pti << ", J= "<< def->JacobienInitial();};
// les coordonnées des noeuds en suivant l'ordre de la table de connection locale
int nbnoe = tab_noeud.Taille();
cout << "\n les coordonnees de l'element \n";
for (int ino=1;ino<=nbnoe;ino++)
{ cout << "\n "<< ino << " " << tab_noeud(ino)->Coord0();
};
};
}
else if (cas == 0)
{ switch (cas_orientation)
{ case 1: // on change l'orientation de l'élément
{retour = true;
int nbnoe = tab_noeud.Taille();
Tableau<Noeud *> tab_inter(tab_noeud); // on crée un tableau intermédiaire
// on récupère la numérotation locale inverse
const Tableau<int> tabi = ElementGeometrique().InvConnec();
for (int n=1;n<=nbnoe;n++)
tab_noeud(n)=tab_inter(tabi(n));
break;
}
case -1: // on sort une information à l'écran
{ cout << "\n element nb= "<< this->Num_elt() << "jacobien negatif " ;
break;;
}
default:
cout << "\n erreur le cas : " << cas_orientation
<< " n'est pas actuellement pris en compte"
<< "\n ElemThermi::Modif_orient_elem(...";
Sortie(1);
};
};
};
// retour
return retour;
};