rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/comportement/Loi_comp_abstraite.cc

4665 lines
244 KiB
C++
Raw Permalink Normal View History

ajout des lois de comportement
2021-09-23 11:21:15 +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)
ajout des lois de comportement
2021-09-23 11:21:15 +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 "Loi_comp_abstraite.h"
#include "Def_Umat.h"
#include "ExceptionsLoiComp.h"
#include "Util.h"
#include "TypeQuelconqueParticulier.h"
#include "CharUtil.h"
#include "TypeConsTens.h"
// CONSTRUCTEURS :
Loi_comp_abstraite::Loi_comp_abstraite () : // Constructeur par defaut
LoiAbstraiteGeneral(),saveResul(NULL),comp_tangent_simplifie(false)
,utilise_vitesse_deformation(false),type_de_deformation(DEFORMATION_STANDART)
,thermo_dependant(false),temperature_tdt(-1.),temperature_0(-1.),temperature_t(-1.)
,temperature(NULL),dilatation(false),def_en_cours(NULL)
,temps_loi()
,epsBB_totale(NULL),epsBB_therm(NULL)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
,delta_epsBB_totale(NULL),delta_epsBB_therm(NULL)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
,DepsBB_totale(NULL),DepsBB_therm(NULL),DepsBB_umat(NULL)
,listQuelc_mis_en_acces_localement(),listdeTouslesQuelc_dispo_localement()
,ptintmeca_en_cours(NULL),permet_affich_loi(0)
,permet_affich_loi_nD(NULL),li_quelconque(),tab_pt_li_quelconque()
{ };
// Constructeur utile si l'identificateur du nom de la loi
// de comportement et la dimension sont connus
Loi_comp_abstraite::Loi_comp_abstraite (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp
,int dimension,bool vit_def) :
LoiAbstraiteGeneral(id_compor,dimension,categorie_comp)
,saveResul(NULL),comp_tangent_simplifie(false)
,utilise_vitesse_deformation(vit_def),type_de_deformation(DEFORMATION_STANDART)
,thermo_dependant(false),temperature_tdt(-1.),temperature_0(-1.),temperature_t(-1.)
,temperature(NULL),dilatation(false),def_en_cours(NULL)
,temps_loi()
,epsBB_totale(NULL),epsBB_therm(NULL)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
,delta_epsBB_totale(NULL),delta_epsBB_therm(NULL)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
,DepsBB_totale(NULL),DepsBB_therm(NULL),DepsBB_umat(NULL)
,listQuelc_mis_en_acces_localement(),listdeTouslesQuelc_dispo_localement()
,ptintmeca_en_cours(NULL),permet_affich_loi(0)
,permet_affich_loi_nD(NULL),li_quelconque(),tab_pt_li_quelconque()
{ };
// Constructeur utile si l'identificateur du nom de la loi
// de comportement et la dimension sont connus
Loi_comp_abstraite::Loi_comp_abstraite (char* nom,Enum_categorie_loi_comp categorie_comp
,int dimension,bool vit_def) :
LoiAbstraiteGeneral(nom,dimension,categorie_comp)
,saveResul(NULL),comp_tangent_simplifie(false)
,utilise_vitesse_deformation(vit_def),type_de_deformation(DEFORMATION_STANDART)
,thermo_dependant(false),temperature_tdt(-1.),temperature_0(-1.),temperature_t(-1.)
,temperature(NULL),dilatation(false),def_en_cours(NULL)
,temps_loi()
,epsBB_totale(NULL),epsBB_therm(NULL)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
,delta_epsBB_totale(NULL),delta_epsBB_therm(NULL)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
,DepsBB_totale(NULL),DepsBB_therm(NULL),DepsBB_umat(NULL)
,listQuelc_mis_en_acces_localement(),listdeTouslesQuelc_dispo_localement()
,ptintmeca_en_cours(NULL),permet_affich_loi(0)
,permet_affich_loi_nD(NULL),li_quelconque(),tab_pt_li_quelconque()
{ };
// Constructeur de copie
Loi_comp_abstraite::Loi_comp_abstraite (const Loi_comp_abstraite & a ) :
LoiAbstraiteGeneral(a),comp_tangent_simplifie(false)
,utilise_vitesse_deformation(a.utilise_vitesse_deformation)
// ,saveResul(a.saveResul->Nevez_SaveResul()),type_de_deformation(a.type_de_deformation)
,saveResul(a.saveResul),type_de_deformation(a.type_de_deformation)
,thermo_dependant(a.thermo_dependant),dilatation(a.dilatation)
,def_en_cours(a.def_en_cours)
,temps_loi(a.temps_loi)
,temperature_tdt(-1.),temperature_0(-1.),temperature_t(-1.),temperature(NULL)
,epsBB_totale(NULL),epsBB_therm(NULL)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
,delta_epsBB_totale(NULL),delta_epsBB_therm(NULL)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
,DepsBB_totale(NULL),DepsBB_therm(NULL),DepsBB_umat(NULL)
,listQuelc_mis_en_acces_localement(a.listQuelc_mis_en_acces_localement)
,listdeTouslesQuelc_dispo_localement(a.listdeTouslesQuelc_dispo_localement)
,ptintmeca_en_cours(NULL),permet_affich_loi(a.permet_affich_loi)
,permet_affich_loi_nD(a.permet_affich_loi_nD)
,li_quelconque(a.li_quelconque),tab_pt_li_quelconque()
{ // idem pour les fonctions nD
if (permet_affich_loi_nD != NULL)
{if (permet_affich_loi_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
permet_affich_loi_nD = Fonction_nD::New_Fonction_nD(*a.permet_affich_loi_nD);
};
// on s'occupe maintenant des stockages des grandeurs complètement quelconques
// qui seront nécessaire pour utiliser la fonction nD
// !! on considère que a.li_quelconque a bien été dimensionné en fonction de
// permet_affich_loi_nD, donc on ne fait que bien définir les pointeurs
if (li_quelconque.size())
{tab_pt_li_quelconque.Change_taille(li_quelconque.size());
List_io <TypeQuelconque >::iterator ili,ilifin = li_quelconque.end();
int i=1; // init
for (ili= li_quelconque.begin();ili != ilifin; ili++,i++)
tab_pt_li_quelconque(i) = &(*ili);
};
};
} ;
// DESTRUCTEUR VIRTUEL :
Loi_comp_abstraite::~Loi_comp_abstraite ()
{ if (epsBB_totale != NULL) delete epsBB_totale;
if (epsBB_therm != NULL) delete epsBB_therm;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (delta_epsBB_totale != NULL) delete delta_epsBB_totale;
if (delta_epsBB_therm != NULL) delete delta_epsBB_therm;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (DepsBB_totale != NULL) delete DepsBB_totale;
if (DepsBB_therm != NULL) delete DepsBB_therm;
if (DepsBB_umat != NULL) delete DepsBB_umat;
if (permet_affich_loi_nD != NULL)
if (permet_affich_loi_nD->NomFonction() == "_") delete permet_affich_loi_nD;
};
// =============== methode de calcul ======================
// définition du type de calcul de déformation sur une instance de déformation passée en paramètre
void Loi_comp_abstraite::Def_type_deformation(Deformation & def)
{ def.Change_type_de_deformation(type_de_deformation);
};
// schema de calcul explicite à t
const Met_abstraite::Expli& Loi_comp_abstraite::Cal_explicit_t
(Loi_comp_abstraite::SaveResul * saveDon
,Deformation & def,DdlElement & tab_ddl,PtIntegMecaInterne& ptintmeca,Tableau <TenseurBB *> & d_epsBB
,double& jacobien,CompThermoPhysiqueAbstraite::SaveResul * saveTP,CompThermoPhysiqueAbstraite* loiTP
,bool dilat,EnergieMeca & energ,const EnergieMeca & energ_t,bool premier_calcul
)
{
Temps_CPU_HZpp& temps_cpu_loi = ptintmeca.Tps_cpu_loi_comp();
temps_cpu_loi.Mise_en_route_du_comptage(); // spécifique pti
temps_loi.Mise_en_route_du_comptage(); // spécifique loi
// passage des infos specifiques aux classes derivantes
saveResul = saveDon;
def_en_cours = &def;
// récup des contraintes
TenseurHH & sigHH_t = *(ptintmeca.SigHH_t());
TenseurHH & sigHH = *(ptintmeca.SigHH());
TenseurBB & epsBB_t = *(ptintmeca.EpsBB());
TenseurBB & DepsBB_ = *(ptintmeca.DepsBB());
TenseurBB & DeltaEpsBB_ = *(ptintmeca.DeltaEpsBB());
double& module_compressibilite = ptintmeca.ModuleCompressibilite();
double& module_cisaillement = ptintmeca.ModuleCisaillement();
Tableau <double>& def_equi = ptintmeca.Deformation_equi();
const Tableau <double>& def_equi_t = ptintmeca.Deformation_equi_t_const();
dilatation=dilat;
ptintmeca_en_cours = &ptintmeca; // pour les méthodes internes
// calcul de : epsBB_t,gijBB_t, gijHH_t,d_epsBB,jacobien
Temps_CPU_HZpp& temps_cpu_metrique = ptintmeca.TpsMetrique();
temps_cpu_metrique.Mise_en_route_du_comptage(); // cpu
const Met_abstraite::Expli& ex = def.Cal_explicit_t(def_equi_t,epsBB_t,d_epsBB,def_equi,DepsBB_,DeltaEpsBB_,premier_calcul);
// on récupère les positions des points
if (premier_calcul)
{ptintmeca.M_0() = def.Position_0();
ptintmeca.M_t() = def.Position_t();
};
ptintmeca.M_tdt() = def.Position_tdt();
temps_cpu_metrique.Arret_du_comptage(); // cpu
// dans le cas d'une loi thermo dépendante, calcul de la température
if ((thermo_dependant) || dilatation)
{temperature_0 = def.DonneeInterpoleeScalaire(TEMP,TEMPS_0);
temperature_t = def.DonneeInterpoleeScalaire(TEMP,TEMPS_t);
temperature = &temperature_t;
};
// dans le cas ou l'on tiens compte de la dilatation on modifie les différents tenseurs de déformation
ThermoDonnee dTP; // init à 0 par défaut
if (dilatation)
{ // attention ici ça peut donner 0 car le pas précédent c'est peut-être 0
// comme on n'en sait rien on fait comme si il y a une valeur déjà enregistrée: pas terrible !!!
// vu que sigma suivant dépend de epsméca et que epsmeca dépend de sigma !!!
double P = -(sigHH_t * (*(ex.gijBB_t))).Trace()/Abs(sigHH_t.Dimension()); // calcul de la pression
loiTP->Cal_donnees_thermiques(P,saveTP,def,P,TEMPS_t,dTP);
// modification de la déformation mécanique et calcul de la déformation thermique
// a- dimensionnement des tenseurs intermédiaires qui pour l'instant ne servent qu'ici
int dim_tens = epsBB_t.Dimension();
// -- cas de la déformation
if (epsBB_totale == NULL) { epsBB_totale = NevezTenseurBB(dim_tens);}
else if (epsBB_totale->Dimension() != dim_tens) { delete epsBB_totale;epsBB_totale = NevezTenseurBB(dim_tens);};
*epsBB_totale = epsBB_t;
if (epsBB_therm == NULL) { epsBB_therm = NevezTenseurBB(dim_tens);}
else if (epsBB_therm->Dimension() != dim_tens) { delete epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// -- cas de l'incrément de déformation
if (delta_epsBB_totale == NULL) { delta_epsBB_totale = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_totale->Dimension() != dim_tens) { delete delta_epsBB_totale;delta_epsBB_totale = NevezTenseurBB(dim_tens);};
*delta_epsBB_totale = DeltaEpsBB_;
if (delta_epsBB_therm == NULL) { delta_epsBB_therm = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_therm->Dimension() != dim_tens) { delete delta_epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// -- cas de la vitesse de déformation
if (DepsBB_totale == NULL) { DepsBB_totale = NevezTenseurBB(dim_tens);}
else if (DepsBB_totale->Dimension() != dim_tens) { delete DepsBB_totale;DepsBB_totale = NevezTenseurBB(dim_tens);};
*DepsBB_totale = DepsBB_;
if (DepsBB_therm == NULL) { DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (DepsBB_therm->Dimension() != dim_tens) { delete DepsBB_therm;DepsBB_totale = NevezTenseurBB(dim_tens);};
// -- calcul du partage entre thermique et mécanique
bool avec_repercution_sur_def_meca = true;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// des pointeurs pour des sorties en absolue
const Met_abstraite::Impli* ex_impli= NULL;
const Met_abstraite::Expli_t_tdt ex_temp = ex.T_dans_tdt();
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &ex_temp;
const Met_abstraite::Umat_cont* ex_umat = NULL;
def.DeformationThermoMecanique(temperature_0,dTP,*epsBB_totale
,*epsBB_therm,temperature_t,epsBB_t
,*delta_epsBB_totale,*delta_epsBB_therm,DeltaEpsBB_
,temperature_t
,*DepsBB_totale,*DepsBB_therm,DepsBB_
,ex_impli,ex_expli_tdt,ex_umat
,false,avec_repercution_sur_def_meca);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
};
if ((thermo_dependant) || dilatation) // répercussion éventuelle du changement de température dans les classes dérivées
{ RepercuteChangeTemperature(TEMPS_t); };
// calcul éventuel des invariants liés à la cinématique
CalculInvariants_cinematique(ptintmeca,*(ex.gijBB_t),*(ex.gijHH_t));
// demande dans les classes dérivées du calcul de grandeurs spécifiques si besoin est
// opération de transmission de la métrique
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt inter_ex_tdt = ex.T_dans_tdt();
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &inter_ex_tdt;
const Met_abstraite::Umat_cont* ex_expli = NULL;
CalculGrandeurTravail(ptintmeca,def,TEMPS_t,dTP,ex_impli,ex_expli_tdt,ex_expli,NULL,NULL);
// calcul de : sigHH
// on utilise la fonction qui calcul à t+dt en mettant les grandeur à t de l'appel idem ceux à t+dt
// pour l'instant (si pb on verra)
// on utilise lafonction T_dans_tdt pour passer avec une grandeur Expli_t_tdt au lieu de Expli
// c'est clair que ce n'est pas très cohérent, pour l'instant on laisse car cette fonction ne sert pas
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
try
{ Calcul_SigmaHH (sigHH_t,DepsBB_,tab_ddl,*(ex.gijBB_t),*(ex.gijHH_t),*(ex.giB_t),*(ex.giH_t),
epsBB_t,DeltaEpsBB_,*(ex.gijBB_t),*(ex.gijHH_t),*(ex.d_gijBB_t),
(*ex.jacobien_0),(*ex.jacobien_t),sigHH,energ,energ_t
,module_compressibilite,module_cisaillement,ex.T_dans_tdt());
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... ) // cas d'une erreur dans la recherche de la contrainte
{ throw ErrNonConvergence_loiDeComportement();
if (ParaGlob::NiveauImpression() >= 1)
cout << "\n warning: exception generee par la loi de comportement ";
if (ParaGlob::NiveauImpression() >= 4) cout << "\n Loi_comp_abstraite::Cal_explicit_t(..";
};
jacobien = (*ex.jacobien_t);
// calcul éventuel des invariants de contraintes
CalculInvariants_contraintes(ptintmeca,*(ex.gijBB_t),*(ex.gijHH_t));
// en sortie il nous faut définir d_epsBB qui doit être égale à D*, donc pas forcément la dérivée
// de epsBB (par exemple dans le cas de la déformation logarithmique ou cumulée)
// variation de la déformation / au ddl
if (type_de_deformation != DEFORMATION_STANDART)
{ int d_epsBBTaille = d_epsBB.Taille();
for (int i=1; i<= d_epsBBTaille; i++)
*(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_t))(i)));
};
def_en_cours=NULL; // plus d'accès possible à la sortie
ptintmeca_en_cours = NULL; // plus d'accès possible à la sortie
temps_cpu_loi.Arret_du_comptage(); // cpu spécifique pti
temps_loi.Arret_du_comptage(); // cpu spécifique loi
return ex; // retour de la métrique
};
// schema de calcul explicite à tdt
const Met_abstraite::Expli_t_tdt& Loi_comp_abstraite::Cal_explicit_tdt
(Loi_comp_abstraite::SaveResul * saveDon,Deformation & def
,DdlElement & tab_ddl
,PtIntegMecaInterne& ptintmeca,Tableau <TenseurBB *> & d_epsBB,double& jacobien
,CompThermoPhysiqueAbstraite::SaveResul * saveTP,CompThermoPhysiqueAbstraite* loiTP
,bool dilat,EnergieMeca & energ,const EnergieMeca & energ_t,bool premier_calcul
)
{
Temps_CPU_HZpp& temps_cpu_loi = ptintmeca.Tps_cpu_loi_comp();
temps_cpu_loi.Mise_en_route_du_comptage(); // spécifique pti
temps_loi.Mise_en_route_du_comptage(); // spécifique loi
// passage des infos specifiques aux classes derivantes
saveResul = saveDon;
def_en_cours = &def;
// récup des contraintes
TenseurHH & sigHH = *(ptintmeca.SigHH()); // ici _tdt est la grandeur finale
TenseurHH & sigHH_t = *(ptintmeca.SigHH_t()); // ici _tdt est la grandeur finale
TenseurBB & epsBB_tdt = *(ptintmeca.EpsBB());
TenseurBB & DepsBB_ = *(ptintmeca.DepsBB());
TenseurBB & delta_epsBB = *(ptintmeca.DeltaEpsBB());
double& module_compressibilite = ptintmeca.ModuleCompressibilite();
double& module_cisaillement = ptintmeca.ModuleCisaillement();
Tableau <double>& def_equi = ptintmeca.Deformation_equi();
Tableau <double>& def_equi_t = ptintmeca.Deformation_equi_t();
dilatation=dilat;
ptintmeca_en_cours = &ptintmeca; // pour les méthodes internes
// calcul de : epsBB_t,gijBB_tdt, gijHH_tdt,d_epsBB,jacobien
Temps_CPU_HZpp& temps_cpu_metrique = ptintmeca.TpsMetrique();
temps_cpu_metrique.Mise_en_route_du_comptage(); // cpu
const Met_abstraite::Expli_t_tdt& ex = def.Cal_explicit_tdt(def_equi_t,epsBB_tdt,d_epsBB,def_equi,DepsBB_,delta_epsBB,premier_calcul);
// on récupère les positions des points
if (premier_calcul)
{ptintmeca.M_0() = def.Position_0();
ptintmeca.M_t() = def.Position_t();
};
ptintmeca.M_tdt() = def.Position_tdt();
temps_cpu_metrique.Arret_du_comptage(); // cpu
// dans le cas d'une loi thermo dépendante, calcul de la température
if ((thermo_dependant) || dilatation)
{temperature_0 = def.DonneeInterpoleeScalaire(TEMP,TEMPS_0);
temperature_t = def.DonneeInterpoleeScalaire(TEMP,TEMPS_t);
temperature_tdt = def.DonneeInterpoleeScalaire(TEMP,TEMPS_tdt);
temperature = &temperature_tdt;
};
// dans le cas ou l'on tiens compte de la dilatation on modifie les différents tenseurs de déformation
ThermoDonnee dTP; // init à 0 par défaut
if (dilatation)
{ // attention ici la seule valeur disponible est la contrainte du temps précédent
// donc pour des questions d'erreur de transport de la contrainte, on utilise également la métrique au pas précédent
double P = -(sigHH_t * (*(ex.gijBB_t))).Trace()/Abs(sigHH_t.Dimension()); // calcul de la pression
loiTP->Cal_donnees_thermiques(P,saveTP,def,P,TEMPS_tdt,dTP);
// modification de la déformation mécanique et calcul de la déformation thermique
// a- dimensionnement des tenseurs intermédiaires qui pour l'instant ne servent qu'ici
int dim_tens = epsBB_tdt.Dimension();
// -- cas de la déformation
if (epsBB_totale == NULL) { epsBB_totale = NevezTenseurBB(dim_tens);}
else if (epsBB_totale->Dimension() != dim_tens) { delete epsBB_totale;epsBB_totale = NevezTenseurBB(dim_tens);};
*epsBB_totale = epsBB_tdt;
if (epsBB_therm == NULL) { epsBB_therm = NevezTenseurBB(dim_tens);}
else if (epsBB_therm->Dimension() != dim_tens) { delete epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// -- cas de l'incrément de déformation
if (delta_epsBB_totale == NULL) { delta_epsBB_totale = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_totale->Dimension() != dim_tens) { delete delta_epsBB_totale;delta_epsBB_totale = NevezTenseurBB(dim_tens);};
*delta_epsBB_totale = delta_epsBB;
if (delta_epsBB_therm == NULL) { delta_epsBB_therm = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_therm->Dimension() != dim_tens) { delete delta_epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
// -- cas de la vitesse de déformation
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (DepsBB_totale == NULL) { DepsBB_totale = NevezTenseurBB(dim_tens);}
else if (DepsBB_totale->Dimension() != dim_tens) { delete DepsBB_totale;DepsBB_totale = NevezTenseurBB(dim_tens);};
*DepsBB_totale = DepsBB_;
if (DepsBB_therm == NULL) { DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (DepsBB_therm->Dimension() != dim_tens) { delete DepsBB_therm;DepsBB_totale = NevezTenseurBB(dim_tens);};
// -- calcul du partage entre thermique et mécanique
bool avec_repercution_sur_def_meca = true;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// des pointeurs pour des sorties en absolue
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt =&ex;
const Met_abstraite::Umat_cont* ex_umat = NULL;
def.DeformationThermoMecanique(temperature_0,dTP,*epsBB_totale
,*epsBB_therm,temperature_tdt,epsBB_tdt
,*delta_epsBB_totale,*delta_epsBB_therm,delta_epsBB
,temperature_t
,*DepsBB_totale,*DepsBB_therm,DepsBB_
,ex_impli,ex_expli_tdt,ex_umat
,true,avec_repercution_sur_def_meca);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
};
if ((thermo_dependant) || dilatation) // répercussion éventuelle du changement de température dans les classes dérivées
{ RepercuteChangeTemperature(TEMPS_tdt); };
// calcul éventuel des invariants liés à la cinématique
CalculInvariants_cinematique(ptintmeca,*(ex.gijBB_t),*(ex.gijHH_t));
// demande dans les classes dérivées du calcul de grandeurs spécifiques si besoin est
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &ex;
const Met_abstraite::Umat_cont* ex_expli = NULL;
CalculGrandeurTravail(ptintmeca,def,TEMPS_tdt,dTP,ex_impli,ex_expli_tdt,ex_expli,NULL,NULL);
// calcul de : sigHH
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
try
{ Calcul_SigmaHH (sigHH_t,DepsBB_,tab_ddl,*(ex.gijBB_t),*(ex.gijHH_t),*(ex.giB_tdt),*(ex.giH_tdt),
epsBB_tdt,delta_epsBB,*(ex.gijBB_tdt),*(ex.gijHH_tdt),*(ex.d_gijBB_tdt),
(*ex.jacobien_0),(*ex.jacobien_tdt),sigHH,energ,energ_t
,module_compressibilite,module_cisaillement,ex);
//--- debug
//cout << "\n sigHH.MaxiComposante() "<<sigHH.MaxiComposante();
//if (sigHH.MaxiComposante() > 20.)
// { cout << "\n debug .... Loi_comp_abstraite::Cal_explicit_tdt "
// << "\n **** sigHH.MaxiComposante() > 20. ";
//
// TenseurBB* ptBB = NevezTenseurBB(epsBB_tdt);
// epsBB_tdt.BaseAbsolue(*ptBB,*(ex.giH_tdt));
// cout << "\n eps: ";
// ptBB->Ecriture(cout);
// delta_epsBB.BaseAbsolue(*ptBB,*(ex.giH_tdt));
// cout << "\n delta eps: ";
// ptBB->Ecriture(cout);
// TenseurHH* ptHH = NevezTenseurHH(sigHH);
// sigHH.BaseAbsolue(*ptHH,*(ex.giB_tdt));
// cout << "\n sigma: ";
// ptHH->Ecriture(cout);
// delete ptBB; delete ptHH;
//
//
// // on recalcule pour le debug le vecteur résidu
// Calcul_SigmaHH (sigHH_t,DepsBB_,tab_ddl,*(ex.gijBB_t),*(ex.gijHH_t),*(ex.giB_tdt),*(ex.giH_tdt),
// epsBB_tdt,delta_epsBB,*(ex.gijBB_tdt),*(ex.gijHH_tdt),*(ex.d_gijBB_tdt),
// (*ex.jacobien_0),(*ex.jacobien_tdt),sigHH,energ,energ_t
// ,module_compressibilite,module_cisaillement,ex);
// };
//--- fin debug
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... ) // cas d'une erreur dans la recherche de la contrainte
{ throw ErrNonConvergence_loiDeComportement();
if (ParaGlob::NiveauImpression() >= 1)
cout << "\n warning: exception generee par la loi de comportement ";
if (ParaGlob::NiveauImpression() >= 4) cout << "\n Loi_comp_abstraite::Cal_explicit_tdt(..";
};
jacobien = (*ex.jacobien_tdt);
// calcul éventuel des invariants de contraintes
CalculInvariants_contraintes(ptintmeca,*(ex.gijBB_tdt),*(ex.gijHH_tdt));
// en sortie il nous faut définir d_epsBB qui doit être égale à D*, donc pas forcément la dérivée
// de epsBB (par exemple dans le cas de la déformation logarithmique ou cumulée)
// variation de la déformation / au ddl
if (type_de_deformation != DEFORMATION_STANDART)
{ int d_epsBBTaille = d_epsBB.Taille();
for (int i=1; i<= d_epsBBTaille; i++)
*(d_epsBB(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
};
def_en_cours=NULL; // plus d'accès possible à la sortie
ptintmeca_en_cours = NULL; // plus d'accès possible à la sortie
temps_cpu_loi.Arret_du_comptage(); // cpu spécifique pti
temps_loi.Arret_du_comptage(); // cpu spécifique loi
// retour
return ex; // retour de la métrique
};
// schema implicit
const Met_abstraite::Impli& Loi_comp_abstraite::Cal_implicit
(Loi_comp_abstraite::SaveResul * saveDon
,Deformation & def,DdlElement & tab_ddl,PtIntegMecaInterne & ptintmeca
, Tableau <TenseurBB *>& d_epsBB_tdt,double& jacobien,Vecteur& d_jacobien_tdt
,Tableau <TenseurHH *>& d_sigHH,const ParaAlgoControle &
,CompThermoPhysiqueAbstraite::SaveResul * saveTP,CompThermoPhysiqueAbstraite* loiTP
,bool dilat,EnergieMeca & energ,const EnergieMeca & energ_t,bool premier_calcul
)
{
Temps_CPU_HZpp& temps_cpu_loi = ptintmeca.Tps_cpu_loi_comp();
temps_cpu_loi.Mise_en_route_du_comptage(); // spécifique pti
temps_loi.Mise_en_route_du_comptage(); // spécifique loi
// passage des infos specifiques aux classes derivantes
saveResul = saveDon;
def_en_cours=&def;
// récup des contraintes
TenseurHH & sigHH = *(ptintmeca.SigHH()); // ici _tdt est la grandeur finale
TenseurHH & sigHH_t = *(ptintmeca.SigHH_t()); // ici _tdt est la grandeur finale
TenseurBB & epsBB_tdt = *(ptintmeca.EpsBB());
TenseurBB & DepsBB_ = *(ptintmeca.DepsBB());
TenseurBB & delta_epsBB = *(ptintmeca.DeltaEpsBB());
double& module_compressibilite = ptintmeca.ModuleCompressibilite();
double& module_cisaillement = ptintmeca.ModuleCisaillement();
Tableau <double>& def_equi = ptintmeca.Deformation_equi();
const Tableau <double>& def_equi_t = ptintmeca.Deformation_equi_t_const();
dilatation=dilat;
ptintmeca_en_cours = &ptintmeca; // pour les méthodes internes
// calcul de : epsBB_tdt, delta_epsBB,gijBB_tdt, gijHH_tdt, d_gijBB_tdt,
// d_gijHH_tdt,jacobien,d_jacobien_tdt
Temps_CPU_HZpp& temps_cpu_metrique = ptintmeca.TpsMetrique();
temps_cpu_metrique.Mise_en_route_du_comptage(); // cpu
bool calcul_varD_ddl = utilise_vitesse_deformation;
const Met_abstraite::Impli& ex = def.Cal_implicit(def_equi_t,epsBB_tdt,d_epsBB_tdt,def_equi,DepsBB_
,delta_epsBB,premier_calcul);
//------- affichage d'erreurs éventuelles ----
bool sortie_metrique = false; // pour ne pas afficher 2 fois
if (Permet_affichage() > 2)
{ // on regarde la taille des déformations
double maxgijBB = ex.gijBB_tdt->MaxiComposante();
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (( maxgijBB > ConstMath::pasmalgrand)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
|| (!isfinite(maxgijBB)) || (isnan(maxgijBB))
)
sortie_metrique = true;
double maxepsBB_tdt = epsBB_tdt.MaxiComposante();
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (( maxepsBB_tdt > ConstMath::pasmalgrand)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
|| (!isfinite(maxepsBB_tdt)) || (isnan(maxepsBB_tdt))
)
sortie_metrique = true;
double maxgijHH = ex.gijHH_tdt->MaxiComposante();
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (( maxgijHH > ConstMath::pasmalgrand)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
|| (!isfinite(maxgijHH)) || (isnan(maxgijHH))
)
sortie_metrique = true;
if (!sortie_metrique) // on peut donc tester une division
{ double maxgijBB_0 = ex.gijBB_0->MaxiComposante();
double ratio1 = maxgijBB/maxgijBB_0;
if ( (ratio1 > 1000.) || (ratio1 < 0.001))
sortie_metrique = true;
double maxgijHH_0 = ex.gijHH_0->MaxiComposante();
double ratio2 = maxgijHH/maxgijHH_0;
if ( (ratio2 > 1000.) || (ratio2 < 0.001))
sortie_metrique = true;
};
if ((sortie_metrique) || (Permet_affichage() > 8))
def.Affiche();
};
// on récupère les positions des points
if (premier_calcul)
{ptintmeca.M_0() = def.Position_0();
ptintmeca.M_t() = def.Position_t();
};
ptintmeca.M_tdt() = def.Position_tdt();
temps_cpu_metrique.Arret_du_comptage(); // cpu
// dans le cas d'une loi thermo dépendante, calcul de la température
if ((thermo_dependant) || dilatation)
{temperature_0 = def.DonneeInterpoleeScalaire(TEMP,TEMPS_0);
temperature_t = def.DonneeInterpoleeScalaire(TEMP,TEMPS_t);
temperature_tdt = def.DonneeInterpoleeScalaire(TEMP,TEMPS_tdt);
temperature = &temperature_tdt;
};
// demande dans les classes dérivées du calcul de grandeurs spécifiques si besoin est
// dans le cas ou l'on tiens compte de la dilatation on modifie les différents tenseurs de déformation
ThermoDonnee dTP; // init à 0 par défaut
if (dilatation)
{ // attention ici on utilise la contrainte du temps précédent, il pourrait y avoir des oscillations
// vu que sigma suivant dépend de epsméca et que epsmeca dépend de sigma !!!
// double P = -(sigHH_t * (*(ex.gijBB_tdt))).Trace()/Abs(sigHH_t.Dimension()); // calcul de la pression
// on utilise la valeur à tdt, qui est donc celle calculée à l'itération précédente en implicite
double P = -(sigHH * (*(ex.gijBB_tdt))).Trace()/Abs(sigHH.Dimension()); // calcul de la pression
// et celle du pas précédent
double P_t = -(sigHH_t * (*(ex.gijBB_t))).Trace()/Abs(sigHH_t.Dimension()); // calcul de la pression à t
loiTP->Cal_donnees_thermiques(P_t,saveTP,def,P,TEMPS_tdt,dTP);
// modification de la déformation mécanique et calcul de la déformation thermique
// a- dimensionnement des tenseurs intermédiaires qui pour l'instant ne servent qu'ici
int dim_tens = epsBB_tdt.Dimension();
// -- cas de la déformation
if (epsBB_totale == NULL) { epsBB_totale = NevezTenseurBB(dim_tens);}
else if (epsBB_totale->Dimension() != dim_tens) { delete epsBB_totale;epsBB_totale = NevezTenseurBB(dim_tens);};
*epsBB_totale = epsBB_tdt;
if (epsBB_therm == NULL) { epsBB_therm = NevezTenseurBB(dim_tens);}
else if (epsBB_therm->Dimension() != dim_tens) { delete epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// -- cas de l'incrément de déformation
if (delta_epsBB_totale == NULL) { delta_epsBB_totale = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_totale->Dimension() != dim_tens) { delete delta_epsBB_totale;delta_epsBB_totale = NevezTenseurBB(dim_tens);};
*delta_epsBB_totale = delta_epsBB;
if (delta_epsBB_therm == NULL) { delta_epsBB_therm = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_therm->Dimension() != dim_tens) { delete delta_epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
// -- cas de la vitesse de déformation
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (DepsBB_totale == NULL) { DepsBB_totale = NevezTenseurBB(dim_tens);}
else if (DepsBB_totale->Dimension() != dim_tens) { delete DepsBB_totale;DepsBB_totale = NevezTenseurBB(dim_tens);};
*DepsBB_totale = DepsBB_;
if (DepsBB_therm == NULL) { DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (DepsBB_therm->Dimension() != dim_tens) { delete DepsBB_therm;DepsBB_totale = NevezTenseurBB(dim_tens);};
// -- calcul du partage entre thermique et mécanique
bool avec_repercution_sur_def_meca = true;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// des pointeurs pour des sorties en absolue
const Met_abstraite::Impli* ex_impli=&ex;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_umat = NULL;
def.DeformationThermoMecanique(temperature_0,dTP,*epsBB_totale
,*epsBB_therm,temperature_tdt,epsBB_tdt
,*delta_epsBB_totale,*delta_epsBB_therm,delta_epsBB
,temperature_t
,*DepsBB_totale,*DepsBB_therm,DepsBB_
,ex_impli,ex_expli_tdt,ex_umat
,true,avec_repercution_sur_def_meca);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
};
if ((thermo_dependant) || dilatation) // répercussion éventuelle du changement de température dans les classes dérivées
{ RepercuteChangeTemperature(TEMPS_tdt); };
// calcul éventuel des invariants liés à la cinématique
CalculInvariants_cinematique(ptintmeca,*(ex.gijBB_t),*(ex.gijHH_t));
// demande dans les classes dérivées du calcul de grandeurs spécifiques si besoin est
const Met_abstraite::Impli* ex_impli = &ex;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_expli = NULL;
CalculGrandeurTravail(ptintmeca,def,TEMPS_tdt,dTP,ex_impli,ex_expli_tdt,ex_expli,NULL,NULL);
// calcul de : sigHH, d_sigHH
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
try
{ Calcul_DsigmaHH_tdt (sigHH_t,DepsBB_,tab_ddl,*(ex.giB_t),*(ex.gijBB_t), *(ex.gijHH_t),
*(ex.giB_tdt),*(ex.d_giB_tdt),*(ex.giH_tdt),*(ex.d_giH_tdt),
epsBB_tdt,d_epsBB_tdt,delta_epsBB, *(ex.gijBB_tdt), *(ex.gijHH_tdt),
*(ex.d_gijBB_tdt), *(ex.d_gijHH_tdt),(*ex.jacobien_0), (*ex.jacobien_tdt),
*(ex.d_jacobien_tdt),sigHH, d_sigHH,energ,energ_t,module_compressibilite,module_cisaillement,ex);
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... ) // cas d'une erreur dans la recherche de la contrainte
{ throw ErrNonConvergence_loiDeComportement();
if (ParaGlob::NiveauImpression() >= 1)
{cout << "\n warning: exception generee par la loi de comportement ";
if (ParaGlob::NiveauImpression() >= 4) cout << "\n Loi_comp_abstraite::Cal_implicit(..";
};
};
jacobien = (*ex.jacobien_tdt);
d_jacobien_tdt = *(ex.d_jacobien_tdt);
// calcul éventuel des invariants de contraintes
CalculInvariants_contraintes(ptintmeca,*(ex.gijBB_tdt),*(ex.gijHH_tdt));
// en sortie il nous faut définir d_epsBB qui doit être égale à D*, donc pas forcément la dérivée
// de epsBB (par exemple dans le cas de la déformation logarithmique ou cumulée)
// variation de la déformation / au ddl
if (type_de_deformation != DEFORMATION_STANDART)
{ int d_epsBB_tdtTaille = d_epsBB_tdt.Taille();
for (int i=1; i<= d_epsBB_tdtTaille; i++)
*(d_epsBB_tdt(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
};
//------- affichage d'erreurs éventuelles ----
if (Permet_affichage() > 2)
{ // on regarde la taille des déformations
double maxsig = sigHH.MaxiComposante();
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (( maxsig > ConstMath::pasmalgrand)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
|| (!isfinite(maxsig)) || (isnan(maxsig))
)
if (!sortie_metrique)
{cout << "\n *** pb sur le calcul de la contrainte: le resultat est "
<< " soit infini soit un nan ";
Signature_pti_encours(cout);
mise à jour 7.010
2023-05-03 17:23:49 +02:00
cout << "\n Loi_comp_abstraite::Cal_implicit: sigHH_tdt ";
sigHH.Ecriture(cout);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
def.Affiche();
};
};
#ifdef MISE_AU_POINT
if (Permet_affichage() > 5)
{cout << "\n Loi_comp_abstraite::Cal_implicit: sigHH_tdt ";
sigHH.Ecriture(cout);
cout << "\n jacobien= "<< jacobien << ", d_jacobien_tdt:";
d_jacobien_tdt.Affiche();
if (Permet_affichage() > 4)
for (int i = 1; i<= d_epsBB_tdt.Taille(); i++)
{ cout << "\n d_sigHH("<<i<<") ";
d_sigHH(i)->Ecriture(cout);
cout << "\n d_epsBB_tdt("<<i<<") ";
d_epsBB_tdt(i)->Ecriture(cout);
};
};
#endif
def_en_cours=NULL; // plus d'accès possible à la sortie
ptintmeca_en_cours = NULL; // plus d'accès possible à la sortie
temps_cpu_loi.Arret_du_comptage(); // cpu spécifique pti
temps_loi.Arret_du_comptage(); // cpu spécifique loi
////debug
//if (temps_cpu_loi.Temps_CPU_User() == 0)
// { cout << "\n Loi_comp_abstraite::Cal_implicit: sigHH_tdt ";
// cout << " *** étrange, le temps loi de comp est nul ! " << endl;
// };
//// fin debug
// retour
return ex; // retour de la métrique
};
// schema pour le flambage linéaire
void Loi_comp_abstraite::Cal_flamb_lin
(SaveResul * saveDon,Deformation & def,DdlElement & tab_ddl
,PtIntegMecaInterne & ptintmeca, Tableau <TenseurBB *>& d_epsBB_tdt
,double& jacobien,Vecteur& d_jacobien_tdt,Tableau <TenseurHH *>& d_sigHH,const ParaAlgoControle &
,CompThermoPhysiqueAbstraite::SaveResul * saveTP,CompThermoPhysiqueAbstraite* loiTP
,bool dilat,EnergieMeca & energ,const EnergieMeca & energ_t,bool premier_calcul
)
{
Temps_CPU_HZpp& temps_cpu_loi = ptintmeca.Tps_cpu_loi_comp();
temps_cpu_loi.Mise_en_route_du_comptage(); // spécifique pti
temps_loi.Mise_en_route_du_comptage(); // spécifique loi
// passage des infos specifiques aux classes derivantes
saveResul = saveDon;
def_en_cours=&def;
// récup des contraintes
TenseurHH & sigHH = *(ptintmeca.SigHH()); // ici _tdt est la grandeur finale
TenseurHH & sigHH_t = *(ptintmeca.SigHH_t()); // ici _tdt est la grandeur finale
TenseurBB & epsBB_tdt = *(ptintmeca.EpsBB());
TenseurBB & DepsBB_ = *(ptintmeca.DepsBB());
TenseurBB & delta_epsBB = *(ptintmeca.DeltaEpsBB());
double& module_compressibilite = ptintmeca.ModuleCompressibilite();
double& module_cisaillement = ptintmeca.ModuleCisaillement();
Tableau <double>& def_equi = ptintmeca.Deformation_equi();
const Tableau <double>& def_equi_t = ptintmeca.Deformation_equi_t_const();
dilatation=dilat;
ptintmeca_en_cours = &ptintmeca; // pour les méthodes internes
// calcul des éléments de métrique
Temps_CPU_HZpp& temps_cpu_metrique = ptintmeca.TpsMetrique();
temps_cpu_metrique.Mise_en_route_du_comptage(); // cpu
const Met_abstraite::flambe_lin& ex = def.Cal_flambe_lin
(def_equi_t,epsBB_tdt,d_epsBB_tdt,def_equi,DepsBB_,delta_epsBB,premier_calcul);
// on récupère les positions des points
if (premier_calcul)
{ptintmeca.M_0() = def.Position_0();
ptintmeca.M_t() = def.Position_t();
};
ptintmeca.M_tdt() = def.Position_tdt();
temps_cpu_metrique.Arret_du_comptage(); // cpu
// dans le cas d'une loi thermo dépendante, calcul de la température
if ((thermo_dependant) || dilatation)
{temperature_0 = def.DonneeInterpoleeScalaire(TEMP,TEMPS_0);
temperature_t = def.DonneeInterpoleeScalaire(TEMP,TEMPS_t);
temperature_tdt = def.DonneeInterpoleeScalaire(TEMP,TEMPS_tdt);
temperature=&temperature_tdt;
};
// dans le cas ou l'on tiens compte de la dilatation on modifie les différents tenseurs de déformation
ThermoDonnee dTP; // init à 0 par défaut
if (dilatation)
{ // attention ici on utilise la contrainte du temps précédent, il pourrait y avoir des oscillations
// vu que sigma suivant dépend de epsméca et que epsmeca dépend de sigma !!!
// double P = -(sigHH_t * (*(ex.gijBB_tdt))).Trace()/Abs(sigHH_t.Dimension()); // calcul de la pression
// on utilise la valeur à tdt, qui est donc celle calculée à l'itération précédente en implicite
// en espérant que cela veuille dire quelque chose ici !!
double P = -(sigHH * (*(ex.gijBB_tdt))).Trace()/Abs(sigHH.Dimension()); // calcul de la pression
// et celle du pas précédent
double P_t = -(sigHH_t * (*(ex.gijBB_t))).Trace()/Abs(sigHH_t.Dimension()); // calcul de la pression à t
loiTP->Cal_donnees_thermiques(P_t,saveTP,def,P,TEMPS_tdt,dTP);
// modification de la déformation mécanique et calcul de la déformation thermique
// a- dimensionnement des tenseurs intermédiaires qui pour l'instant ne servent qu'ici
int dim_tens = epsBB_tdt.Dimension();
// -- cas de la déformation
if (epsBB_totale == NULL) { epsBB_totale = NevezTenseurBB(dim_tens);}
else if (epsBB_totale->Dimension() != dim_tens) { delete epsBB_totale;epsBB_totale = NevezTenseurBB(dim_tens);};
*epsBB_totale = epsBB_tdt;
if (epsBB_therm == NULL) { epsBB_therm = NevezTenseurBB(dim_tens);}
else if (epsBB_therm->Dimension() != dim_tens) { delete epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// -- cas de l'incrément de déformation
if (delta_epsBB_totale == NULL) { delta_epsBB_totale = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_totale->Dimension() != dim_tens) { delete delta_epsBB_totale;delta_epsBB_totale = NevezTenseurBB(dim_tens);};
*delta_epsBB_totale = delta_epsBB;
if (delta_epsBB_therm == NULL) { delta_epsBB_therm = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_therm->Dimension() != dim_tens) { delete delta_epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// -- cas de la vitesse de déformation
if (DepsBB_totale == NULL) { DepsBB_totale = NevezTenseurBB(dim_tens);}
else if (DepsBB_totale->Dimension() != dim_tens) { delete DepsBB_totale;DepsBB_totale = NevezTenseurBB(dim_tens);};
*DepsBB_totale = DepsBB_;
if (DepsBB_therm == NULL) { DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (DepsBB_therm->Dimension() != dim_tens) { delete DepsBB_therm;DepsBB_totale = NevezTenseurBB(dim_tens);};
// -- calcul du partage entre thermique et mécanique
bool avec_repercution_sur_def_meca = true;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// des pointeurs pour des sorties en absolue
const Met_abstraite::Impli* ex_impli=&ex;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_umat = NULL;
def.DeformationThermoMecanique(temperature_0,dTP,*epsBB_totale
,*epsBB_therm,temperature_tdt,epsBB_tdt
,*delta_epsBB_totale,*delta_epsBB_therm,delta_epsBB
,temperature_tdt
,*DepsBB_totale,*DepsBB_therm,DepsBB_
,ex_impli,ex_expli_tdt,ex_umat
,true,avec_repercution_sur_def_meca);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
};
if ((thermo_dependant) || dilatation) // répercussion éventuelle du changement de température dans les classes dérivées
{ RepercuteChangeTemperature(TEMPS_tdt); };
// calcul éventuel des invariants liés à la cinématique
CalculInvariants_cinematique(ptintmeca,*(ex.gijBB_t),*(ex.gijHH_t));
// demande dans les classes dérivées du calcul de grandeurs spécifiques si besoin est
const Met_abstraite::Impli* ex_impli = &ex;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_expli = NULL;
CalculGrandeurTravail(ptintmeca,def,TEMPS_tdt,dTP,ex_impli,ex_expli_tdt,ex_expli,NULL,NULL);
// calcul de : sigHH, d_sigHH
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
try
{ Calcul_DsigmaHH_tdt (sigHH_t,DepsBB_,tab_ddl,*(ex.giB_t),*(ex.gijBB_t), *(ex.gijHH_t),*(ex.giB_tdt)
,*(ex.d_giB_tdt),*(ex.giH_tdt),*(ex.d_giH_tdt)
,epsBB_tdt,d_epsBB_tdt,delta_epsBB, *(ex.gijBB_tdt), *(ex.gijHH_tdt)
,*(ex.d_gijBB_tdt), *(ex.d_gijHH_tdt),(*ex.jacobien_0), (*ex.jacobien_tdt)
,*(ex.d_jacobien_tdt),sigHH, d_sigHH,energ,energ_t,module_compressibilite,module_cisaillement,
ex);
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... ) // cas d'une erreur dans la recherche de la contrainte
{ throw ErrNonConvergence_loiDeComportement();
if (ParaGlob::NiveauImpression() >= 1)
cout << "\n warning: exception generee par la loi de comportement ";
if (ParaGlob::NiveauImpression() >= 4) cout << "\n Loi_comp_abstraite::Cal_flamb_lin(..";
};
jacobien = (*ex.jacobien_tdt);
d_jacobien_tdt = *(ex.d_jacobien_tdt);
// calcul éventuel des invariants de contraintes
CalculInvariants_contraintes(ptintmeca,*(ex.gijBB_tdt),*(ex.gijHH_tdt));
// en sortie il nous faut définir d_epsBB qui doit être égale à D*, donc pas forcément la dérivée
// de epsBB (par exemple dans le cas de la déformation logarithmique ou cumulée)
// variation de la déformation / au ddl
if (type_de_deformation != DEFORMATION_STANDART)
{ int d_epsBB_tdtTaille = d_epsBB_tdt.Taille();
for (int i=1; i<= d_epsBB_tdtTaille; i++)
*(d_epsBB_tdt(i)) = 0.5 * (*((*(ex.d_gijBB_tdt))(i)));
};
def_en_cours=NULL; // plus d'accès possible à la sortie
ptintmeca_en_cours = NULL; // plus d'accès possible à la sortie
temps_cpu_loi.Arret_du_comptage(); // cpu spécifique pti
temps_loi.Arret_du_comptage(); // cpu spécifique loi
// retour
};
// schema pour le calcul de la loi de comportement dans le cas de l'umat
void Loi_comp_abstraite::ComportementUmat
(Loi_comp_abstraite::SaveResul * saveDon
,Deformation & def,PtIntegMecaInterne& ptintmeca,ParaAlgoControle & para
,CompThermoPhysiqueAbstraite::SaveResul * saveTP,CompThermoPhysiqueAbstraite* loiTP
,bool dilat,UmatAbaqus& umatAbaqusqus,bool premier_calcul)
{
Temps_CPU_HZpp& temps_cpu_loi = ptintmeca.Tps_cpu_loi_comp();
temps_cpu_loi.Mise_en_route_du_comptage(); // spécifique pti
temps_loi.Mise_en_route_du_comptage(); // spécifique loi
// passage des infos specifiques aux classes derivantes
saveResul = saveDon;
dilatation=dilat;
def_en_cours=&def;
ptintmeca_en_cours = &ptintmeca; // pour les méthodes internes
// on calcul les éléments de la métrique dans le cas particulier de l'umat d'où d'une métrique
// Met_ElemPoint et d'une déformation particulière Def_Umat
Def_Umat& defUmat = *((Def_Umat*) &def); // récup
// on regarde si c'est le premier incrément ou pas (d'où le calcul du premier pas)
// bool premier_calcul=false;
// if (umatAbaqusqus.nb_increment <= 1) premier_calcul = true;
Temps_CPU_HZpp& temps_cpu_metrique = ptintmeca.TpsMetrique();
temps_cpu_metrique.Mise_en_route_du_comptage(); // cpu
const Met_abstraite::Umat_cont& ex= defUmat.Cal_defUmat(premier_calcul);
// on récupère les positions des points
if (premier_calcul)
{ptintmeca.M_0() = def.Position_0();
ptintmeca.M_t() = def.Position_t();
};
ptintmeca.M_tdt() = def.Position_tdt();
umatAbaqusqus.coor_pt = ptintmeca.M_tdt();// on met à jour les coordonnées absolues du pti temps_cpu_metrique.Arret_du_comptage(); // cpu
// -- modification du pas de temps: on part du principe que le pas de temps est géré par le programme appelant donc ici on ne fait que répercuter
// -- le pas de temps et le temps courant au niveau des paramètres globaux du calcul, qui normalement ne sont (et ne doivent être)
// -- utilisées que pour la loi de comportement Umat. Les autres cas ne sont pas prévus
para.Modif_temps(umatAbaqusqus.temps_tdt,umatAbaqusqus.delta_t);
// // modif le 21 avril 2016 -> on passe les infos via des paramètres de passage
// // double deltat = umatAbaqusqus.delta_t;
// // VariablesTemps tempo_umat(para.Variables_de_temps());
// // ici on accède directement car c'est loi_comp_abstraite est friend de Variables Temps
// tempo_umat.temps = umatAbaqusqus.temps_tdt;
// tempo_umat.deltat = umatAbaqusqus.delta_t;
// dans le cas d'une loi thermo dépendante, calcul de la température
if ((thermo_dependant) || dilatation)
{temperature_0 = 0.; // a priori pour l'instant
temperature_t = umatAbaqusqus.temper_t;
temperature_tdt = umatAbaqusqus.temper_t+umatAbaqusqus.delta_temper;
temperature=&temperature_tdt;
};
#ifdef MISE_AU_POINT
if ( Permet_affichage() > 5)
{cout << "\n temperature_0= "<< temperature_0
<< " temperature_t= "<< temperature_t
<< " temperature= "<< temperature
<< "\n thermo_dependant= "<<thermo_dependant
<< " dilatation= "<< dilatation
<< flush;
};
#endif
// -- cas de la vitesse de déformation à la umat
// construction du tenseur vitesse de déformation, constant sur le pas de temps
// = delta epsilonBB /delta t
if (DepsBB_umat == NULL) // cas du premier passage, on cré le tenseur
{ DepsBB_umat = NevezTenseurBB(*(umatAbaqusqus.delta_eps_meca));}
else if (DepsBB_umat->Dimension() != (umatAbaqusqus.delta_eps_meca)->Dimension()) // tenseur inadapté, on change
{ delete DepsBB_umat;DepsBB_umat = NevezTenseurBB(*(umatAbaqusqus.delta_eps_meca));}
else // cas courant on récupère la valeur de delta_epsilon_BB
{ *DepsBB_umat = *(umatAbaqusqus.delta_eps_meca); };
// puis on divise par l'incrément de temsp
if (Abs(umatAbaqusqus.delta_t) >= ConstMath::trespetit)
{(*DepsBB_umat) /= umatAbaqusqus.delta_t;}
else {DepsBB_umat->Inita(0.);}; // si trop petit on met à 0
ThermoDonnee dTP; // init à 0 par défaut
if (dilatation)
{ // attention ici on utilise la contrainte du temps précédent, car a priori on ne dispose pas de la nouvelle valeur
// il pourrait y avoir des oscillations
// dans le cas de l'umat, le tenseur initial sigma est en orthonormée
double P=0.;
for (int i=1;i<=umatAbaqusqus.dim_tens;i++) P += (*umatAbaqusqus.t_sigma)(i,i);
P /= umatAbaqusqus.dim_tens; // pression
loiTP->Cal_donnees_thermiques(P,saveTP,def,P,TEMPS_tdt,dTP);
// Dans le cas de la procédure Umat, la déformation récupérée est directement la déformation mécanique
// l'influence de la dilatation thermique s'effectue éventullement au sein de la loi de comportement
// comme par exemple Hysteresis_bulk, quand on utilise pas directement les deformations méca
// a- dimensionnement des tenseurs intermédiaires qui pour l'instant ne servent qu'ici
// epsBB_totale et DepsBB_totale ne servent à rien en fait ... l'objectif est uniquement de satifaire
// les variables de passages dans DeformationThermoMecanique, car comme avec_repercution_sur_def_meca = false
// on ne change pas les tenseurs de déformation et de vitesse de déformation
// par contre : epsBB_therm et DepsBB_therm sont calculées dans DeformationThermoMecanique
int dim_tens = 3; // a priori on est en dimension 3
// -- cas de la déformation
if (epsBB_totale == NULL) { epsBB_totale = NevezTenseurBB(dim_tens);}
else if (epsBB_totale->Dimension() != dim_tens) { delete epsBB_totale;epsBB_totale = NevezTenseurBB(dim_tens);};
if (epsBB_therm == NULL) { epsBB_therm = NevezTenseurBB(dim_tens);}
else if (epsBB_therm->Dimension() != dim_tens) { delete epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// -- cas de l'incrément de déformation
if (delta_epsBB_totale == NULL) { delta_epsBB_totale = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_totale->Dimension() != dim_tens) { delete delta_epsBB_totale;delta_epsBB_totale = NevezTenseurBB(dim_tens);};
if (delta_epsBB_therm == NULL) { delta_epsBB_therm = NevezTenseurBB(dim_tens);}
else if (delta_epsBB_therm->Dimension() != dim_tens) { delete delta_epsBB_therm;epsBB_therm = NevezTenseurBB(dim_tens);};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// -- cas de la vitesse de déformation
if (DepsBB_totale == NULL) { DepsBB_totale = NevezTenseurBB(dim_tens);}
else if (DepsBB_totale->Dimension() != dim_tens) { delete DepsBB_totale;DepsBB_totale = NevezTenseurBB(dim_tens);};
if (DepsBB_therm == NULL) { DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (DepsBB_therm->Dimension() != dim_tens) { delete DepsBB_therm;DepsBB_totale = NevezTenseurBB(dim_tens);};
// -- calcul du partage entre thermique et mécanique
bool avec_repercution_sur_def_meca = false;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// des pointeurs pour des sorties en absolue
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_umat =&ex;
def.DeformationThermoMecanique(temperature_0,dTP,*epsBB_totale
,*epsBB_therm,temperature_tdt,*(umatAbaqusqus.eps_meca)
,*delta_epsBB_totale,*delta_epsBB_therm,*(umatAbaqusqus.delta_eps_meca)
,temperature_t
,*DepsBB_totale,*DepsBB_therm,*DepsBB_umat
,ex_impli,ex_expli_tdt,ex_umat
,true,avec_repercution_sur_def_meca);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
};
// cas des énergie méca
// tout d'abord celles précédentes
EnergieMeca ener_t(umatAbaqusqus.energie_elastique,umatAbaqusqus.dissipation_plastique
,umatAbaqusqus.dissipation_visqueuse);
EnergieMeca ener; // la sortie initialisée à 0
// transfert de l'Umat vers les pt d'integ avant le calcul de la contrainte donc à t
umatAbaqusqus.Transfert_Umat_ptInteg_t(ptintmeca);
*(ptintmeca.DepsBB())=*DepsBB_umat; // stockage de vitesse de déformation
// calcul des def équivalentes
// - récup des infos
TenseurBB & epsBB_tdt = *(ptintmeca.EpsBB());
TenseurBB & DepsBB_ = *(ptintmeca.DepsBB());
TenseurBB & delta_epsBB = *(ptintmeca.DeltaEpsBB());
Tableau <double>& def_equi = ptintmeca.Deformation_equi();
const Tableau <double>& def_equi_t = ptintmeca.Deformation_equi_t_const();
// - calcul
defUmat.CalDefEqui(def_equi_t,epsBB_tdt,def_equi,DepsBB_,delta_epsBB,ex,premier_calcul);
if ((thermo_dependant) || dilatation) // répercussion éventuelle du changement de température dans les classes dérivées
{ RepercuteChangeTemperature(TEMPS_tdt); };
// calcul éventuel des invariants liés à la cinématique
CalculInvariants_cinematique(ptintmeca,*(ex.gijBB_t),*(ex.gijHH_t));
// calcul éventuel préalable avant le calcul de la contrainte
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_expli = &ex;
CalculGrandeurTravail(ptintmeca,defUmat,TEMPS_tdt,dTP,ex_impli,ex_expli_tdt,ex_expli,NULL,NULL);
// affichage de certaines données de def d'entrée du calcul Calcul_dsigma_deps
// car elles sont différentes (car calculées) en fonction des entrées
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 4)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\n affichage des def avant appel de Calcul_dsigma_deps \n ";
int dim = ParaGlob::Dimension();
// def en orthonormee
TenseurBB* epsBB_inter= (NevezTenseurBB(dim)) ;
(umatAbaqusqus.eps_meca)->BaseAbsolue(*epsBB_inter,*(ex.giH_tdt));
cout << "\n epsBB meca d'entree =en orthonormee ";epsBB_inter->Ecriture(cout);
delete epsBB_inter;
// Deps
TenseurBB* DepsBB_inter= (NevezTenseurBB(dim)) ;
DepsBB_umat->BaseAbsolue(*DepsBB_inter,*(ex.giH_tdt));
cout << "\n DepsBB meca d'entree =en orthonormee "; DepsBB_inter->Ecriture(cout);
delete DepsBB_inter;
// Delta eps
TenseurBB* delta_epsBB_inter= (NevezTenseurBB(dim)) ;
(umatAbaqusqus.delta_eps_meca)->BaseAbsolue(*delta_epsBB_inter,*(ex.giH_tdt));
cout << "\n delta eps meca d'entree =en orthonormee "; delta_epsBB_inter->Ecriture(cout);
delete delta_epsBB_inter;
cout << endl;
};
#endif
// calcul de : sigHH, d_sigHH
bool en_base_orthonormee=true; // ici les tenseurs sont en orthonormee a priori
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
try
{ Calcul_dsigma_deps(en_base_orthonormee, *(umatAbaqusqus.t_sigma),*DepsBB_umat
,*(umatAbaqusqus.eps_meca),*(umatAbaqusqus.delta_eps_meca)
,(*ex.jacobien_0),(*ex.jacobien_tdt)
,*(umatAbaqusqus.t_sigma),*(umatAbaqusqus.d_sigma_deps)
,ener,ener_t,ptintmeca.ModuleCompressibilite(),ptintmeca.ModuleCisaillement(),ex
) ;
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... ) // cas d'une erreur dans la recherche de la contrainte
{ throw ErrNonConvergence_loiDeComportement();
if (ParaGlob::NiveauImpression() >= 1)
cout << "\n warning: exception generee par la loi de comportement ";
if (ParaGlob::NiveauImpression() >= 4) cout << "\n Loi_comp_abstraite::ComportementUmat(..";
};
umatAbaqusqus.energie_elastique = ener.EnergieElastique();
umatAbaqusqus.dissipation_plastique = ener.DissipationPlastique();
umatAbaqusqus.dissipation_visqueuse= ener.DissipationVisqueuse();
// dans le cas où il s'agit d'une loi de contrainte plane, on met à jour les normales finales
// ce qui permettra de transmettre la variation d'épaisseur
if ((Comp_3D_CP_DP_1D(this->Id_comport()) == COMP_CONTRAINTES_PLANES)
&& (ParaGlob::Dimension() == 3))
{ umatAbaqusqus.N_tdt(3) = HsurH0(saveDon);
};
// transfert de l'Umat vers les pt d'integ
// umatAbaqusqus.Transfert_Umat_ptInteg_t(ptintmeca); // correction erreur le 7/sept/ 2013
umatAbaqusqus.Transfert_Umat_ptInteg_tdt(ptintmeca); // correction erreur le 7/sept/ 2013
// calcul éventuel des invariants de contraintes
CalculInvariants_contraintes(ptintmeca,*(ex.gijBB_tdt),*(ex.gijHH_tdt));
def_en_cours=NULL; // plus d'accès possible à la sortie
ptintmeca_en_cours = NULL; // plus d'accès possible à la sortie
temps_cpu_loi.Arret_du_comptage(); // cpu spécifique pti
temps_loi.Arret_du_comptage(); // cpu spécifique loi
};
// calcul d'un module d'young équivalent à la loi, ceci pour un
// chargement nul
double Loi_comp_abstraite::Module_young_equivalent(Enum_dure ,const Deformation & ,SaveResul * )
{ cout << "\n erreur , aucun module d'young equivalent n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// récupération d'un module de compressibilité équivalent à la loi pour un chargement nul
// il s'agit ici de la relation -pression = sigma_trace/3. = module de compressibilité * I_eps
double Loi_comp_abstraite::Module_compressibilite_equivalent(Enum_dure ,const Deformation & ,SaveResul * )
{ cout << "\n erreur , aucun module de compressibilite equivalent n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// récupération de la dernière déformation d'épaisseur calculée: cette déformaion n'est utile que pour des lois en contraintes planes ou doublement planes
// - pour les lois 3D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// - pour les lois 2D def planes: retour de 0
// les infos nécessaires à la récupération de la def, sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double Loi_comp_abstraite::Eps33BH(SaveResul * saveResul) const
{ cout << "\n erreur , la methode Eps33BH() n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// récupération de la dernière déformation de largeur calculée: cette déformaion n'est utile que pour des lois en contraintes doublement planes
// - pour les lois 3D et 2D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// les infos nécessaires à la récupération de la def, sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double Loi_comp_abstraite::Eps22BH(SaveResul * saveResul) const
{ cout << "\n erreur , la methode Eps22BH() n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// récupération de la variation relative d'épaisseur calculée: h/h0
// et de sa variation par rapport aux ddls la concernant: d_hsurh0
// cette variation n'est utile que pour des lois en contraintes planes
// - pour les lois 3D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// - pour les lois 2D def planes: retour de 0
// les infos nécessaires à la récupération , sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double Loi_comp_abstraite::d_HsurH0(SaveResul * saveResul,Vecteur & d_hsurh0) const
{ cout << "\n erreur , la methode d_HsurH0() n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// récupération de la variation relative de largeur calculée: b/b0
// cette variation n'est utile que pour des lois en contraintes planes double
// - pour les lois 3D et 2D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// les infos nécessaires à la récupération , sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double Loi_comp_abstraite::BsurB0(SaveResul * saveResul) const
{ cout << "\n erreur , la methode BsurB0() n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// récupération de la variation relative de largeur calculée: b/b0
// et de sa variation par rapport aux ddls la concernant: d_bsurb0
// cette variation n'est utile que pour des lois en contraintes planes
// - pour les lois 3D et 2D : retour d'un nombre très grand, indiquant que cette fonction est invalide
// les infos nécessaires à la récupération , sont stockées dans saveResul
// qui est le conteneur spécifique au point où a été calculé la loi
double Loi_comp_abstraite::d_BsurB0(SaveResul * saveResul,Vecteur & d_bsurb0) const
{ cout << "\n erreur , la methode d_BsurB0() n'est defini pour cette loi ";
this->Affiche();
Sortie(1);
return 0.; // pour éviter le warning à la compil
};
// mise à jour des températures d'une manière unilatérale: sert par exemple pour
// le calcul du module d'young équivalent, quand il est calculé en dehors du calcul des contraintes
void Loi_comp_abstraite::Mise_a_jour_temperature(Enum_dure temps,Deformation & def)
{
if (thermo_dependant)
{switch (temps)
{ case TEMPS_0:
temperature_0 = def.DonneeInterpoleeScalaire(TEMP,TEMPS_0);
temperature = &temperature_0;
break;
case TEMPS_t:
temperature_t = def.DonneeInterpoleeScalaire(TEMP,TEMPS_t);
temperature = &temperature_t;
break;
case TEMPS_tdt:
temperature_tdt = def.DonneeInterpoleeScalaire(TEMP,TEMPS_tdt);
temperature = &temperature_tdt;
break;
default :
cout << "\n cas du temps non implante temps= " << Nom_dure(temps)
<< "\n Loi_comp_abstraite::Mise_a_jour_temperature(Enum_dure temps...";
Sortie(1);
};
// puis on répercute la température
RepercuteChangeTemperature(temps);
};
};
// lecture éventuelle du type de déformation et du niveau de commentaire
void Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(UtilLecture& entreePrinc,LesFonctions_nD& lesFonctionsnD
,bool avec_passage_nouvelle_donnee)
{
string nom_class_methode("Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire(..");
// recherche du mot clé éventuel indiquant le type de déformation
if (avec_passage_nouvelle_donnee)
entreePrinc.NouvelleDonnee();
// on note ce qu'il faut lire
int a_lire_type_deformation = 0;
if (strstr(entreePrinc.tablcar,"type_de_deformation")!= NULL)
a_lire_type_deformation = 1;
int a_lire_permet_affichage = 0;
if (strstr(entreePrinc.tablcar,"permet_affichage_")!= NULL)
a_lire_permet_affichage = 1;
// on initialise par défaut: la déformation standart
type_de_deformation = DEFORMATION_STANDART;
// le niveau de commentaire local
// on ne l'initialise pas car il peut déjà avoir une valeur
// via une lecture de paramètres particuliers permet_affichage = 0;
if ( a_lire_type_deformation || a_lire_permet_affichage)
{while ( a_lire_type_deformation || a_lire_permet_affichage)
{ // on lit le mot clé
string nom;
*(entreePrinc.entree) >> nom;
if (nom == "type_de_deformation")
{ *(entreePrinc.entree) >> nom ;
type_de_deformation=Id_nom_type_deformation(nom.c_str());
a_lire_type_deformation = 0;
entreePrinc.NouvelleDonnee();
if (strstr(entreePrinc.tablcar,"permet_affichage_")!= NULL)
a_lire_permet_affichage = 1;
}
else if (nom == "permet_affichage_")
{Loi_comp_abstraite::Lecture_permet_affichage(&entreePrinc,lesFonctionsnD);
a_lire_permet_affichage = 0;
entreePrinc.NouvelleDonnee();
if (strstr(entreePrinc.tablcar,"type_de_deformation")!= NULL)
a_lire_type_deformation = 1;
}
else
{cout << "\n erreur*** : on a lue " << nom
<< " et on attendait un des deux mots clefs suivant: "
<< " type_de_deformation ou permet_affichage_ "
<< "\n Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire(.. "
<< endl;
Sortie(1);
};
};
};
};
// affichage des données de la classe comp_abstraite
void Loi_comp_abstraite::Affiche_don_classe_abstraite() const
{ cout << "\n type_de_deformation: " << Nom_type_deformation(type_de_deformation);
{cout << "\n list_grandeur_particuliere_dispo_localement:taille= "
<< listdeTouslesQuelc_dispo_localement.size()<<" ";
list <EnumTypeQuelconque >::const_iterator il,ilfin=listdeTouslesQuelc_dispo_localement.end();
for (il=listdeTouslesQuelc_dispo_localement.begin();il != ilfin;il++)
cout << NomTypeQuelconque(*il) << " ";
};
{cout << "\n list_grandeur_particuliere_demandee_localement:taille= "
<< listQuelc_mis_en_acces_localement.size()<<" ";
list <EnumTypeQuelconque >::const_iterator il,ilfin=listQuelc_mis_en_acces_localement.end();
for (il=listQuelc_mis_en_acces_localement.begin();il != ilfin;il++)
cout << NomTypeQuelconque(*il) << " ";
};
cout << endl;
};
// affichage et definition interactive des commandes particulières à la classe loi_comp_abstraite
void Loi_comp_abstraite::Info_commande_don_LoisDeComp(UtilLecture& entreePrinc) const
{ ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
cout << "\n definition standart de la deformation (rep o (defaut) ou n) ? ";
string rep = "_";
// procédure de lecture avec prise en charge d'un retour chariot
rep = lect_return_defaut(true,"o");
if ((Minuscules(rep) != "f")&&(Minuscules(rep) != "0")&&(Minuscules(rep) != "o"))
{rep = "_";
while ((Minuscules(rep) != "f")&&(Minuscules(rep) != "0")&&(Minuscules(rep) != "o")
&&(Minuscules(rep) != "1")&&(Minuscules(rep) != "2"))
{try
{
cout << "\n -- definition du type de mesure de deformation : --- "
<< "\n (0 ou f) (fin) "
<< "\n (1) Almansi (defaut) "
<< "\n (2) Logarithmique (experimental) "
<< "\n ";
// procédure de lecture avec prise en charge d'un retour chariot
rep = lect_return_defaut(true,"1");
if ((Minuscules(rep) != "f")&&(Minuscules(rep) != "0")&&(Minuscules(rep) != "o")
&&(Minuscules(rep) != "1"))
{ if (Minuscules(rep) == "2")
{cout << "\n choix de la deformation Logarithmique ";
sort << "\n# -- definition du type de deformation (par defaut: DEFORMATION_STANDART) -- "
<< "\n type_de_deformation DEFORMATION_LOGARITHMIQUE ";
}
else { cout << "\n Erreur on attendait un entier entre 0 et 2 !!, "
<< "\n redonnez une bonne valeur"
<< "\n ou taper f ou 0 pour stoper la definition de la deformation ";
};
}
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch (...)//(UtilLecture::ErrNouvelleDonnee erreur)
{ cout << "\n Erreur on attendait un des mots cles proposes !!, "
<< "\n redonnez une bonne valeur"
<< "\n ou taper f ou 0 pour sortir ";
};
}; //-- fin du while
}; // fin du cas ou on veut une définition non standart
// sort << "\n# -- definition du type de deformation (par defaut: DEFORMATION_STANDART) -- "
// << "\n type_de_deformation DEFORMATION_STANDART ";
};
//----- lecture écriture de restart spécifique aux données de la classe -----
// 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 Loi_comp_abstraite::Lecture_don_base_info(ifstream& ent,const int cas,LesReferences& lesRef
,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string toto;
if (cas == 1)
{ // tout d'abord appel de la classe mère
LoiAbstraiteGeneral::Lect_base_info_loi(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
// puis les infos propres
// le type de déformation
ent >> toto ;
if (toto != "type_deformation")
{ cout << "\n erreur en lecture du type de deformation, nom lu: " << toto
<< "\n Loi_comp_abstraite::Lecture_don_base_info(...";
Sortie(1);
};
// lecture des infos
ent >> toto; type_de_deformation=Id_nom_type_deformation(toto.c_str());
// lecture concernant les grandeurs particulières de liaison
{int taille=0;string nom;
ent >> toto >> taille;
listdeTouslesQuelc_dispo_localement.clear();
for (int i=1;i<=taille;i++)
{ ent >> nom;listdeTouslesQuelc_dispo_localement.push_back(Id_nomTypeQuelconque(nom));};
};
{int taille=0;string nom;
listQuelc_mis_en_acces_localement.clear();
ent >> toto >> taille;
for (int i=1;i<=taille;i++)
{ ent >> nom;listQuelc_mis_en_acces_localement.push_back(Id_nomTypeQuelconque(nom));};
};
};
ent >> toto >> temps_loi;
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Loi_comp_abstraite::Ecriture_don_base_info(ofstream& sort,const int cas) const
{ if (cas == 1)
{ // tout d'abord appel de la classe mère
LoiAbstraiteGeneral::Ecrit_base_info_loi(sort,cas);
// puis les infos propres
// le type de déformation
sort << " \n type_deformation " << Nom_type_deformation(type_de_deformation) << " ";
{sort << "\n list_grandeur_particuliere_dispo_localement:taille= "
<< listdeTouslesQuelc_dispo_localement.size() << " ";
list <EnumTypeQuelconque >::const_iterator il,ilfin=listdeTouslesQuelc_dispo_localement.end();
for (il=listdeTouslesQuelc_dispo_localement.begin();il != ilfin;il++)
sort << NomTypeQuelconque(*il) << " ";
};
{sort << "\n list_grandeur_particuliere_demandee_localement:taille= "
<< listQuelc_mis_en_acces_localement.size() << " ";
list <EnumTypeQuelconque >::const_iterator il,ilfin=listQuelc_mis_en_acces_localement.end();
for (il=listQuelc_mis_en_acces_localement.begin();il != ilfin;il++)
sort << NomTypeQuelconque(*il) << " ";
};
};
sort << "\n tps_cumule_loi: "<<temps_loi;
};
// activation des données des noeuds et/ou elements nécessaires au fonctionnement de la loi
// exemple: mise en service des ddl de température aux noeuds
// méthode appelée par Activation_donnees principal, ou des classes dérivées
void Loi_comp_abstraite::Activ_donnees(Tableau<Noeud *>& tabnoeud,bool dilatation,LesPtIntegMecaInterne& lesPtMecaInt)
{ if ((thermo_dependant)|| (dilatation))
{ // dans le cas où la loi est thermo dépendante on active les ddl de thermique
int nbnoeud = tabnoeud.Taille();
for (int i=1;i<=nbnoeud;i++)
{ // on vérifie que la variable TEMP existe sinon erreur
if (tabnoeud(i)->Existe_ici(TEMP))
{tabnoeud(i)->Met_en_service(TEMP);}
else
{ cout << "\n erreur: la variable temperature n'existe pas "
<< " il n'est pas possible d'utiliser une loi thermodependante "
<< " il manque sans doute des donnees !!! "
<< "\n Loi_comp_abstraite::Activ_donnees(...";
Sortie(1);
}
}
};
// cas de l'utilisation d'une fonction nD pour l'affichage
if (permet_affich_loi_nD != NULL)
{// l'objectif ici est de définir les conteneurs de grandeurs totalement quelconques
// qui seront ensuite nécessaire pour appeler la fonction
// on ne fait qu'une fois l'opération, car la fonction est appelée à chaque pti
// 1) on récupère le tableau des énumérés de quelconques de la fonction
const Tableau <EnumTypeQuelconque>& tqi = permet_affich_loi_nD->Tab_enu_quelconque();
if (tqi.Taille() != li_quelconque.size())
{// on initialise la liste
li_quelconque.clear();
// on va boucler sur les énumérés quelconque et
// créer les conteneurs de retour
int nb_tqi = tqi.Taille();
for (int i_tqi = 1;i_tqi <= nb_tqi;i_tqi++ )
// on crée les conteneurs
{EnumTypeQuelconque& enuconq = tqi(i_tqi);
switch (enuconq)
{
case CONTRAINTE_MISES :
{Grandeur_scalaire_double grand_courant(0.);
TypeQuelconque typQ1(CONTRAINTE_MISES,SIG11,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case CONTRAINTE_TRESCA :
{Grandeur_scalaire_double grand_courant(0.);
TypeQuelconque typQ1(CONTRAINTE_TRESCA,SIG11,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case CONTRAINTE_MISES_T :
{Grandeur_scalaire_double grand_courant(0.);
TypeQuelconque typQ1(CONTRAINTE_MISES_T,SIG11,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case CONTRAINTE_TRESCA_T :
{Grandeur_scalaire_double grand_courant(0.);
TypeQuelconque typQ1(CONTRAINTE_TRESCA_T,SIG11,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case NUM_ELEMENT:
{Grandeur_scalaire_entier grand_courant(0);
TypeQuelconque typQ1(NUM_ELEMENT,NU_DDL,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case NUM_MAIL_ELEM:
{Grandeur_scalaire_entier grand_courant(0);
TypeQuelconque typQ1(NUM_MAIL_ELEM,NU_DDL,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case NUM_PTI:
{Grandeur_scalaire_entier grand_courant(0);
TypeQuelconque typQ1(NUM_PTI,NU_DDL,grand_courant);
li_quelconque.push_back(typQ1);
break;
}
case REPERE_LOCAL_ORTHO :
{Coordonnee v_rep(ParaGlob::Dimension()); // un vecteur intermédiaire
Tab_Grandeur_Coordonnee gr(v_rep,3); // def d'une grandeur intermédiaire
// stockage
TypeQuelconque typQ6(REPERE_LOCAL_ORTHO,SIG11,gr);
li_quelconque.push_back(typQ6);
break;
}
case REPERE_LOCAL_H :
{Coordonnee v_rep(ParaGlob::Dimension()); // un vecteur intermédiaire
Tab_Grandeur_Coordonnee gr(v_rep,3); // def d'une grandeur intermédiaire
// stockage
TypeQuelconque typQ6(REPERE_LOCAL_H,SIG11,gr);
li_quelconque.push_back(typQ6);
break;
}
case REPERE_LOCAL_B :
{Coordonnee v_rep(ParaGlob::Dimension()); // un vecteur intermédiaire
Tab_Grandeur_Coordonnee gr(v_rep,3); // def d'une grandeur intermédiaire
// stockage
TypeQuelconque typQ6(REPERE_LOCAL_B,SIG11,gr);
li_quelconque.push_back(typQ6);
break;
}
default :
{ cout << "\n erreur dans la definition du stockage inter pour la fonction nD "
<< " d'affichage : "<< permet_affich_loi_nD->NomFonction()
<< " cas non traite: "
<< NomTypeQuelconque(enuconq) << "!\n";
cout << "\n Loi_comp_abstraite::Activ_donnees(....";
Sortie(1);
};
};
};
// on s'occupe maintenant des pointeurs
// qui seront nécessaire pour utiliser la fonction nD
if (li_quelconque.size())
{tab_pt_li_quelconque.Change_taille(li_quelconque.size());
List_io <TypeQuelconque >::iterator ili,ilifin = li_quelconque.end();
int i=1; // init
for (ili= li_quelconque.begin();ili != ilifin; ili++,i++)
tab_pt_li_quelconque(i) = &(*ili);
};
};
};
};
// ***** a supprimer à terme, -> écriture d'un message d'erreur ******
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormeee
// le tenseur de déformation et son incrémentsont également en orthonormees
// si = false: les bases transmises sont utilisées
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void Loi_comp_abstraite::Calcul_dsigma_deps
(bool en_base_orthonormee, TenseurHH & ,TenseurBB&
,TenseurBB & ,TenseurBB & ,double& ,double&
,TenseurHH& ,TenseurHHHH&
,EnergieMeca & ,const EnergieMeca & ,double& , double& ,const Met_abstraite::Umat_cont&
)
{ cout << "\n ERREUR: methode non encore implante !!!!"
<< "\n methode Loi_comp_abstraite::Calcul_dsigma_deps (...";
Affiche();
Sortie(1);
};
// calcul éventuel des invariants liés aux contraintes
void Loi_comp_abstraite::CalculInvariants_contraintes
(PtIntegMecaInterne& ptIntegMeca, TenseurBB & gijBB,TenseurHH & gijHH)
{ // def de la dimension des tenseurs
int dimT = ptIntegMeca.EpsBB()->Dimension();
int caas=0; // pour la gestion d'erreur éventuelle
// on calcul maintenant systématiquement les contraintes de Mises et Tresca
TenseurHB& sigHB = *(NevezTenseurHB(dimT)) ;
sigHB = (*(ptIntegMeca.SigHH()) * gijBB);
// calcul de Mises et Tresca
Coordonnee valPropreSig(sigHB.ValPropre(caas));
if (caas == -1)
{cout << "\n *** erreur dans le calcul des valeurs propres de sigma "
<< " pb eventuel si on doit ensuite les utiliser (dans une fonction nD par exemple) "
<< " on continue neanmoins ...";
cout << "\n Loi_comp_abstraite::CalculInvariants_contraintes(..." << flush;
};
Tableau <double>& tab_sig_equi = ptIntegMeca.Sig_equi();
Coordonnee& vv = valPropreSig; int dimvec=vv.Dimension();// pour condenser l'écriture
double& sig_mises = tab_sig_equi(1); // pour être plus explicite
double& Tresca = tab_sig_equi(2); // "
switch (dimvec) // dans le cas où dimvec=0 on ne fait rien, cas ou on n'a pas besoin de mises
{ case 1: sig_mises = Dabs(vv(1)); Tresca=0.5 * vv(1);
break;
case 2: sig_mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + vv(1) * vv(1)
+ vv(2) * vv(2)) * 0.5);
Tresca=0.5 * (vv(1)-vv(2));
break;
case 3: sig_mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + (vv(1)-vv(3))*(vv(1)-vv(3))
+ (vv(3)-vv(2))*(vv(3)-vv(2))) * 0.5);
Tresca=0.5 * (vv(1)-vv(3));
break;
default: cout << "\n erreur de dimension, ne doit pas arriver ?? "
<< "\n Loi_comp_abstraite::CalculInvariants_contraintes(..." << flush;
Sortie(1);
};
// calcul éventuel des invariants
if (ptIntegMeca.Statut_Invariants_contrainte())
{ Vecteur & invariant = ptIntegMeca.SigInvar();
invariant.Change_taille(Abs(dimT));
switch (abs(dimT))
{case 3: invariant(3)= sigHB.Det();
case 2: invariant(2)= sigHB.II();
case 1: invariant(1)= sigHB.Trace();
};
////----- debug
////if ((invariant(1) > 0.) && ((*(ptIntegMeca.SigHH()))(1,1) < 0.))
// { cout << "\n Loi_comp_abstraite::CalculInvariants "
// << " invariant(1)= "<<invariant(1)<< ", (*(ptIntegMeca.SigHH()))(1,1)= "<< (*(ptIntegMeca.SigHH()))(1,1)
// << endl;
// };
////--- fin debug
};
TenseurHB * ptHB = &sigHB; delete ptHB;
};
// calcul éventuel des invariants liés uniquement à la cinématique
void Loi_comp_abstraite::CalculInvariants_cinematique
(PtIntegMecaInterne& ptIntegMeca, TenseurBB & gijBB,TenseurHH & gijHH)
{ // def de la dimension des tenseurs
int dimT = ptIntegMeca.EpsBB()->Dimension();
int caas=0; // pour la gestion d'erreur éventuelle
if (ptIntegMeca.Statut_Invariants_deformation())
{ TenseurHB& epsHB = *(NevezTenseurHB(dimT)) ;
epsHB = gijHH * (*(ptIntegMeca.EpsBB()));
Vecteur & invariant = ptIntegMeca.EpsInvar();
invariant.Change_taille(Abs(dimT));
switch (abs(dimT))
{case 3: invariant(3)= epsHB.Det();
case 2: invariant(2)= epsHB.II();
case 1: invariant(1)= epsHB.Trace();
};
// ptIntegMeca.EpsInvar() = epsHB.ValPropre(caas);
// if (caas == -1)
// { cout << "\n warning *** erreur dans le calcul des valeurs propres de la deformation";
// if (ParaGlob::NiveauImpression() >= 7) {epsHB.Ecriture(cout); cout << "\nLoi_comp_abstraite::CalculInvariants(...";};
// cout << endl;
// };
TenseurHB * ptHB = &epsHB; delete ptHB;
};
if (ptIntegMeca.Statut_Invariants_vitesseDeformation())
{ TenseurHB& DepsHB = *(NevezTenseurHB(dimT)) ;
DepsHB = gijHH * (*(ptIntegMeca.DepsBB()));
Vecteur & invariant = ptIntegMeca.DepsInvar();
invariant.Change_taille(Abs(dimT));
switch (abs(dimT))
{case 3: invariant(3)= DepsHB.Det();
case 2: invariant(2)= DepsHB.II();
case 1: invariant(1)= DepsHB.Trace();
};
TenseurHB * ptHB = &DepsHB; delete ptHB;
// ptIntegMeca.DepsInvar() = DepsHB.ValPropre(caas);
// if (caas == -1)
// { cout << "\n warning *** erreur dans le calcul des valeurs propres de la vitesse de deformation";
// if (ParaGlob::NiveauImpression() >= 7) {DepsHB.Ecriture(cout); cout << "\nLoi_comp_abstraite::CalculInvariants(...";};
// cout << endl;
// };
};
};
// récupération de valeurs interpolées pour les grandeur enu ou directement calculées
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
// une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul
// exclure_dd_etend: donne une liste de Ddl_enum_etendu à exclure de la recherche
// parce que par exemple, ils sont calculés par ailleurs
// on peut également ne pas définir de métrique, dans ce cas on ne peut pas calculer certaines grandeurs
// -> il y a vérification
Tableau <double> Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer
(bool absolue, Enum_dure temps,const List_io<Ddl_enum_etendu>& enu
,const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat
,const List_io<Ddl_enum_etendu>* exclure_dd_etend
)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ // def du tableau de retour: init 0.
Tableau <double> tab_ret (enu.size(),0.);
// si on détecte qu'il n'y a pas de valeur à calculer, on retourne directement
if (enu.size() == 0)
return tab_ret;
// sinon on continue
// on a besoin a priori de ptintmeca_en_cours, donc s'il n'est pas définit -> erreur
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
#ifdef MISE_AU_POINT
if (ptintmeca_en_cours == NULL)
{cout << "\n *** cas non prevu : aucun conteneur de point d'integration transmis "
<< "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(..." << endl;
Sortie(1);
};
if (def_en_cours == NULL)
{cout << "\n *** cas non prevu : la deformation n'est pas transmise "
<< "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(..." << endl;
Sortie(1);
};
#endif
// 2) le fait que l'on veut une sortie dans une base ad hoc ou pas
int dim = ptintmeca_en_cours->EpsBB_const().Dimension();
int dim_sortie_tenseur = dim;
// dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue
mise à jour 7.010
2023-05-03 17:23:49 +02:00
int dim_espace = ParaGlob::Dimension();
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (absolue)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
dim_sortie_tenseur = dim_espace;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// --- pour ne faire qu'un seul test ensuite
bool prevoir_change_dim_tenseur = false;
if (absolue)
// mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire
// car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue
// il faut un tenseur intermédiaire qui a ou non une dimension différente
prevoir_change_dim_tenseur = true;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// si on n'est pas en absolue, donc on est en ad hoc (c'est l'un ou l'autre)
// et si la dim ad hoc est différente de dim_espace, alors il faut aussi
// un tenseur intermédiaire mais de dim ad hoc
bool transfert_tenseur = false;
if ((!absolue)&&(dim != dim_espace))
{dim_sortie_tenseur = dim; // la dimension des tenseurs intermédiaires créés
transfert_tenseur = true;
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat=0;
if (Abs(deltat) >= ConstMath::trespetit)
{unSurDeltat = 1./deltat;}
else // si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
if (unSurDeltat < 0)
{ cout << "\n le pas de temps est négatif !! "; };
unSurDeltat = ConstMath::tresgrand;
};
// -- def des tenseurs locaux
mise à jour 7.010
2023-05-03 17:23:49 +02:00
TenseurHH* sigHH=NULL;TenseurBB* epsBB=NULL;TenseurBB* epslogBB=NULL;
TenseurBB* epsAlmBB=NULL;TenseurBB* eps0BB=NULL;TenseurBB* DepsBB = NULL;
TenseurBB* epsAlmTotalBB=NULL;TenseurBB* epsGLTotalBB=NULL;TenseurBB* epsLogTotalBB=NULL;
TenseurBB* DeltaEpsBB = NULL;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
Coordonnee* Mtdt = NULL; // coordonnées finales éventuelles du point d'intégration considéré
Coordonnee* Mt=NULL; // coordonnées à t
Coordonnee* M0 = NULL; // coordonnées initiales éventuelles du point d'intégration considéré
Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal actuel si c'est adéquate
Coordonnee* N_surf_t = NULL; // coordonnée d'un vecteur normal à t si c'est adéquate
Coordonnee* N_surf_t0 = NULL; // coordonnée d'un vecteur normal à t0 si c'est adéquate
Coordonnee* Vitesse = NULL; // cas des vitesses
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Coordonnee* Deplacement = NULL; // cas du déplacement
// pour les valeurs propres
Coordonnee* valPropreSig=NULL;Coordonnee* valPropreEps=NULL;
Coordonnee* valPropreDeps=NULL;
// pour les vecteurs propres
Tableau <Coordonnee *> base_propre_sigma , base_propre_eps , base_propre_D;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// pour les bases
Tableau <Coordonnee *> base_ad_hoc , base_giH , base_giB;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// grandeurs scalaires
double* Mises = NULL; double* defDualMises = NULL; double* Tresca = NULL;
double* erreur = NULL;double* erreur_rel = NULL;
double* spherique_eps=NULL,* Q_eps=NULL,* cos3phi_eps=NULL;
double* spherique_sig=NULL,* Q_sig=NULL,* cos3phi_sig=NULL;
double* spherique_Deps=NULL,* Q_Deps=NULL,* cos3phi_Deps=NULL;
double* def_equivalente=NULL, * defDualMisesMaxi=NULL;
double* sig_mises= NULL; // *** peut-être en trop
double* sig_tresca = NULL; // *** peut-être en trop
// --- dev d'un ensemble de variable booléenne pour gérer les sorties en une passe -----
// on se réfère au informations définit dans la méthode: Les_type_evolues_internes()
bool deformationCourante=false; bool vitesseDeformationCourante=false;
bool almansi=false; bool greenLagrange=false;
bool logarithmique=false; bool deltaDef=false;
bool almansiTotal=false; bool greenLagrangeTotale=false;
bool logarithmiqueTotale=false; bool contrainteCourante=false;
bool defPrincipales=false; bool sigmaPrincipales=false;
bool vitPrincipales=false; bool contrainteMises=false;
// bool contraintesTresca=false; bool erreurQ=false;
// bool erreurRel = false;
bool defPlastiqueCumulee=false; bool def_duale_mises=false;
bool besoin_des_contraintes=false; bool besoin_des_deformation=false;
bool besoin_des_contraintes_barre=false; bool besoin_des_deformation_barre=false;
bool besoin_des_vitesses_deformation=false; bool besoin_des_vitesses_deformation_barre=false;
bool besoin_des_valpropre_sigma=false;
bool besoin_des_valpropre_deformation = false; bool besoin_des_valpropre_vitdef = false;
bool besoin_coordonnees = false;
bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false;
bool besoin_dir_princ_sig = false; bool besoin_dir_princ_eps = false;
bool besoin_dir_princ_D = false;
bool besoin_rep_local_ortho=false;
bool besoin_rep_giH = false; bool besoin_rep_giB = false;
bool def_SPHERIQUE_EPS = false; bool def_Q_EPS = false; bool def_COS3PHI_EPS = false;
bool def_SPHERIQUE_SIG = false; bool def_Q_SIG = false; bool def_COS3PHI_SIG = false;
bool def_SPHERIQUE_DEPS = false; bool def_Q_DEPS = false; bool def_COS3PHI_DEPS = false;
bool def_def_equivalente = false; bool def_duale_mises_maxi = false;
bool def_sig_mises= false; // *** peut-être en trop
bool def_sig_tresca= false; // *** peut-être en trop
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// éléments de métrique et matrices de passage
TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH;
BaseB* giB_0;BaseB* giB_t;
Mat_pleine jB0(dim,dim),jBfin(dim,dim);
bool pas_de_metrique_dispo = false; // init
if (ex_impli != NULL)
{ gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt;
giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt;
giB_0 = ex_impli->giB_0;giB_t = ex_impli->giB_t;
}
else if (ex_expli_tdt != NULL)
{gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt;
giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt;
giB_0 = ex_expli_tdt->giB_0; giB_t = ex_expli_tdt->giB_t;
}
else if (ex_umat != NULL)
{gijHH = ex_umat->gijHH_t;gijBB = ex_umat->gijBB_t;
giB = giB_t = ex_umat->giB_t; giH_0 = ex_umat->giH_0;giH = ex_umat->giH_t;
giB_0 = ex_umat->giB_0;
}
else
{ pas_de_metrique_dispo = true;
// cout << "\n *** cas non prevu : aucune metrique transmise "
// << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(..." << endl;
// Sortie(1);
};
// on définie des indicateurs pour ne pas faire plusieurs fois le même calcul
List_io<Ddl_enum_etendu>::const_iterator ie,iefin=enu.end();
mise à jour 7.010
2023-05-03 17:23:49 +02:00
bool besoin_des_vitesse_deformation_barre=false;
bool besoin_deformation_greenlagrange = false;
bool besoin_deformation_logarithmique = false;
bool besoin_deformation_almansi = false;
bool besoin_deplacements = false;
// init pour l'utilisation de la liste d'exclusion
List_io<Ddl_enum_etendu>::const_iterator itt_deb;
List_io<Ddl_enum_etendu>::const_iterator itt_fin;
if (exclure_dd_etend != NULL)
{ itt_deb = exclure_dd_etend->begin();
itt_fin = exclure_dd_etend->end();
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
for (ie=enu.begin(); ie!=iefin;ie++)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{// on commence par regarder s'il faut traiter ou pas
bool a_atraiter = true;
if (exclure_dd_etend != NULL)
if (find(itt_deb,itt_fin,(*ie)) != itt_fin)
a_atraiter = false;
if (a_atraiter)
{ if (Meme_famille((*ie).Enum(),SIG11)) besoin_des_contraintes=true;
if (Meme_famille((*ie).Enum(),EPS11)) besoin_des_deformation=true;
if (Meme_famille((*ie).Enum(),DEPS11)) besoin_des_vitesses_deformation=true;
if (Meme_famille((*ie).Enum(),X1)) besoin_coordonnees=true;
if (Meme_famille((*ie).Enum(),UX)) {besoin_deplacements=true;besoin_coordonnees=true;};
int posi = (*ie).Position()-NbEnum_ddl();
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
switch (posi)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{case 1: case 2: case 3: case 4: case 5: case 6:
{besoin_deformation_greenlagrange=true;
eps0BB = NevezTenseurBB(dim_sortie_tenseur) ;break;}
case 7: case 8: case 9: case 10: case 11: case 12:
{besoin_deformation_almansi=true;
epsAlmBB =NevezTenseurBB(dim_sortie_tenseur) ;break;}
case 28: case 29: case 30: case 31: case 32:
{besoin_des_valpropre_sigma=true;
valPropreSig = new Coordonnee(dim_sortie_tenseur);break;}
case 25: case 26: case 27: case 77:
{besoin_des_valpropre_deformation=true;
valPropreEps = new Coordonnee(dim_sortie_tenseur);break;}
case 40: case 41: case 42:
{besoin_des_valpropre_vitdef=true;
valPropreDeps = new Coordonnee(dim_sortie_tenseur);break;}
case 49: case 50: case 51: case 52: case 53: case 54:
{besoin_deformation_logarithmique=true;
epslogBB =NevezTenseurBB(dim_sortie_tenseur) ;break;}
case 55: case 56: case 57: case 58: case 59: case 60:
{if ((epsAlmTotalBB == NULL) && (dilatation))
{epsAlmTotalBB = (NevezTenseurBB(dim_sortie_tenseur));};
break;
}
case 61: case 62: case 63: case 64: case 65: case 66:
{ if ((epsGLTotalBB == NULL) && (dilatation))
{epsGLTotalBB = (NevezTenseurBB(dim_sortie_tenseur));};
break;
}
case 67: case 68: case 69: case 70: case 71: case 72:
{if ((epsLogTotalBB == NULL) && (dilatation))
{epsLogTotalBB = (NevezTenseurBB(dim_sortie_tenseur));};
break;
}
case 78: case 79: case 80:
{besoin_des_deformation_barre=true;break;}
case 81: case 82: case 83:
{besoin_des_contraintes_barre=true;break;}
case 84: case 85: case 86:
{besoin_des_vitesse_deformation_barre=true;break;}
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
case 114: case 115: case 116: // le vecteur normal
{ N_surf = new Coordonnee(ParaGlob::Dimension()); break;}
case 117: case 118: case 119: // le vecteur normal à t
{ N_surf_t = new Coordonnee(ParaGlob::Dimension()); break;}
case 120: case 121: case 122: // le vecteur normal à t0
{ N_surf_t0 = new Coordonnee(ParaGlob::Dimension()); break;}
case 123: case 124: case 125: // la position à t
{ Mt = new Coordonnee(ParaGlob::Dimension());
besoin_coordonnees_t = true;
break;
}
case 126: case 127: case 128: // la position à t0
{ M0 = new Coordonnee(ParaGlob::Dimension());
besoin_coordonnees_t0 = true;
break;
}
default:
break;
};
};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// on complète les définitions de tenseurs si besoin
// les tenseurs restants en locale
TenseurHB* sigHB = NULL ;TenseurHB* sig_barreHB = NULL ;
TenseurHB* epsHB = NULL ;TenseurHB* eps_barreHB = NULL ;
TenseurHB* DepsHB = NULL ; TenseurHB* Deps_barreHB = NULL ;
// -- def de tenseurs pour la sortie
if ((besoin_des_contraintes) || (besoin_des_valpropre_sigma))
{sigHH = NevezTenseurHH(dim_sortie_tenseur) ;
sigHB = (NevezTenseurHB(dim)) ;
sig_barreHB = (NevezTenseurHB(dim)) ;
}
if ((besoin_des_deformation) || (besoin_des_valpropre_deformation))
{epsBB = NevezTenseurBB(dim_sortie_tenseur) ;
epsHB = (NevezTenseurHB(dim)) ;
eps_barreHB = (NevezTenseurHB(dim)) ;
}
if ((besoin_des_vitesses_deformation) || (besoin_des_valpropre_vitdef))
{DepsBB = NevezTenseurBB(dim_sortie_tenseur);
DepsHB = (NevezTenseurHB(dim));
Deps_barreHB = (NevezTenseurHB(dim)) ;
}
// on statut sur le fait d'avoir besoin ou non des chgt de base
bool besoin_matrice_chg_base = false;
if ( besoin_des_contraintes
|| besoin_des_deformation
|| besoin_des_vitesses_deformation
|| besoin_rep_local_ortho
)
{besoin_matrice_chg_base = true;};
// on ne change pas le numéro de point d'intégration courant
// on considère ici que c'est déjà le bon, puisque c'est appelé par un élément
// a priori ...
// on recupère le tableau pour la lecture des coordonnées des tenseurs
int nbcompo = ParaGlob::NbCompTens();
// définition des grandeurs qui sont indépendante de la boucle sur les ddl_enum_etendue
// matrices de passage
int dim_effective = dim; // init
if (absolue) dim_effective = ParaGlob::Dimension();
Mat_pleine* Aa0 = NULL;Mat_pleine* Aafin = NULL;
Mat_pleine* gamma0 = NULL;Mat_pleine* gammafin = NULL;
Mat_pleine* beta0 = NULL;Mat_pleine* betafin = NULL;
if (besoin_matrice_chg_base)
// dans le cas où on n'est pas en absolue => on sort dans un repère ad hoc donc
// il a la dimension locale
// sinon on sort dans le repère globale => il a la dimension globale
{Aa0 = new Mat_pleine(dim_effective,dim_effective);
Aafin = new Mat_pleine(dim_effective,dim_effective);
gamma0 = new Mat_pleine(dim_effective,dim_effective);
gammafin = new Mat_pleine(dim_effective,dim_effective);
beta0 = new Mat_pleine(dim_effective,dim_effective);
betafin = new Mat_pleine(dim_effective,dim_effective);
};
// on considère que les métriques sont directement utilisables
// si elles sont disponibles
if (besoin_matrice_chg_base)
{if (pas_de_metrique_dispo)
{cout << "\n *** erreur : on a besoin de matrices de changement de base et il n'y a pas "
<< " de metrique disponible: impossible de continuer"
<< "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer..." << endl;
Sortie(1);
}
else if (def_en_cours== NULL)
{cout << "\n *** erreur : on a besoin de l'objet Deformation et il n'y a pas "
<< " de pointeur de deformation disponible: impossible de continuer"
<< "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer..." << endl;
Sortie(1);
}
else
// on peut calculer
{Mat_pleine Aat(dim_effective,dim_effective);
// a priori Aat ne sert pas par la suite, mais est nécessaire pour le passage de par
const Met_abstraite::Info_et_metrique_0_t_tdt ex
= def_en_cours->Remont_et_metrique_0_t_tdtSansCalMet(absolue,*Aa0,Aat,*Aafin);
// pour les formules de passage de repère il nous faut :
// Ip_a = beta_a^{.j} g_j et Ip^b = gamma^b_{.j} g^j
// on a: [beta_a^{.j}] = [Aa^j_{.a}]^T
// et [gamma^b_{.j}] = [beta_a^{.j}]^{-1T} = [Aa^j_{.a}]^{-1}
(*gamma0) = (Aa0->Inverse());
(*gammafin) = (Aafin->Inverse());
// on détermine également les matrices beta
(*beta0) = (Aa0->Transpose());
(*betafin) = (Aafin->Transpose());
};
};
// récup des bases si besoin
if (besoin_rep_local_ortho)
{if ((!absolue)&&(dim_espace==3)&&(dim==2))
// cas d'éléments 2D, pour lesquels on veut un repère local ad hoc
// on ramène une base à 3 vecteurs
{ *(base_ad_hoc(1)) = jBfin(1,1) * giB->Coordo(1) + jBfin(2,1) * giB->Coordo(2);
*(base_ad_hoc(2)) = jBfin(1,2) * giB->Coordo(1) + jBfin(2,2) * giB->Coordo(2);
*(base_ad_hoc(3)) = Util::ProdVec_coor(*(base_ad_hoc(1)),*(base_ad_hoc(2)));
}
else if((!absolue)&&(dim_espace>1)&&(dim==1))
// cas d'éléments 1D, dans un espace 2D ou 3D
// on ramène un seul vecteur non nul, les autres ne peuvent être calculé sans info supplémentaire
{ *(base_ad_hoc(1)) = jBfin(1,1) * giB->Coordo(1);
(base_ad_hoc(2))->Zero(); (base_ad_hoc(3))->Zero(); // init à 0 par défaut
}
else // dans tous les autres cas
{ switch (dim_espace) // on se contente de ramener le repère identité
{case 3: (base_ad_hoc(3))->Zero();(*(base_ad_hoc(3)))(3)=1.;
case 2: (base_ad_hoc(2))->Zero();(*(base_ad_hoc(2)))(2)=1.;
case 1: (base_ad_hoc(1))->Zero();(*(base_ad_hoc(1)))(1)=1.;
default:break;
};
};
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (besoin_rep_giH)
{switch (dim) {case 3: *(base_giH(3)) = giH->Coordo(3);
case 2: *(base_giH(2)) = giH->Coordo(2);
case 1: *(base_giH(1)) = giH->Coordo(1);
default:break;
};
};
if (besoin_rep_giB)
{switch (dim) {case 3: *(base_giB(3)) = giB->Coordo(3);
case 2: *(base_giB(2)) = giB->Coordo(2);
case 1: *(base_giB(1)) = giB->Coordo(1);
default:break;
};
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// ----- maintenant on calcule les grandeurs nécessaires -----
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// calcul des tenseurs initiaux
bool plusZero = true; // s'il faut rajouter des termes, on met des 0
if (besoin_des_contraintes)
{(*sigHB) = (ptintmeca_en_cours->SigHH_const()) * (*gijBB);
if (contrainteCourante) // on n'intervient ici que si on veut une sortie des contraintes courantes
{if (absolue) {(ptintmeca_en_cours->SigHH_const()).BaseAbsolue(*sigHH,*giB);}// changement de base finale
else if (transfert_tenseur)
{ TenseurHH* sigHH_inter = NevezTenseurHH(dim) ;
*sigHH_inter = (ptintmeca_en_cours->SigHH_const()); sigHH_inter->ChBase(*gammafin);
sigHH->Affectation_trans_dimension(*sigHH_inter,false);
delete sigHH_inter;
}
else
{ *sigHH = (ptintmeca_en_cours->SigHH_const());sigHH->ChBase(*gammafin);};
};
};
if (besoin_des_deformation)
{// cas de delta_eps
if(DeltaEpsBB != NULL)
{if (absolue)// changement de base finale
{(ptintmeca_en_cours->DeltaEpsBB_const()).BaseAbsolue(*DeltaEpsBB,*giH);}
else if (transfert_tenseur)
{ TenseurBB* DeltaEpsBB_inter = NevezTenseurBB(dim) ;
*DeltaEpsBB_inter = (ptintmeca_en_cours->DeltaEpsBB_const()); DeltaEpsBB_inter->ChBase(*betafin);
DeltaEpsBB->Affectation_trans_dimension(*DeltaEpsBB_inter,false);
delete DeltaEpsBB_inter;
}
else {*DeltaEpsBB = (ptintmeca_en_cours->DeltaEpsBB_const());DeltaEpsBB->ChBase(*betafin);};
};
// cas de la déformation
(*epsHB) = (*gijHH) * ((ptintmeca_en_cours->EpsBB_const()));
if (deformationCourante)// on n'intervient ici que si on veut une sortie des déformations courantes
{if (absolue)// changement de base finale
{(ptintmeca_en_cours->EpsBB_const()).BaseAbsolue(*epsBB,*giH);}
else if (transfert_tenseur)
{ TenseurBB* epsBB_inter = NevezTenseurBB(dim) ;
*epsBB_inter = (ptintmeca_en_cours->EpsBB_const()); epsBB_inter->ChBase(*betafin);
epsBB->Affectation_trans_dimension(*epsBB_inter,false);
delete epsBB_inter;
}
else {*epsBB = (ptintmeca_en_cours->EpsBB_const());epsBB->ChBase(*betafin);};
};
switch (def_en_cours->Type_de_deformation())
{case DEFORMATION_STANDART : // c'est à dire almansi
{if (greenLagrange)
{if (absolue)
{(ptintmeca_en_cours->EpsBB_const()).BaseAbsolue(*eps0BB,*giH_0);}// changement de base finale
else if (transfert_tenseur)
{ TenseurBB* eps0BB_inter = NevezTenseurBB(dim) ;
*eps0BB_inter = (ptintmeca_en_cours->EpsBB_const()); eps0BB_inter->ChBase(*beta0);
eps0BB->Affectation_trans_dimension(*eps0BB_inter,false);
delete eps0BB_inter;
}
else {eps0BB->ChBase(*beta0);};
};
if (almansi) *epsAlmBB = *epsBB;// le changement de base a été fait juste plus haut
if (almansiTotal) // cas avec dilatation et demande de def Almansi totale
{TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsAlmTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsAlmTotal_local_BB);
if (absolue)// changement de base finale
{epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);}
else if (transfert_tenseur)
{ epsAlmTotal_local_BB->ChBase(*betafin);
epsAlmTotalBB->Affectation_trans_dimension(*epsAlmTotal_local_BB,false);
}
else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite
};
if (greenLagrangeTotale) // cas avec dilatation et demande de def Green_Lagrange totale
{TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsGLTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsGLTotal_local_BB);
if (absolue)// changement de base finale
{epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);}
else if (transfert_tenseur)
{ epsGLTotal_local_BB->ChBase(*beta0);
epsGLTotalBB->Affectation_trans_dimension(*epsGLTotal_local_BB,false);
}
else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite
};
// si l'on veut sortir la déformation logarithmique le plus simple est de la calculer
if (logarithmiqueTotale || logarithmique)
{def_en_cours->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE);
if (logarithmique) // cas du calcul de la def logarithmique
{TenseurBB* epslog_local_BB = epslogBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epslog_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epslog_local_BB);
if (absolue)// changement de base finale
{epslog_local_BB->BaseAbsolue(*epslogBB,*giH);}
else if (transfert_tenseur)
{ epslog_local_BB->ChBase(*betafin);
epslogBB->Affectation_trans_dimension(*epslog_local_BB,false);
}
else {epslogBB->ChBase(*betafin);}; // ici epslog_local_BB == epslogBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epslog_local_BB; // car pas utilisé ensuite
};
if (logarithmiqueTotale) // cas avec dilatation et demande de def log totale
{TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsLogTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsLogTotal_local_BB);
if (absolue)// changement de base finale
{epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);}
else if (transfert_tenseur)
{ epsLogTotal_local_BB->ChBase(*betafin);
epsLogTotalBB->Affectation_trans_dimension(*epsLogTotal_local_BB,false);
}
else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite
};
def_en_cours->Change_type_de_deformation(DEFORMATION_STANDART); // on revient au type initial
};
break;
}
case DEFORMATION_LOGARITHMIQUE :
{ if (logarithmique) *epslogBB=*epsBB;
if (logarithmiqueTotale) // cas avec dilatation et demande de def log totale
{TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsLogTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsLogTotal_local_BB);
if (absolue)// changement de base finale
{epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);}
else if (transfert_tenseur)
{ epsLogTotal_local_BB->ChBase(*betafin);
epsLogTotalBB->Affectation_trans_dimension(*epsLogTotal_local_BB,false);
}
else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite
};
// si l'on veut sortir la déformation d'Almansi ou de green-lagrange le plus simple est de les calculer
if (almansi || greenLagrangeTotale || almansiTotal)
{def_en_cours->Change_type_de_deformation(DEFORMATION_STANDART);
if (almansi) // cas de la def d'almansi
{ TenseurBB* eps_local_BB = epsAlmBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
eps_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*eps_local_BB);
if (absolue)// changement de base finale
{eps_local_BB->BaseAbsolue(*epsAlmBB,*giH);}
else if (transfert_tenseur)
{ eps_local_BB->ChBase(*betafin);
epsAlmBB->Affectation_trans_dimension(*eps_local_BB,false);
}
else {epsAlmBB->ChBase(*betafin);};// ici eps_local_BB == epsAlmBB
if(greenLagrange)
{if (absolue)// changement de base finale
{eps_local_BB->BaseAbsolue(*eps0BB,*giH_0);}
else if (transfert_tenseur)
{ eps_local_BB->ChBase(*beta0);
eps0BB->Affectation_trans_dimension(*eps_local_BB,false);
}
else {eps0BB->ChBase(*beta0);}; // ici eps_local_BB == eps0BB
};
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete eps_local_BB; // car pas utilisé ensuite
};
if (almansiTotal) // cas avec dilatation et demande de def Almansi totale
{TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsAlmTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsAlmTotal_local_BB);
if (absolue)// changement de base finale
{epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);}
else if (transfert_tenseur)
{ epsAlmTotal_local_BB->ChBase(*betafin);
epsAlmTotalBB->Affectation_trans_dimension(*epsAlmTotal_local_BB,false);
}
else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite
};
if (greenLagrangeTotale) // cas avec dilatation et demande de def Green_Lagrange totale
{TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsGLTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsGLTotal_local_BB);
if (absolue)// changement de base finale
{epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);}
else if (transfert_tenseur)
{ epsGLTotal_local_BB->ChBase(*beta0);
epsGLTotalBB->Affectation_trans_dimension(*epsGLTotal_local_BB,false);
}
else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite
};
def_en_cours->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); // on revient au type initial
};
break;
}
default:
cout << "\n cas de deformation non encore implante en sortie de visualisation "
<< Nom_type_deformation(def_en_cours->Type_de_deformation())
<< " affichage donc errone des valeurs !!!";
};
}; //-- fin de if (besoin_des_deformation)
if (besoin_des_vitesses_deformation)
{ *DepsHB = (*gijHH) * ((ptintmeca_en_cours->DepsBB_const()));
if (absolue)// changement de base finale
{(ptintmeca_en_cours->DepsBB_const()).BaseAbsolue(*DepsBB,*giH);}
else if (transfert_tenseur)
{ TenseurBB* DepsBB_inter = NevezTenseurBB(dim) ;
*DepsBB_inter = (ptintmeca_en_cours->DepsBB_const()); DepsBB_inter->ChBase(*betafin);
DepsBB->Affectation_trans_dimension(*DepsBB_inter,false);
delete DepsBB_inter;
}
else {*DepsBB = (ptintmeca_en_cours->DepsBB_const());DepsBB->ChBase(*betafin);};
};
if (besoin_des_contraintes_barre)
{double Isig = sigHB->Trace(); // trace de la déformation
*sig_barreHB = (*sigHB) - (Isig/dim_espace) * (*Id_dim_HB(dim));
};
if (besoin_des_deformation_barre)
{double Ieps = epsHB->Trace(); // trace de la déformation
*eps_barreHB = (*epsHB) - (Ieps/dim_espace) * (*Id_dim_HB(dim));
};
if (besoin_des_vitesses_deformation_barre)
{double IDeps = DepsHB->Trace(); // trace de la déformation
*Deps_barreHB = (*DepsHB) - (IDeps/dim_espace) * (*Id_dim_HB(dim));
};
// cas des valeurs propres et éventuellement des vecteurs propres
{int caas=0;
////----- debug
//cout << "\n besoin_des_valpropre_sigma= "<< besoin_des_valpropre_sigma
// << " besoin_dir_princ_sig= "<< besoin_dir_princ_sig << endl;
////--- fin debug
if (besoin_des_valpropre_sigma && besoin_dir_princ_sig)// on veut les directions et les valeurs propres
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(sigHB->ValPropre(caas,mat));
inter.Change_dim(dim_sortie_tenseur); // on étant la taille éventuellement
*valPropreSig = inter; // sauvegarde des valeurs propres
// puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_sigma(3)),tab(3));(base_propre_sigma(3))->Normer();
giB->BaseAbsolue( *(base_propre_sigma(2)),tab(2));(base_propre_sigma(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
break;
case 2: // en 2D le premier vecteur en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
giH->BaseAbsolue( *(base_propre_sigma(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_sigma(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_dir_princ_sig) // cas où on ne veut que les directions principales
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(sigHB->ValPropre(caas,mat));
// sauvegarde des directions principales
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_sigma(3)),tab(3));(base_propre_sigma(3))->Normer();
giB->BaseAbsolue( *(base_propre_sigma(2)),tab(2));(base_propre_sigma(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
giH->BaseAbsolue( *(base_propre_sigma(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_sigma(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_des_valpropre_sigma)// on ne veut que les valeurs propres
{ Coordonnee inter(sigHB->ValPropre(caas));
inter.Change_dim(dim_sortie_tenseur); // on étend la taille éventuellement
*valPropreSig = inter; // sauvegarde des valeurs propres
};
// gestion d'une erreur éventuelle
if (caas == -1)
{ cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la contrainte ";
if (ParaGlob::NiveauImpression() > 5) {sigHB->Ecriture(cout);cout << "\n cas = "<< caas ;
cout << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(...";};
cout << endl;
};
}; // -- fin de l'encapsulation du cas des contraintes
{int caas=0;
if (besoin_des_valpropre_deformation && besoin_dir_princ_eps)// on veut les directions et les valeurs propres
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
// 1) calcul des valeurs propres dans un vecteur de travail inter, et des vecteurs propres stockés
Coordonnee inter(epsHB->ValPropre(caas,mat)); // dans la matrice mat
inter.Change_dim(dim_sortie_tenseur); // on étend la taille, dans le cas ou le nb de valeurs
// propre est inférieur à la dimension de l'espace
*valPropreEps = inter; // sauvegarde des valeurs propres
// 2) puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_eps(3)),tab(3));(base_propre_eps(3))->Normer();
giB->BaseAbsolue( *(base_propre_eps(2)),tab(2));(base_propre_eps(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
giH->BaseAbsolue( *(base_propre_eps(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_eps(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_dir_princ_eps) // cas où on ne veut que les directions principales
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(epsHB->ValPropre(caas,mat)); // dans la matrice mat
// puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_eps(3)),tab(3));(base_propre_eps(3))->Normer();
giB->BaseAbsolue( *(base_propre_eps(2)),tab(2));(base_propre_eps(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
giH->BaseAbsolue( *(base_propre_eps(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_eps(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_des_valpropre_deformation)// on ne veut que les valeurs propres
{ Coordonnee inter(epsHB->ValPropre(caas));
inter.Change_dim(dim_sortie_tenseur); // on étant la taille éventuellement
*valPropreEps = inter; // sauvegarde des valeurs propres
};
// gestion d'une erreur éventuelle
if (caas == -1)
{ cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la deformation";
if (ParaGlob::NiveauImpression() >= 7) {epsHB->Ecriture(cout);cout << "\n cas = "<< caas ;
cout << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(...";};
cout << endl;
};
}; // -- fin de l'encapsulation du cas des déformations
{int caas=0; // il faut l'initialiser sinon il peut prendre la valeur du cas précédant
if (besoin_des_valpropre_vitdef && besoin_dir_princ_D)// on veut les directions et les valeurs propres
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
// 1) calcul des valeurs propres dans un vecteur de travail inter, et des vecteurs propres stockés
Coordonnee inter(DepsHB->ValPropre(caas,mat)); // dans la matrice mat
inter.Change_dim(dim_sortie_tenseur); // on étend la taille, dans le cas ou le nb de valeurs
// propre est inférieur à la dimension de l'espace
*valPropreDeps = inter; // sauvegarde des valeurs propres
// 2) puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_D(3)),tab(3));(base_propre_D(3))->Normer();
giB->BaseAbsolue( *(base_propre_D(2)),tab(2));(base_propre_D(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
giH->BaseAbsolue( *(base_propre_D(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_D(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_dir_princ_D) // cas où on ne veut que les directions principales
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(DepsHB->ValPropre(caas,mat)); // dans la matrice mat
// puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_D(3)),tab(3));(base_propre_D(3))->Normer();
giB->BaseAbsolue( *(base_propre_D(2)),tab(2));(base_propre_D(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
giH->BaseAbsolue( *(base_propre_D(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_D(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_des_valpropre_vitdef)// on ne veut que les valeurs propres
{ Coordonnee inter(DepsHB->ValPropre(caas));
inter.Change_dim(dim_sortie_tenseur); // on étant la taille
*valPropreDeps = inter; // sauvegarde des valeurs propres
};
// gestion d'une erreur éventuelle
if (caas == -1)
{ cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la vitesse de deformation";
if (ParaGlob::NiveauImpression() >= 7) {DepsHB->Ecriture(cout);cout << "\n cas = "<< caas ;
cout << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(...";};
cout << endl;
};
}; // -- fin de l'encapsulation du cas des vitesses de déformations
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (besoin_coordonnees)
{Mtdt = new Coordonnee(ParaGlob::Dimension());
*Mtdt = def_en_cours->Position_tdt();
};
if (besoin_coordonnees_t )
{*Mt = def_en_cours->Position_tdt();
};
if (besoin_deplacements || besoin_coordonnees_t0)
{if (M0 == NULL)
M0 = new Coordonnee(ParaGlob::Dimension());
(*M0) = def_en_cours->Position_0();
};
if (Vitesse != NULL)
{Vitesse = new Coordonnee(ParaGlob::Dimension());
(*Vitesse) = def_en_cours->VitesseM_tdt();
};
// def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D
if (N_surf != NULL)
{ if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur N_surf necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// on vérifie que la métrique est correcte
if (giB->NbVecteur() != 2)
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
<< " le vecteur normal ne sera pas correctement calcule";
if (ParaGlob::NiveauImpression() > 2)
cout << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(... ";
};
cout << endl;
}
else // sinon c'est ok
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
(*N_surf) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2));
N_surf->Normer(); // que l'on norme
};
};
// idem à l'instant t
if (N_surf_t != NULL)
{ if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur N_surf_t necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// on vérifie que la métrique est correcte
if (giB_t->NbVecteur() != 2)
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
<< " le vecteur normal ne sera pas correctement calcule";
if (ParaGlob::NiveauImpression() > 2)
cout << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(... ";
};
cout << endl;
}
else // sinon c'est ok
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
(*N_surf_t) = Util::ProdVec_coorBN( (*giB_t)(1), (*giB_t)(2));
N_surf_t->Normer(); // que l'on norme
};
};
// idem à l'instant t0
if (N_surf_t0 != NULL)
{ if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur N_surf_0 necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// on vérifie que la métrique est correcte
if (giB_0->NbVecteur() != 2)
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
<< " le vecteur normal ne sera pas correctement calcule";
if (ParaGlob::NiveauImpression() > 2)
cout << "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(... ";
};
cout << endl;
}
else // sinon c'est ok
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
(*N_surf_t0) = Util::ProdVec_coorBN( (*giB_0)(1), (*giB_0)(2));
N_surf_t0->Normer(); // que l'on norme
};
// on n'arrête pas l'exécution, car on pourrait vouloir sortir les normales pour un ensemble
// d'éléments contenant des volumes, des surfaces, des lignes: bon... il y aura quand même des
// pb au niveau des iso par exemple, du au fait que l'on va faire des moyennes sur des éléments
// de type différents (à moins de grouper par type du coup on n'aura pas le warning
};
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// cas du déplacement et de la vitesse
if (Deplacement != NULL)
(*Deplacement) = (*Mtdt) - def_en_cours->Position_0();
if (Vitesse != NULL)
(*Vitesse) = def_en_cours->VitesseM_tdt();
// cas de la deformation équivalente cumulée
if (def_def_equivalente) {*def_equivalente = ptintmeca_en_cours->Deformation_equi_const()(1);};
// cas de la deformation duale au sens de mises,
if (def_duale_mises) {*defDualMises = ptintmeca_en_cours->Deformation_equi_const()(2);};
// cas de la deformation maxi duale au sens de mises,
if (def_duale_mises_maxi) {*defDualMisesMaxi = ptintmeca_en_cours->Deformation_equi_const()(3);};
// cas des contraintes: mises et tresca
if (def_sig_mises) {*sig_mises = ptintmeca_en_cours->Sig_equi_const()(1);};
if (def_sig_tresca) {*sig_tresca = ptintmeca_en_cours->Sig_equi_const()(2);};
// //contrainte_tresca
// if (contraintesTresca)
// { switch (dim) {case 1: *Tresca=0.5 * (*valPropreSig)(1);break;
// case 2: *Tresca=0.5 * ((*valPropreSig)(1)-(*valPropreSig)(2));break;
// case 3: *Tresca=0.5 * ((*valPropreSig)(1)-(*valPropreSig)(3));break;
// };
// };
// --- cas des grandeurs de la décomposition polaire
// cas de la déformation
if (def_SPHERIQUE_EPS) {*spherique_eps = epsHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012
double mini_Q = 5.e-5;
if (def_Q_EPS) {*Q_eps = sqrt(eps_barreHB->II());};
if (cos3phi_eps)
{ double Qepsilon = ( def_Q_EPS ? *Q_eps : sqrt(eps_barreHB->II()));
double Qepsilon3 = Qepsilon * Qepsilon * Qepsilon;
if (Qepsilon > mini_Q )
{ // on peut calculer un cos3phi pas débile
double bIIIb = eps_barreHB->III() / 3.;
*cos3phi_eps = 3. * sqrt(6.) * bIIIb/ Qepsilon3;
}
else *cos3phi_eps=0.; // sinon on le met à 0
};
// cas de la contrainte
if (def_SPHERIQUE_SIG) {*spherique_sig = sigHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012
if (def_Q_SIG) {*Q_sig = sqrt(sig_barreHB->II());};
if (cos3phi_sig)
{ double Qsig = ( def_Q_SIG ? *Q_sig : sqrt(sig_barreHB->II()));
double Qsig3 = Qsig * Qsig * Qsig;
if (Qsig > mini_Q )
{ // on peut calculer un cos3phi pas débile
double bIIIb = sig_barreHB->III() / 3.;
*cos3phi_sig = 3. * sqrt(6.) * bIIIb/ Qsig3;
}
else *cos3phi_sig=0.; // sinon on le met à 0
};
// cas de la vitesse de déformation
if (def_SPHERIQUE_DEPS) {*spherique_Deps = DepsHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012
if (def_Q_DEPS) {*Q_Deps = sqrt(Deps_barreHB->II());};
if (cos3phi_Deps)
{ double QDepsilon = ( def_Q_DEPS ? *Q_Deps : sqrt(Deps_barreHB->II()));
double QDepsilon3 = QDepsilon * QDepsilon * QDepsilon;
if (QDepsilon > mini_Q )
{ // on peut calculer un cos3phi pas débile
double bIIIb = Deps_barreHB->III() / 3.;
*cos3phi_Deps = 3. * sqrt(6.) * bIIIb/ QDepsilon3;
}
else *cos3phi_Deps=0.; // sinon on le met à 0
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
//----- fin du calcul des grandeurs nécessaires -----
mise à jour 7.010
2023-05-03 17:23:49 +02:00
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// on balaie maintenant la liste des grandeurs à sortir
int it; // it est l'indice dans le tableau de retour
for (it=1,ie=enu.begin(); ie!=iefin;ie++,it++)
{ if ((Meme_famille((*ie).Enum(),SIG11)) || (Meme_famille((*ie).Enum(),EPS11))
|| (Meme_famille((*ie).Enum(),DEPS11)) || (Meme_famille((*ie).Enum(),X1))
|| (Meme_famille((*ie).Enum(),UX)) )
{ // on recherche en générale une interpolation en fonction des noeuds: il faut donc que la grandeur soit
// présente aux noeuds !!
// def du numéro de référence du ddl_enum_etendue
int posi = (*ie).Position()-NbEnum_ddl();
// récupération des informations en fonction des différents cas
// **** 1 >>>>> -- cas des ddl pur, que l'on sort dans le repère global par défaut
// cas des contraintes
if ((Meme_famille((*ie).Enum(),SIG11)) && ((*ie).Nom_vide()))
{ tab_ret(it)= def_en_cours->DonneeInterpoleeScalaire((*ie).Enum(),temps);
}
else if ((Meme_famille((*ie).Enum(),EPS11)) && ((*ie).Nom_vide()))
{ tab_ret(it)= def_en_cours->DonneeInterpoleeScalaire((*ie).Enum(),temps);
}
else if ((Meme_famille((*ie).Enum(),DEPS11)) && ((*ie).Nom_vide()))
{ tab_ret(it)= def_en_cours->DonneeInterpoleeScalaire((*ie).Enum(),temps);
}
else if ((Meme_famille((*ie).Enum(),X1)) && ((*ie).Nom_vide()))
{ tab_ret(it)= (*Mtdt)((*ie).Enum() - X1 +1);
}
else if ((Meme_famille((*ie).Enum(),UX)) && ((*ie).Nom_vide()))
{ int i_cor = (*ie).Enum() - UX +1; // l'indice de coordonnée
tab_ret(it)= (*Mtdt)(i_cor) - (*M0)(i_cor);
}
else if ((Meme_famille((*ie).Enum(),V1)) && ((*ie).Nom_vide()))
{ int i_cor = (*ie).Enum() - V1 +1; // l'indice de coordonnée
tab_ret(it)= (*Vitesse)(i_cor);
}
// --- a complèter ----
else
{// **** 2 >>>>> -- cas des grandeurs déduites des ddl pures
switch (posi)
{ case 31: /*contrainte_mises*/
tab_ret(it)=ptintmeca_en_cours->Sig_equi_const()(1);break;
case 32: /* contrainte_tresca */
tab_ret(it)=ptintmeca_en_cours->Sig_equi_const()(2);break;
case 77: /*def_duale_mises*/ tab_ret(it)=ptintmeca_en_cours->Deformation_equi_const()(2);break;
case 87: /*def_equivalente*/ tab_ret(it)=ptintmeca_en_cours->Deformation_equi_const()(1);break;
case 88: /*def_duale_mises_maxi*/ tab_ret(it)=ptintmeca_en_cours->Deformation_equi_const()(3);break;
case 89: /*vitesse_def_equivalente*/ tab_ret(it)=ptintmeca_en_cours->Deformation_equi_const()(4) * unSurDeltat;break;
case 114: // le vecteur normal N_surf_1
{tab_ret(it)= (*N_surf)(1);break;}
case 115: // le vecteur normal N_surf_2
{tab_ret(it)= (*N_surf)(2);break;}
case 116: // le vecteur normal N_surf_3
{tab_ret(it)= (*N_surf)(3);break;}
case 117: // le vecteur normal N_surf_1_t
{tab_ret(it)= (*N_surf_t)(1);break;}
case 118: // le vecteur normal N_surf_2_t
{tab_ret(it)= (*N_surf_t)(2);break;}
case 119: // le vecteur normal N_surf_3_t
{tab_ret(it)= (*N_surf_t)(3);break;}
case 120: // le vecteur normal N_surf_1_t0
{tab_ret(it)= (*N_surf_t0)(1);break;}
case 121: // le vecteur normal N_surf_2_t0
{tab_ret(it)= (*N_surf_t0)(2);break;}
case 122: // le vecteur normal N_surf_3_t0
{tab_ret(it)= (*N_surf_t0)(3);break;}
case 123: // la position géométrique Mt
{tab_ret(it)= (*Mt)(1);break;}
case 124: // la position géométrique Mt
{tab_ret(it)= (*Mt)(2);break;}
case 125: // la position géométrique Mt
{tab_ret(it)= (*Mt)(3);break;}
case 126: // la position géométrique M0
{tab_ret(it)= (*M0)(1);break;}
case 127: // la position géométrique M0
{tab_ret(it)= (*M0)(2);break;}
case 128: // la position géométrique M0
{tab_ret(it)= (*M0)(3);break;}
default :
{// on regarde si la grandeur en fait n'est pas à calculer
bool a_atraiter = true;
if (exclure_dd_etend != NULL)
if (find(itt_deb,itt_fin,(*ie)) != itt_fin)
a_atraiter = false;
if (a_atraiter)
{cout << "\n cas de ddl actuellement non traite "
<< "\n pas de ddl " << (*ie).Nom() << " dans l'element "
<< "\n ou cas non implante pour l'instant"
<< "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(...";
tab_ret(it) = 0.;
};
};
} // fin cas **** 2 >>>>>
} // " " "
} // -- fin
else if (( (*ie).Enum() == TEMP) && ((*ie).Nom_vide()))
{tab_ret(it)= def_en_cours->DonneeInterpoleeScalaire(TEMP,temps);
}
else
{ tab_ret(it) = 0.;
cout << "\n cas de ddl actuellement non traite "
<< "\n pas de ddl " << (*ie).Nom() << " dans l'element "
<< "\n ou cas non implante pour l'instant, on retourne 0"
<< "\n Loi_comp_abstraite::Valeur_multi_interpoler_ou_calculer(...";
};
};// -- fin de la boucle sur la liste de Ddl_enum_etendu
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// delete des tenseurs
if (sigHH != NULL); delete sigHH;
if (eps0BB != NULL); delete eps0BB;
if (epsBB != NULL); delete epsBB;
if (epslogBB != NULL); delete epslogBB;
if (epsAlmBB != NULL); delete epsAlmBB;
if (sigHB != NULL); delete sigHB;
if (epsHB != NULL); delete epsHB;
if (DepsBB != NULL); delete DepsBB;
if (DepsHB != NULL); delete DepsHB;
if (DeltaEpsBB != NULL); delete DeltaEpsBB;
if (eps_barreHB != NULL); delete eps_barreHB;
if (Deps_barreHB != NULL); delete Deps_barreHB;
if (sig_barreHB != NULL); delete sig_barreHB;
// cas des pointeurs
if (epsAlmTotalBB!=NULL) delete epsAlmTotalBB; // pour la déformation totale d'almansi
if (epsGLTotalBB!=NULL) delete epsGLTotalBB; // pour la déformation totale de green_lagrange
if (epsLogTotalBB!=NULL) delete epsLogTotalBB; // pour la déformation totale logarithmique
if (Q_sig != NULL) delete Q_sig; // grandeurs polaires
if (Q_eps != NULL) delete Q_eps; // grandeurs polaires
if (Mtdt != NULL) delete Mtdt; // coordonnée du point à t
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (Mt != NULL ) delete Mt; // la position à t
if (M0 != NULL ) delete M0; // coordonnée du point à 0
if (N_surf != NULL) delete N_surf; // vecteur normal à la surface
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (N_surf_t != NULL) delete N_surf_t; // vecteur normal à la surface à t
if (N_surf_t0 != NULL) delete N_surf_t0; // vecteur normal à la surface à t0
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
if (Vitesse != NULL) delete Vitesse; // vitesse
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// pointeurs de matrice
if (Aa0 != NULL) delete Aa0;
if (Aafin != NULL) delete Aafin;
if (gamma0 != NULL) delete gamma0;
if (gammafin != NULL) delete gammafin;
if (beta0 != NULL) delete beta0;
if (betafin != NULL) delete betafin;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// liberation des tenseurs intermediaires
LibereTenseur();
return tab_ret;
};
// affichage de la liste des grandeurs possible à calculer avec Valeur_multi_interpoler_ou_calculer
void Loi_comp_abstraite::Affichage_grandeurs_Valeur_multi_interpoler_ou_calculer()
{ cout << "\n grandeurs scalaires disponibles :";
cout << "\n Mtdt(i) -> X1, X2 ..., Mt(i) -> X1_t, X2_t ..., M0(i) -> X1_t0, X2_t0..,"
<< "\n si pertinent: N_surf(i) -> N_surf_1.., N_surf_t(i) -> N_surf_1_t.., N_surf_t0(i) -> N_surf_1_t0..."
<< "\n SIG11 SIG12 ..., EPS11, EPS12 ..., UX, UY, UZ , V1,V2,V3"
<< "\n def_duale_mises, def_equivalente, def_duale_mises_maxi, vitesse_def_equivalente"
<< "\n contrainte de Mises, contrainte de Tresca"
<< "\n temperature -> TEMP "
;
};
// récupération de valeurs interpolées pour les grandeur enu
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
// une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul
// exclure_Q: donne une liste de grandeur quelconque à exclure de la recherche
// parce que par exemple, ils sont calculés par ailleurs
// on peut également ne pas définir de métrique, dans ce cas on ne peut pas calculer certaines grandeurs
// -> il y a vérification
void Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer
(bool absolue, Enum_dure temps,List_io<TypeQuelconque>& enu
,const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat
,const List_io<EnumTypeQuelconque>* exclure_Q
)
// la dimension des tenseurs de retour est obligatoirement celle de l'espace
// Dans le cas ou on veut une base ad hoc, on calcule dans la base ad hoc et on
// transfert dans le conteneur: ainsi dans ce cas seule les composantes de la base
// ad hoc sont alimentées:
// exemple: base ad hoc : dim 1, espace 3 -> eps11 sera stocké dans le conteneur
// en 11, les autres seront nulles
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ // on ne continue que s'il y a des grandeurs à calculer
if (enu.size() == 0)
return;
// on a besoin a priori de ptintmeca_en_cours, donc s'il n'est pas définit -> erreur
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
#ifdef MISE_AU_POINT
if (ptintmeca_en_cours == NULL)
{cout << "\n *** cas non prevu : aucun conteneur de point d'integration transmis "
<< "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(..." << endl;
Sortie(1);
};
if (def_en_cours == NULL)
{cout << "\n *** cas non prevu : la deformation n'est pas transmise "
<< "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(..." << endl;
Sortie(1);
};
#endif
// ----- def de grandeurs de travail
// def de la dimension des tenseurs
// il y a deux pb a gérer: 1) le fait que la dimension absolue peut-être différente de la dimension actuelle des tenseurs calculés
// 2) le fait que l'on veut une sortie dans une base ad hoc ou pas
// 3) le fait que les conteneurs de retour peuvent ne pas avoir la dimension voulue
int dim = ptintmeca_en_cours->EpsBB_const().Dimension(); // == a dimension actuelle des tenseurs calculés
int dim_sortie_tenseur = dim; // init
// dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue
int dim_espace = ParaGlob::Dimension();
if (absolue)
dim_sortie_tenseur = dim_espace;
// pour ne faire qu'un seul test ensuite
bool prevoir_change_dim_tenseur = false;
// initialement on faisait le test suivant,
// if ((absolue) && (dim != dim_sortie_tenseur))
if (absolue)
// mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire
// car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue
// il faut un tenseur intermédiaire qui a ou non une dimension différente
prevoir_change_dim_tenseur = true;
// si on n'est pas en absolue, donc on est en ad hoc (c'est l'un ou l'autre)
// et si la dim ad hoc est différente de dim_espace, alors il faut aussi
// un tenseur intermédiaire mais de dim ad hoc
bool transfert_tenseur = false;
if ((!absolue)&&(dim != dim_espace))
{dim_sortie_tenseur = dim; // la dimension des tenseurs intermédiaires créés
transfert_tenseur = true;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// éléments de métrique et matrices de passage
TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH;
Mat_pleine jB0(dim,dim),jBfin(dim,dim);
BaseB* giB_0;BaseB* giB_t;
bool pas_de_metrique_dispo = false; // init
if (ex_impli != NULL)
{ gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt;
giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt;
giB_t = ex_impli->giB_t; giB_0 = ex_impli->giB_0;
}
else if (ex_expli_tdt != NULL)
{gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt;
giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt;
giB_t = ex_expli_tdt->giB_t; giB_0 = ex_expli_tdt->giB_0;
}
else if (ex_umat != NULL)
{gijHH = ex_umat->gijHH_t;gijBB = ex_umat->gijBB_t;
giB = giB_t = ex_umat->giB_t; giH_0 = ex_umat->giH_0;giH = ex_umat->giH_t;
giB_0 = ex_umat->giB_0;
}
else
{ pas_de_metrique_dispo = true;
};
// def de tenseurs pour la sortie
// les tenseurs restants en locale
TenseurHH* sigHH=NULL;TenseurBB* epsBB=NULL;TenseurBB* epslogBB=NULL;
TenseurBB* epsAlmBB=NULL;TenseurBB* eps0BB=NULL;TenseurBB* DepsBB = NULL;
TenseurBB* epsAlmTotalBB=NULL;TenseurBB* epsGLTotalBB=NULL;TenseurBB* epsLogTotalBB=NULL;
TenseurBB* DeltaEpsBB = NULL;
Coordonnee* Mtdt=NULL; // coordonnées éventuelles du point d'intégration considéré
Coordonnee* Mt=NULL; // coordonnées à t
Coordonnee* M0=NULL; // coordonnées à t0
Coordonnee* N_surf = NULL; // coordonnée d'un vecteur normal si c'est adéquate
Coordonnee* N_surf_t = NULL; // coordonnée d'un vecteur normal à t si c'est adéquate
Coordonnee* N_surf_t0 = NULL; // coordonnée d'un vecteur normal à t0 si c'est adéquate
Coordonnee* Vitesse = NULL; // cas des vitesses
Coordonnee* Deplacement = NULL; // cas du déplacement
// pour les valeurs propres
Coordonnee* valPropreSig=NULL;Coordonnee* valPropreEps=NULL;
Coordonnee* valPropreDeps=NULL;
// pour les vecteurs propres
Tableau <Coordonnee *> base_propre_sigma , base_propre_eps , base_propre_D;
// pour les bases
Tableau <Coordonnee *> base_ad_hoc , base_giH , base_giB;
// grandeurs scalaires
double* Mises = NULL; double* defDualMises = NULL; double* Tresca = NULL;
double* erreur = NULL;double* erreur_rel = NULL;
double* spherique_eps=NULL,* Q_eps=NULL,* cos3phi_eps=NULL;
double* spherique_sig=NULL,* Q_sig=NULL,* cos3phi_sig=NULL;
double* spherique_Deps=NULL,* Q_Deps=NULL,* cos3phi_Deps=NULL;
double* def_equivalente=NULL, * defDualMisesMaxi=NULL;
double* sig_mises= NULL; // *** peut-être en trop
double* sig_tresca = NULL; // *** peut-être en trop
// --- dev d'un ensemble de variable booléenne pour gérer les sorties en une passe -----
// on se réfère au informations définit dans la méthode: Les_type_evolues_internes()
bool deformationCourante=false; bool vitesseDeformationCourante=false;
bool almansi=false; bool greenLagrange=false;
bool logarithmique=false; bool deltaDef=false;
bool almansiTotal=false; bool greenLagrangeTotale=false;
bool logarithmiqueTotale=false; bool contrainteCourante=false;
bool defPrincipales=false; bool sigmaPrincipales=false;
bool vitPrincipales=false; bool contrainteMises=false;
// bool contraintesTresca=false; bool erreurQ=false;
// bool erreurRel = false;
bool defPlastiqueCumulee=false; bool def_duale_mises=false;
bool besoin_des_contraintes=false; bool besoin_des_deformation=false;
bool besoin_des_contraintes_barre=false; bool besoin_des_deformation_barre=false;
bool besoin_des_vitesses_deformation=false; bool besoin_des_vitesses_deformation_barre=false;
bool besoin_des_valpropre_sigma=false;
bool besoin_des_valpropre_deformation = false; bool besoin_des_valpropre_vitdef = false;
bool besoin_coordonnees = false;
bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false;
bool besoin_dir_princ_sig = false; bool besoin_dir_princ_eps = false;
bool besoin_dir_princ_D = false;
bool besoin_rep_local_ortho=false;
bool besoin_rep_giH = false; bool besoin_rep_giB = false;
bool def_SPHERIQUE_EPS = false; bool def_Q_EPS = false; bool def_COS3PHI_EPS = false;
bool def_SPHERIQUE_SIG = false; bool def_Q_SIG = false; bool def_COS3PHI_SIG = false;
bool def_SPHERIQUE_DEPS = false; bool def_Q_DEPS = false; bool def_COS3PHI_DEPS = false;
bool def_def_equivalente = false; bool def_duale_mises_maxi = false;
// --- dev d'un ensemble de variable booléenne pour gérer les sorties en une passe -----
// on se réfère au informations définit dans la méthode: Les_type_evolues_internes()
bool def_sig_mises= false; // *** peut-être en trop
bool def_sig_tresca= false; // *** peut-être en trop
// init pour l'utilisation de la liste d'exclusion
List_io<EnumTypeQuelconque>::const_iterator itt_deb;
List_io<EnumTypeQuelconque>::const_iterator itt_fin;
if (exclure_Q != NULL)
{ itt_deb = exclure_Q->begin();
itt_fin = exclure_Q->end();
};
// on initialise ces variables booléennes et les conteneurs
List_io<TypeQuelconque>::iterator ipq,ipqfin=enu.end();
for (ipq=enu.begin();ipq!=ipqfin;ipq++)
{EnumTypeQuelconque enuconq = (*ipq).EnuTypeQuelconque().EnumTQ();
switch (enuconq)
{ case CONTRAINTE_COURANTE : {contrainteCourante=true; besoin_des_contraintes=true;
Grandeur_TenseurHH& gr= *((Grandeur_TenseurHH*) ((*ipq).Grandeur_pointee()));
sigHH = gr.ConteneurTenseur();
break;}
case DEFORMATION_COURANTE : {deformationCourante=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
epsBB = gr.ConteneurTenseur(); break;}
case VITESSE_DEFORMATION_COURANTE : {vitesseDeformationCourante=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
DepsBB = gr.ConteneurTenseur(); besoin_des_vitesses_deformation=true; break;}
case ALMANSI : {almansi=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
epsAlmBB = gr.ConteneurTenseur(); break;}
case GREEN_LAGRANGE : {greenLagrange=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
eps0BB = gr.ConteneurTenseur(); break;}
case LOGARITHMIQUE : {logarithmique=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
epslogBB = gr.ConteneurTenseur(); break;}
case DELTA_DEF : {deltaDef=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
DeltaEpsBB = gr.ConteneurTenseur(); break;}
case ALMANSI_TOTAL : {almansiTotal=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
epsAlmTotalBB = gr.ConteneurTenseur(); break;}
case GREEN_LAGRANGE_TOTAL : {greenLagrangeTotale=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
epsGLTotalBB = gr.ConteneurTenseur(); break;}
case LOGARITHMIQUE_TOTALE : {logarithmiqueTotale=true; besoin_des_deformation=true;
Grandeur_TenseurBB& gr= *((Grandeur_TenseurBB*) ((*ipq).Grandeur_pointee()));
epsLogTotalBB = gr.ConteneurTenseur(); break;}
case DEF_PRINCIPALES : {defPrincipales=true; besoin_des_deformation=true;
besoin_des_valpropre_deformation=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
valPropreEps = gr.ConteneurCoordonnee(); break;}
case SIGMA_PRINCIPALES : {sigmaPrincipales=true; besoin_des_contraintes=true;
besoin_des_valpropre_sigma=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
valPropreSig = gr.ConteneurCoordonnee(); break;}
case VIT_PRINCIPALES : {vitPrincipales=true; besoin_des_vitesses_deformation=true;
besoin_des_valpropre_vitdef=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
valPropreDeps = gr.ConteneurCoordonnee(); break;}
case DEF_DUALE_MISES : {def_duale_mises=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
defDualMises = gr.ConteneurDouble();
//besoin_des_valpropre_deformation=true; // non car maintenant on utilise directement ce qui est stocké au pt d'integ
break;}
case DEF_EQUIVALENTE : {def_def_equivalente=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
def_equivalente = gr.ConteneurDouble(); break;}
case DEF_DUALE_MISES_MAXI : {def_duale_mises_maxi=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
defDualMisesMaxi = gr.ConteneurDouble(); break;}
case CONTRAINTE_MISES : {def_sig_mises=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
sig_mises = gr.ConteneurDouble(); break;}
case CONTRAINTE_TRESCA : {def_sig_tresca=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
sig_tresca = gr.ConteneurDouble(); break;}
case POSITION_GEOMETRIQUE : {besoin_coordonnees=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
Mtdt = gr.ConteneurCoordonnee(); break;}
case POSITION_GEOMETRIQUE_t : {besoin_coordonnees_t=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
Mt = gr.ConteneurCoordonnee(); break;}
case POSITION_GEOMETRIQUE_t0 : {besoin_coordonnees_t0=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
M0 = gr.ConteneurCoordonnee(); break;}
case SPHERIQUE_EPS : {def_SPHERIQUE_EPS=true; besoin_des_deformation=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
spherique_eps = gr.ConteneurDouble(); break;}
case Q_EPS : {def_Q_EPS=true; besoin_des_deformation=true;besoin_des_deformation_barre=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
Q_eps = gr.ConteneurDouble(); break;}
case COS3PHI_EPS : {def_COS3PHI_EPS=true; besoin_des_deformation=true;besoin_des_deformation_barre=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
cos3phi_eps = gr.ConteneurDouble(); besoin_des_valpropre_deformation=true; break;}
case SPHERIQUE_SIG : {def_SPHERIQUE_SIG=true; besoin_des_contraintes=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
spherique_sig = gr.ConteneurDouble(); besoin_des_valpropre_sigma=true; break;}
case Q_SIG : {def_Q_SIG=true; besoin_des_contraintes=true;besoin_des_contraintes_barre=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
Q_sig = gr.ConteneurDouble(); break;}
case COS3PHI_SIG : {def_COS3PHI_SIG=true; besoin_des_contraintes=true;besoin_des_contraintes_barre=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
cos3phi_sig = gr.ConteneurDouble(); break;}
case SPHERIQUE_DEPS : {def_SPHERIQUE_DEPS=true; besoin_des_vitesses_deformation=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
spherique_Deps = gr.ConteneurDouble(); break;}
case Q_DEPS : {def_Q_DEPS=true; besoin_des_vitesses_deformation=true;besoin_des_vitesses_deformation_barre=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
Q_Deps = gr.ConteneurDouble(); break;}
case COS3PHI_DEPS : {def_COS3PHI_DEPS=true; besoin_des_vitesses_deformation=true;besoin_des_vitesses_deformation_barre=true;
Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
cos3phi_Deps = gr.ConteneurDouble(); break;}
case REPERE_LOCAL_ORTHO : {besoin_rep_local_ortho=true;
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur REPERE_LOCAL_ORTHO necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
base_ad_hoc.Change_taille(dim_espace);
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
switch (dim_espace) {case 3: base_ad_hoc(3) = &(gr(3));
case 2: base_ad_hoc(2) = &(gr(2));
case 1: base_ad_hoc(1) = &(gr(1));
default:break;
};
break;}
case REPERE_LOCAL_H :
{besoin_rep_giH=true;
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur REPERE_LOCAL_H necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
base_giH.Change_taille(dim);
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
switch (dim) {case 3: base_giH(3) = &(gr(3));
case 2: base_giH(2) = &(gr(2));
case 1: base_giH(1) = &(gr(1));
default:break;
};
break;
}
case REPERE_LOCAL_B :
{besoin_rep_giB=true;
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur REPERE_LOCAL_B necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
base_giB.Change_taille(dim);
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
switch (dim) {case 3: base_giB(3) = &(gr(3));
case 2: base_giB(2) = &(gr(2));
case 1: base_giB(1) = &(gr(1));
default:break;
};
break;
}
case DIRECTIONS_PRINC_SIGMA : {besoin_des_contraintes=true;
besoin_des_valpropre_sigma=true;besoin_dir_princ_sig=true;
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur DIRECTIONS_PRINC_SIGMA necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
base_propre_sigma.Change_taille(dim);
switch (dim) {case 3: base_propre_sigma(3) = &(gr(3));
case 2: base_propre_sigma(2) = &(gr(2));
case 1: base_propre_sigma(1) = &(gr(1));
default:break;
};
break;}
case DIRECTIONS_PRINC_DEF : {besoin_des_deformation=true;
besoin_des_valpropre_deformation=true;besoin_dir_princ_eps=true;
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur DIRECTIONS_PRINC_DEF necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
base_propre_eps.Change_taille(dim);
switch (dim) {case 3: base_propre_eps(3) = &(gr(3));
case 2: base_propre_eps(2) = &(gr(2));
case 1: base_propre_eps(1) = &(gr(1));
default:break;
};
break;}
case DIRECTIONS_PRINC_D : {vitesseDeformationCourante=true;besoin_des_vitesses_deformation=true;
besoin_des_valpropre_deformation=true;besoin_dir_princ_D=true;
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur DIRECTIONS_PRINC_D necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
base_propre_D.Change_taille(dim);
switch (dim) {case 3: base_propre_D(3) = &(gr(3));
case 2: base_propre_D(2) = &(gr(2));
case 1: base_propre_D(1) = &(gr(1));
default:break;
};
break;}
case NN_SURF:{
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
N_surf = gr.ConteneurCoordonnee(); break;}
case NN_SURF_t:{
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
N_surf_t = gr.ConteneurCoordonnee(); break;}
case NN_SURF_t0:{
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
N_surf_t0 = gr.ConteneurCoordonnee(); break;}
case DEPLACEMENT:{besoin_coordonnees=true;
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
Deplacement = gr.ConteneurCoordonnee(); break;}
case VITESSE:{
Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee()));
Vitesse = gr.ConteneurCoordonnee(); break;}
// dans le cas des numéros, traitement direct ici
case NUM_ELEMENT:
{if (ptintmeca_en_cours != NULL)
{*((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= ptintmeca_en_cours->Nb_ele();}
// sinon on le laisse à 0
break;
}
case NUM_MAIL_ELEM:
{ if (ptintmeca_en_cours != NULL)
*((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= ptintmeca_en_cours->Nb_mail();
break;
}
case NUM_PTI:
{ if (ptintmeca_en_cours != NULL)
*((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= ptintmeca_en_cours->Nb_pti();
break;
}
default :
{ // on regarde si la grandeur en fait n'est pas à calculer
bool a_atraiter = true;
if (exclure_Q != NULL)
if (find(itt_deb,itt_fin,enuconq) != itt_fin)
a_atraiter = false;
if (a_atraiter)
// cas où ce n'est pas à exclure, donc cela manque
{// on initialise la grandeur pour éviter d'avoir des valeurs aléatoires
((*ipq).Grandeur_pointee())->InitParDefaut();
if (ParaGlob::NiveauImpression() > 0)
{cout << "\nWarning : attention cas non traite: "
<< (*ipq).EnuTypeQuelconque().NomPlein() << "!\n";
if (ParaGlob::NiveauImpression() > 5)
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(....";
};
}; // sinon on ne fait rien, car les grandeurs exclues sont sensées
// être calculés par ailleurs
}
};
};
// on complète les définitions de tenseurs si besoin
// les tenseurs restants en locale
TenseurHB* sigHB = NULL ;TenseurHB* sig_barreHB = NULL ;
TenseurHB* epsHB = NULL ;TenseurHB* eps_barreHB = NULL ;
TenseurHB* DepsHB = NULL ; TenseurHB* Deps_barreHB = NULL ;
if ((besoin_des_contraintes && !contrainteCourante)
|| (besoin_des_valpropre_sigma && !sigmaPrincipales)
)
{sigHH = NevezTenseurHH(dim_sortie_tenseur) ;
sigHB = (NevezTenseurHB(dim)) ;
sig_barreHB = (NevezTenseurHB(dim)) ;
}
else if (besoin_des_contraintes
|| (besoin_des_valpropre_sigma && !sigmaPrincipales)
)
// là il s'agit uniquement des tenseurs locaux
{sigHB = (NevezTenseurHB(dim)) ;
sig_barreHB = (NevezTenseurHB(dim)) ;
};
if ((besoin_des_deformation && !deformationCourante)
|| (besoin_des_valpropre_deformation && !defPrincipales)
)
{epsBB = NevezTenseurBB(dim_sortie_tenseur) ;
epsHB = (NevezTenseurHB(dim)) ;
eps_barreHB = (NevezTenseurHB(dim)) ;
}
else if (besoin_des_deformation
|| (besoin_des_valpropre_deformation && !defPrincipales)
)
// là il s'agit uniquement des tenseurs locaux
{epsHB = (NevezTenseurHB(dim)) ;
eps_barreHB = (NevezTenseurHB(dim)) ;
};
if ((besoin_des_vitesses_deformation && !vitesseDeformationCourante)
|| (besoin_des_valpropre_vitdef && !vitPrincipales)
)
{DepsBB = NevezTenseurBB(dim_sortie_tenseur);
DepsHB = (NevezTenseurHB(dim));
Deps_barreHB = (NevezTenseurHB(dim)) ;
}
else if(besoin_des_vitesses_deformation
|| (besoin_des_valpropre_vitdef && !vitPrincipales)
)
// là il s'agit uniquement des tenseurs locaux
{DepsHB = (NevezTenseurHB(dim));
Deps_barreHB = (NevezTenseurHB(dim)) ;
};
if (besoin_des_valpropre_deformation && !defPrincipales) valPropreEps = new Coordonnee(dim_sortie_tenseur);
if (besoin_des_valpropre_vitdef && !vitPrincipales) valPropreDeps = new Coordonnee(dim_sortie_tenseur);
if (besoin_des_valpropre_sigma && !sigmaPrincipales) valPropreSig = new Coordonnee(dim_sortie_tenseur);
// on statut sur le fait d'avoir besoin ou non des chgt de base
bool besoin_matrice_chg_base = false;
if ( besoin_des_contraintes
|| besoin_des_deformation
|| besoin_des_vitesses_deformation
|| besoin_rep_local_ortho
)
{besoin_matrice_chg_base = true;};
// on ne change pas le numéro de point d'intégration courant
// on considère ici que c'est déjà le bon, puisque c'est appelé par un élément
// a priori ...
// matrices de passage
int dim_effective = dim; // init
if (absolue) dim_effective = ParaGlob::Dimension();
Mat_pleine* Aa0 = NULL;Mat_pleine* Aafin = NULL;
Mat_pleine* gamma0 = NULL;Mat_pleine* gammafin = NULL;
Mat_pleine* beta0 = NULL;Mat_pleine* betafin = NULL;
if (besoin_matrice_chg_base)
// dans le cas où on n'est pas en absolue => on sort dans un repère ad hoc donc
// il a la dimension locale
// sinon on sort dans le repère globale => il a la dimension globale
{Aa0 = new Mat_pleine(dim_effective,dim_effective);
Aafin = new Mat_pleine(dim_effective,dim_effective);
gamma0 = new Mat_pleine(dim_effective,dim_effective);
gammafin = new Mat_pleine(dim_effective,dim_effective);
beta0 = new Mat_pleine(dim_effective,dim_effective);
betafin = new Mat_pleine(dim_effective,dim_effective);
};
// on considère que les métriques sont directement utilisables
// si elles sont disponibles
if (besoin_matrice_chg_base)
{if (pas_de_metrique_dispo)
{cout << "\n *** erreur : on a besoin de matrices de changement de base et il n'y a pas "
<< " de metrique disponible: impossible de continuer"
<< "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer..." << endl;
Sortie(1);
}
else if (def_en_cours== NULL)
{cout << "\n *** erreur : on a besoin de l'objet Deformation et il n'y a pas "
<< " de pointeur de deformation disponible: impossible de continuer"
<< "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer..." << endl;
Sortie(1);
}
else
// on peut calculer
{Mat_pleine Aat(dim_effective,dim_effective);
// a priori Aat ne sert pas par la suite, mais est nécessaire pour le passage de par
const Met_abstraite::Info_et_metrique_0_t_tdt ex
= def_en_cours->Remont_et_metrique_0_t_tdtSansCalMet(absolue,*Aa0,Aat,*Aafin);
// pour les formules de passage de repère il nous faut :
// Ip_a = beta_a^{.j} g_j et Ip^b = gamma^b_{.j} g^j
// on a: [beta_a^{.j}] = [Aa^j_{.a}]^T
// et [gamma^b_{.j}] = [beta_a^{.j}]^{-1T} = [Aa^j_{.a}]^{-1}
(*gamma0) = (Aa0->Inverse());
(*gammafin) = (Aafin->Inverse());
// on détermine également les matrices beta
(*beta0) = (Aa0->Transpose());
(*betafin) = (Aafin->Transpose());
};
};
// récup des bases si besoin
if (besoin_rep_local_ortho)
{if ((!absolue)&&(dim_espace==3)&&(dim==2))
// cas d'éléments 2D, pour lesquels on veut un repère local ad hoc
// on ramène une base à 3 vecteurs
{ *(base_ad_hoc(1)) = jBfin(1,1) * giB->Coordo(1) + jBfin(2,1) * giB->Coordo(2);
*(base_ad_hoc(2)) = jBfin(1,2) * giB->Coordo(1) + jBfin(2,2) * giB->Coordo(2);
*(base_ad_hoc(3)) = Util::ProdVec_coor(*(base_ad_hoc(1)),*(base_ad_hoc(2)));
}
else if((!absolue)&&(dim_espace>1)&&(dim==1))
// cas d'éléments 1D, dans un espace 2D ou 3D
// on ramène un seul vecteur non nul, les autres ne peuvent être calculé sans info supplémentaire
{ *(base_ad_hoc(1)) = jBfin(1,1) * giB->Coordo(1);
(base_ad_hoc(2))->Zero(); (base_ad_hoc(3))->Zero(); // init à 0 par défaut
}
else // dans tous les autres cas
{ switch (dim_espace) // on se contente de ramener le repère identité
{case 3: (base_ad_hoc(3))->Zero();(*(base_ad_hoc(3)))(3)=1.;
case 2: (base_ad_hoc(2))->Zero();(*(base_ad_hoc(2)))(2)=1.;
case 1: (base_ad_hoc(1))->Zero();(*(base_ad_hoc(1)))(1)=1.;
default:break;
};
};
};
if (besoin_rep_giH)
{switch (dim) {case 3: *(base_giH(3)) = giH->Coordo(3);
case 2: *(base_giH(2)) = giH->Coordo(2);
case 1: *(base_giH(1)) = giH->Coordo(1);
default:break;
};
};
if (besoin_rep_giB)
{switch (dim) {case 3: *(base_giB(3)) = giB->Coordo(3);
case 2: *(base_giB(2)) = giB->Coordo(2);
case 1: *(base_giB(1)) = giB->Coordo(1);
default:break;
};
};
// ----- calcul des grandeurs à sortir
// calcul des tenseurs initiaux
bool plusZero = true; // s'il faut rajouter des termes, on met des 0
if (besoin_des_contraintes)
{(*sigHB) = (ptintmeca_en_cours->SigHH_const()) * (*gijBB);
if (contrainteCourante) // on n'intervient ici que si on veut une sortie des contraintes courantes
{if (absolue) {(ptintmeca_en_cours->SigHH_const()).BaseAbsolue(*sigHH,*giB);}// changement de base finale
else if (transfert_tenseur)
{ TenseurHH* sigHH_inter = NevezTenseurHH(dim) ;
*sigHH_inter = (ptintmeca_en_cours->SigHH_const()); sigHH_inter->ChBase(*gammafin);
sigHH->Affectation_trans_dimension(*sigHH_inter,false);
delete sigHH_inter;
}
else
{ *sigHH = (ptintmeca_en_cours->SigHH_const());sigHH->ChBase(*gammafin);};
};
};
if (besoin_des_deformation)
{// cas de delta_eps
if(DeltaEpsBB != NULL)
{if (absolue)// changement de base finale
{(ptintmeca_en_cours->DeltaEpsBB_const()).BaseAbsolue(*DeltaEpsBB,*giH);}
else if (transfert_tenseur)
{ TenseurBB* DeltaEpsBB_inter = NevezTenseurBB(dim) ;
*DeltaEpsBB_inter = (ptintmeca_en_cours->DeltaEpsBB_const()); DeltaEpsBB_inter->ChBase(*betafin);
DeltaEpsBB->Affectation_trans_dimension(*DeltaEpsBB_inter,false);
delete DeltaEpsBB_inter;
}
else {*DeltaEpsBB = (ptintmeca_en_cours->DeltaEpsBB_const());DeltaEpsBB->ChBase(*betafin);};
};
// cas de la déformation
(*epsHB) = (*gijHH) * ((ptintmeca_en_cours->EpsBB_const()));
if (deformationCourante)// on n'intervient ici que si on veut une sortie des déformations courantes
{if (absolue)// changement de base finale
{(ptintmeca_en_cours->EpsBB_const()).BaseAbsolue(*epsBB,*giH);}
else if (transfert_tenseur)
{ TenseurBB* epsBB_inter = NevezTenseurBB(dim) ;
*epsBB_inter = (ptintmeca_en_cours->EpsBB_const()); epsBB_inter->ChBase(*betafin);
epsBB->Affectation_trans_dimension(*epsBB_inter,false);
delete epsBB_inter;
}
else {*epsBB = (ptintmeca_en_cours->EpsBB_const());epsBB->ChBase(*betafin);};
};
switch (def_en_cours->Type_de_deformation())
{case DEFORMATION_STANDART : // c'est à dire almansi
{if (greenLagrange)
{if (absolue)
{(ptintmeca_en_cours->EpsBB_const()).BaseAbsolue(*eps0BB,*giH_0);}// changement de base finale
else if (transfert_tenseur)
{ TenseurBB* eps0BB_inter = NevezTenseurBB(dim) ;
*eps0BB_inter = (ptintmeca_en_cours->EpsBB_const()); eps0BB_inter->ChBase(*beta0);
eps0BB->Affectation_trans_dimension(*eps0BB_inter,false);
delete eps0BB_inter;
}
else {eps0BB->ChBase(*beta0);};
};
if (almansi) *epsAlmBB = *epsBB;// le changement de base a été fait juste plus haut
if (almansiTotal) // cas avec dilatation et demande de def Almansi totale
{TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsAlmTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsAlmTotal_local_BB);
if (absolue)// changement de base finale
{epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);}
else if (transfert_tenseur)
{ epsAlmTotal_local_BB->ChBase(*betafin);
epsAlmTotalBB->Affectation_trans_dimension(*epsAlmTotal_local_BB,false);
}
else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite
};
if (greenLagrangeTotale) // cas avec dilatation et demande de def Green_Lagrange totale
{TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsGLTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsGLTotal_local_BB);
if (absolue)// changement de base finale
{epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);}
else if (transfert_tenseur)
{ epsGLTotal_local_BB->ChBase(*beta0);
epsGLTotalBB->Affectation_trans_dimension(*epsGLTotal_local_BB,false);
}
else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite
};
// si l'on veut sortir la déformation logarithmique le plus simple est de la calculer
if (logarithmiqueTotale || logarithmique)
{def_en_cours->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE);
if (logarithmique) // cas du calcul de la def logarithmique
{TenseurBB* epslog_local_BB = epslogBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epslog_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epslog_local_BB);
if (absolue)// changement de base finale
{epslog_local_BB->BaseAbsolue(*epslogBB,*giH);}
else if (transfert_tenseur)
{ epslog_local_BB->ChBase(*betafin);
epslogBB->Affectation_trans_dimension(*epslog_local_BB,false);
}
else {epslogBB->ChBase(*betafin);}; // ici epslog_local_BB == epslogBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epslog_local_BB; // car pas utilisé ensuite
};
if (logarithmiqueTotale) // cas avec dilatation et demande de def log totale
{TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsLogTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsLogTotal_local_BB);
if (absolue)// changement de base finale
{epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);}
else if (transfert_tenseur)
{ epsLogTotal_local_BB->ChBase(*betafin);
epsLogTotalBB->Affectation_trans_dimension(*epsLogTotal_local_BB,false);
}
else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite
};
def_en_cours->Change_type_de_deformation(DEFORMATION_STANDART); // on revient au type initial
};
break;
}
case DEFORMATION_LOGARITHMIQUE :
{ if (logarithmique) *epslogBB=*epsBB;
if (logarithmiqueTotale) // cas avec dilatation et demande de def log totale
{TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsLogTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsLogTotal_local_BB);
if (absolue)// changement de base finale
{epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);}
else if (transfert_tenseur)
{ epsLogTotal_local_BB->ChBase(*betafin);
epsLogTotalBB->Affectation_trans_dimension(*epsLogTotal_local_BB,false);
}
else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite
};
// si l'on veut sortir la déformation d'Almansi ou de green-lagrange le plus simple est de les calculer
if (almansi || greenLagrangeTotale || almansiTotal)
{def_en_cours->Change_type_de_deformation(DEFORMATION_STANDART);
if (almansi) // cas de la def d'almansi
{ TenseurBB* eps_local_BB = epsAlmBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
eps_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*eps_local_BB);
if (absolue)// changement de base finale
{eps_local_BB->BaseAbsolue(*epsAlmBB,*giH);}
else if (transfert_tenseur)
{ eps_local_BB->ChBase(*betafin);
epsAlmBB->Affectation_trans_dimension(*eps_local_BB,false);
}
else {epsAlmBB->ChBase(*betafin);};// ici eps_local_BB == epsAlmBB
if(greenLagrange)
{if (absolue)// changement de base finale
{eps_local_BB->BaseAbsolue(*eps0BB,*giH_0);}
else if (transfert_tenseur)
{ eps_local_BB->ChBase(*beta0);
eps0BB->Affectation_trans_dimension(*eps_local_BB,false);
}
else {eps0BB->ChBase(*beta0);}; // ici eps_local_BB == eps0BB
};
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete eps_local_BB; // car pas utilisé ensuite
};
if (almansiTotal) // cas avec dilatation et demande de def Almansi totale
{TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsAlmTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsAlmTotal_local_BB);
if (absolue)// changement de base finale
{epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);}
else if (transfert_tenseur)
{ epsAlmTotal_local_BB->ChBase(*betafin);
epsAlmTotalBB->Affectation_trans_dimension(*epsAlmTotal_local_BB,false);
}
else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite
};
if (greenLagrangeTotale) // cas avec dilatation et demande de def Green_Lagrange totale
{TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut
if ((prevoir_change_dim_tenseur)|| transfert_tenseur)
epsGLTotal_local_BB = NevezTenseurBB(dim);
def_en_cours->Cal_deformation (temps,*epsGLTotal_local_BB);
if (absolue)// changement de base finale
{epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);}
else if (transfert_tenseur)
{ epsGLTotal_local_BB->ChBase(*beta0);
epsGLTotalBB->Affectation_trans_dimension(*epsGLTotal_local_BB,false);
}
else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB
if ((prevoir_change_dim_tenseur)|| transfert_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite
};
def_en_cours->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); // on revient au type initial
};
break;
}
default:
cout << "\n cas de deformation non encore implante en sortie de visualisation "
<< Nom_type_deformation(def_en_cours->Type_de_deformation())
<< " affichage donc errone des valeurs !!!";
};
}; //-- fin de if (besoin_des_deformation)
if (besoin_des_vitesses_deformation)
{ *DepsHB = (*gijHH) * ((ptintmeca_en_cours->DepsBB_const()));
if (absolue)// changement de base finale
{(ptintmeca_en_cours->DepsBB_const()).BaseAbsolue(*DepsBB,*giH);}
else if (transfert_tenseur)
{ TenseurBB* DepsBB_inter = NevezTenseurBB(dim) ;
*DepsBB_inter = (ptintmeca_en_cours->DepsBB_const()); DepsBB_inter->ChBase(*betafin);
DepsBB->Affectation_trans_dimension(*DepsBB_inter,false);
delete DepsBB_inter;
}
else {*DepsBB = (ptintmeca_en_cours->DepsBB_const());DepsBB->ChBase(*betafin);};
};
if (besoin_des_contraintes_barre)
{double Isig = sigHB->Trace(); // trace de la déformation
*sig_barreHB = (*sigHB) - (Isig/dim_espace) * (*Id_dim_HB(dim));
};
if (besoin_des_deformation_barre)
{double Ieps = epsHB->Trace(); // trace de la déformation
*eps_barreHB = (*epsHB) - (Ieps/dim_espace) * (*Id_dim_HB(dim));
};
if (besoin_des_vitesses_deformation_barre)
{double IDeps = DepsHB->Trace(); // trace de la déformation
*Deps_barreHB = (*DepsHB) - (IDeps/dim_espace) * (*Id_dim_HB(dim));
};
// cas des valeurs propres et éventuellement des vecteurs propres
{int caas=0;
////----- debug
//cout << "\n besoin_des_valpropre_sigma= "<< besoin_des_valpropre_sigma
// << " besoin_dir_princ_sig= "<< besoin_dir_princ_sig << endl;
////--- fin debug
if (besoin_des_valpropre_sigma && besoin_dir_princ_sig)// on veut les directions et les valeurs propres
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(sigHB->ValPropre(caas,mat));
inter.Change_dim(dim_sortie_tenseur); // on étant la taille éventuellement
*valPropreSig = inter; // sauvegarde des valeurs propres
// puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_sigma(3)),tab(3));(base_propre_sigma(3))->Normer();
giB->BaseAbsolue( *(base_propre_sigma(2)),tab(2));(base_propre_sigma(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
break;
case 2: // en 2D le premier vecteur en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
giH->BaseAbsolue( *(base_propre_sigma(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_sigma(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_dir_princ_sig) // cas où on ne veut que les directions principales
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(sigHB->ValPropre(caas,mat));
// sauvegarde des directions principales
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_sigma(3)),tab(3));(base_propre_sigma(3))->Normer();
giB->BaseAbsolue( *(base_propre_sigma(2)),tab(2));(base_propre_sigma(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_sigma(1)),tab(1));(base_propre_sigma(1))->Normer();
giH->BaseAbsolue( *(base_propre_sigma(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_sigma(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_des_valpropre_sigma)// on ne veut que les valeurs propres
{ Coordonnee inter(sigHB->ValPropre(caas));
inter.Change_dim(dim_sortie_tenseur); // on étend la taille éventuellement
*valPropreSig = inter; // sauvegarde des valeurs propres
};
// gestion d'une erreur éventuelle
if (caas == -1)
{ cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la contrainte ";
if (ParaGlob::NiveauImpression() > 5) {sigHB->Ecriture(cout);cout << "\n cas = "<< caas ;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(...";};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << endl;
};
}; // -- fin de l'encapsulation du cas des contraintes
{int caas=0;
if (besoin_des_valpropre_deformation && besoin_dir_princ_eps)// on veut les directions et les valeurs propres
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
// 1) calcul des valeurs propres dans un vecteur de travail inter, et des vecteurs propres stockés
Coordonnee inter(epsHB->ValPropre(caas,mat)); // dans la matrice mat
inter.Change_dim(dim_sortie_tenseur); // on étend la taille, dans le cas ou le nb de valeurs
// propre est inférieur à la dimension de l'espace
*valPropreEps = inter; // sauvegarde des valeurs propres
// 2) puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_eps(3)),tab(3));(base_propre_eps(3))->Normer();
giB->BaseAbsolue( *(base_propre_eps(2)),tab(2));(base_propre_eps(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
giH->BaseAbsolue( *(base_propre_eps(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_eps(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_dir_princ_eps) // cas où on ne veut que les directions principales
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(epsHB->ValPropre(caas,mat)); // dans la matrice mat
// puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_eps(3)),tab(3));(base_propre_eps(3))->Normer();
giB->BaseAbsolue( *(base_propre_eps(2)),tab(2));(base_propre_eps(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_eps(1)),tab(1));(base_propre_eps(1))->Normer();
giH->BaseAbsolue( *(base_propre_eps(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_eps(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_des_valpropre_deformation)// on ne veut que les valeurs propres
{ Coordonnee inter(epsHB->ValPropre(caas));
inter.Change_dim(dim_sortie_tenseur); // on étant la taille éventuellement
*valPropreEps = inter; // sauvegarde des valeurs propres
};
// gestion d'une erreur éventuelle
if (caas == -1)
{ cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la deformation";
if (ParaGlob::NiveauImpression() >= 7) {epsHB->Ecriture(cout);cout << "\n cas = "<< caas ;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(...";};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << endl;
};
}; // -- fin de l'encapsulation du cas des déformations
{int caas=0; // il faut l'initialiser sinon il peut prendre la valeur du cas précédant
if (besoin_des_valpropre_vitdef && besoin_dir_princ_D)// on veut les directions et les valeurs propres
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
// 1) calcul des valeurs propres dans un vecteur de travail inter, et des vecteurs propres stockés
Coordonnee inter(DepsHB->ValPropre(caas,mat)); // dans la matrice mat
inter.Change_dim(dim_sortie_tenseur); // on étend la taille, dans le cas ou le nb de valeurs
// propre est inférieur à la dimension de l'espace
*valPropreDeps = inter; // sauvegarde des valeurs propres
// 2) puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_D(3)),tab(3));(base_propre_D(3))->Normer();
giB->BaseAbsolue( *(base_propre_D(2)),tab(2));(base_propre_D(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
giH->BaseAbsolue( *(base_propre_D(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_D(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_dir_princ_D) // cas où on ne veut que les directions principales
{ Mat_pleine mat(dim,dim); // contiendra par colonne les vecteurs principaux, exprimé en giH
Coordonnee inter(DepsHB->ValPropre(caas,mat)); // dans la matrice mat
// puis des directions principales, que l'on doit exprimer dans le repère absolue
Tableau <CoordonneeH > tab = mat.CoordonneeH_Base_associee(); // récup en tab de coor
// on regarde le cas particulier ou les valeurs propres sont toutes nulles
// si c'est le cas, on ne fait rien, en sortie on aura des vecteurs propres nulles
if (inter.Max_val_abs() > ConstMath::petit)
// et s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
{if (caas != -1)
{switch (dim)
{case 3: giB->BaseAbsolue( *(base_propre_D(3)),tab(3));(base_propre_D(3))->Normer();
giB->BaseAbsolue( *(base_propre_D(2)),tab(2));(base_propre_D(2))->Normer();
case 1: giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
break;
case 2: // en 2D le premier vecteur en en giB et le second en giH !!
giB->BaseAbsolue( *(base_propre_D(1)),tab(1));(base_propre_D(1))->Normer();
giH->BaseAbsolue( *(base_propre_D(2)),mat.CoordonneeB_Base_associee(2));
(base_propre_D(2))->Normer();
break;
default:break;
};
};
};
}
else if (besoin_des_valpropre_vitdef)// on ne veut que les valeurs propres
{ Coordonnee inter(DepsHB->ValPropre(caas));
inter.Change_dim(dim_sortie_tenseur); // on étant la taille
*valPropreDeps = inter; // sauvegarde des valeurs propres
};
// gestion d'une erreur éventuelle
if (caas == -1)
{ cout << "\n warning *** erreur dans le calcul des valeurs propres et-ou vecteurs propres de la vitesse de deformation";
if (ParaGlob::NiveauImpression() >= 7) {DepsHB->Ecriture(cout);cout << "\n cas = "<< caas ;
mise à jour 7.010
2023-05-03 17:23:49 +02:00
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(...";};
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << endl;
};
}; // -- fin de l'encapsulation du cas des vitesses de déformations
if (besoin_coordonnees)
(*Mtdt) = def_en_cours->Position_tdt();
if (besoin_coordonnees_t)
(*Mt) = def_en_cours->Position_t();
if (besoin_coordonnees_t0)
(*M0) = def_en_cours->Position_0();
// def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D
if (N_surf != NULL)
{ if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur N_surf necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// on vérifie que la métrique est correcte
if (giB->NbVecteur() != 2)
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n *** attention il ne s'agit pas d'un element 2D,"
<< " le vecteur normal ne sera pas disponible";
if (ParaGlob::NiveauImpression() > 2)
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(... ";
};
cout << endl;
}
else // sinon c'est ok
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
(*N_surf) = Util::ProdVec_coorBN( (*giB)(1), (*giB)(2));
N_surf->Normer(); // que l'on norme
};
// on n'arrête pas l'exécution, car on pourrait vouloir sortir les normales pour un ensemble
// d'éléments contenant des volumes, des surfaces, des lignes: bon... il y aura quand même des
// pb au niveau des iso par exemple, du au fait que l'on va faire des moyennes sur des éléments
// de type différents (à moins de grouper par type du coup on n'aura pas le warning
};
// idem à l'instant t
if (N_surf_t != NULL)
{ if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur N_surf_t necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// on vérifie que la métrique est correcte
if (giB_t->NbVecteur() != 2)
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
<< " le vecteur normal ne sera pas correctement calcule";
if (ParaGlob::NiveauImpression() > 2)
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(... ";
};
cout << endl;
}
else // sinon c'est ok
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
(*N_surf_t) = Util::ProdVec_coorBN( (*giB_t)(1), (*giB_t)(2));
N_surf_t->Normer(); // que l'on norme
};
};
// idem à l'instant t0
if (N_surf_t0 != NULL)
{ if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur N_surf_0 necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// on vérifie que la métrique est correcte
if (giB_0->NbVecteur() != 2)
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n *** attention il ne s'agit pas d'une metrique 2D 2D,"
<< " le vecteur normal ne sera pas correctement calcule";
if (ParaGlob::NiveauImpression() > 2)
cout << "\n Loi_comp_abstraite::Valeurs_Tensorielles_interpoler_ou_calculer(... ";
};
cout << endl;
}
else // sinon c'est ok
{// la normale vaut le produit vectoriel des 2 premiers vecteurs
(*N_surf_t0) = Util::ProdVec_coorBN( (*giB_0)(1), (*giB_0)(2));
N_surf_t0->Normer(); // que l'on norme
};
};
// cas du déplacement et de la vitesse
if (Deplacement != NULL)
(*Deplacement) = (*Mtdt) - def_en_cours->Position_0();
if (Vitesse != NULL)
(*Vitesse) = def_en_cours->VitesseM_tdt();
// cas de la deformation équivalente cumulée
if (def_def_equivalente) {*def_equivalente = ptintmeca_en_cours->Deformation_equi_const()(1);};
// cas de la deformation duale au sens de mises,
if (def_duale_mises) {*defDualMises = ptintmeca_en_cours->Deformation_equi_const()(2);};
// cas de la deformation maxi duale au sens de mises,
if (def_duale_mises_maxi) {*defDualMisesMaxi = ptintmeca_en_cours->Deformation_equi_const()(3);};
// cas des contraintes: mises et tresca
if (def_sig_mises) {*sig_mises = ptintmeca_en_cours->Sig_equi_const()(1);};
if (def_sig_tresca) {*sig_tresca = ptintmeca_en_cours->Sig_equi_const()(2);};
// //contrainte_tresca
// if (contraintesTresca)
// { switch (dim) {case 1: *Tresca=0.5 * (*valPropreSig)(1);break;
// case 2: *Tresca=0.5 * ((*valPropreSig)(1)-(*valPropreSig)(2));break;
// case 3: *Tresca=0.5 * ((*valPropreSig)(1)-(*valPropreSig)(3));break;
// };
// };
// --- cas des grandeurs de la décomposition polaire
// cas de la déformation
if (def_SPHERIQUE_EPS) {*spherique_eps = epsHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012
double mini_Q = 5.e-5;
if (def_Q_EPS) {*Q_eps = sqrt(eps_barreHB->II());};
if (cos3phi_eps)
{ double Qepsilon = ( def_Q_EPS ? *Q_eps : sqrt(eps_barreHB->II()));
double Qepsilon3 = Qepsilon * Qepsilon * Qepsilon;
if (Qepsilon > mini_Q )
{ // on peut calculer un cos3phi pas débile
double bIIIb = eps_barreHB->III() / 3.;
*cos3phi_eps = 3. * sqrt(6.) * bIIIb/ Qepsilon3;
}
else *cos3phi_eps=0.; // sinon on le met à 0
};
// cas de la contrainte
if (def_SPHERIQUE_SIG) {*spherique_sig = sigHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012
if (def_Q_SIG) {*Q_sig = sqrt(sig_barreHB->II());};
if (cos3phi_sig)
{ double Qsig = ( def_Q_SIG ? *Q_sig : sqrt(sig_barreHB->II()));
double Qsig3 = Qsig * Qsig * Qsig;
if (Qsig > mini_Q )
{ // on peut calculer un cos3phi pas débile
double bIIIb = sig_barreHB->III() / 3.;
*cos3phi_sig = 3. * sqrt(6.) * bIIIb/ Qsig3;
}
else *cos3phi_sig=0.; // sinon on le met à 0
};
// cas de la vitesse de déformation
if (def_SPHERIQUE_DEPS) {*spherique_Deps = DepsHB->Trace()/3.;}; //ParaGlob::Dimension();}; modif 5/2/2012
if (def_Q_DEPS) {*Q_Deps = sqrt(Deps_barreHB->II());};
if (cos3phi_Deps)
{ double QDepsilon = ( def_Q_DEPS ? *Q_Deps : sqrt(Deps_barreHB->II()));
double QDepsilon3 = QDepsilon * QDepsilon * QDepsilon;
if (QDepsilon > mini_Q )
{ // on peut calculer un cos3phi pas débile
double bIIIb = Deps_barreHB->III() / 3.;
*cos3phi_Deps = 3. * sqrt(6.) * bIIIb/ QDepsilon3;
}
else *cos3phi_Deps=0.; // sinon on le met à 0
};
// delete des tenseurs
if (sigHB != NULL) delete sigHB;
if (epsHB != NULL) delete epsHB;
if (DepsHB != NULL) delete DepsHB;
if (eps_barreHB != NULL) delete eps_barreHB;
if (Deps_barreHB != NULL) delete Deps_barreHB;
if (sig_barreHB != NULL) delete sig_barreHB;
// effacement conditionnel
if (besoin_des_contraintes && !contrainteCourante) delete sigHH;
if (besoin_des_deformation && !deformationCourante) delete epsBB;
if (vitesseDeformationCourante && !vitesseDeformationCourante) delete DepsBB;
if (besoin_des_valpropre_deformation && !defPrincipales) delete valPropreEps;
if (besoin_des_valpropre_vitdef && !vitPrincipales) delete valPropreDeps ;
if (besoin_des_valpropre_sigma && !sigmaPrincipales) delete valPropreSig;
// pointeurs de matrice
if (Aa0 != NULL) delete Aa0;
if (Aafin != NULL) delete Aafin;
if (gamma0 != NULL) delete gamma0;
if (gammafin != NULL) delete gammafin;
if (beta0 != NULL) delete beta0;
if (betafin != NULL) delete betafin;
// liberation des tenseurs intermediaires
LibereTenseur();
};
// affichage de la liste des grandeurs possible à calculer avec Valeur_multi_interpoler_ou_calculer
void Loi_comp_abstraite::Affichage_grandeurs_Valeurs_Tensorielles_interpoler_ou_calculer()
{ cout << "\n grandeurs tensorielles ou evoluees disponibles :";
cout << "\n CONTRAINTE_COURANTE, DEFORMATION_COURANTE, VITESSE_DEFORMATION_COURANTE "
<< "\n ALMANSI, GREEN_LAGRANGE, LOGARITHMIQUE, DELTA_DEF, ALMANSI_TOTAL"
<< "\n GREEN_LAGRANGE_TOTAL, LOGARITHMIQUE_TOTALE, DEF_PRINCIPALES"
<< "\n SIGMA_PRINCIPALES, VIT_PRINCIPALES, SPHERIQUE_EPS, Q_EPS, COS3PHI_EPS"
<< "\n SPHERIQUE_SIG, Q_SIG, COS3PHI_SIG, SPHERIQUE_DEPS, Q_DEPS"
<< "\n COS3PHI_DEPS, DIRECTIONS_PRINC_SIGMA, DIRECTIONS_PRINC_DEF"
<< "\n DIRECTIONS_PRINC_D, NUM_ELEMENT, NUM_MAIL_ELEM, NUM_PTI"
<< "\n POSITION_GEOMETRIQUE, POSITION_GEOMETRIQUE_t, POSITION_GEOMETRIQUE_t0, "
<< "\n si pertinent: REPERE_LOCAL_ORTHO, REPERE_LOCAL_H, REPERE_LOCAL_B"
<< "\n si pertinent: NN_SURF, NN_SURF_t, NN_SURF_t0, DEPLACEMENT, VITESSE"
<< "\n DEF_DUALE_MISES, DEF_EQUIVALENTE, DEF_DUALE_MISES_MAXI "
<< "\n ,CONTRAINTE_MISES,CONTRAINTE_TRESCA"
;
};
// récupération de valeurs interpolées pour les grandeur ici considéré quelconque enu
// ces grandeurs ne sont pas définies dans la liste des Ddl_enum_etendu : ex mises à t
// ou le numéro de l'élément etc.
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
// une seule des 3 métriques doit-être renseigné, les autres doivent être un pointeur nul
// exclure_Q: donne une liste de grandeur quelconque à exclure de la recherche
// parce que par exemple, ils sont calculés par ailleurs
// on peut également ne pas définir de métrique, dans ce cas on ne peut pas calculer certaines grandeurs
// -> il y a vérification
// retour: la list li_quelc qui en entrée doit avoir la même dimension que tqi
void Loi_comp_abstraite::Valeurs_quelconque_interpoler_ou_calculer
(bool absolue, Enum_dure temps
,const Tableau <EnumTypeQuelconque>& tqi
,List_io<TypeQuelconque>& li_quelc
,const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat
,const List_io<EnumTypeQuelconque>* exclure_Q
)
{ // on a besoin a priori de ptintmeca_en_cours, donc s'il n'est pas définit -> erreur
#ifdef MISE_AU_POINT
if (ptintmeca_en_cours == NULL)
{cout << "\n *** cas non prevu : aucun conteneur de point d'integration transmis "
<< "\n Loi_comp_abstraite::Valeurs_quelconque_interpoler_ou_calculer(..." << endl;
Sortie(1);
};
if (def_en_cours == NULL)
{cout << "\n *** cas non prevu : la deformation n'est pas transmise "
<< "\n Loi_comp_abstraite::Valeurs_quelconque_interpoler_ou_calculer(..." << endl;
Sortie(1);
};
#endif
// ----- def de grandeurs de travail
// def de la dimension des tenseurs
// il y a deux pb a gérer: 1) le fait que la dimension absolue peut-être différente de la dimension des tenseurs
// 2) le fait que l'on veut une sortie dans une base ad hoc ou pas
int dim = ptintmeca_en_cours->EpsBB_const().Dimension();
int dim_sortie_tenseur = dim;
// dans le cas ou l'on veut une sortie en base absolue, il faut que dim_sortie_tenseur = la dimension de la base absolue
int dim_espace = ParaGlob::Dimension();
if (absolue)
dim_sortie_tenseur = dim_espace;
// pour ne faire qu'un seul test ensuite
bool prevoir_change_dim_tenseur = false;
// initialement on faisait le test suivant,
// if ((absolue) && (dim != dim_sortie_tenseur))
if (absolue)
// mais en fait, quand on sort en absolu on est obligé d'utiliser un tenseur intermédiaire
// car BaseAbsolue(.. modifie tenseur passé en paramètre, donc dans tous les cas de sortie absolue
// il faut un tenseur intermédiaire qui a ou non une dimension différente
prevoir_change_dim_tenseur = true;
// éléments de métrique et matrices de passage
TenseurHH* gijHH;TenseurBB* gijBB;BaseB* giB; BaseH* giH_0;BaseH* giH;
Mat_pleine jB0(dim,dim),jBfin(dim,dim);
BaseB* giB_0;BaseB* giB_t;
bool pas_de_metrique_dispo = false; // init
if (ex_impli != NULL)
{ gijHH = ex_impli->gijHH_tdt;gijBB = ex_impli->gijBB_tdt;
giB = ex_impli->giB_tdt; giH_0 = ex_impli->giH_0;giH = ex_impli->giH_tdt;
giB_t = ex_impli->giB_t; giB_0 = ex_impli->giB_0;
}
else if (ex_expli_tdt != NULL)
{gijHH = ex_expli_tdt->gijHH_tdt;gijBB = ex_expli_tdt->gijBB_tdt;
giB = ex_expli_tdt->giB_tdt; giH_0 = ex_expli_tdt->giH_0;giH = ex_expli_tdt->giH_tdt;
giB_t = ex_expli_tdt->giB_t; giB_0 = ex_expli_tdt->giB_0;
}
else if (ex_umat != NULL)
{gijHH = ex_umat->gijHH_t;gijBB = ex_umat->gijBB_t;
giB = giB_t = ex_umat->giB_t; giH_0 = ex_umat->giH_0;giH = ex_umat->giH_t;
giB_0 = ex_umat->giB_0;
}
else
{ pas_de_metrique_dispo = true;
// cout << "\n *** cas non prevu : aucune metrique transmise "
// << "\n Loi_comp_abstraite::Valeurs_quelconque_interpoler_ou_calculer(..." << endl;
// Sortie(1);
};
// def de tenseurs pour la sortie
// les tenseurs restants en locale
// pour les valeurs propres
// pour les vecteurs propres
// pour les bases
Tableau <Coordonnee *> base_ad_hoc , base_giH , base_giB;
// --- dev d'un ensemble de variable booléenne pour gérer les sorties en une passe -----
// on se réfère au informations définit dans la méthode: Les_type_evolues_internes()
bool besoin_rep_local_ortho=false;
bool besoin_rep_giH = false; bool besoin_rep_giB = false;
// init pour l'utilisation de la liste d'exclusion
List_io<EnumTypeQuelconque>::const_iterator itt_deb;
List_io<EnumTypeQuelconque>::const_iterator itt_fin;
if (exclure_Q != NULL)
{ itt_deb = exclure_Q->begin();
itt_fin = exclure_Q->end();
};
// on va boucler sur les énumérés quelconque
int nb_tqi = tqi.Taille();
List_io <TypeQuelconque >::iterator ipq=li_quelc.begin();
for (int i_tqi = 1;i_tqi <= nb_tqi;i_tqi++,ipq++ )
// on initialise ces variables booléennes et les conteneurs
{EnumTypeQuelconque& enuconq = tqi(i_tqi);
switch (enuconq)
{
case CONTRAINTE_MISES :
{Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurDouble()) = ptintmeca_en_cours->Sig_equi_const()(1);
break;
}
case CONTRAINTE_TRESCA :
{Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurDouble()) = ptintmeca_en_cours->Sig_equi_const()(2);
break;
}
case CONTRAINTE_MISES_T :
{Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurDouble()) = ptintmeca_en_cours->Sig_equi_t_const()(1);
break;
}
case CONTRAINTE_TRESCA_T :
{Grandeur_scalaire_double& gr= *((Grandeur_scalaire_double*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurDouble()) = ptintmeca_en_cours->Sig_equi_t_const()(2);
break;
}
case NUM_ELEMENT:
{Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurEntier())=ptintmeca_en_cours->Nb_ele();
break;
}
case NUM_MAIL_ELEM:
{Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurEntier())=ptintmeca_en_cours->Nb_mail();
break;
}
case NUM_PTI:
{Grandeur_scalaire_entier& gr= *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()));
(*gr.ConteneurEntier())=ptintmeca_en_cours->Nb_pti();
break;
}
case REPERE_LOCAL_ORTHO :
{base_ad_hoc.Change_taille(dim_espace);
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur REPERE_LOCAL_ORTHO necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
// calcul de la base
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
switch (dim_espace) {case 3: base_ad_hoc(3) = &(gr(3));
case 2: base_ad_hoc(2) = &(gr(2));
case 1: base_ad_hoc(1) = &(gr(1));
default:break;
};
if ((!absolue)&&(dim_espace==3)&&(dim==2))
// cas d'éléments 2D, pour lesquels on veut un repère local ad hoc
// on ramène une base à 3 vecteurs
{ *(base_ad_hoc(1)) = jBfin(1,1) * giB->Coordo(1) + jBfin(2,1) * giB->Coordo(2);
*(base_ad_hoc(2)) = jBfin(1,2) * giB->Coordo(1) + jBfin(2,2) * giB->Coordo(2);
*(base_ad_hoc(3)) = Util::ProdVec_coor(*(base_ad_hoc(1)),*(base_ad_hoc(2)));
}
else if((!absolue)&&(dim_espace>1)&&(dim==1))
// cas d'éléments 1D, dans un espace 2D ou 3D
// on ramène un seul vecteur non nul, les autres ne peuvent être calculé sans info supplémentaire
{ *(base_ad_hoc(1)) = jBfin(1,1) * giB->Coordo(1);
(base_ad_hoc(2))->Zero(); (base_ad_hoc(3))->Zero(); // init à 0 par défaut
}
else // dans tous les autres cas
{ switch (dim_espace) // on se contente de ramener le repère identité
{case 3: (base_ad_hoc(3))->Zero();(*(base_ad_hoc(3)))(3)=1.;
case 2: (base_ad_hoc(2))->Zero();(*(base_ad_hoc(2)))(2)=1.;
case 1: (base_ad_hoc(1))->Zero();(*(base_ad_hoc(1)))(1)=1.;
default:break;
};
};
break;
}
case REPERE_LOCAL_H :
{base_giH.Change_taille(dim);
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur REPERE_LOCAL_H necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
switch (dim) {case 3: base_giH(3) = &(gr(3));
case 2: base_giH(2) = &(gr(2));
case 1: base_giH(1) = &(gr(1));
default:break;
};
// calcul de la base
switch (dim) {case 3: *(base_giH(3)) = giH->Coordo(3);
case 2: *(base_giH(2)) = giH->Coordo(2);
case 1: *(base_giH(1)) = giH->Coordo(1);
default:break;
};
// stockage
TypeQuelconque typQ6(REPERE_LOCAL_H,SIG11,gr);
li_quelc.push_back(typQ6);
break;
}
case REPERE_LOCAL_B :
{base_giB.Change_taille(dim);
Tab_Grandeur_Coordonnee& gr= *((Tab_Grandeur_Coordonnee*) ((*ipq).Grandeur_pointee()));
if (pas_de_metrique_dispo)
{ cout << "\n ***erreur, la grandeur REPERE_LOCAL_B necessite la donnee de la metrique "
<< " or elle n'est pas ici disponible !! ";
Sortie(1);
};
switch (dim_espace) {case 3: gr(3).Zero(); case 2: gr(2).Zero(); default:break;};//init
switch (dim) {case 3: base_giB(3) = &(gr(3));
case 2: base_giB(2) = &(gr(2));
case 1: base_giB(1) = &(gr(1));
default:break;
};
// calcul de la base
switch (dim) {case 3: *(base_giB(3)) = giB->Coordo(3);
case 2: *(base_giB(2)) = giB->Coordo(2);
case 1: *(base_giB(1)) = giB->Coordo(1);
default:break;
};
// stockage
TypeQuelconque typQ6(REPERE_LOCAL_B,SIG11,gr);
li_quelc.push_back(typQ6);
break;
}
default :
{ // on regarde si la grandeur en fait n'est pas à calculer
bool a_atraiter = true;
if (exclure_Q != NULL)
if (find(itt_deb,itt_fin,enuconq) != itt_fin)
a_atraiter = false;
if (a_atraiter)
// cas où ce n'est pas à exclure, donc cela manque
{if (ParaGlob::NiveauImpression() > 0)
{cout << "\n warning : attention cas non traite: "
<< NomTypeQuelconque(enuconq) << "!\n";
if (ParaGlob::NiveauImpression() > 5)
cout << "\n Loi_comp_abstraite::Valeurs_quelconque_interpoler_ou_calculer(....";
Sortie(1);
};
};
}
};
};
// delete des tenseurs
// liberation des tenseurs intermediaires
LibereTenseur();
};
// affichage de la liste des grandeurs possible à calculer avec Valeur_multi_interpoler_ou_calculer
void Loi_comp_abstraite::Affichage_grandeurs_Valeurs_quelconque_interpoler_ou_calculer()
{ cout << "\n grandeurs quelconques disponibles :";
cout << "\n si pertinent: REPERE_LOCAL_ORTHO, REPERE_LOCAL_H, REPERE_LOCAL_B"
<< "\n ,CONTRAINTE_MISES_T,CONTRAINTE_TRESCA_T"
;
};
// calcul de la valeur et retour dans tab_ret d'une fonction nD
// à l'aide des grandeurs disponibles pour la loi de comportement
// nb_retour: nombre de composantes attendues en retour, si > 0
// deja_calculer: donne une liste de grandeur Ddl_enum_etendu et quelconque
// à exclure de la recherche car ils doivent avoir été calculés par ailleurs
// c'est une donnée d'entrée qui peut-être utilisée pour l'appel de la fonction nD
// list_save : est censé contenir la ou les save_result à consulter pour avoir
// des infos supplémentaires
Tableau <double> & Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Fonction_nD* fct,int nb_retour
,const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat
,const List_io<Ddl_etendu>* deja_calculer_etend
,const List_io<const TypeQuelconque *>* deja_calculer_Q
,list <SaveResul*>* list_save
)
{ // ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
// 1) la partie grandeurs évoluées
List_io <Ddl_enum_etendu>& li_enu_scal = fct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct->Li_equi_Quel_evolue();
// on récupère le type d'expression des tenseurs
bool absolue = fct->Absolue();
// --- on s'occupe des Ddl_enum_etendu nécessaires pour l'appel de la fonction nD
List_io<Ddl_enum_etendu> exclure_dd_etend;
// --> on rempli la liste à exclure
if (deja_calculer_etend != NULL)
{ List_io<Ddl_etendu>::const_iterator ik,ikfin=deja_calculer_etend->end();
for (ik= deja_calculer_etend->begin(); ik != ikfin;ik++)
exclure_dd_etend.push_back((*ik).Const_DdlEnumEtendu());
};
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_umat,&exclure_dd_etend));
// on va récupérer les infos qui ont été exclus
int index = 1;
if (deja_calculer_etend != NULL)
{ List_io<Ddl_etendu>::const_iterator ik,ikfin=deja_calculer_etend->end();
List_io <Ddl_enum_etendu>::iterator il,ilfin=li_enu_scal.end();
for (il = li_enu_scal.begin(); il != ilfin; il++,index++)
{// on parcours la liste: normalement il n'y a pas beaucoup d'éléments ...
// donc ce n'est "peut-être" pas trop couteux
for (ik= deja_calculer_etend->begin(); ik != ikfin;ik++)
{if ((*ik).Const_DdlEnumEtendu() == (*il))
// on a trouvé
{val_ddl_enum(index) = (*ik).ConstValeur();
break;
};
};
};
};
// pour les grandeurs quelconque, on va exclure certaines grandeurs de l'appel de
// Valeur_multi_interpoler_ou_calculer(..)
// pour cela on se sert d'une liste locale que l'on abonde ensuite avec la liste passée
// en paramètre si elle existe
List_io<EnumTypeQuelconque> exclure_local;
if (list_save != NULL)
{// on essaie de récupérer les types quelconques
// récupération des type quelconque sous forme d'un arbre pour faciliter la recherche
list <SaveResul*>::iterator isave=list_save->begin(); // pour les saveResul des lois
list <SaveResul*>::iterator isavedebut=list_save->begin(); // pour les saveResul des lois
list <SaveResul*>::iterator isavefin=list_save->end(); // pour les saveResul des lois
List_io <TypeQuelconque>::iterator it,itfin = li_quelc.end();
for (it = li_quelc.begin(); it != itfin; it++)
{bool valeur_recuperer = false; // pour savoir s'il y a eu récupération
EnumTypeQuelconque enu=(*it).EnuTypeQuelconque().EnumTQ();
// on boucle sur les save_results pour récupérer les types quelconques
for (isave = isavedebut; isave != isavefin;isave++)
// l'idée ici est de regarder si on ne trouve pas une des grandeurs demandées dans la map
{const map < EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >*
map_quelcon = (*isave)->Map_type_quelconque();
map < EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >::const_iterator imap;
if (map_quelcon != NULL) // sinon il n'y a pas de map
{imap = map_quelcon->find(enu);
if (imap != map_quelcon->end())
{ // on a trouvé la grandeur
switch (Type_de_grandeur_associee(enu))
{ case SCALAIRE_DOUBLE:
{ const Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*)
(*imap).second.Const_Grandeur_pointee());
// on renseigne l'argument
(*(*it).Grandeur_pointee()) = tyTQ;
valeur_recuperer = true;
break;
}
case SCALAIRE_ENTIER:
{ const Grandeur_scalaire_entier& tyTQ= *((Grandeur_scalaire_entier*)
(*imap).second.Const_Grandeur_pointee());
// on renseigne l'argument
(*(*it).Grandeur_pointee()) = tyTQ;
valeur_recuperer = true;
break;
}
default:
{cout << "\n erreur fatale dans l'utilisation de la grandeur locale,"
<< NomTypeQuelconque(enu)
<< " necessaire pour le calcul de la ponderation "
<< " cette grandeur n'est pas un scalaire, pour l'instant "
<< " on ne sait pas comment l'utiliser !! "
<< "\n Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee(... ";
Sortie(1);
};
break;
};
// on sort de la boucle
break;
};
// si la grandeur a été récupérée, on l'exclue de la liste à transmettre
// à la méthode Valeurs_Tensorielles_interpoler_ou_calculer
if (valeur_recuperer)
exclure_local.push_back(enu);
};
// on merge les des deux exclure
if (deja_calculer_Q != NULL)
{List_io <const TypeQuelconque *>::const_iterator iv,ivfin = deja_calculer_Q->end();
for (iv=deja_calculer_Q->begin();iv != ivfin;iv++)
{ if (*iv != NULL)
{EnumTypeQuelconque enutq = (*iv)->EnuTypeQuelconque().EnumTQ();
exclure_local.push_back(enutq);
}
}
};
exclure_local.sort();
exclure_local.unique();
};
};
};
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_umat,&exclure_local);
// arrivée ici toutes les grandeurs évoluées de la fonction sont connues
// 2) on s'occupe maintenant des grandeurs quelconques autres qu'évoluées
// celles-ci , si elles ne sont pas déjà connues via la liste d'exclusion
// sont supposées être uniquement dans la map
// -> si on ne les trouve pas ==> pb !
const Tableau <EnumTypeQuelconque>& tqi = fct->Tab_enu_quelconque();
int taille_tqi = tqi.Taille();
Tableau < const TypeQuelconque * > tquelc(taille_tqi);
// on réutilise la map qui est aussi sensé contenir les quelconques autres
if (taille_tqi != 0)
{if(list_save == NULL)
{cout << "\n erreur fatale dans l'utilisation de grandeurs quelconque,"
<< " necessaire pour le calcul de la fonction nD "
<< " les grandeurs ne sont pas disponibles, a priori ??? "
<< "\n Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee(... ";
Sortie(1);
}
else
{ // on essaie de récupérer les types quelconques
// --- préparation pour les grandeurs déjà calculées et passée en paramètre
List_io<const TypeQuelconque *>::const_iterator ik,ikfin,ikdebut;
if (deja_calculer_Q != NULL)
{ ikfin = deja_calculer_Q->end();
ikdebut = deja_calculer_Q->begin();
};
// récupération des type quelconque sous forme d'un arbre pour faciliter la recherche
list <SaveResul*>::iterator isave=list_save->begin(); // pour les saveResul des lois
list <SaveResul*>::iterator isavedebut=list_save->begin(); // pour les saveResul des lois
list <SaveResul*>::iterator isavefin=list_save->end(); // pour les saveResul des lois
for (int i = 1; i <= taille_tqi; i++)
{bool valeur_recuperer = false; // pour savoir s'il y a eu récupération
EnumTypeQuelconque enu= tqi(i);
//a) on regarde si la grandeur ne fait pas partie des grandeurs déjà connues
// via l'exclusion passée en paramètre
if (deja_calculer_Q != NULL)
{// on parcours la liste: normalement il n'y a pas beaucoup d'éléments ...
// donc ce n'est "peut-être" pas trop couteux
for (ik= ikdebut; ik != ikfin;ik++)
{if ((*ik)->EnuTypeQuelconque().EnumTQ() == enu)
// on a trouvé
{tquelc(i) = (*ik);
valeur_recuperer = true;
break;
};
};
};
// on continue si on n'a pas récupéré
// on boucle sur les save_results pour récupérer les types quelconques
if (!valeur_recuperer)
for (isave = isavedebut; isave != isavefin;isave++)
// l'idée ici est de regarder si on ne trouve pas une des grandeurs demandées dans la map
{const map < EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >*
map_quelcon = (*isave)->Map_type_quelconque();
map < EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >::const_iterator imap;
if (map_quelcon != NULL) // sinon il n'y a pas de map
{imap = map_quelcon->find(enu);
if (imap != map_quelcon->end())
{ // on a trouvé la grandeur
tquelc(i) = &(*imap).second;
valeur_recuperer = true;
break;
};
};
};
if (!valeur_recuperer)// si on n'a pas trouvé, pb car on ne pourra pas fournir l'info nécessaire
// pour la fonction nD
{cout << "\n erreur fatale dans l'utilisation de la grandeur locale,"
<< NomTypeQuelconque(enu)
<< " necessaire pour le calcul de la fonction nD "
<< " cette grandeur n'est pas disponible, a priori ??? "
<< "\n Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee(... ";
Sortie(1);
};
};
};
};
try
{// maintenant calcul de la valeur et retour dans tab_ret
// sauf les valeurs à exclure
Tableau <double> & tab_val = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,&tquelc,NULL);
#ifdef MISE_AU_POINT
if (nb_retour > 0)
{if (tab_val.Taille() != nb_retour)
{ cout << "\nErreur : la fonction nD relative " << fct->NomFonction()
<< " doit calculer "<< nb_retour << " or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << "\n **** erreur: Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee \n";
Sortie(1);
};
};
#endif
// on retourne les infos
return tab_val;
}
catch(ErrCalculFct_nD &e)
{ cout << "\n Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee \n";
cout << "\n ** erreur dans l'appel de la fonction " << fct->NomFonction();
//cout << "\n pour debug Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee \n";
//Tableau <double> & tab_val = fct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,&tquelc,NULL);
throw e;
Sortie(1);
}
};
// calcul l'affichage, si celui-ci dépend d'une fonction nD
// ne doit être appelé que si la fonction nD existe (pas de vérification)
// en fait est utilisé uniquement par la méthode inline: Permet_affichage()
int Loi_comp_abstraite::Cal_permet_affichage()
{
// il faut calculer la fonction nD
// ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = permet_affich_loi_nD->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = permet_affich_loi_nD->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
// ici toutes les métriques sont NULL (les 3 premiers NULL)
// pas d'exclusion de Ddl_enum_etendu (le dernier NULL)
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,NULL,NULL,NULL,NULL));
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,NULL,NULL,NULL,NULL);
// arrivée ici toutes les grandeurs évoluées de la fonction sont connues
// 2) on s'occupe maintenant des grandeurs quelconques autres qu'évoluées
const Tableau <EnumTypeQuelconque>& tqi = permet_affich_loi_nD->Tab_enu_quelconque();
Valeurs_quelconque_interpoler_ou_calculer
(absolue,TEMPS_tdt,tqi,li_quelconque,NULL,NULL,NULL,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val
= permet_affich_loi_nD->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,&tab_pt_li_quelconque,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 1)
{ cout << "\nErreur : la fonction nD relative a l'affichage dans la loi de comportement "
<< " doit calculer un vecteur de dimention 1 or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " Loi_comp_abstraite::Cal_permet_affichage()\n";
Sortie(1);
};
#endif
// on récupère les valeurs du tableau
permet_affich_loi= tab_val(1);
// retour
return permet_affich_loi;
};
// sortie du niveau d'affichage
void Loi_comp_abstraite::Affiche_niveau_affichage()const
{ cout << " niveau_affichage_local= ";
if (permet_affich_loi_nD == NULL)
{ cout << permet_affich_loi; }
else
{ cout << " controle via la fonction nD: "<< permet_affich_loi_nD->NomFonction()
<< " (derniere valeur: "<<permet_affich_loi << ") ";
};
};
void Loi_comp_abstraite::Affiche_niveau_affichage(ofstream& sort,const int cas)
{ // on traite en fonction de l'existance ou non d'une fonction nD
if (permet_affich_loi_nD != NULL)
{sort << " 1 permet_affich_loi_nD: ";
LesFonctions_nD::Ecriture_pour_base_info(sort,cas,permet_affich_loi_nD);
}
else
{sort << " 0 " << permet_affich_loi ; };
};
void Loi_comp_abstraite::Lecture_permet_affichage
(ifstream& ent,const int cas,LesFonctions_nD& lesFonctionsnD)
{
// on traite en fonction de l'existance ou non d'une fonction nD
int test =0;string nom;
ent >> test;
if (test)
{ent >> nom;
if (nom != " permet_affich_loi_nD: ")
{ cout << "\n erreur en lecture de la fonction nD, on attendait "
<< " permet_affich_loi_nD: et on a lue " << nom
<< "\n Loi_comp_abstraite::Lecture_permet_affichage(...";
Sortie(1);
};
permet_affich_loi_nD = lesFonctionsnD.Lecture_pour_base_info
( ent,cas,permet_affich_loi_nD);
}
else
{ent >> permet_affich_loi;
if (permet_affich_loi_nD != NULL)
if (permet_affich_loi_nD->NomFonction() == "_") delete permet_affich_loi_nD;
permet_affich_loi_nD = NULL;
};
};
// lecture de l'affichage avec éventuellement une fonction nD
void Loi_comp_abstraite::Lecture_permet_affichage(UtilLecture * entreePrinc,LesFonctions_nD& lesFonctionsnD)
{// il faut voir s'il s'agit d'un niveau paramétrée par une fonction nD ou pas
string nom_class_methode("Loi_comp_abstraite::Lecture_permet_affichage(..");
if (strstr(entreePrinc->tablcar,"affichage_fonction_nD:")!= NULL)
{ string nom;
string mot_cle2="affichage_fonction_nD:";
string nom_fonct;
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer("**erreur en lecture** "+mot_cle2);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{permet_affich_loi_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
permet_affich_loi_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
permet_affich_loi_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (permet_affich_loi_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << permet_affich_loi_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur03** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde si la fonction nD intègre la température
const Tableau <Ddl_enum_etendu>& tab_enu = permet_affich_loi_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
};
}
else // sinon il s'agit d'un simple entier
{*(entreePrinc->entree) >> permet_affich_loi;
// gestion du pointeur de fonction
if (permet_affich_loi_nD != NULL)
if (permet_affich_loi_nD->NomFonction() == "_")
delete permet_affich_loi_nD;
};
};
Reference in a new issue Copy permalink