From 9c402fb5d8ee6013d363a709dd47c85c3f8f87b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?G=C3=A9rard=20Rio?= Date: Sat, 27 May 2023 11:33:09 +0200 Subject: [PATCH] =?UTF-8?q?mise=20=C3=A0=20jour=20version=20inter=20V=207.?= =?UTF-8?q?011?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../Geometrie/Frontiere/Ligne/FrontSegCub.cc | 8 +- Elements/Mecanique/ElemMeca2.cc | 5 +- Elements/Mecanique/ElemMeca4.cc | 5 +- Elements/ResRaid_MPI.cc | 165 ++ Elements/ResRaid_MPI.h | 183 ++ Elements/Trans_val_multi_tensoriel.cc | 928 ++++++++ Elements/Trans_val_multi_tensoriel.h | 212 ++ Enumeration/Enum_TypeQuelconque.cc | 15 +- Enumeration/Enum_TypeQuelconque.h | 2 +- Enumeration/Enum_ddl.cc | 4 +- Maillage/LesCondLim2.cc | 50 + Maillage/Maillage.cc | 367 +-- Parametres/EnteteParaGlob.h | 2 +- comportement/Hypo_elastique/Hypo_hooke1D.cc | 2045 +++++++++++++++++ comportement/Hypo_elastique/Hypo_hooke1D.h | 367 +++ contact/ElContact.cc | 84 +- contact/ElContact_2.cc | 383 +-- contact/LesContacts_3.cc | 201 +- 18 files changed, 4528 insertions(+), 498 deletions(-) create mode 100755 Elements/ResRaid_MPI.cc create mode 100755 Elements/ResRaid_MPI.h create mode 100755 Elements/Trans_val_multi_tensoriel.cc create mode 100755 Elements/Trans_val_multi_tensoriel.h create mode 100755 comportement/Hypo_elastique/Hypo_hooke1D.cc create mode 100755 comportement/Hypo_elastique/Hypo_hooke1D.h diff --git a/Elements/Geometrie/Frontiere/Ligne/FrontSegCub.cc b/Elements/Geometrie/Frontiere/Ligne/FrontSegCub.cc index 42163f4..2cc9ff2 100644 --- a/Elements/Geometrie/Frontiere/Ligne/FrontSegCub.cc +++ b/Elements/Geometrie/Frontiere/Ligne/FrontSegCub.cc @@ -62,10 +62,10 @@ FrontSegCub::FrontSegCub ( const Tableau & tab, const DdlElement& ddlEl { // au premier appel on construit la metrique associee if ( met == NULL) DefMetrique(); - // définition de d_T - int nb_ddl = 4 * tab(1)->Coord0().Dimension(); - d_T.Change_taille(nb_ddl); - }; + // définition de d_T + int nb_ddl = 4 * tab(1)->Coord0().Dimension(); + d_T.Change_taille(nb_ddl); + }; // de copie FrontSegCub::FrontSegCub( const FrontSegCub& a) : ElFrontiere(a),refP(a.refP),droite(a.droite),theta(a.theta),theta_repere(a.theta_repere) diff --git a/Elements/Mecanique/ElemMeca2.cc b/Elements/Mecanique/ElemMeca2.cc index 5e6fde5..061a5f6 100644 --- a/Elements/Mecanique/ElemMeca2.cc +++ b/Elements/Mecanique/ElemMeca2.cc @@ -1774,7 +1774,10 @@ void ElemMeca::Valeurs_Tensorielles_interpoler_ou_calculer case VITESSE:{ Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); Vitesse = gr.ConteneurCoordonnee(); break;} - + case ACCELERATION:{ + Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); + Vitesse = gr.ConteneurCoordonnee(); break;} + // dans le cas des numéros, traitement direct ici case NUM_ELEMENT: { *((Grandeur_scalaire_entier*) ((*ipq).Grandeur_pointee()))= num_elt; break;} diff --git a/Elements/Mecanique/ElemMeca4.cc b/Elements/Mecanique/ElemMeca4.cc index 701f722..1a99826 100644 --- a/Elements/Mecanique/ElemMeca4.cc +++ b/Elements/Mecanique/ElemMeca4.cc @@ -1779,7 +1779,10 @@ void ElemMeca::Valeurs_Tensorielles(bool absolue, Enum_dure temps,List_io. +// +// Herezh++ is distributed under GPL 3 license ou ultérieure. +// +// Copyright (C) 1997-2022 Université Bretagne Sud (France) +// AUTHOR : Gérard Rio +// E-MAIL : gerardrio56@free.fr +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty +// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// For more information, please consult: . + +#include "ResRaid_MPI.h" +#include +#include +#include +#include +namespace mpi = boost::mpi; + +// constructeur par défaut +ResRaid_MPI::ResRaid_MPI(): + res(NULL),raid(NULL),val_de_base() + {}; + +// constructeur de copie +ResRaid_MPI::ResRaid_MPI (const ResRaid_MPI& a): + res(NULL),raid(NULL) + ,val_de_base(a.val_de_base) + { // création et affectation des vecteurs + Creation_res_raid(true,a.res,a.raid); + }; + +// fonction d'un resraid et d'un élément particulier +// avec_recopie : indique si on veut une recopie ou pas +// si false: c'est seulement les tailles qui sont utilisées +ResRaid_MPI::ResRaid_MPI(const DeuxEntiers& elem, const Vecteur* ress, const Mat_pleine* raidd,bool avec_recopie): + res(NULL),raid(NULL),val_de_base() + { // on remplit val_de_base + // 1) on dimensionne + int taille = 2; // init avec les numéros + if (ress != NULL) + taille +=ress->Taille(); + if (raidd != NULL) + taille +=raidd->Nb_ligne() * raidd->Nb_colonne(); + // mise à jour de la taille globale + val_de_base.Change_taille(taille); + // affectation des numéros + val_de_base(1) = elem.un; + val_de_base(2) = elem.deux; + // création et affectation des vecteurs + Creation_res_raid(avec_recopie,ress,raidd); + }; + +// DESTRUCTEUR : +ResRaid_MPI::~ResRaid_MPI() + { // comme on a fait des new, il faut faire des delete + delete res; delete raid; + } ; + +// affectation à partir d'un élément particulier et d'un resraid +// opération analogue à la construction, mais dédié à une instance déjà existante +// correspond à la surcharge d'affectation, mais je préfère une fonction plus explicite ! +// pour éviter des affectations non voulues +void ResRaid_MPI::Affectation(const DeuxEntiers& elem, const Vecteur* ress, const Mat_pleine* raidd) + {// tout d'abord les num +// num_elt=elem.un; +// num_maillage=elem.deux; +// // puis les conteneurs +// if (ress != NULL) +// res = *(ress); +// if (raidd != NULL) +// raid = *(raidd); + }; + + +// 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 ResRaid_MPI::Lecture_base_info(ifstream& ent,const int cas) + { // on suit exactement la même procédure que pour archive +// load(ent,cas); + }; +// cas donne le niveau de sauvegarde +// = 1 : on sauvegarde tout +// = 2 : on sauvegarde uniquement les données variables (supposées comme telles) +void ResRaid_MPI::Ecriture_base_info(ofstream& sort,const int cas) + { // on suit exactement la même procédure que pour archive +// save(sort,cas); + }; + +// méthode interne pour créer les vecteurs tranches +// ce qui crée une liaison entre les deux stockages +// avec_recopie : indique si on veut une recopie ou pas +void ResRaid_MPI::Creation_res_raid(bool avec_recopie,const Vecteur* ress, const Mat_pleine* raidd) + { + double * pt = val_de_base.Pointeur_vect(); // le début du tableau + pt++; pt++; // on passe les numéros + bool memoire = false; // on va faire des tranches + if (ress != NULL) + {int taille1 = ress->Taille(); + // on crée un vecteur lié + res = new Vecteur(taille1,pt,memoire); + if (avec_recopie) // on recopie si c'est demandé + *res = *ress; + pt = &pt[2+taille1]; // on pointe sur le début de la suite + }; + if (raidd != NULL) + {// on crée une matrice liée + raid = new Mat_pleine(avec_recopie,pt,*raidd); + }; + + }; +// idem pour modifier éventuellement les tailles uniquement +void ResRaid_MPI::Change_tailles_res_raid(bool avec_recopie,const Vecteur* ress, const Mat_pleine* raidd) + { + double * pt = val_de_base.Pointeur_vect(); // le début du tableau + pt++; pt++; // on passe les numéros + bool memoire = false; // on travaille avec des tranches + + #ifdef MISE_AU_POINT + int taille_globale = 2; // init + if (ress != NULL) + taille_globale += ress->Taille(); + if (raidd != NULL) + taille_globale += raidd->Nb_ligne() * raidd->Nb_colonne(); + if (taille_globale != val_de_base.Taille()) + { cout << "\n*** Erreur : la taille de val_de_base: " << val_de_base.Taille() + << " est differente de 2 + la taille du vecteur: " << ress->Taille() + << " + la taille de la matrice : nb_lig*nb_col: " << raidd->Nb_ligne() * raidd->Nb_colonne(); + cout << "ResRaid_MPI::Change_tailles_res_raid(.. \n" << endl; + Sortie(1); + } + #endif +//à continuer à modifier ... +// if (ress != NULL) +// {int taille1 = ress->Taille(); +// // on crée un vecteur lié +// res = new Vecteur(taille1,pt,memoire); +// if (avec_recopie) // on recopie si c'est demandé +// *res = *ress; +// pt = &pt[2+taille1]; // on pointe sur le début de la suite +// }; +// if (raidd != NULL) +// {// on crée une matrice liée +// raid = new Mat_pleine(avec_recopie,pt,*raidd); +// }; + + }; + diff --git a/Elements/ResRaid_MPI.h b/Elements/ResRaid_MPI.h new file mode 100755 index 0000000..c6b8be8 --- /dev/null +++ b/Elements/ResRaid_MPI.h @@ -0,0 +1,183 @@ +// 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) . +// +// Herezh++ is distributed under GPL 3 license ou ultérieure. +// +// Copyright (C) 1997-2022 Université Bretagne Sud (France) +// AUTHOR : Gérard Rio +// E-MAIL : gerardrio56@free.fr +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty +// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// For more information, please consult: . + +/************************************************************************ + * DATE: 03/01/2022 * + * $ * + * AUTEUR: G RIO (mailto:gerardrio56@free.fr) * + * $ * + * PROJET: Herezh++ * + * $ * + ************************************************************************ + * BUT: une classe qui a pour objectif de servir, * + * d'interface entre un conteneur Element::ResRaid et * + * un conteneur qui est transmis entre CPU via MPI * + * $ * + * '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * + ************************************************************************/ + +#ifndef RESRAID_MPI_H +#define RESRAID_MPI_H + +#include +#include + +#include "Vecteur.h" +#include "Mat_pleine.h" +#include "Basiques.h" + +/** +* +* BUT: une classe qui a pour objectif de servir, +* d'interface entre un conteneur Element::ResRaid et +* un conteneur qui est transmis entre CPU via MPI +* +* +* \author Gérard Rio +* \version 1.0 +* \date 03/01/2022 +* \brief classe d'interface entre Element::ResRaid et un conteneur transmis entre CPU via MPI +* +*/ + +class ResRaid_MPI +{ + public : + // CONSTRUCTEURS : + // par défaut + ResRaid_MPI(); + + // fonction d'un resraid et d'un élément particulier + // avec_recopie : indique si on veut une recopie ou pas + // si false: c'est seulement les tailles qui sont utilisées + ResRaid_MPI(const DeuxEntiers& elem, const Vecteur* res, const Mat_pleine* raid,bool avec_recopie); + + // constructeur de copie + ResRaid_MPI (const ResRaid_MPI& a); + + // DESTRUCTEUR : + ~ResRaid_MPI() ; + + // METHODES PUBLIQUES : + + // affectation à partir d'un élément particulier et d'un resraid + // opération analogue à la construction, mais dédié à une instance déjà existante + // correspond à la surcharge d'affectation, mais je préfère une fonction plus explicite ! + // pour éviter des affectations non voulues + void Affectation(const DeuxEntiers& elem, const Vecteur* res, const Mat_pleine* raid); + + //============= lecture écriture dans base info ========== + // 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 Lecture_base_info(ifstream& ent,const int cas); + // cas donne le niveau de sauvegarde + // = 1 : on sauvegarde tout + // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) + void Ecriture_base_info(ofstream& sort,const int cas); + + // récupération des infos + int Num_elt() const {return val_de_base(1);}; // numero d'identification de l'element + int Num_mail() const {return val_de_base(2);}; // numéro de maillage + + const Vecteur & Const_res() const {return *res;}; // vecteur résidu + const Mat_pleine & Const_raid() const {return *raid ;}; // vecteur raideur + + + private : + // VARIABLES PROTEGEES : + // l'idée c'est d'avoir un seul gros conteneur qui englobe toutes les infos + // pour pouvoir faire un passage MPI natif: c-a-d un pointeur + un nombre + + // --1) le vecteur qui contient tout + Vecteur val_de_base; + // val_de_base(1) contient le num de l'element + // val_de_base(2) contient le num du maillage + + // --2) puis les grandeurs qui sont une tranche de val_de_base (vecteur et mat_pleine) + + Vecteur* res; // vecteur résidu + Mat_pleine* raid ; // vecteur raideur + + // méthode interne pour créer les vecteurs tranches + // ce qui crée une liaison entre les deux stockages + // avec_recopie : indique si on veut une recopie ou pas + void Creation_res_raid(bool avec_recopie,const Vecteur* ress, const Mat_pleine* raidd); + // idem pour modifier éventuellement les tailles uniquement + void Change_tailles_res_raid(bool avec_recopie,const Vecteur* ress, const Mat_pleine* raidd); + + +// on supprime la sérialisation, qui ne devrait plus être utilisée + +// // METHODES PROTEGEES : +// // -- serialisation --- +// // NB: pas trouvé comment déclarer l'implantation des méthodes template +// // en dehors du .h, donc on décrit le fct ici en inline +// // déclaration en friend pour l'acces direct de boost +// friend class boost::serialization::access; +// // on spécialise la sauvegarde et la restitution +// // version == 0 pour la première sauvegarde et ensuite > 0 +// // NB: c'est toujours la version en cours au moment de la sauvegarde +// // ==> dans notre cas, on ne sent sert pas pour l'instant: supposé tjs == 0 +// template +// void save(Archive & ar, unsigned int version) const +// { // on commence par les num mail et élément +// ar << std::string(" num_maillage ")<< num_maillage +// << std::string(" num_elem ") << num_elt; +// // puis le vecteur résidu +// ar << res; +// // puis la raideur +// ar << raid; +// } +// +// // en lecture, le num de version permet de ce positionner sur une version particulière +// template +// void load(Archive & ar, const unsigned int version) +// { // opération inverse de save +// std::string toto; +// // on commence par les num mail et élément +// ar >> toto >> num_maillage +// >> toto >> num_elt; +// // puis le vecteur résidu +// ar >> res; +// // puis la raideur +// ar >> raid; +// } +// +// // la macro suivante va définir automatiquement la méthode : "serialize" +// BOOST_SERIALIZATION_SPLIT_MEMBER() +// // pour mémoire on indique l'entête de la méthode "serialize" +//// // la méthode serialize fonctionne dans les deux sens: lecture et écriture dans ar +//// // << et >> est remplacé par & +//// // le choix dépend du fait que ar est un flux entrant ou sortant +//// template +//// void serialize(Archive & ar, const unsigned int version); + + }; + +#endif diff --git a/Elements/Trans_val_multi_tensoriel.cc b/Elements/Trans_val_multi_tensoriel.cc new file mode 100755 index 0000000..43cec07 --- /dev/null +++ b/Elements/Trans_val_multi_tensoriel.cc @@ -0,0 +1,928 @@ +// FICHIER : Trans_val_multi_tensoriel.cc +// CLASSE : Trans_val_multi_tensoriel + + +// 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) . +// +// Herezh++ is distributed under GPL 3 license ou ultérieure. +// +// Copyright (C) 1997-2022 Université Bretagne Sud (France) +// AUTHOR : Gérard Rio +// E-MAIL : gerardrio56@free.fr +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty +// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// For more information, please consult: . + + +// récupération des valeurs au numéro d'ordre = iteg pour +// les grandeur enu +// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière +// cas = -1 ou -2: cela signifie que la métrique vient juste d'être calculée au pti iteg +// sinon il faut recalculer qq éléments de la métrique +// NB Importante: il faut faire attention à ce que ces métriques soient identiques à celles qui ont servit +// pour le calcul des tenseurs: en particulier si c'est utilisé pour calculer les grandeurs pour le chargement +// il faut s'assurer que ce sont les "mêmes pti" qui servent pour la charge et pour la raideur !! +//Tableau ElemMeca::Valeur_multi +// (bool absolue, Enum_dure temps,const List_io& enu,int iteg,int cas +// ) +// +// { // ----- def de grandeurs de travail +// // 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 = lesPtIntegMecaInterne->DimTens();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 +// if (absolue) +// dim_sortie_tenseur = ParaGlob::Dimension(); +// // --- 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; +// +// #ifdef MISE_AU_POINT +// if(iteg > lesPtIntegMecaInterne->NbPti()) +// { cout << "\n erreur, les informations demandees au pti "< 3 ) +// cout << "\n ElemMeca::Valeur_multi(.."; +// cout << endl; +// Sortie (1); +// }; +// #endif +// +// PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(iteg); +// // 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 +// TenseurHB* sigHB = NULL ; TenseurHB* sig_barreHB = NULL ; +// TenseurHB* epsHB = NULL ; TenseurHB* eps_barreHB = NULL ; +// TenseurBB* DepsBB = NULL ; TenseurHB* Deps_barreHB = NULL ; +// TenseurHB* DepsHB = NULL ; +// +// TenseurBB* epsAlmTotalBB=NULL; // pour la déformation totale d'almansi +// TenseurBB* epsGLTotalBB=NULL; // pour la déformation totale de green_lagrange +// TenseurBB* epsLogTotalBB=NULL; // pour la déformation totale logarithmique +// Coordonnee* Mtdt = NULL; // coordonnées finales éventuelles du point d'intégration considéré +// Coordonnee* Mt = NULL; // coordonnées à t éventuelles du point d'intégration considéré +// Coordonnee* M0 = NULL; // coordonnées initiales éventuelles du point d'intégration considéré +// 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 +// +// // on définie des indicateurs pour ne pas faire plusieurs fois le même calcul +// List_io::const_iterator ie,iefin=enu.end(); +// 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_vitesse_deformation_barre=false; +// bool besoin_des_valpropre_sigma=false; +// bool besoin_des_valpropre_deformation = false; bool besoin_des_valpropre_vitdef = false; +// bool besoin_deformation_logarithmique = false; bool besoin_deformation_greenlagrange = false; +// bool besoin_deformation_almansi = false; +// bool besoin_coordonnees = false; bool besoin_deplacements = false; +// bool besoin_coordonnees_t = false;bool besoin_coordonnees_t0 = false; +// for (ie=enu.begin(); ie!=iefin;ie++) +// { 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(); +// +// switch (posi) +// { case 1: case 2: case 3: case 4: case 5: case 6: +// {besoin_deformation_greenlagrange=true; break;} +// case 7: case 8: case 9: case 10: case 11: case 12: +// {besoin_deformation_almansi=true; break;} +// case 28: case 29: case 30: case 31: case 32: +// {besoin_des_valpropre_sigma=true; break;} +// case 25: case 26: case 27: case 77: +// {besoin_des_valpropre_deformation=true;break;} +// case 40: case 41: case 42: +// {besoin_des_valpropre_vitdef=true;break;} +// case 49: case 50: case 51: case 52: case 53: case 54: +// {besoin_deformation_logarithmique=true;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;} +// 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; +// }; +// +// }; +// +// +// // -- def de tenseurs pour la sortie +// TenseurHH* sigHH = NULL ; +// TenseurBB* eps0BB = NULL ; // pour def green-lagrange +// TenseurBB* epsBB = NULL ; // pour def courante +// TenseurBB* epslogBB = NULL ; // pour def logarithmique +// TenseurBB* epsAlmBB = NULL ; // pour def d'almansi +// TenseurBB* DeltaEpsBB = NULL ; +// bool besoin_matrice_chg_base = false; +// if (besoin_des_contraintes || besoin_des_valpropre_sigma || besoin_des_contraintes_barre) +// {sigHH = (NevezTenseurHH(dim_sortie_tenseur)) ; +// sigHB = (NevezTenseurHB(dim)) ; +// sig_barreHB = (NevezTenseurHB(dim)) ;besoin_matrice_chg_base=true; +// }; +// if (besoin_des_deformation || besoin_deformation_almansi || besoin_des_deformation_barre +// || besoin_deformation_greenlagrange || besoin_deformation_logarithmique ) +// {eps0BB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def green-lagrange +// epsBB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def courante +// epslogBB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def logarithmique +// epsAlmBB = (NevezTenseurBB(dim_sortie_tenseur)) ; // pour def d'almansi +// epsHB = (NevezTenseurHB(dim)) ; eps_barreHB = (NevezTenseurHB(dim)) ; +// DeltaEpsBB = (NevezTenseurBB(dim_sortie_tenseur)) ; +// besoin_matrice_chg_base=true; +// }; +// if (besoin_des_valpropre_vitdef || besoin_des_vitesses_deformation ) +// {DepsBB = (NevezTenseurBB(dim)) ; Deps_barreHB = (NevezTenseurHB(dim)) ; +// DepsHB = (NevezTenseurHB(dim)) ; +// besoin_matrice_chg_base=true; +// }; +// +// +// Tableau tab_ret (enu.size()); +// // 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 +// +// def->ChangeNumInteg(iteg); // on change le numéro de point d'intégration courant +// if (((cas == 1) || (cas == 2))) +// { // cas d'une premiere initialisation +// Tableau tab(5); +// tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt; +// tab(4) = igiH_0; tab(5) = igijHH_0; +// met->PlusInitVariables(tab) ; +// }; +// // éléments de métrique et matrices de passage +// TenseurHH* gijHH=NULL;TenseurBB* gijBB=NULL;BaseB* giB=NULL; +// BaseH* giH_0=NULL;BaseH* giH=NULL; +// BaseB* giB_0=NULL; +// BaseB* giB_t=NULL; // n'est définit "que" pour certains cas +// 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 +// {int dim_effective = dim; // init +// if (absolue) dim_effective = ParaGlob::Dimension(); +// 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); +// }; +// +// switch (cas) +// { case 1: case 11: +// // calcul d'une ortho interessance de visualisation des tenseurs +// // cas de tenseur 3D -> Ia, cas 1D on prend un vecteur norme collineaire +// // avec g1, dans le cas 2D +// // la nouvelle base gamma est calculee dans def par projection de "Ipa" sur Galpha +// // le resultat est une matrice de passage utilisable pour le changement de base +// // Aa(i,a) = Aa^i_{.a}, avec g^i = Aa^i_{.a} * Ip^a +// // tout ce passe comme si Ip^a est la nouvelle base vers laquelle on veut évoluer +// {if (besoin_matrice_chg_base) +// {const Met_abstraite::InfoImp& ex = def->RemontImp(absolue,*Aa0,*Aafin); +// gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; +// giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; +// giB_0 = ex.giB_0; +// } +// else +// {const Met_abstraite::InfoImp& ex = def->RemontImp(); +// gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; +// giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; +// giB_0 = ex.giB_0; +// }; +// break; +// } +// case -1: case -2: +// // identique au cas 1 mais avec le fait que la métrique est directement disponible +// // car le calcul est supposé suivre un calcul implicit au bon pti +// {if (besoin_matrice_chg_base) +// {Mat_pleine Aat(dim,dim); +// // 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->Remont_et_metrique_0_t_tdtSansCalMet(absolue,*Aa0,Aat,*Aafin); +// gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; +// giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; +// giB_0 = ex.giB_0;giB_t = ex.giB_t; +// } +// else +// {const Met_abstraite::Info_et_metrique_0_t_tdt ex = def->Remont_et_metrique_0_t_tdtSansCalMet(); +// gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; +// giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; +// giB_0 = ex.giB_0;giB_t = ex.giB_t; +// } +// break; +// } +// case 2: case 12: +// // resultat a t +// {if (besoin_matrice_chg_base) +// {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(absolue,*Aa0,*Aafin); +// gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; +// giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; +// giB_0 = ex.giB_0; +// } +// else +// {const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(); +// gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt; +// giB = ex.giB_tdt; giH_0 = ex.giH_0;giH = ex.giH_tdt; +// giB_0 = ex.giB_0; +// } +// break; +// }; +// default: +// cout << "\n *** cas non prevu cas = "<< cas +// << "\n ElemMeca::Valeur_multi(..." << endl; +// Sortie(1); +// }; +// if (besoin_matrice_chg_base) +// {// 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()); +// }; +// +// // définition des tenseurs si nécessaire +// +// +// // ----- maintenant on calcule les grandeurs nécessaires ----- +// // calcul des tenseurs +// bool plusZero = true; // s'il faut rajouter des termes, on met des 0 +// if (besoin_des_contraintes) +// { if (absolue) {(ptIntegMeca.SigHH())->BaseAbsolue(*sigHH,*giB);}// changement de base finale +// else {(*sigHH) = *(ptIntegMeca.SigHH());sigHH->ChBase(*gammafin);}; +// (*sigHB) = *(ptIntegMeca.SigHH()) * (*gijBB); +// }; +// if (besoin_des_deformation) +// {(*epsHB) = (*gijHH) * (*(ptIntegMeca.EpsBB())); +// if (absolue)// changement de base finale +// {(ptIntegMeca.DeltaEpsBB())->BaseAbsolue((*DeltaEpsBB),*giH); +// (ptIntegMeca.EpsBB())->BaseAbsolue((*epsBB),*giH);} +// else {(*DeltaEpsBB) = *(ptIntegMeca.DeltaEpsBB());DeltaEpsBB->ChBase(*betafin); +// (*epsBB) = *(ptIntegMeca.EpsBB());epsBB->ChBase(*betafin);}; +// switch (def->Type_de_deformation()) +// { case DEFORMATION_STANDART : // c'est à dire almansi +// { // Green-Lagrange +// if ( besoin_deformation_greenlagrange) +// { if (absolue) {(ptIntegMeca.EpsBB())->BaseAbsolue((*eps0BB),*giH_0);}// changement de base finale +// else {(*eps0BB) = *(ptIntegMeca.EpsBB());eps0BB->ChBase(*beta0);}; +// }; +// (*epsAlmBB) = (*epsBB); +// if (epsAlmTotalBB!=NULL) // cas avec dilatation et demande de def Almansi totale +// {TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epsAlmTotal_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epsAlmTotal_local_BB); +// if (absolue)// changement de base finale +// {epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);} +// else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB +// if (prevoir_change_dim_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite +// }; +// +// if (epsGLTotalBB!=NULL) // cas avec dilatation et demande de def Green_Lagrange totale +// {TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epsGLTotal_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epsGLTotal_local_BB); +// if (absolue)// changement de base finale +// {epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);} +// else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB +// if (prevoir_change_dim_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 ((besoin_deformation_logarithmique) || (epsLogTotalBB!=NULL)) +// {def->Change_type_de_deformation(DEFORMATION_LOGARITHMIQUE); +// if (besoin_deformation_logarithmique) // cas du calcul de la def logarithmique +// {TenseurBB* epslog_local_BB = epslogBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epslog_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epslog_local_BB); +// if (absolue)// changement de base finale +// {epslog_local_BB->BaseAbsolue((*epslogBB),*giH);} +// else {epslogBB->ChBase(*betafin);}; // ici epslog_local_BB == epslogBB +// if (prevoir_change_dim_tenseur) +// delete epslog_local_BB; // car pas utilisé ensuite +// }; +// if (epsLogTotalBB!=NULL) // cas avec dilatation et demande de def log totale +// {TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epsLogTotal_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epsLogTotal_local_BB); +// if (absolue)// changement de base finale +// {epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);} +// else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB +// if (prevoir_change_dim_tenseur) delete epsLogTotal_local_BB; // car pas utilisé ensuite +// }; +// def->Change_type_de_deformation(DEFORMATION_STANDART); // on revient au type initial +// }; +// break; +// } +// case DEFORMATION_LOGARITHMIQUE : +// { (*epslogBB) = (*epsBB); +// if (epsLogTotalBB!=NULL) // cas avec dilatation et demande de def log totale +// {TenseurBB* epsLogTotal_local_BB = epsLogTotalBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epsLogTotal_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epsLogTotal_local_BB); +// if (absolue)// changement de base finale +// {epsLogTotal_local_BB->BaseAbsolue(*epsLogTotalBB,*giH);} +// else {epsLogTotalBB->ChBase(*betafin);}; // ici epsLogTotal_local_BB == epsLogTotalBB +// if (prevoir_change_dim_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 (( besoin_deformation_greenlagrange || besoin_deformation_almansi) +// || (epsAlmTotalBB!=NULL) || (epsGLTotalBB!=NULL)) +// {def->Change_type_de_deformation(DEFORMATION_STANDART); +// if ( ( besoin_deformation_greenlagrange || besoin_deformation_almansi) ) // cas de la def d'almansi +// { TenseurBB* eps_local_BB = epsAlmBB; // par défaut +// if (prevoir_change_dim_tenseur) +// eps_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*eps_local_BB); +// if (absolue)// changement de base finale +// {eps_local_BB->BaseAbsolue(*epsAlmBB,*giH); +// eps_local_BB->BaseAbsolue(*eps0BB,*giH_0);} +// else {epsAlmBB->ChBase(*betafin); +// eps0BB->ChBase(*beta0);};// ici eps_local_BB == epsAlmBB +// }; +// if (epsAlmTotalBB!=NULL) // cas avec dilatation et demande de def Almansi totale +// {TenseurBB* epsAlmTotal_local_BB = epsAlmTotalBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epsAlmTotal_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epsAlmTotal_local_BB); +// if (absolue)// changement de base finale +// {epsAlmTotal_local_BB->BaseAbsolue(*epsAlmTotalBB,*giH);} +// else {epsAlmTotalBB->ChBase(*betafin);}; // ici epsAlmTotal_local_BB == epsAlmTotalBB +// if (prevoir_change_dim_tenseur) delete epsAlmTotal_local_BB; // car pas utilisé ensuite +// }; +// if (epsGLTotalBB!=NULL) // cas avec dilatation et demande de def Green_Lagrange totale +// {TenseurBB* epsGLTotal_local_BB = epsGLTotalBB; // par défaut +// if (prevoir_change_dim_tenseur) +// epsGLTotal_local_BB = NevezTenseurBB(dim); +// def->Cal_deformation (temps,*epsGLTotal_local_BB); +// if (absolue)// changement de base finale +// {epsGLTotal_local_BB->BaseAbsolue(*epsGLTotalBB,*giH_0);} +// else {epsGLTotalBB->ChBase(*beta0);}; // ici epsGLTotal_local_BB == epsGLTotalBB +// if (prevoir_change_dim_tenseur) delete epsGLTotal_local_BB; // car pas utilisé ensuite +// }; +// def->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->Type_de_deformation()) +// << " affichage donc errone des valeurs !!!"; +// }; +// }; +// if (besoin_des_vitesses_deformation) +// { if (absolue) {(ptIntegMeca.DepsBB())->BaseAbsolue(*DepsBB,*giH);}// changement de base finale +// else {(*DepsBB) = *(ptIntegMeca.DepsBB());DepsBB->ChBase(*betafin);}; +// (*DepsHB) = (*gijHH) * (*(ptIntegMeca.DepsBB())); +// }; +// +// int caas=0; Coordonnee valPropreSig,valPropreEps,valPropreDeps; // init a dim=0 +// if (besoin_des_valpropre_sigma) +// {valPropreSig = sigHB->ValPropre(caas); +// if (caas == -1) +// { cout << "\n warning *** erreur dans le calcul des valeurs propres de la contrainte (Valeur_multi)"; +// if (ParaGlob::NiveauImpression() >= 7) {sigHB->Ecriture(cout); cout << "\nElemMeca::Valeur_multi(...";}; +// cout << endl; +// // valPropreSig = sigHB.ValPropre(caas); // !!!!!!! pour débug +// }; +// }; +// if (besoin_des_valpropre_deformation) +// {valPropreEps = 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 << "\nElemMeca::Valeur_multi(...";}; +// cout << endl;}; +// }; +// if (besoin_des_valpropre_vitdef) +// {valPropreDeps = 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 << "\nElemMeca::Valeur_multi(...";}; +// cout << endl;}; +// }; +// if (besoin_coordonnees) +// {Mtdt = new Coordonnee(ParaGlob::Dimension()); +// *Mtdt = def->Position_tdt(); +// } +// if (besoin_coordonnees_t ) +// {*Mt = def->Position_tdt(); +// }; +// if (besoin_deplacements || besoin_coordonnees_t0) +// {if (M0 == NULL) +// M0 = new Coordonnee(ParaGlob::Dimension()); +// (*M0) = def->Position_0(); +// }; +// if (Vitesse != NULL) +// {Vitesse = new Coordonnee(ParaGlob::Dimension()); +// (*Vitesse) = def->VitesseM_tdt(); +// } +// if (besoin_des_contraintes_barre) +// {double Isig = sigHB->Trace(); // trace de la déformation +// (*sig_barreHB) = (*sigHB) - (Isig/ParaGlob::Dimension()) * (*Id_dim_HB(dim)); +// }; +// if (besoin_des_deformation_barre) +// {double Ieps = epsHB->Trace(); // trace de la déformation +// (*eps_barreHB) = (*epsHB) - (Ieps/ParaGlob::Dimension()) * (*Id_dim_HB(dim)); +// }; +// if (besoin_des_vitesse_deformation_barre) +// {double IDeps = DepsHB->Trace(); // trace de la déformation +// (*Deps_barreHB) = (*DepsHB) - (IDeps/ParaGlob::Dimension()) * (*Id_dim_HB(dim)); +// }; +// +// // def éventuelle de la contrainte de Mises +// double Mises = 0.; Coordonnee& vv = valPropreSig; int dimvec=vv.Dimension();// pour condenser l'écriture +// switch (dimvec) // dans le cas où dimvec=0 on ne fait rien, cas ou on n'a pas besoin de mises +// { case 1: Mises = Dabs(vv(1)); break; +// case 2: Mises = sqrt( ((vv(1)-vv(2))*(vv(1)-vv(2)) + vv(1) * vv(1) +// + vv(2) * vv(2)) * 0.5); break; +// case 3: 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); break; +// }; +// +// // def éventuelle du vecteur normal: ceci n'est correct qu'avec une métrique 2D +// if (N_surf != NULL) +// { // 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 element: " << Num_elt() << " pti "<< iteg +// << "\n ElemMeca::Valeur_multi(... "; +// }; +// 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) +// { // 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 ElemMeca::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) +// { // 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 ElemMeca::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 +// }; +// }; +// +// // def éventuelle de la déformation duale de mises = sqrt(2/3 * epsB:epsB) +////--- on utilise directement la grandeur stockée au pt d'integ, mais ici ne sert à rien, puis cette grandeur n'est pas utilisée par la suite ! +// +// +//// // l'expression est la même que celle de mises, ormis un coeff 4/9 qui permet de passer de 3/2 à 2/3 +//// double defDualMises = 0.; Coordonnee& vdef = valPropreEps; int dimvdef=vdef.Dimension();// pour condenser l'écriture +//// switch (dimvdef) // dans le cas où dimvdef=0 on ne fait rien, cas ou on n'a pas besoin de defDualMises +//// { case 1: defDualMises = Dabs(vdef(1)); break; +//// case 2: defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + vdef(1) * vdef(1) +//// + vdef(2) * vdef(2)) * 2./9.); break; +//// case 3: defDualMises = sqrt( ((vdef(1)-vdef(2))*(vdef(1)-vdef(2)) + (vdef(1)-vdef(3))*(vdef(1)-vdef(3)) +//// + (vdef(3)-vdef(2))*(vdef(3)-vdef(2))) * 2./9.); break; +//// }; +// +// // donnees propre a la loi mécanique au pt d'integ +// Loi_comp_abstraite::SaveResul * sDon = tabSaveDon(iteg); +// // donnees propre a la loi thermo physique au pt d'integ +// CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques +// if (loiTP != NULL) {sTP = tabSaveTP(iteg);}; // au pt d'integ si elles existes +// // donnees propre a la déformation mécanique au pt d'integ +// Deformation::SaveDefResul * sDefDon = tabSaveDefDon(iteg); +// +// // pour la sortie des grandeurs polaires (def et contrainte) +// double mini_Q = 5.e-5; +// double * Q_eps=NULL,* Q_sig=NULL,* Q_Deps; +// +// //----- fin du calcul des grandeurs nécessaires ----- +// +// // 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++) +// { // dans le cas où c'est une contrainte, une déformation ou d'une vitesse de déformation +// // il y a préparation des grandeurs à sortir +// 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)) ) +// {// 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())) +// { // récup de l'ordre +// Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*sigHH)(ij.i,ij.j); +// } +// else if ((Meme_famille((*ie).Enum(),EPS11)) && ((*ie).Nom_vide())) +// { // récup de l'ordre +// Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*epsBB)(ij.i,ij.j); +// } +// else if ((Meme_famille((*ie).Enum(),DEPS11)) && ((*ie).Nom_vide())) +// { // récup de l'ordre +// Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*DepsBB)(ij.i,ij.j); +// } +// 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 1: case 2: case 3: case 4: case 5: case 6: +// /*Green-Lagrange */ +// { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*eps0BB)(ij.i,ij.j);break;} +// case 7: case 8: case 9: case 10: case 11: case 12: +// /*Almansi */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*epsAlmBB)(ij.i,ij.j);break;} +// case 49: case 50: case 51: case 52: case 53: case 54: +// /*logarithmique */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*epslogBB)(ij.i,ij.j);break;} +// case 55: case 56: case 57: case 58: case 59: case 60: +// /*Almansi totale */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*epsAlmTotalBB)(ij.i,ij.j);break;} +// case 61: case 62: case 63: case 64: case 65: case 66: +// /*Green_Lagrange totale */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*epsGLTotalBB)(ij.i,ij.j);break;} +// case 67: case 68: case 69: case 70: case 71: case 72: +// /*Log totale */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*epsLogTotalBB)(ij.i,ij.j);break;} +// case 13: case 14: case 15: case 16: case 17: case 18: +// /*Cauchy_local */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*ptIntegMeca.SigHH())(ij.i,ij.j);break;} +// case 19: case 20: case 21: case 22: case 23: case 24: +// /*Almansi_local */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*ptIntegMeca.EpsBB())(ij.i,ij.j);break; } +// case 25: /*Def_principaleI*/ tab_ret(it)=valPropreEps(1);break; +// case 26: /*Def_principaleII*/ +// if (valPropreEps.Dimension() > 1) {tab_ret(it)=valPropreEps(2);} else {tab_ret(it)= 0.;};break; +// case 27: /*Def_principaleIII*/ +// if (valPropreEps.Dimension() > 2) {tab_ret(it)=valPropreEps(3);} else {tab_ret(it)= 0.;};break; +// case 28: /*Sigma_principaleI*/ tab_ret(it)=valPropreSig(1);break; +// case 29: /*Sigma_principaleII*/ +// if (valPropreSig.Dimension() > 1) {tab_ret(it)=valPropreSig(2);} else {tab_ret(it)= 0.;};break; +// case 30: /*Sigma_principaleIII*/ +// if (valPropreSig.Dimension() > 2) {tab_ret(it)=valPropreSig(3);} else {tab_ret(it)= 0.;};break; +// case 31: /*contrainte_mises*/ tab_ret(it)=Mises;break; +//// case 77: /*def_duale_mises*/ tab_ret(it)=defDualMises;break; +// case 77: /*def_duale_mises*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(2);break; +// case 87: /*def_equivalente*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(1);break; +// case 88: /*def_duale_mises_maxi*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(3);break; +// case 89: /*vitesse_def_equivalente*/ tab_ret(it)=ptIntegMeca.Deformation_equi_const()(4) * unSurDeltat;break; +// case 32: /*contrainte_tresca*/ +// { switch (dim) +// {case 1: tab_ret(it)=0.5 * valPropreSig(1);break; +// case 2: tab_ret(it)=0.5 * (valPropreSig(1)-valPropreSig(2));break; +// case 3: tab_ret(it)=0.5 * (valPropreSig(1)-valPropreSig(3));break; +// } +// break; +// } +// case 33: /*def_plastique_cumulee*/ +// { string nom_comport(loiComp->Nom_comport()); +// if ((strstr(nom_comport.c_str(),"PRANDTL_REUSS")!=NULL)) +// { tab_ret(it) = tabSaveDon(iteg)->Deformation_plastique(); +// } +// else +// { cout << "\n erreur, la déformation plastique n'est pas disponible" +// << "dans l element " << this->Num_elt() +// << "\n ElemMeca::Valeur_multi(...."; +// Sortie(1); +// } +// break; +// } +// case 34: case 35: case 36: case 37: case 38: case 39: +// /*Vit_def11 et suite */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*DepsBB)(ij.i,ij.j);break; } +// case 40: /*Vit_principaleI*/ tab_ret(it)=valPropreDeps(1);break; +// case 41: /*Vit_principaleII*/ +// if (valPropreDeps.Dimension() > 1) {tab_ret(it)=valPropreDeps(2);} else {tab_ret(it)= 0.;};break; +// case 42: /*Vit_principaleIII*/ tab_ret(it)=valPropreDeps(3);break; +// if (valPropreDeps.Dimension() > 2) {tab_ret(it)=valPropreDeps(3);} else {tab_ret(it)= 0.;};break; +// case 43: case 44: case 45: case 46: case 47: case 48: +// /*Delta_def11 et suite */ { Deuxentiers_enu ij = IJind((*ie).Enum(),nbcompo); +// tab_ret(it)= (*DeltaEpsBB)(ij.i,ij.j);break; } +// case 73: /*energie_elastique*/ +// if (temps==TEMPS_tdt) {tab_ret(it)=tab_energ(iteg).EnergieElastique();} +// else {tab_ret(it)=tab_energ_t(iteg).EnergieElastique();};break; +// case 74: /*dissipation_plastique*/ +// if (temps==TEMPS_tdt) {tab_ret(it)=tab_energ(iteg).DissipationPlastique();} +// else {tab_ret(it)=tab_energ_t(iteg).DissipationPlastique();};break; +// case 75: /*dissipation_visqueuse*/ +// if (temps==TEMPS_tdt) {tab_ret(it)=tab_energ(iteg).DissipationVisqueuse();} +// else {tab_ret(it)=tab_energ_t(iteg).DissipationVisqueuse();};break; +// +// case 78: /*Spherique_eps*/ tab_ret(it)=epsHB->Trace()/3.; break; //ParaGlob::Dimension();break; modi du 5/2/2012 +// case 79: /*Q_eps*/ {Q_eps = new double; tab_ret(it)=*Q_eps = sqrt(eps_barreHB->II());break;} +// case 80: /*Cos3phi_eps*/ +// { double Qepsilon = ( (Q_eps!=NULL) ? *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.; +// tab_ret(it) = 3. * sqrt(6.) * bIIIb/ Qepsilon3; +//////------ debug +////cout << "\n debug: ElemMeca::Valeur_multi " +//// << "\n Qepsilon3= "<< Qepsilon3 << " bIIIb= "<< bIIIb +//// << " Cos3phi_eps= " << tab_ret(it) << " "; +//// +//////--- fin debug +// } +// else tab_ret(it)=0.; // sinon on le met à 0 +// break; +// } +// case 81: /*Spherique_sig*/ tab_ret(it)=sigHB->Trace()/ParaGlob::Dimension();;break; +// case 82: /*Q_sig*/ {Q_sig = new double; tab_ret(it)=*Q_sig = sqrt(sig_barreHB->II());break;} +// case 83: /*Cos3phi_sig*/ +// { double Qsig = ( (Q_sig!=NULL) ? *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.; +// tab_ret(it) = 3. * sqrt(6.) * bIIIb/ Qsig3; +// } +// else tab_ret(it)=0.; // sinon on le met à 0 +// break; +// } +// case 84: /*Spherique_Deps*/ tab_ret(it)=DepsHB->Trace()/ParaGlob::Dimension();break; +// case 85: /*Q_Deps*/ {Q_Deps = new double; tab_ret(it)=*Q_Deps = sqrt(Deps_barreHB->II());break;} +// case 86: /*Cos3phi_Deps*/ +// { double QDepsilon = ( (Q_Deps!=NULL) ? *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.; +// tab_ret(it) = 3. * sqrt(6.) * bIIIb/ QDepsilon3; +// } +// else tab_ret(it)=0.; // sinon on le met à 0 +// break; +// } +//// le cas 94 est a priori ok, mais il y a un pb de logique: en fait les pti d'efforts externes ne sont pas +//// a priori les mêmes que ceux de la raideur, donc cela va entrainer des pb +//// il faudrait une autre méthode spécifique aux efforts externes. +//// je laisse le code en prévision, mais je commente pour éviter les pb +// +//// case 94: /* "pression_ext */ +//// { if (lesChargeExtSurEle != NULL) +//// {if (lesChargeExtSurEle->LesPressionsExternes() != NULL) // cas où il existe des pressions sauvegardées +//// { Tableau >& lesPressionsExternes = *lesChargeExtSurEle->LesPressionsExternes(); +//// int nb_face = lesPressionsExternes.Taille(); +//// if (nb_face != 1) +//// {cout << "\n pas de sortie de pression possible en dehors d'element face "<< endl; +//// tab_ret(it)= 0.; +//// } +//// else +//// {int n_face=1; +//// int t_tail = lesPressionsExternes(n_face).Taille(); +//// if (t_tail != 0) +//// { Tableau & tab_press_appliquee = (lesPressionsExternes(n_face)); // pour simplifier +//// if (t_tail != 0) // cas où on a une face chargée +//// { tab_ret(it)=tab_press_appliquee(it).press;} +//// else {tab_ret(it)= 0.;}; // sinon rien +//// } +//// else {tab_ret(it)= 0.;}; // sinon rien +//// }; +//// } +//// else {tab_ret(it)= 0.;}; // sinon rien +//// } +//// else {tab_ret(it)= 0.;}; // sinon rien +//// 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;} +// case 137: // l'erreur relative +// {if (sigErreur_relative != NULL) +// {tab_ret(it)= *sigErreur_relative;} +// else {tab_ret(it)= 0.;}; +// break; +// } +// default : +// {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 ElemMeca::Valeur_multi(...."; +// tab_ret(it) = 0.; +// }; +// } // fin cas **** 2 >>>>> +// } // " " " +// } // -- fin du cas ou c'est une grandeur liée aux contraintes ou déformations +// // cas de l'erreur +// else if (( (*ie).Enum() == ERREUR) && ((*ie).Nom_vide())) +// {if (sigErreur != NULL) +// tab_ret(it)= *sigErreur; +// else if (ParaGlob::NiveauImpression()>4) +// {cout << "\n pas encore de ddl erreur dans l'element " +// << "\n ElemMeca::Valeur_multi(...."; +// // this->Affiche(); +// }; +// } +// else if (( (*ie).Enum() == TEMP) && ((*ie).Nom_vide())) +// {// on vérifie que le ddl existe au premier noeud +// if (tab_noeud(1)->Existe_ici(TEMP)) +// {tab_ret(it)= def->DonneeInterpoleeScalaire(TEMP,temps);} +// else if (ParaGlob::NiveauImpression()>3) +// {cout << "\n pas de ddl temperature disponible au noeud " +// << "\n ElemMeca::Valeur_multi(...."; +// // this->Affiche(); +// }; +// } +// 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 ElemMeca::Valeur_multi(...."; +// }; +// };// -- fin de la boucle sur la liste de Ddl_enum_etendu +// +// // 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 +// 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 +// 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 +// if (Vitesse != NULL) delete Vitesse; // vitesse +// // 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(); +// def->Retour_pti_precedant(); // on revient au pti précédent +// +// return tab_ret; +// }; +// + diff --git a/Elements/Trans_val_multi_tensoriel.h b/Elements/Trans_val_multi_tensoriel.h new file mode 100755 index 0000000..c132fba --- /dev/null +++ b/Elements/Trans_val_multi_tensoriel.h @@ -0,0 +1,212 @@ + +// 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) . +// +// Herezh++ is distributed under GPL 3 license ou ultérieure. +// +// Copyright (C) 1997-2022 Université Bretagne Sud (France) +// AUTHOR : Gérard Rio +// E-MAIL : gerardrio56@free.fr +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty +// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// For more information, please consult: . + +/************************************************************************ + * DATE: 27/07/2022 * + * $ * + * AUTEUR: G RIO (mailto:gerardrio56@free.fr) * + * $ * + * PROJET: Herezh++ * + * $ * + ************************************************************************ + * BUT: Méthodes génériques de transformations * + * de tenseurs. Utilisées en particulier par ElemMeca et les + * lois de comportement + * $ * + * '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * + * VERIFICATION: * + * * + * ! date ! auteur ! but ! * + * ------------------------------------------------------------ * + * ! ! ! ! * + * $ * + * '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * + * MODIFICATIONS: * + * ! date ! auteur ! but ! * + * ------------------------------------------------------------ * + * $ * + ************************************************************************/ +#ifndef TRANS_VAL_MULTI_TENSORIEL_H +#define TRANS_VAL_MULTI_TENSORIEL_H + +#include "Tenseur.h" +#include "NevezTenseur.h" +#include "Enum_calcul_masse.h" +#include "Basiques.h" +#include "Enum_dure.h" +#include "LesPtIntegMecaInterne.h" +#include "Enum_StabHourglass.h" +#include "Temps_CPU_HZpp.h" + +/** +* +* BUT: Méthodes génériques de transformations +* de tenseurs. Utilisées en particulier par ElemMeca et les +* lois de comportement +* +* \author Gérard Rio +* \version 1.0 +* \date 27/07/2022 +* \brief Méthodes génériques de transformations de tenseurs +* +*/ + + +class Trans_val_multi_tensoriel : +{ + public : + // VARIABLES PUBLIQUES : + +/* // CONSTRUCTEURS par défaut: + Trans_val_multi_tensoriel(); + // Constructeur de copie + Trans_val_multi_tensoriel(const Trans_val_multi_tensoriel& a); + + // DESTRUCTEUR : + ~Trans_val_multi_tensoriel(); + + // METHODES PUBLIQUES : + + + // récupération des valeurs au numéro d'ordre = iteg pour + // les grandeur enu + // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière + // NB Importante: il faut faire attention à ce que ces métriques soient identiques à celles qui ont servit + // pour le calcul des tenseurs: en particulier si c'est utilisé pour calculer les grandeurs pour le chargement + // il faut s'assurer que ce sont les "mêmes pti" qui servent pour la charge et pour la raideur !! + Tableau Valeur_multi(bool absolue,Enum_dure enu_t,const List_io& enu + ,int iteg,int cas + ) ; + + // récupération des valeurs Tensorielles (et non scalaire comme avec Valeur_multi) + // au numéro d'ordre = iteg pour les grandeur enu + // enu contient les grandeurs de retour + // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière + // NB Importante: il faut faire attention à ce que ces métriques soient identiques à celles qui ont servit + // pour le calcul des tenseurs: en particulier si c'est utilisé pour calculer les grandeurs pour le chargement + // il faut s'assurer que ce sont les "mêmes pti" qui servent pour la charge et pour la raideur !! + void Valeurs_Tensorielles(bool absolue, Enum_dure enu_t,List_io& enu + ,int iteg,int cas + ) ; + + + + + + // 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 + void Valeurs_Tensorielles_interpoler_ou_calculer + (bool absolue, Enum_dure temps,List_io& enu + ,Deformation & defor + ,const Met_abstraite::Impli* ex_impli + ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt + ,const Met_abstraite::Expli* ex_expli + ); + + + // 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 + // 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 et matrice de passage (c-a-d gijHH == NULL) + // ==> tous les autres pointeurs d'éléments de métrique et matrices de passages sont alors réputés NULL + // donc inutilisable + // , dans ce cas on ne peut pas calculer certaines grandeurs + // -> il y a vérification + Tableau Valeur_multi_interpoler_ou_calculer + (bool absolue, Enum_dure temps,const List_io& enu + // é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 // n'est définit "que" pour certains cas + ,Mat_pleine* Aa0,Mat_pleine* Aafin + ,Mat_pleine* gamma0,Mat_pleine* gammafin + ,Mat_pleine* beta0,Mat_pleine* betafin + // exclusion du calcul + ,const List_io* exclure_dd_etend + ); + + // 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 + // 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 et matrice de passage (c-a-d gijHH == NULL) + // ==> tous les autres pointeurs d'éléments de métrique et matrices de passages sont alors réputés NULL + // donc inutilisable + // , dans ce cas on ne peut pas calculer certaines grandeurs + // -> il y a vérification + void Valeurs_Tensorielles_interpoler_ou_calculer + (bool absolue, Enum_dure temps,List_io& enu + // é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 // n'est définit "que" pour certains cas + ,Mat_pleine* Aa0,Mat_pleine* Aafin + ,Mat_pleine* gamma0,Mat_pleine* gammafin + ,Mat_pleine* beta0,Mat_pleine* betafin + // exclusion du calcul + ,const List_io* exclure_Q + ); + + // 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 + // 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 et matrice de passage (c-a-d gijHH == NULL) + // ==> tous les autres pointeurs d'éléments de métrique et matrices de passages sont alors réputés NULL + // donc inutilisable + // , dans ce cas on ne peut pas calculer certaines grandeurs + // -> il y a vérification + // retour: la list li_quelc + void Valeurs_quelconque_interpoler_ou_calculer + (bool absolue, Enum_dure temps + ,const Tableau & tqi + ,List_io& li_quelc + // é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 // n'est définit "que" pour certains cas + ,Mat_pleine* Aa0,Mat_pleine* Aafin + ,Mat_pleine* gamma0,Mat_pleine* gammafin + ,Mat_pleine* beta0,Mat_pleine* betafin + // exclusion du calcul + ,const List_io* exclure_Q + ); +*/ + }; + + + +#endif diff --git a/Enumeration/Enum_TypeQuelconque.cc b/Enumeration/Enum_TypeQuelconque.cc index 97fdfe6..f8366e7 100644 --- a/Enumeration/Enum_TypeQuelconque.cc +++ b/Enumeration/Enum_TypeQuelconque.cc @@ -66,7 +66,7 @@ using namespace std; //introduces namespace std map_EnumTypeQuelconque["CONTRAINTE_COURANTE"]=CONTRAINTE_COURANTE; map_EnumTypeQuelconque["DEFORMATION_COURANTE"]=DEFORMATION_COURANTE; map_EnumTypeQuelconque["VITESSE_DEFORMATION_COURANTE"]=VITESSE_DEFORMATION_COURANTE; - map_EnumTypeQuelconque["ALMANSI"]=ALMANSI; + map_EnumTypeQuelconque["ALMANSI"]=ALMANSI; map_EnumTypeQuelconque["GREEN_LAGRANGE"]=GREEN_LAGRANGE; map_EnumTypeQuelconque["LOGARITHMIQUE"]=LOGARITHMIQUE; map_EnumTypeQuelconque["DELTA_DEF"]=DELTA_DEF; @@ -277,6 +277,7 @@ using namespace std; //introduces namespace std map_EnumTypeQuelconque["GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT"]=GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT; map_EnumTypeQuelconque["DEPLACEMENT"]=DEPLACEMENT; map_EnumTypeQuelconque["VITESSE"]=VITESSE; + map_EnumTypeQuelconque["ACCELERATION"]=ACCELERATION; map_EnumTypeQuelconque["DELTA_XI"]=DELTA_XI; map_EnumTypeQuelconque["XI_ITER_0"]=XI_ITER_0; map_EnumTypeQuelconque["MASSE_RELAX_DYN"]=MASSE_RELAX_DYN; @@ -290,7 +291,7 @@ using namespace std; //introduces namespace std map_EnumTypeQuelconque["NUM_ARETE"]=NUM_ARETE; // remplissage de tt_GLOB - tt_GLOB.Change_taille(223); // ****** ne pas oublier de redéfinir la taille si l'on ajoute des termes enum + tt_GLOB.Change_taille(224); // ****** ne pas oublier de redéfinir la taille si l'on ajoute des termes enum tt_GLOB(RIEN_TYPEQUELCONQUE)=1; tt_GLOB(SIGMA_BARRE_BH_T)=0; tt_GLOB(CONTRAINTE_INDIVIDUELLE_A_CHAQUE_LOI_A_T_SANS_PROPORTION)=0; @@ -298,7 +299,7 @@ using namespace std; //introduces namespace std tt_GLOB(CONTRAINTE_COURANTE)=1; tt_GLOB(DEFORMATION_COURANTE)=1; tt_GLOB(VITESSE_DEFORMATION_COURANTE)=1; - tt_GLOB(ALMANSI)=1; + tt_GLOB(ALMANSI)=1; tt_GLOB(GREEN_LAGRANGE)=1; tt_GLOB(LOGARITHMIQUE)=1; tt_GLOB(DELTA_DEF)=1; @@ -503,6 +504,7 @@ using namespace std; //introduces namespace std tt_GLOB(GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT)=1; tt_GLOB(DEPLACEMENT)=1; tt_GLOB(VITESSE)=1; + tt_GLOB(ACCELERATION)=1; tt_GLOB(DELTA_XI)=1; tt_GLOB(XI_ITER_0)=1; tt_GLOB(MASSE_RELAX_DYN)=1; @@ -516,7 +518,7 @@ using namespace std; //introduces namespace std tt_GLOB(NUM_ARETE)=1; // remplissage de tt_TQ_temps - tt_TQ_temps.Change_taille(223); // ****** ne pas oublier de redéfinir la taille si l'on ajoute des termes enum + tt_TQ_temps.Change_taille(224); // ****** ne pas oublier de redéfinir la taille si l'on ajoute des termes enum tt_TQ_temps(RIEN_TYPEQUELCONQUE)=TEMPS_0; tt_TQ_temps(SIGMA_BARRE_BH_T)=TEMPS_t; tt_TQ_temps(CONTRAINTE_INDIVIDUELLE_A_CHAQUE_LOI_A_T_SANS_PROPORTION)=TEMPS_tdt; @@ -729,6 +731,7 @@ using namespace std; //introduces namespace std tt_TQ_temps(GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT)=TEMPS_tdt; tt_TQ_temps(DEPLACEMENT)=TEMPS_tdt; tt_TQ_temps(VITESSE)=TEMPS_tdt; + tt_TQ_temps(ACCELERATION)=TEMPS_tdt; tt_TQ_temps(DELTA_XI)=TEMPS_tdt; tt_TQ_temps(XI_ITER_0)=TEMPS_tdt; tt_TQ_temps(MASSE_RELAX_DYN)=TEMPS_tdt; @@ -963,6 +966,7 @@ string NomTypeQuelconque(EnumTypeQuelconque id_TypeQuelconque) case GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT : result="GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT"; break; case DEPLACEMENT : result="DEPLACEMENT"; break; case VITESSE : result="VITESSE"; break; + case ACCELERATION : result="ACCELERATION"; break; case DELTA_XI : result="DELTA_XI"; break; case XI_ITER_0 : result="XI_ITER_0"; break; case MASSE_RELAX_DYN : result="MASSE_RELAX_DYN"; break; @@ -1206,6 +1210,7 @@ string NomTypeQuelconque_court (EnumTypeQuelconque id_TypeQuelconque) case GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT : result="GEN_VAR_GLOB_DOUBLE_USER_TDT"; break; case DEPLACEMENT : result="DEPLACEMENT"; break; case VITESSE : result="VITESSE"; break; + case ACCELERATION : result="ACCELERATION"; break; case DELTA_XI : result="DELTA_XI"; break; case XI_ITER_0 : result="XI_ITER_0"; break; case MASSE_RELAX_DYN : result="MASSE_RELAX_DYN"; break; @@ -1440,6 +1445,7 @@ string NomGeneriqueTypeQuelconque(EnumTypeQuelconque id_TypeQuelconque) case GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT : result="GEN_VAR_GLOB_DOUBLE_USER_TDT"; break; case DEPLACEMENT : result="DEPLACEMENT"; break; case VITESSE : result="VITESSE"; break; + case ACCELERATION : result="ACCELERATION"; break; case DELTA_XI : result="DELTA_XI"; break; case XI_ITER_0 : result="XI_ITER_0"; break; case MASSE_RELAX_DYN : result="MASSE_RELAX_DYN"; break; @@ -1727,6 +1733,7 @@ EnumTypeGrandeur Type_de_grandeur_associee(EnumTypeQuelconque typa) case GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT : return SCALAIRE; break; case DEPLACEMENT : return COORDONNEE; break; case VITESSE : return COORDONNEE; break; + case ACCELERATION : return COORDONNEE; break; case DELTA_XI : return COORDONNEE; break; case XI_ITER_0 : return COORDONNEE; break; case MASSE_RELAX_DYN : return COORDONNEE; break; diff --git a/Enumeration/Enum_TypeQuelconque.h b/Enumeration/Enum_TypeQuelconque.h index 3a419a8..f07d6e1 100644 --- a/Enumeration/Enum_TypeQuelconque.h +++ b/Enumeration/Enum_TypeQuelconque.h @@ -124,7 +124,7 @@ enum EnumTypeQuelconque { RIEN_TYPEQUELCONQUE = 1, SIGMA_BARRE_BH_T ,GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_0 ,GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_T ,GENERIQUE_UNE_VARIABLE_GLOB_DOUBLE_UTILISATEUR_TDT - ,DEPLACEMENT, VITESSE + ,DEPLACEMENT, VITESSE, ACCELERATION ,DELTA_XI,XI_ITER_0,MASSE_RELAX_DYN ,COMP_TORSEUR_REACTION ,NUM_NOEUD,NUM_MAIL_NOEUD,NUM_ELEMENT,NUM_MAIL_ELEM,NUM_PTI diff --git a/Enumeration/Enum_ddl.cc b/Enumeration/Enum_ddl.cc index ceb9cef..ba2c80f 100644 --- a/Enumeration/Enum_ddl.cc +++ b/Enumeration/Enum_ddl.cc @@ -953,9 +953,7 @@ EnumTypeQuelconque Equi_Enum_ddl_en_enu_quelconque(Enum_ddl a) case UX : case UY : case UZ : result=DEPLACEMENT; break; case V1 : case V2 : case V3 : result=VITESSE; break; case PR : result=UN_DDL_ENUM_ETENDUE; break; - case GAMMA1 : result=UN_DDL_ENUM_ETENDUE; break; - case GAMMA2 : result=UN_DDL_ENUM_ETENDUE; break; - case GAMMA3 : result=UN_DDL_ENUM_ETENDUE; break; + case GAMMA1 : case GAMMA2 : case GAMMA3 :result=ACCELERATION; break; case SIG11 : case SIG22 : case SIG33 : case SIG12 : case SIG23 : case SIG13 : result=CONTRAINTE_COURANTE; break; case ERREUR : result=UN_DDL_ENUM_ETENDUE; break; diff --git a/Maillage/LesCondLim2.cc b/Maillage/LesCondLim2.cc index 80b3f49..fdacaf9 100644 --- a/Maillage/LesCondLim2.cc +++ b/Maillage/LesCondLim2.cc @@ -1950,6 +1950,56 @@ void LesCondLim::Valeurs_Tensorielles_interpoler_ou_calculer }; break; } + case ACCELERATION : + {Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) ((*ipq).Grandeur_pointee())); + Coordonnee* gamma = gr.ConteneurCoordonnee(); + // on vérifie l'existante + #ifdef MISE_AU_POINT + if (!(noe.Existe_ici(GAMMA1))) + {cout << "\nErreur : l'acceleration n'existe pas au noeud " + << noe.Num_noeud() << " du maillage "<< noe.Num_Mail() + << "\nLesCondLim::Valeurs_Tensorielles_interpoler_ou_calculer(.. \n"; + if (tempsCL.Comptage_en_cours()) tempsCL.Arret_du_comptage(); + if (tempsCLL.Comptage_en_cours()) tempsCLL.Arret_du_comptage(); + Sortie(1); + }; + #endif + // récupération de la vitesse + switch (temps) + { case TEMPS_0 : + {switch (dim_espace) + {case 3: (*gamma)(3) = noe.Valeur_0(GAMMA3); + case 2: (*gamma)(2) = noe.Valeur_0(GAMMA2); + case 1: (*gamma)(1) = noe.Valeur_0(GAMMA1); + break; + default: break; + }; + break; + } + case TEMPS_t : + {switch (dim_espace) + {case 3: (*gamma)(3) = noe.Valeur_t(GAMMA3); + case 2: (*gamma)(2) = noe.Valeur_t(GAMMA2); + case 1: (*gamma)(1) = noe.Valeur_t(GAMMA1); + break; + default: break; + }; + break; + } + case TEMPS_tdt : + {switch (dim_espace) + {case 3: (*gamma)(3) = noe.Valeur_tdt(GAMMA3); + case 2: (*gamma)(2) = noe.Valeur_tdt(GAMMA2); + case 1: (*gamma)(1) = noe.Valeur_tdt(GAMMA1); + break; + default: break; + }; + break; + } + default: break; + }; + break; + } default : {// on initialise la grandeur pour éviter d'avoir des valeurs aléatoires ((*ipq).Grandeur_pointee())->InitParDefaut(); diff --git a/Maillage/Maillage.cc b/Maillage/Maillage.cc index 9ae2b34..66a8599 100755 --- a/Maillage/Maillage.cc +++ b/Maillage/Maillage.cc @@ -1388,70 +1388,77 @@ void Maillage::InitNormaleAuxNoeuds() Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Grandeur_pointee())); Coordonnee& normale = *gr.ConteneurCoordonnee(); normale.Zero(); // init - - List_io & li_elem = indice_NFr(ine);// la liste des Front qui contiennent le noeud - List_io ::iterator il,ilfin = li_elem.end(); - // on balaie les éléments - int nb_normale=0; // le nombre de normales trouvées - for (il = li_elem.begin(); il != ilfin; il++) - {Front& elem = *(*il); - Enum_type_geom enutygeom = elem.Eleme_const()->Type_geom_front(); - int cas = elem.Eleme()->CalculNormale_noeud(TEMPS_0,*tab_noeud(ine),coor_inter); -// //------ debug -// if ((noe.Num_noeud() == 6) && (noe.Num_Mail()==2)) -// { cout << "\n debug Maillage::InitNormaleAuxNoeuds():" -// << " noe "<< noe.Num_noeud() << " N= "; coor_inter.Affiche(); -// }; -// // ----- end debug - - if (cas == 0) - {cout << "\n *** erreur, le calcul de la normale n'est pas possible " - << " pour le noeud "<Num_noeud() - << " du maillage " << tab_noeud(ine)->Num_Mail() - << " par rapport a l'element de frontiere "; - elem.Affiche(1); - cout << "\n Maillage::InitNormaleAuxNoeuds() "; - Sortie(1); - } - else if (cas == 2) - {// le calcul n'est pas licite, mais ce n'est pas une erreur - // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien - } - else // sinon c'est ok - {// on peut avoir des directions très proches mais de sens inverse ... ce qui va - // conduire à une somme nulle, -// double intens_normale = normale.Norme(); // l'intensité actuelle -// Coordonnee inter_co = normale + coor_inter; -// if ((intens_normale > ConstMath::unpeupetit) && (inter_co.Norme() < ConstMath::petit)) -// // normalement ne devrait pas arriver -// { -// -// } - normale += coor_inter; // ici cela me semble plus juste: on ne doit pas avoir de pb d'inversion de normale + + // si on est en dimension 1, la normale est directement déterminé + if (dima == 1) + {normale(1) = 1.; + } + else + {// ici dima > 1, les éléments ont une métrique car ils sont différents des éléments points + List_io & li_elem = indice_NFr(ine);// la liste des Front qui contiennent le noeud + List_io ::iterator il,ilfin = li_elem.end(); + // on balaie les éléments + int nb_normale=0; // le nombre de normales trouvées + for (il = li_elem.begin(); il != ilfin; il++) + {Front& elem = *(*il); + Enum_type_geom enutygeom = elem.Eleme_const()->Type_geom_front(); + int cas = elem.Eleme()->CalculNormale_noeud(TEMPS_0,*tab_noeud(ine),coor_inter); + // //------ debug + // if ((noe.Num_noeud() == 6) && (noe.Num_Mail()==2)) + // { cout << "\n debug Maillage::InitNormaleAuxNoeuds():" + // << " noe "<< noe.Num_noeud() << " N= "; coor_inter.Affiche(); + // }; + // // ----- end debug - - -// //pour éviter cela, on regarde le produit scalaire -// // s'il est négatif on utilise l'inverse de la normale -// // donc en définitif on gardera globalement la direction précédente -// if ((normale * coor_inter) > 0.) -// {normale += coor_inter;} -// else -// {normale -= coor_inter;}; - nb_normale++; + if (cas == 0) + {cout << "\n *** erreur, le calcul de la normale n'est pas possible " + << " pour le noeud "<Num_noeud() + << " du maillage " << tab_noeud(ine)->Num_Mail() + << " par rapport a l'element de frontiere "; + elem.Affiche(1); + cout << "\n Maillage::InitNormaleAuxNoeuds() "; + Sortie(1); + } + else if (cas == 2) + {// le calcul n'est pas licite, mais ce n'est pas une erreur + // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien + } + else // sinon c'est ok + {// on peut avoir des directions très proches mais de sens inverse ... ce qui va + // conduire à une somme nulle, + // double intens_normale = normale.Norme(); // l'intensité actuelle + // Coordonnee inter_co = normale + coor_inter; + // if ((intens_normale > ConstMath::unpeupetit) && (inter_co.Norme() < ConstMath::petit)) + // // normalement ne devrait pas arriver + // { + // + // } + normale += coor_inter; // ici cela me semble plus juste: on ne doit pas avoir de pb d'inversion de normale + + + + // //pour éviter cela, on regarde le produit scalaire + // // s'il est négatif on utilise l'inverse de la normale + // // donc en définitif on gardera globalement la direction précédente + // if ((normale * coor_inter) > 0.) + // {normale += coor_inter;} + // else + // {normale -= coor_inter;}; + nb_normale++; + }; + }; + if (nb_normale != 0) + {normale /= nb_normale; + // enfin on normalise la normale + normale.Normer(); + // ce qui fini la mise à jour de la normale au noeud + // //------ debug + // if ((noe.Num_noeud() == 6) && (noe.Num_Mail()==2)) + // { cout << "\n debug Maillage::InitNormaleAuxNoeuds():" + // << " noe "<< noe.Num_noeud() << " N= "; normale.Affiche(); + // }; + // // ----- end debug }; - }; - if (nb_normale != 0) - {normale /= nb_normale; - // enfin on normalise la normale - normale.Normer(); - // ce qui fini la mise à jour de la normale au noeud -// //------ debug -// if ((noe.Num_noeud() == 6) && (noe.Num_Mail()==2)) -// { cout << "\n debug Maillage::InitNormaleAuxNoeuds():" -// << " noe "<< noe.Num_noeud() << " N= "; normale.Affiche(); -// }; -// // ----- end debug }; }; #ifdef MISE_AU_POINT @@ -1529,40 +1536,48 @@ void Maillage::MiseAjourNormaleAuxNoeuds() normale.Zero(); // init // calcul éventuel du tableau Indice this->Indice(); + int dima = ParaGlob::Dimension(); - List_io < Element*>& li_elem = indice(ine);// la liste des éléments qui contiennent le noeud - List_io < Element*>::iterator il,ilfin = li_elem.end(); - // on balaie les éléments - int nb_normale=0; // le nombre de normales trouvées - for (il = li_elem.begin(); il != ilfin; il++) - {Element& elem = *(*il); - Enum_type_geom enutygeom = Type_geom_generique(elem.Id_geometrie()); - // if ((enutygeom == LIGNE) || (enutygeom == SURFACE)) - if (enutygeom == SURFACE) // pour l'intant on ne traite que les surfaces - {int cas = elem.CalculNormale_noeud(TEMPS_t,*tab_noeud(ine),coor_inter); - if (cas == 0) - {cout << "\n *** erreur, le calcul de la normale n'est pas possible " - << " pour le noeud "<Num_noeud() - << " du maillage " << tab_noeud(ine)->Num_Mail() - << " par rapport a l'element "<< elem.Num_elt() - << " du maillage " << elem.Num_maillage() - << "\n Maillage::MiseAjourNormaleAuxNoeuds() "; - Sortie(1); - } - else if (cas == 2) - {// le calcul n'est pas licite, mais ce n'est pas une erreur - // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien - } - else // sinon c'est ok - {normale += coor_inter; - nb_normale++; + // si on est en dimension 1, la normale est directement déterminé + if (dima == 1) + {normale(1) = 1.; + } + else + {// ici dima > 1, les éléments ont une métrique car ils sont différents des éléments points + List_io < Element*>& li_elem = indice(ine);// la liste des éléments qui contiennent le noeud + List_io < Element*>::iterator il,ilfin = li_elem.end(); + // on balaie les éléments + int nb_normale=0; // le nombre de normales trouvées + for (il = li_elem.begin(); il != ilfin; il++) + {Element& elem = *(*il); + Enum_type_geom enutygeom = Type_geom_generique(elem.Id_geometrie()); + // if ((enutygeom == LIGNE) || (enutygeom == SURFACE)) + if (enutygeom == SURFACE) // pour l'intant on ne traite que les surfaces + {int cas = elem.CalculNormale_noeud(TEMPS_t,*tab_noeud(ine),coor_inter); + if (cas == 0) + {cout << "\n *** erreur, le calcul de la normale n'est pas possible " + << " pour le noeud "<Num_noeud() + << " du maillage " << tab_noeud(ine)->Num_Mail() + << " par rapport a l'element "<< elem.Num_elt() + << " du maillage " << elem.Num_maillage() + << "\n Maillage::MiseAjourNormaleAuxNoeuds() "; + Sortie(1); + } + else if (cas == 2) + {// le calcul n'est pas licite, mais ce n'est pas une erreur + // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien + } + else // sinon c'est ok + {normale += coor_inter; + nb_normale++; + }; }; }; - }; - if (nb_normale != 0) - {normale /= nb_normale; - // enfin on normalise la normale - normale.Normer(); + if (nb_normale != 0) + {normale /= nb_normale; + // enfin on normalise la normale + normale.Normer(); + }; }; // ce qui fini la mise à jour de la normale au noeud }; @@ -1578,44 +1593,52 @@ void Maillage::MiseAjourNormaleAuxNoeuds() Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Grandeur_pointee())); Coordonnee& normale = *gr.ConteneurCoordonnee(); normale.Zero(); // init - List_io & li_elem = indice_NFr(ine);// la liste des Front qui contiennent le noeud - List_io ::iterator il,ilfin = li_elem.end(); - // on balaie les éléments - int nb_normale=0; // le nombre de normales trouvées - for (il = li_elem.begin(); il != ilfin; il++) - {Front& elem = *(*il); - Enum_type_geom enutygeom = elem.Eleme_const()->Type_geom_front(); - int cas = elem.Eleme()->CalculNormale_noeud(TEMPS_t,*tab_noeud(ine),coor_inter); - if (cas == 0) - {cout << "\n *** erreur, le calcul de la normale n'est pas possible " - << " pour le noeud "<Num_noeud() - << " du maillage " << tab_noeud(ine)->Num_Mail() - << " par rapport a l'element de frontiere "; - elem.Affiche(1); - cout << "\n Maillage::InitNormaleAuxNoeuds() "; - Sortie(1); - } - else if (cas == 2) - {// le calcul n'est pas licite, mais ce n'est pas une erreur - // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien - } - else // sinon c'est ok - {// on peut avoir des directions très proches mais de sens inverse ... ce qui va - // conduire à une somme nulle, pour éviter cela, on regarde le produit scalaire - // s'il est négatif on utilise l'inverse de la normale - // donc en définitif on gardera globalement la direction précédente -// if ((normale * coor_inter) > 0.) - {normale += coor_inter;} -// else -// {normale -= coor_inter;}; - nb_normale++; + + // si on est en dimension 1, la normale est directement déterminé + if (dima == 1) + {normale(1) = 1.; + } + else + {// ici dima > 1, les éléments ont une métrique car ils sont différents des éléments points + List_io & li_elem = indice_NFr(ine);// la liste des Front qui contiennent le noeud + List_io ::iterator il,ilfin = li_elem.end(); + // on balaie les éléments + int nb_normale=0; // le nombre de normales trouvées + for (il = li_elem.begin(); il != ilfin; il++) + {Front& elem = *(*il); + Enum_type_geom enutygeom = elem.Eleme_const()->Type_geom_front(); + int cas = elem.Eleme()->CalculNormale_noeud(TEMPS_t,*tab_noeud(ine),coor_inter); + if (cas == 0) + {cout << "\n *** erreur, le calcul de la normale n'est pas possible " + << " pour le noeud "<Num_noeud() + << " du maillage " << tab_noeud(ine)->Num_Mail() + << " par rapport a l'element de frontiere "; + elem.Affiche(1); + cout << "\n Maillage::InitNormaleAuxNoeuds() "; + Sortie(1); + } + else if (cas == 2) + {// le calcul n'est pas licite, mais ce n'est pas une erreur + // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien + } + else // sinon c'est ok + {// on peut avoir des directions très proches mais de sens inverse ... ce qui va + // conduire à une somme nulle, pour éviter cela, on regarde le produit scalaire + // s'il est négatif on utilise l'inverse de la normale + // donc en définitif on gardera globalement la direction précédente + // if ((normale * coor_inter) > 0.) + {normale += coor_inter;} + // else + // {normale -= coor_inter;}; + nb_normale++; + }; + }; + if (nb_normale != 0) + {normale /= nb_normale; + // enfin on normalise la normale + normale.Normer(); }; }; - if (nb_normale != 0) - {normale /= nb_normale; - // enfin on normalise la normale - normale.Normer(); - }; // ce qui fini la mise à jour de la normale au noeud }; #ifdef MISE_AU_POINT @@ -1717,47 +1740,55 @@ void Maillage::MiseAjourNormaleAuxNoeuds_de_tdt_vers_T() Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Grandeur_pointee())); Coordonnee& normale = *gr.ConteneurCoordonnee(); normale.Zero(); // init - // la norme est nulle, on continue - List_io & li_elem = indice_NFr(ine);// la liste des Front qui contiennent le noeud - List_io ::iterator il,ilfin = li_elem.end(); - // on balaie les éléments - int nb_normale=0; // le nombre de normales trouvées - for (il = li_elem.begin(); il != ilfin; il++) - {Front& elem = *(*il); - Enum_type_geom enutygeom = elem.Eleme_const()->Type_geom_front(); - // ici contrairement à la méthode: MiseAjourNormaleAuxNoeuds - // on demande à l'élément de calculer la normale avec la géométrie à tdt - int cas = elem.Eleme()->CalculNormale_noeud(TEMPS_tdt,*tab_noeud(ine),coor_inter); - if (cas == 0) - {cout << "\n *** erreur, le calcul de la normale n'est pas possible " - << " pour le noeud "<Num_noeud() - << " du maillage " << tab_noeud(ine)->Num_Mail() - << " par rapport a l'element de frontiere "; - elem.Affiche(1); - cout << "\n Maillage::InitNormaleAuxNoeuds() "; - Sortie(1); - } - else if (cas == 2) - {// le calcul n'est pas licite, mais ce n'est pas une erreur - // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien - } - else // sinon c'est ok - {// on peut avoir des directions très proches mais de sens inverse ... ce qui va - // conduire à une somme nulle, pour éviter cela, on regarde le produit scalaire - // s'il est négatif on utilise l'inverse de la normale - // donc en définitif on gardera globalement la direction précédente -// if ((normale * coor_inter) > 0.) - {normale += coor_inter;} -// else -// {normale -= coor_inter;}; - nb_normale++; + + // si on est en dimension 1, la normale est directement déterminé + if (dima == 1) + {normale(1) = 1.; + } + else + {// ici dima > 1, les éléments ont une métrique car ils sont différents des éléments points + // la norme est nulle, on continue + List_io & li_elem = indice_NFr(ine);// la liste des Front qui contiennent le noeud + List_io ::iterator il,ilfin = li_elem.end(); + // on balaie les éléments + int nb_normale=0; // le nombre de normales trouvées + for (il = li_elem.begin(); il != ilfin; il++) + {Front& elem = *(*il); + Enum_type_geom enutygeom = elem.Eleme_const()->Type_geom_front(); + // ici contrairement à la méthode: MiseAjourNormaleAuxNoeuds + // on demande à l'élément de calculer la normale avec la géométrie à tdt + int cas = elem.Eleme()->CalculNormale_noeud(TEMPS_tdt,*tab_noeud(ine),coor_inter); + if (cas == 0) + {cout << "\n *** erreur, le calcul de la normale n'est pas possible " + << " pour le noeud "<Num_noeud() + << " du maillage " << tab_noeud(ine)->Num_Mail() + << " par rapport a l'element de frontiere "; + elem.Affiche(1); + cout << "\n Maillage::InitNormaleAuxNoeuds() "; + Sortie(1); + } + else if (cas == 2) + {// le calcul n'est pas licite, mais ce n'est pas une erreur + // simplement, l'élément n'est pas rusé pour ce noeud, on ne fait rien + } + else // sinon c'est ok + {// on peut avoir des directions très proches mais de sens inverse ... ce qui va + // conduire à une somme nulle, pour éviter cela, on regarde le produit scalaire + // s'il est négatif on utilise l'inverse de la normale + // donc en définitif on gardera globalement la direction précédente + // if ((normale * coor_inter) > 0.) + {normale += coor_inter;} + // else + // {normale -= coor_inter;}; + nb_normale++; + }; + }; + if (nb_normale != 0) + {normale /= nb_normale; + // enfin on normalise la normale + normale.Normer(); }; }; - if (nb_normale != 0) - {normale /= nb_normale; - // enfin on normalise la normale - normale.Normer(); - }; // ce qui fini la mise à jour de la normale au noeud }; #ifdef MISE_AU_POINT diff --git a/Parametres/EnteteParaGlob.h b/Parametres/EnteteParaGlob.h index 5677d97..cc1096e 100644 --- a/Parametres/EnteteParaGlob.h +++ b/Parametres/EnteteParaGlob.h @@ -41,7 +41,7 @@ EnumLangue ParaGlob::langueHZ = FRANCAIS; // langue utilisée pour les entrées sorties int ParaGlob::nbComposantesTenseur = 1; // nombre de composantes par defaut a 1 int ParaGlob::nivImpression = 2; // niveau d'impression - string ParaGlob::nbVersion = "7.010" ; // numéro de version du logiciel + string ParaGlob::nbVersion = "7.011" ; // numéro de version du logiciel string ParaGlob::NbVersionsurfichier = ""; // numéro de version lue en entrée fichier int ParaGlob::nb_diggit_double_calcul= 17; // nombre de chiffre significatifs utilisé pour // l'affichage des double précision pour l'archivage diff --git a/comportement/Hypo_elastique/Hypo_hooke1D.cc b/comportement/Hypo_elastique/Hypo_hooke1D.cc new file mode 100755 index 0000000..ea163ad --- /dev/null +++ b/comportement/Hypo_elastique/Hypo_hooke1D.cc @@ -0,0 +1,2045 @@ +// FICHIER : Hypo_hooke1D.cc +// CLASSE : Hypo_hooke1D + + +// 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) . +// +// Herezh++ is distributed under GPL 3 license ou ultérieure. +// +// Copyright (C) 1997-2022 Université Bretagne Sud (France) +// AUTHOR : Gérard Rio +// E-MAIL : gerardrio56@free.fr +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, +// or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty +// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +// See the GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +// For more information, please consult: . + + +//#include "Debug.h" + +# include +using namespace std; //introduces namespace std +#include +#include +#include "Sortie.h" +#include "TypeConsTens.h" +#include "ParaGlob.h" +#include "ConstMath.h" +#include "TypeQuelconqueParticulier.h" +#include "CharUtil.h" + + +#include "Hypo_hooke1D.h" + + +// constructeur par défaut à ne pas utiliser +Hypo_hooke1D::SaveResulLoi_Hypo1D::SaveResulLoi_Hypo1D() : + Kc(0.),Kc_t(0.),f(0.),f_t(0.) + ,eps22(0.),eps22_t(0.),eps33(0.),eps33_t(0.) + ,eps_cumulBB(),eps_cumulBB_t() + { cout << "\n erreur, le constructeur par defaut ne doit pas etre utilise !" + << "\n Hypo_hooke1D::SaveResulLoi_Hypo1D::SaveResulLoi_Hypo1D()"; + Sortie(1); + }; + +// le constructeur courant +Hypo_hooke1D::SaveResulLoi_Hypo1D::SaveResulLoi_Hypo1D + (SaveResul* ): + Kc(0.),Kc_t(0.),f(0.),f_t(0.) + ,eps22(0.),eps22_t(0.),eps33(0.),eps33_t(0.) + ,eps_cumulBB(),eps_cumulBB_t() + + { }; + +// 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 Hypo_hooke1D::SaveResulLoi_Hypo1D::Lecture_base_info + (ifstream& ent,const int cas) + { // ici toutes les données sont toujours a priori variables + // ou en tout cas pour les méthodes appelées, elles sont gérées par le paramètre: cas + string toto; ent >> toto; + #ifdef MISE_AU_POINT + if (toto != "HYPO_E_ISO_1D") + { cout << "\n erreur en lecture du conteneur pour la loi 1D, on a lue: " + <> toto >> Kc >> toto >> f + >> toto >> eps22 >> toto >> eps33 + >> toto >> eps_cumulBB; + Kc_t = Kc; f_t=f; eps22_t=eps22;eps33_t=eps33; + eps_cumulBB_t=eps_cumulBB; + }; + + +// cas donne le niveau de sauvegarde +// = 1 : on sauvegarde tout +// = 2 : on sauvegarde uniquement les données variables (supposées comme telles) +void Hypo_hooke1D::SaveResulLoi_Hypo1D::Ecriture_base_info + (ofstream& sort,const int cas) + { // ici toutes les données sont toujours a priori variables + // ou en tout cas pour les méthodes appelées, elles sont gérées par le paramètre: cas + sort << "\n HYPO_E_ISO_1D "; + // données de la loi + sort << "\n Kc= "<< Kc << " f= " << f + << " eps22= "<< eps22 << " eps33= "<< eps33 + << " eps_cumulBB= "<NomCourbe() == "_") + f_temperature = Courbe1D::New_Courbe1D(*(loi.f_temperature)); + if (Kc_temperature != NULL) + if (Kc_temperature->NomCourbe() == "_") + Kc_temperature = Courbe1D::New_Courbe1D(*(loi.Kc_temperature));; + if (f_IIeps != NULL) + if (f_IIeps->NomCourbe() == "_") + f_IIeps = Courbe1D::New_Courbe1D(*(loi.f_IIeps));; + if (Kc_IIeps != NULL) + if (Kc_IIeps->NomCourbe() == "_") + Kc_IIeps = Courbe1D::New_Courbe1D(*(loi.Kc_IIeps));; + // idem pour les fonctions nD + if (f_nD != NULL) + if (f_nD->NomFonction() == "_") + {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi) + string non_courbe("_"); + f_nD = Fonction_nD::New_Fonction_nD(*loi.f_nD); + }; + if (Kc_nD != NULL) + if (Kc_nD->NomFonction() == "_") + {// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi) + string non_courbe("_"); + Kc_nD = Fonction_nD::New_Fonction_nD(*loi.Kc_nD); + }; + }; + +Hypo_hooke1D::~Hypo_hooke1D () +// Destructeur +{ if (f_temperature != NULL) + if (f_temperature->NomCourbe() == "_") delete f_temperature; + if (Kc_temperature != NULL) + if (Kc_temperature->NomCourbe() == "_") delete Kc_temperature; + if (f_IIeps != NULL) + if (f_IIeps->NomCourbe() == "_") delete f_IIeps; + if (Kc_IIeps != NULL) + if (Kc_IIeps->NomCourbe() == "_") delete Kc_IIeps; + if (f_nD != NULL) + if (f_nD->NomFonction() == "_") delete f_nD; + if (Kc_nD != NULL) + if (Kc_nD->NomFonction() == "_") delete Kc_nD; +}; + +// Lecture des donnees de la classe sur fichier +void Hypo_hooke1D::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D + ,LesFonctions_nD& lesFonctionsnD) + { string nom; + string nom_class_methode("Hypo_hooke1D::LectureDonneesParticulieres"); + // lecture de la compressibilité instantanée + *(entreePrinc->entree) >> nom ; + #ifdef MISE_AU_POINT + if (nom != "Kc=") + { cout << "\n erreur en lecture de la loi de HYPO_ELAS1D, on attendait Kc= suivi d'un scalaire " + << " ou du mot cle Kc_thermo_dependant_ et ensuite une courbe " + << " ou du mot cle Kc_fonction_nD: et ensuite une fonction nD "; + entreePrinc->MessageBuffer("**erreur1 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + #endif + + // gestion d'erreur de syntaxe non permise + if(((strstr(entreePrinc->tablcar,"Kc_thermo_dependant_")!=0) || + (strstr(entreePrinc->tablcar,"Kc_fonction_nD:")!=0) + ) + && + (strstr(entreePrinc->tablcar,"Kc_avec_compressibilite_externe")!=0) + ) + { cout << "\n erreur en lecture de la loi de HYPO_ELAS1D, il n'est pas possible d'avoir une" + << " compressibilite externe " + << " et une compressibilite thermo_dependante ou via une fonction nD, simultanement: " + << " c'est incoherent !! "; + entreePrinc->MessageBuffer("**erreur1 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + + // on regarde si Kc est thermo dépendant + compress_thermophysique=false; // init par défaut + if(strstr(entreePrinc->tablcar,"Kc_thermo_dependant_")!=0) + { thermo_dependant=true; *(entreePrinc->entree) >> nom; + if (nom != "Kc_thermo_dependant_") + { cout << "\n erreur en lecture de la thermodependance de Kc, on aurait du lire le mot" + << " cle Kc_thermo_dependant_" + << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; + entreePrinc->MessageBuffer("**erreur2 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // lecture de la loi d'évolution de Kc en fonction de la température + *(entreePrinc->entree) >> nom; + // on regarde si la courbe existe, si oui on récupère la référence + if (lesCourbes1D.Existe(nom)) { Kc_temperature = lesCourbes1D.Trouve(nom);} + else + { // sinon il faut la lire maintenant + string non_courbe("_"); + Kc_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); + // lecture de la courbe + Kc_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc); + }; + // prepa du flot de lecture + entreePrinc->NouvelleDonnee(); + } + // sinon on regarde si la compressibilité sera fourni en externe + else if (strstr(entreePrinc->tablcar,"Kc_avec_compressibilite_externe")!=0) + {*(entreePrinc->entree) >> nom; + compress_thermophysique=true; + } + // sinon on regarde si le module dépend d'une fonction nD + else if(strstr(entreePrinc->tablcar,"Kc_fonction_nD:")!=0) + { string nom; + string mot_cle1="Kc="; + string mot_cle2="Kc_fonction_nD:"; + + // on passe le mot clé générique + bool lec = entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle1); + // on lit le nom de la fonction + string nom_fonct; + lec = lec && entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct); + if (!lec ) + { entreePrinc->MessageBuffer("**erreur02 en lecture** "+mot_cle2); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // maintenant on définit la fonction + if (lesFonctionsnD.Existe(nom_fonct)) + {Kc_nD = lesFonctionsnD.Trouve(nom_fonct); + } + else + {// sinon il faut la lire maintenant + string non("_"); + Kc_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct)); + // lecture de la courbe + Kc_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc); + // maintenant on vérifie que la fonction est utilisable + if (Kc_nD->NbComposante() != 1 ) + { cout << "\n erreur en lecture de Kc, la fonction " << nom_fonct + << " est une fonction vectorielle a " << Kc_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 & tab_enu = Kc_nD->Tab_enu_etendu(); + if (tab_enu.Contient(TEMP)) + thermo_dependant=true; + } + // sinon c'est directement le module que l'on lit + else + { // lecture de la valeur de Kc + *(entreePrinc->entree) >> Kc ; + }; + + // on regarde si Kc dépend multiplicativement directement du deuxième invariant + if(strstr(entreePrinc->tablcar,"Kc_IIeps_")!=0) + { // on vérifie que Kc_avec_compressibilite_externe n'est pas valide car dans ce cas + // la dépendance au deuxième invariant de la déformation n'est pas valide + if (compress_thermophysique) + { cout << "\n erreur on ne peut pas avoir Kc_avec_compressibilite_externe " + << " et Kc dependant du deuxieme invariant d'epsilon, "; + entreePrinc->MessageBuffer("**erreur5 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // si tout est ok on lit + *(entreePrinc->entree) >> nom; + if (nom != "Kc_IIeps_") + { cout << "\n erreur en lecture de la dependance de Kc au deuxieme invariant d'epsilon, " + << " on aurait du lire le mot cle Kc_IIeps_" + << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; + entreePrinc->MessageBuffer("**erreur6 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // lecture de la loi d'évolution en fonction du deuxième invariant + *(entreePrinc->entree) >> nom; + // on regarde si la courbe existe, si oui on récupère la référence + if (lesCourbes1D.Existe(nom)) + { Kc_IIeps = lesCourbes1D.Trouve(nom); + } + else + { // sinon il faut la lire maintenant + string non_courbe("_"); + Kc_IIeps = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); + // lecture de la courbe + Kc_IIeps->LectDonnParticulieres_courbes (non_courbe,entreePrinc); + } + // prepa du flot de lecture + if(strstr(entreePrinc->tablcar,"fin_loi_HYPO_ELAS1D")==0) entreePrinc->NouvelleDonnee(); + }; + + + // ---- lecture de la fonction f reliant \dot sig_1^{.1} et D_1^{.1} + *(entreePrinc->entree) >> nom; + if (nom != "f=") + { cout << "\n erreur en lecture de la fonction f , on aurait du lire le mot f="; + entreePrinc->MessageBuffer("**erreur3 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // on regarde si f est thermo dépendante + if(strstr(entreePrinc->tablcar,"f_thermo_dependant_")!=0) + { thermo_dependant=true; + *(entreePrinc->entree) >> nom; + if (nom != "f_thermo_dependant_") + { cout << "\n erreur en lecture de la thermodependance de f, on aurait du lire " + << " le mot cle f_thermo_dependant_" + << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; + entreePrinc->MessageBuffer("**erreur3 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // lecture de la loi d'évolution en fonction de la température + *(entreePrinc->entree) >> nom; + // on regarde si la courbe existe, si oui on récupère la référence + if (lesCourbes1D.Existe(nom)) + { f_temperature = lesCourbes1D.Trouve(nom); + } + else + { // sinon il faut la lire maintenant + string non_courbe("_"); + f_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); + // lecture de la courbe + f_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc); + }; + // prepa du flot de lecture + if(strstr(entreePrinc->tablcar,"fin_loi_HYPO_ELAS1D")==0) entreePrinc->NouvelleDonnee(); + } + // sinon on regarde si la fonction dépend d'une fonction nD + else if(strstr(entreePrinc->tablcar,"f_fonction_nD:")!=0) + { string nom; + string mot_cle1="f="; + string mot_cle2="f_fonction_nD:"; + + // on passe le mot clé générique + bool lec = entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle1); + // on lit le nom de la fonction + string nom_fonct; + lec = lec && entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct); + if (!lec ) + { entreePrinc->MessageBuffer("**erreur02 en lecture** "+mot_cle2); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // maintenant on définit la fonction + if (lesFonctionsnD.Existe(nom_fonct)) + {f_nD = lesFonctionsnD.Trouve(nom_fonct); + } + else + {// sinon il faut la lire maintenant + string non("_"); + f_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct)); + // lecture de la courbe + f_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc); + // maintenant on vérifie que la fonction est utilisable + if (f_nD->NbComposante() != 1 ) + { cout << "\n erreur en lecture de f , la fonction " << nom_fonct + << " est une fonction vectorielle a " << f_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 & tab_enu = f_nD->Tab_enu_etendu(); + if (tab_enu.Contient(TEMP)) + thermo_dependant=true; + } + else + { // lecture de f + *(entreePrinc->entree) >> f ; + }; + // on regarde si la fonction f dépend multiplicativement et directement du deuxième invariant + if(strstr(entreePrinc->tablcar,"f_IIeps_")!=0) + { *(entreePrinc->entree) >> nom; + if (nom != "f_IIeps_") + { cout << "\n erreur en lecture de la dependance de f au deuxieme invariant d'epsilon, " + << " on aurait du lire le mot cle f_IIeps_" + << " suivi du nom d'une courbe de charge ou de la courbe elle meme "; + entreePrinc->MessageBuffer("**erreur4 Hypo_hooke1D::LectureDonneesParticulieres (...**"); + throw (UtilLecture::ErrNouvelleDonnee(-1)); + Sortie(1); + }; + // lecture de la loi d'évolution en fonction du deuxième invariant + *(entreePrinc->entree) >> nom; + // on regarde si la courbe existe, si oui on récupère la référence + if (lesCourbes1D.Existe(nom)) + { f_IIeps = lesCourbes1D.Trouve(nom); + } + else + { // sinon il faut la lire maintenant + string non_courbe("_"); + f_IIeps = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str())); + // lecture de la courbe + f_IIeps->LectDonnParticulieres_courbes (non_courbe,entreePrinc); + } + // prepa du flot de lecture + if(strstr(entreePrinc->tablcar,"fin_loi_HYPO_ELAS1D")==0) entreePrinc->NouvelleDonnee(); + }; + + // ---- on regarde ensuite si le type de dérivée est indiqué + string toto; + if (strstr(entreePrinc->tablcar,"type_derivee")!=NULL) + { // lecture du type + *(entreePrinc->entree) >> toto >> type_derive; + if ((type_derive!=0)&&(type_derive!=-1)&&(type_derive!=1)) + { cout << "\n le type de derivee indique pour la loi de Hypo_hooke1D: "<< type_derive + << " n'est pas acceptable (uniquement -1 ou 0 ou 1), on utilise le type par defaut (-1)" + << " qui correspond a la derivee de Jauman "; + type_derive = -1; + }; + }; + // prepa du flot de lecture + if(strstr(entreePrinc->tablcar,"fin_loi_HYPO_ELAS1D")==0) entreePrinc->NouvelleDonnee(); + + // ---- on regarde ensuite s'il y a une restriction sur la contrainte + if (strstr(entreePrinc->tablcar,"restriction_traction_compression_")!=NULL) + { // lecture du type + *(entreePrinc->entree) >> toto >> restriction_traction_compression; + if ( (restriction_traction_compression!=0) + &&(restriction_traction_compression!=-1) + &&(restriction_traction_compression!=1) + ) + { cout << "\n le type de restriction lu pour la loi de Hypo_hooke1D: "<< restriction_traction_compression + << " n'est pas acceptable (uniquement -1 ou 0 ou 1), on utilise le type par defaut (0)" + << " qui correspond a 'aucune restriction' "; + restriction_traction_compression = 0; + }; + }; + // prepa du flot de lecture + if(strstr(entreePrinc->tablcar,"fin_loi_HYPO_ELAS1D")==0) entreePrinc->NouvelleDonnee(); + + // ---------- appel au niveau de la classe mère + Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire + (*entreePrinc,lesFonctionsnD); + + }; + +// affichage de la loi +void Hypo_hooke1D::Affiche() const + { cout << " \n loi de comportement HYPO_ELAS 1D "; + + if (Kc_temperature != NULL) { cout << " coef Kc thermo dependant explicitement " + << " courbe Kc=f(T): " << Kc_temperature->NomCourbe() <<" ";} + else if(compress_thermophysique) { cout << " compressibilite calculee avec une loi thermophysique " ;} + else if (Kc_nD != NULL) + {cout << "\n compressibilite calculee avec une fonction nD: "; + cout << "Kc_fonction_nD:" << " "; + if (Kc_nD->NomFonction() != "_") + cout << Kc_nD->NomFonction(); + else + Kc_nD->Affiche(); + } + else + { cout << " coef Kc= " << Kc ;} + if ( Kc_IIeps != NULL) { cout << " compressibilite dependant multiplicativement du deuxieme invariant d'epsilon " + << " courbe Kc_final=Kc*f(eps:eps): " << Kc_IIeps->NomCourbe() <<" ";}; + if ( f_temperature != NULL) { cout << " fonction f thermo dependante explicite " + << " courbe f=f(T): " << f_temperature->NomCourbe() <<" ";} + else if (f_nD != NULL) + {cout << "\n fonction f calculee avec une fonction nD: "; + cout << "f_fonction_nD:" << " "; + if (f_nD->NomFonction() != "_") + cout << f_nD->NomFonction(); + else + f_nD->Affiche(); + } + else + { cout << " cisaillement f= " << f ;} + if ( f_IIeps != NULL) { cout << " cisaillement dependant multiplicativement du deuxieme invariant d'epsilon " + << " courbe f_final=f*f(eps:eps): " << f_IIeps->NomCourbe() <<" ";}; + switch (type_derive) + { case -1: cout << ", et derivee de Jauman pour la contrainte" << endl;break; + case 0: cout << ", et derivee de Lie deux fois covariantes pour la contrainte" << endl; break; + case 1: cout << ", et derivee de Lie deux fois contravariantes pour la contrainte" << endl; break; + }; + // =0 -> pas de restriction + // = -1 : traction uniquement autorisée, la compression est mise à 0 + // = 1 : compression uniquement autorisée, la traction est mise à 0 + switch (restriction_traction_compression) + { case -1: cout << ", pas de compression possible " << endl;break; + case 0: cout << ", pas de restriction traction-compression " << endl; break; + case 1: cout << ", pas de traction possible " << endl; break; + }; + cout << endl; + // appel de la classe mère + Loi_comp_abstraite::Affiche_don_classe_abstraite(); + }; + +// affichage et definition interactive des commandes particulières à chaques lois +void Hypo_hooke1D::Info_commande_LoisDeComp(UtilLecture& entreePrinc) + { ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier + cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? "; + string rep = "_"; + // procédure de lecture avec prise en charge d'un retour chariot + rep = lect_return_defaut(true,"o"); + + if (Kc == -ConstMath::tresgrand) + { // on initialise à une valeur arbitraire + Kc = 10.;} + if (f == -ConstMath::trespetit) + { // on initialise à une valeur arbitraire + f = 0.15;} + sort << "\n# ......................... loi de comportement HYPO_ELAS 1D ........................." + << "\n# | coef compressibilite | rapport contrainte | type de derivee objective utilisee |" + << "\n# | instantane | deformation | pour le calcul de la contrainte |" + << "\n# | Kc | f |type_derivee (facultatif) |" + << "\n#......................................................................................" + << "\n Kc= " << setprecision(8) << Kc << " f= " << setprecision(8) << f + << " type_derivee " << type_derive + << "\n fin_loi_HYPO_ELAS1D " << endl; + if ((rep != "o") && (rep != "O" ) && (rep != "0") ) + {sort + << "\n# --------- type de derivee ----------" + << "\n# le type de derivee est optionnel: = -1 -> derivee de Jauman (valeur par defaut), " + << "\n# = 0 -> derivee deux fois covariantes , " + << "\n# = 1 -> derivee deux fois contravariantes" + << "\n# dans le cas ou l'on veut une valeur differente de la valeur par defaut" + << "\n# il faut mettre le mot cle" + << "\n# suivi de la valeur -1 ou 0 ou 1" + << "\n# "; + sort << "\n# --------- Kc calcule via une loi thermophysique ---------" + << "\n# pour le module de compressibilite: Kc, il est possible d'indiquer que" + << "\n# le module sera fourni " + << "\n# par une loi thermophysique (celle associee au meme element), pour ce faire" + << "\n# on indique uniquement: " + << "\n# Kc= Kc_avec_compressibilite_externe " + << "\n# NB: ensuite aucune autre info concernant Kc ne doit etre fournie. " + << "\n# "; + sort << "\n# ------- dependance explicite a la temperature ----------" + << "\n# \n# chacun des 2 parametres materiau Kc et f " + << "\n# peut etre remplace par une fonction dependante de la temperature " + << "\n# (sauf si le cas Kc_avec_compressibilite_externe est actif, si oui " + << "\n# aucune info complementaire ne peut-etre fournie pour Kc) " + << "\n# pour ce faire on utilise un mot cle puis une nom de courbe ou la courbe" + << "\n# directement comme avec les autre loi de comportement " + << "\n# exemple pour Kc: Kc= Kc_thermo_dependant_ courbe3 " + << "\n# exemple pour f: f= f_thermo_dependant_ courbe2 " + << "\n# IMPORTANT: a chaque fois qu'il y a une thermodependance, il faut passer une ligne " + << "\n# apres la description de la grandeur thermodependante, mais pas de passage " + << "\n# a la ligne si ce n'est pas thermo dependant" + << "\n# " + << "\n# -------- dependance a une fonction nD -------------" + << "\n# Kc et f peuvent au choix dependre d'une fonction nD. " + << "\n# Dans ce cas Kc et f ne doivent pas avoir ete declare thermodependant, " + << "\n# mais on peut inclure une thermo-dependance dans la fonction nD, " + << "\n# Exemple pour Kc : " + << "\n# Kc= Kc_fonction_nD: fonction_1 " + << "\n# " + << "\n# Exemple pour f : " + << "\n# f= f_fonction_nD: fonction_2 " + << "\n# " + << "\n# La declaration des fonctions suit la meme logique que pour les courbes 1D " + << "\n# " + << "\n# -------- dependance directe a (Eps:Eps) -------------" + << "\n# il est également possible de definir un facteur multiplicatif pour Kc et f " + << "\n# fonction explicite du deuxieme invariant de la deformation = (Eps:Eps) " + << "\n# = eps_1^1 * eps_1^1 " + << "\n# Ceci permet d'obtenir une evolution finale non lineaire dont l'operateur tangent" + << "\n# prend en compte les variations de la deformation." + << "\n# pour cela on indique le mot cle Kc_IIeps_ suivi du nom d'une courbe" + << "\n# suivi d'un passage a la ligne." + << "\n# Exemple de declaration d'un Kc thermo dependant et fonction " + << "\n# egalement du deuxieme invariant " + << "\n# Kc= Kc_thermo_dependant_ courbe2 " + << "\n# Kc_IIeps_ courbe5 " + << "\n# " + << "\n# Exemple de declaration d'un Kc non thermo-dependant, " + << "\n# mais fonction du deuxieme invariant" + << "\n# Kc= 10000 Kc_IIeps_ courbe4 " + << "\n# " + << "\n# NB: 1) on peut cumuler une dependance a une fonction nD et une dependance " + << "\n# explicite multiplicative a (Eps:Eps)" + << "\n# 2) a la fin de la definition de Kc_IIeps_ , il faut passer une ligne " + << "\n# " + << "\n# de la meme maniere on peut definir un facteur multiplicatif pour f " + << "\n# fonction du deuxieme invariant de la deformation = (Eps:Eps) . " + << "\n# pour cela on indique le mot cle f_IIeps_ suivi du nom d'une " + << "\n# courbe suivi d'un passage a la ligne" + << "\n# Exemple de declaration d'un f thermo dependant et fonction egalement " + << "\n# du deuxieme invariant " + << "\n# f= f_thermo_dependant_ courbe2 " + << "\n# f_IIeps_ courbe4 " + << "\n# " + << "\n# Exemple de declaration d'un f non thermo-dependant, mais fonction " + << "\n# du deuxieme invariant " + << "\n# f= 1200 f_IIeps_ courbe4 " + << "\n# " + << "\n# NB: 1) comme pour Kc, on peut cumuler une dependance a une fonction nD et " + << "\n# une dependance explicite multiplicative a (Eps:Eps)" + << "\n# 2) a la fin de la definition de f_IIeps_ , il faut passer une ligne " + << "\n# " + << "\n# " + << "\n# NB: pour toutes les definition de courbe ou de fct nD il est aussi possible" + << "\n# d'introduire directement la courbe a la place de son nom, comme pour " + << "\n# toutes les autres lois " + << "\n#" + << "\n# --------- restriction traction-compression ----------" + << "\n# a la suite du type de derivee on peut indiquer de maniere optionnelle " + << "\n# une restriction sur la contrainte calculee, via le mot cle: " + << "\n# restriction_traction_compression_ " + << "\n# suivi d'un scalaire tel que : " + << "\n# = 0 -> pas de restriction (valeur par defaut) " + << "\n# = -1 : traction uniquement autorisee, la compression est mise a 0 " + << "\n# = 1 : compression uniquement autorisée, la traction est mise a 0 " + << "\n# " + << "\n# ----------- fin des parametres --------" + << "\n# la derniere ligne doit contenir uniquement le mot cle: " + << "\n# fin_loi_HYPO_ELAS1D " + << "\n# "; + }; + sort << endl; + // appel de la classe mère + Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc); + }; + +// test si la loi est complete +int Hypo_hooke1D::TestComplet() + { int ret = LoiAbstraiteGeneral::TestComplet(); + if (!compress_thermophysique &&((Kc_temperature == NULL) && (Kc == -ConstMath::tresgrand)) + && (Kc_nD == NULL) + ) + { cout << " \n le coef de compressibilite instantane n'est pas defini pour la loi " << Nom_comp(id_comp) + << '\n'; + ret = 0; + } + if ((f_temperature == NULL) && (f == -ConstMath::trespetit) + && (f_nD == NULL) + ) + { cout << " \n la fonction f n'est pas defini pour la loi " << Nom_comp(id_comp) + << '\n'; + ret = 0; + } + return ret; + }; + +//----- lecture écriture de restart ----- + // 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 Hypo_hooke1D::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef + ,LesCourbes1D& lesCourbes1D + ,LesFonctions_nD& lesFonctionsnD) + { string nom; + if (cas == 1) + { ent >> nom; + if (nom != "HYPO_ELAS1D") + { cout << "\n erreur en lecture de la loi : Hypo_hooke1D, on attendait le mot cle : HYPO_ELAS1D " + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + // ensuite normalement il n'y a pas de pb de lecture puisque c'est écrit automatiquement (sauf en debug) + ent >> nom >> type_derive >> nom >> restriction_traction_compression; + // --- compressibilité + ent >> nom >> compress_thermophysique ; + int test; // sert pour le test des courbes + if (!compress_thermophysique) + { ent >> test; + switch (test) + {case 1: // cas d'une valeur numérique + ent >> Kc; break; + case 2: // cas d'une fonction nD + {ent >> nom ; + if (nom != " Kc_fonction_nD: ") + { cout << "\n erreur en lecture de la fonction nD, on attendait " + << " Kc_fonction_nD: et on a lue " << nom + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + Kc_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Kc_nD); + break; + } + case 3: // cas d'une fonction_Kc_temperature + {ent >> nom; + if (nom != " fonction_Kc_temperature ") + { cout << "\n erreur en lecture de la fonction Kc_temperature, on attendait " + << " fonction_Kc_temperature et on a lue " << nom + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + Kc_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Kc_temperature); + break; + } + default: + { cout << "\n erreur en lecture de Kc, on attendait un nombre 1,2 ou 3 " + << " et on a lue " << test + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + }; + // dépendance multiplicative à IIeps éventuelle + ent >> nom >> test; + if (!test) + { if (Kc_IIeps != NULL) {if (Kc_IIeps->NomCourbe() == "_") + delete Kc_IIeps; Kc_IIeps = NULL;}; + } + else + { ent >> nom; Kc_IIeps = lesCourbes1D.Lecture_pour_base_info(ent,cas,Kc_IIeps); }; + }; + // --- fonction f + ent >> nom >> test; + + switch (test) + {case 1: // cas d'une valeur numérique + ent >> f; break; + case 2: // cas d'une fonction nD + {ent >> nom ; + if (nom != " f_fonction_nD: ") + { cout << "\n erreur en lecture de la fonction nD, on attendait " + << " f_fonction_nD: et on a lue " << nom + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + f_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,f_nD); + break; + } + case 3: // cas d'une fonction_f_temperature + {ent >> nom; + if (nom != " fonction_f_temperature ") + { cout << "\n erreur en lecture de la fonction f_temperature, on attendait " + << " fonction_f_temperature et on a lue " << nom + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + f_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,f_temperature); + break; + } + default: + { cout << "\n erreur en lecture de f, on attendait un nombre 1,2 ou 3 " + << " et on a lue " << test + << "\n Hypo_hooke1D::Lecture_base_info_loi(..."; + Sortie(1); + }; + }; + // dépendance multiplicative à IIeps éventuelle + ent >> nom >> test; + if (!test) + { if (f_IIeps != NULL) {if (f_IIeps->NomCourbe() == "_") + delete f_IIeps; f_IIeps = NULL;}; + } + else + { ent >> nom; f_IIeps = lesCourbes1D.Lecture_pour_base_info(ent,cas,f_IIeps); }; + } + Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD); + }; + + // cas donne le niveau de sauvegarde + // = 1 : on sauvegarde tout + // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) +void Hypo_hooke1D::Ecriture_base_info_loi(ofstream& sort,const int cas) +{ if (cas == 1) + { sort << " HYPO_ELAS1D " + << " type_derivee " << type_derive + << " restriction_traction_compression "<& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list& decal ) const + {// tout d'abord on récupère le conteneur + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + // on passe en revue la liste + List_io::iterator itq,itqfin=liTQ.end(); + list::iterator idecal=decal.begin(); + for (itq=liTQ.begin();itq!=itqfin;itq++,idecal++) + {TypeQuelconque& tipParticu = (*itq); // pour simplifier + if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur + { EnumTypeQuelconque enuTQ = tipParticu.EnuTypeQuelconque().EnumTQ(); + switch (enuTQ) + {// 1) -----cas du module de compressibilité dépendant de la température + case COMPRESSIBILITE_TANGENTE: + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + tyTQ(1+(*idecal))=save_resul.Kc/3.;(*idecal)++; + break; + } + case MODULE_TANGENT_1D: + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + tyTQ(1+(*idecal))=save_resul.f;(*idecal)++; + break; + } + case DEF_EPAISSEUR: + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + tyTQ(1+(*idecal))=save_resul.eps33;(*idecal)++; + break; + } + case DEF_LARGEUR: + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + tyTQ(1+(*idecal))=save_resul.eps22;(*idecal)++; + break; + } + case DEF_ASSO_LOI: + { Tab_Grandeur_TenseurBB& tyTQ= *((Tab_Grandeur_TenseurBB*) (*itq).Grandeur_pointee()); // pour simplifier + tyTQ(1+(*idecal))=save_resul.eps_cumulBB;(*idecal)++; + break; + } + default: break; // sinon rien à faire + }; + }; // fin du if + }; // fin de la boucle + }; + +// récupération et création de la liste de tous les grandeurs particulières +// ces grandeurs sont ajoutées à la liste passées en paramètres +// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière +void Hypo_hooke1D::ListeGrandeurs_particulieres(bool absolue,List_io& liTQ) const +{ //on commence par définir une grandeur_scalaire_double + Tableau tab_1(1); + + Tab_Grandeur_scalaire_double grand_courant(tab_1); + // def d'un type quelconque représentatif à chaque grandeur + // a priori ces grandeurs sont défini aux points d'intégration identique à la contrainte par exemple + // enu_ddl_type_pt est définit dans la loi Abtraite générale + //on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée + + // $$$ cas de MODULE_COMPRESSIBILITE: + {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; + for (itq=liTQ.begin();itq!=itqfin;itq++) + if ((*itq).EnuTypeQuelconque() == COMPRESSIBILITE_TANGENTE) + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + int taille = tyTQ.Taille()+1; + tyTQ.Change_taille(taille); nexistePas = false; + }; + if (nexistePas) + {TypeQuelconque typQ1(COMPRESSIBILITE_TANGENTE,enu_ddl_type_pt,grand_courant); + liTQ.push_back(typQ1); + }; + }; + + // $$$ cas du module tangent + {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; + for (itq=liTQ.begin();itq!=itqfin;itq++) + if ((*itq).EnuTypeQuelconque() == MODULE_TANGENT_1D) + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + int taille = tyTQ.Taille()+1; + tyTQ.Change_taille(taille); nexistePas = false; + }; + if (nexistePas) + {TypeQuelconque typQ2(MODULE_TANGENT_1D,enu_ddl_type_pt,grand_courant); + liTQ.push_back(typQ2); + }; + }; + + // $$$ cas de la déformation en épaisseur (suivant 3) + {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; + for (itq=liTQ.begin();itq!=itqfin;itq++) + if ((*itq).EnuTypeQuelconque() == DEF_EPAISSEUR) + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + int taille = tyTQ.Taille()+1; + tyTQ.Change_taille(taille); nexistePas = false; + }; + if (nexistePas) + {TypeQuelconque typQ1(DEF_EPAISSEUR,enu_ddl_type_pt,grand_courant); + liTQ.push_back(typQ1); + }; + }; + + // $$$ cas de la déformation en largeur (suivant 2) + {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; + for (itq=liTQ.begin();itq!=itqfin;itq++) + if ((*itq).EnuTypeQuelconque() == DEF_LARGEUR) + { Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier + int taille = tyTQ.Taille()+1; + tyTQ.Change_taille(taille); nexistePas = false; + }; + if (nexistePas) + {TypeQuelconque typQ1(DEF_LARGEUR,enu_ddl_type_pt,grand_courant); + liTQ.push_back(typQ1); + }; + }; + // $$$ cas de la déformation cumulée associée à la loi, c-a-d au type d'intégration + {List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true; + for (itq=liTQ.begin();itq!=itqfin;itq++) + if ((*itq).EnuTypeQuelconque() == DEF_ASSO_LOI) + {Tab_Grandeur_TenseurBB& tyTQ= *((Tab_Grandeur_TenseurBB*) (*itq).Grandeur_pointee()); // pour simplifier + int taille = tyTQ.Taille()+1; + tyTQ.Change_taille(taille); nexistePas = false; + }; + if (nexistePas) + {TenseurBB* tens = NevezTenseurBB(1); // un tenseur typique + Tab_Grandeur_TenseurBB epsassoBB(*tens,1); + // def d'un type quelconque représentatif + TypeQuelconque typQ(DEF_ASSO_LOI,EPS11,epsassoBB); + liTQ.push_back(typQ); + delete tens; // car on n'en a plus besoin + }; + }; + + }; + +// calcul d'un module d'young équivalent à la loi +double Hypo_hooke1D::Module_young_equivalent(Enum_dure temps,const Deformation & ,SaveResul * saveDon) + { // ici le module d'Young correspond à la fonction f, cependant f est le module tangent et non + // sécant ... peut-être sera changé par la suite + // compte tenu du fait que l'on ne connait pas la métrique etc... on ramène le module en cours + + // on récupère le conteneur + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + // au niveau instantané on a: E=(3*Kc*mu)/(2*Kc+mu); + // et nu = (Kc-mu)/(2*Kc+mu) + double E=0.; + switch (temps) + {case TEMPS_0 : // rien n'a été calculé on utilise les grandeurs initiales + E=f; + break; + case TEMPS_t : // on utilise les grandeurs stockées à t + {E= save_resul.f_t; + break; + } + case TEMPS_tdt : // on utilise les grandeurs stockées à tdt + {E= save_resul.f; + break; + } + }; + return E; + }; + +// récupération d'un module de compressibilité équivalent à la loi, ceci pour un chargement nul +// il s'agit ici de la relation -pression = sigma_trace/3. = module de compressibilité * I_eps +// >>> en fait ici il s'agit du dernier module tangent calculé !! +double Hypo_hooke1D::Module_compressibilite_equivalent + (Enum_dure temps,const Deformation & def,SaveResul * saveDon) + { // ici le module correspond à Kc, cependant Kc est un module tangent et non + // sécant ... peut-être sera changé par la suite + // compte tenu du fait que l'on ne connait pas la métrique etc... on ramène le module en cours + + // on récupère le conteneur + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + double K=0.; + switch (temps) + {case TEMPS_0 : // rien n'a été calculé on utilise les grandeurs initiales + K=Kc/3.; + break; + case TEMPS_t : // on utilise les grandeurs stockées à t + {K=save_resul.Kc_t/3.; + break; + } + case TEMPS_tdt : // on utilise les grandeurs stockées à tdt + {K = save_resul.Kc/3.; + break; + } + }; + + return K; + }; + +// récupération de la variation relative d'épaisseur calculée: h/h0 +// 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 Hypo_hooke1D::HsurH0(SaveResul * saveDon) const + { // on récupère le conteneur + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + + // la déformation d'épaisseur obtenue dépend du type de dérivée retenue + double h_sur_h0 = 0.; // init + switch (type_derive) + {case 0: // cas d'une dérivée de Lie deux fois covariantes + // cas d'une déformation d'Almansi + { // epsBB33 = 1/2 * (1. - (h0/h)^2)= 1/2 * (1. - 1./(var_epai)^2), en orthonormee + h_sur_h0 = sqrt(1./(1.-2.* save_resul.eps33)); + }; + break; + case -1: // cas d'une dérivée de jauman: + // cas d'une def logarithmique cumulée que l'on approxime à une def logarithmique + { // eps_33 = log(var_epai); + h_sur_h0 = exp(save_resul.eps33); + }; + break; +//**** à finir +// case 1: // cas d'une dérivée de Lie deux fois contravariantes +// // cas d'une déformation e_{2} +// { // epsBB33 = 1/2 * (1. - (h0/h)^2)= 1/2 * (1. - 1./(var_epai)^2), en orthonormee +// h_sur_h0 = sqrt(1./(1.-2.* save_resul.eps33)); +// }; +// break; + default : + cout << "\nErreur : type de derivee qui n'est pas actuellement pris en compte, type= " + << type_derive; + cout << "\n Hypo_hooke1D::HsurH0 (... \n"; + Sortie(1); + }; + + + return h_sur_h0; + }; + + // 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 Hypo_hooke1D::BsurB0(SaveResul * saveDon) const + { // on récupère le conteneur + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + // la déformation d'épaisseur obtenue dépend du type de dérivée retenue + double b_sur_b0 = 0.; // init + switch (type_derive) + {case 0: // cas d'une dérivée de Lie deux fois covariantes + // cas d'une déformation d'Almansi + { // epsBB33 = 1/2 * (1. - (h0/h)^2)= 1/2 * (1. - 1./(var_epai)^2), en orthonormee + b_sur_b0 = sqrt(1./(1.-2.* save_resul.eps22)); + }; + break; + case -1: // cas d'une dérivée de jauman: + // cas d'une def logarithmique cumulée que l'on approxime à une def logarithmique + { // eps_33 = log(var_epai); + b_sur_b0 = exp(save_resul.eps22); + }; + break; + //**** à finir + // case 1: // cas d'une dérivée de Lie deux fois contravariantes + // // cas d'une déformation e_{2} + // { // epsBB33 = 1/2 * (1. - (h0/h)^2)= 1/2 * (1. - 1./(var_epai)^2), en orthonormee + // b_sur_b0 = sqrt(1./(1.-2.* save_resul.eps33)); + // }; + // break; + default : + cout << "\nErreur : type de derivee qui n'est pas actuellement pris en compte, type= " + << type_derive; + cout << "\n Hypo_hooke1D::HsurH0 (... \n"; + Sortie(1); + }; + + + return b_sur_b0; + }; + + + // ========== codage des METHODES VIRTUELLES protegees:================ + // calcul des contraintes a t+dt +void Hypo_hooke1D::Calcul_SigmaHH (TenseurHH& sigHH_t,TenseurBB& ,DdlElement & + ,TenseurBB & gijBB_t_,TenseurHH & gijHH_t_,BaseB& giB,BaseH& gi_H,TenseurBB& epsBB_ + ,TenseurBB& delta_epsBB_, TenseurBB& gijBB_ ,TenseurHH & gijHH_,Tableau & + ,double& jacobien_0,double& jacobien,TenseurHH & sigHH_ + ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement + ,const Met_abstraite::Expli_t_tdt& ex) + { + #ifdef MISE_AU_POINT + if (delta_epsBB_.Dimension() != 1) + { cout << "\nErreur : la dimension devrait etre 1 !\n"; + cout << " Hypo_hooke1D::Calcul_SigmaHH\n"; + Sortie(1); + }; + #endif + const Tenseur1BB & epsBB = *((Tenseur1BB*) &epsBB_); // passage en dim 1 + const Tenseur1BB & delta_epsBB = *((Tenseur1BB*) &delta_epsBB_); // passage en dim 1 + const Tenseur1HH & gijHH = *((Tenseur1HH*) &gijHH_); // " " " " + const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_); // " " " " + const Tenseur1HH & gijHH_t = *((Tenseur1HH*) &gijHH_t_); // " " " " + const Tenseur1BB & gijBB_t = *((Tenseur1BB*) &gijBB_t_); // " " " " + Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_); // " " " " + Tenseur1HH & sigHH_nn = *((Tenseur1HH*) &sigHH_t); // " " " " + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveResul); + + // tenseur intermédiaires utilisées selon les cas (par forcément dans tous les cas !!) + Tenseur1BH sigBH_n;Tenseur1HH sigHH_n; Tenseur1BB sigBB_n;Tenseur1BB sig_interBB_n; + switch (type_derive) + {case -1: // cas d'une dérivée de jauman: 1/2 deux fois covariant + deux fois contra + {sig_interBB_n = gijBB_t * sigHH_nn * gijBB_t; + sigBH_n = 0.5*( sig_interBB_n * gijHH + gijBB * sigHH_nn) ; + sigHH_n = gijHH * sigBH_n ; + sigBB_n = sigBH_n * gijBB; + // pour la déformation cumulée associée, + Tenseur1HH eps_interHH_n = gijHH_t * save_resul.eps_cumulBB_t * gijHH_t; + Tenseur1BB epsBB_n = 0.5*( gijBB * eps_interHH_n * gijBB + save_resul.eps_cumulBB_t); + save_resul.eps_cumulBB = delta_epsBB + epsBB_n; + break;} + case 0: //cas d'une dérivée de Lie deux fois covariantes + {sigBB_n = gijBB_t * sigHH_nn * gijBB_t; + sigBH_n = sigBB_n * gijHH ; + sigHH_n = gijHH * sigBH_n ; + // pour la déformation cumulée associée, on sait que cela donne directement la def d'Almansi + // mais la def passée en paramètre pourrait ne pas être d'Almansi, donc on calcule + // quand même la def cumulée + save_resul.eps_cumulBB = save_resul.eps_cumulBB_t + delta_epsBB; + break; + } + + case 1: // cas d'une dérivée de Lie deux fois contravariante + {sigHH_n = sigHH_nn; + sigBH_n = gijBB * sigHH_n; + sigBB_n = sigBH_n * gijBB; + // pour la déformation cumulée associée, + Tenseur1HH eps_interHH_n = gijHH_t * save_resul.eps_cumulBB_t * gijHH_t; + Tenseur1BB epsBB_n = gijBB * eps_interHH_n * gijBB ; + save_resul.eps_cumulBB = delta_epsBB + epsBB_n; + break; + } + }; + + // opération de transmission de la métrique + 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; + + // cas de la thermo dépendance, on calcul les grandeurs + if (f_temperature != NULL) + { f = f_temperature->Valeur(*temperature);} + else if (f_nD != NULL) + // là il faut calculer la fonction nD + { // on utilise la méthode générique de loi abstraite + list list_save; // inter pour l'appel de la fonction + list_save.push_back(saveResul); + Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee + (f_nD,1 // une seule valeur attendue en retour + ,ex_impli,ex_expli_tdt,ex_expli + ,NULL + ,NULL + ,&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 + List_io & li_enu_scal = f_nD->Li_enu_etendu_scalaire(); + List_io & li_quelc = f_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 + Tableau val_ddl_enum(Valeur_multi_interpoler_ou_calculer + (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,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,ex_impli,ex_expli_tdt,ex_expli,NULL); + // calcul de la valeur et retour dans tab_ret + Tableau & tab_val = f_nD->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL); + #ifdef MISE_AU_POINT + if (tab_val.Taille() != 1) + { cout << "\nErreur : la fonction nD relative a la fonction f " + << " doit calculer un scalaire or le tableau de retour est de taille " + << tab_val.Taille() << " ce n'est pas normal !\n"; + cout << " Hypo_hooke1D::Calcul_SigmaHH\n"; + Sortie(1); + }; + #endif + */ + // on récupère le premier élément du tableau uniquement + f = tab_val(1); + }; + + if (!compress_thermophysique) // sinon kc est mis à jour avec la méthode CalculGrandeurTravail(..) + {if (Kc_temperature != NULL) + {Kc = Kc_temperature->Valeur(*temperature);} + else if (Kc_nD != NULL) + // là il faut calculer la fonction nD + { // on utilise la méthode générique de loi abstraite + list list_save; // inter pour l'appel de la fonction + list_save.push_back(saveResul); + Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee + (Kc_nD,1 // une seule valeur attendue en retour + ,ex_impli,ex_expli_tdt,ex_expli + ,NULL + ,NULL + ,&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 + List_io & li_enu_scal = Kc_nD->Li_enu_etendu_scalaire(); + List_io & li_quelc = Kc_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 + Tableau val_ddl_enum(Valeur_multi_interpoler_ou_calculer + (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,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,ex_impli,ex_expli_tdt,ex_expli,NULL); + // calcul de la valeur et retour dans tab_ret + Tableau & tab_val = Kc_nD->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL); + #ifdef MISE_AU_POINT + if (tab_val.Taille() != 1) + { cout << "\nErreur : la fonction nD relative a la fonction Kc " + << " doit calculer un scalaire or le tableau de retour est de taille " + << tab_val.Taille() << " ce n'est pas normal !\n"; + cout << " Hypo_hooke1D::Calcul_SigmaHH\n"; + Sortie(1); + }; + #endif + */ + // on récupère le premier élément du tableau uniquement + Kc = tab_val(1); + }; + }; + // cas d'une non-linéarité multiplicative en fonction de la déformation + double Kc_use=Kc;// pour ne pas changer les valeurs à chaque passage !! + double f_use=f;// " + if ((Kc_IIeps != NULL)||(f_IIeps != NULL)) + { Tenseur1BH epsBH = epsBB * gijHH;double II_eps=epsBH.II(); + if (Kc_IIeps != NULL) + { double coef1 = Kc_IIeps->Valeur(II_eps); + Kc_use *= coef1; + }; + if (f_IIeps != NULL) + { double coef2 = f_IIeps->Valeur(II_eps); + f_use *= coef2; + }; + }; + + // sauvegarde des paramètres matériau + save_resul.Kc = Kc_use; + save_resul.f = f_use; + + // calcul de la contrainte unidirectionnelle + Tenseur1BH delta_epsBH = delta_epsBB * gijHH; + Tenseur1BH delta_sigmaBH(delta_epsBH*f_use); + Tenseur1BH sigBH = sigBH_n + delta_sigmaBH; + sigHH = gijHH * sigBH; + + // on applique éventuellement une restriction sur la contrainte + // =0 -> pas de restriction + // = -1 : traction uniquement autorisée, la compression est mise à 0 + // = 1 : compression uniquement autorisée, la traction est mise à 0 + switch (restriction_traction_compression) + { case 0: break; // cas sans restriction + case -1: // traction uniquement autorisée, la compression est mise à 0 + {if (sigHH(1,1)<0.) + sigHH.Coor(1,1) = 0.; + break; + } + case 1: // compression uniquement autorisée, la traction est mise à 0 + {if (sigHH(1,1)>0.) + sigHH.Coor(1,1) = 0.; + break; + } + default: + cout << "\n erreur Hypo_hooke1D::Calcul_SigmaHH " + << " restriction_traction_compression= "<< restriction_traction_compression + << " ce qui n'est pas exploitable !! "; + Sortie(1); + }; + + // traitement des énergies + // on incrémente l'énergie élastique + energ.ChangeEnergieElastique(energ_t.EnergieElastique() + 0.5 * ((sigHH + sigHH_nn) && delta_epsBB)); + + // récup de la compressibilité (-p_point = compress * I_D, S_point = cisaille * D_barre) + module_compressibilite = Kc_use/3.; // en fait il s'agit ici de la compressibilité tangente + + // en faisant l'analogie : K == Kc/3. = E/3/(1-2nu) et f==E + // on obtient: nu = 1/2 (1-f/Kc) et G = Kc*f/(3*Kc - f) + module_cisaillement = Kc_use * f_use/(3. * Kc_use - f_use); + + // calcul des déformations transversales + // \varepsilon_i^{.i}~(t+\Delta t) = \varepsilon_i^{.i} (t) + + //\frac{1}{2} \left ( \frac{{\Delta_t^{t+\Delta t} \sigma_1^{.1}}}{3~K_c - \Delta_t^{t+\Delta t} \varepsilon_1^{.1} } \right ) + save_resul.eps22 = save_resul.eps22_t + 0.5 *(delta_sigmaBH(1,1)/(Kc_use) - delta_epsBH(1,1) ); + save_resul.eps33 = save_resul.eps22; + + +// ******** arrête ici ******** +// il faut également traiter les déformations eps22 et 33 qui sont stockées et qui dépendent du type de +// transport (contrairement à ce qui a été défini dans Hypo_hooke1D::HsurH0 qu'il faut revoir +// du coup il faut également revoir la def du type de def, ou plutôt sa cohérence avec le type de transport +// ici, c'est le type de transport qui va définir la def utilisée et non l'utilisateur +// donc au moment de la lecture il faut mettre un message d'erreur si l'utilisateur choisit une def différente du type de transport -> donc modifier : Hypo_hooke1D::LectureDonneesParticulieres +// ... à faire + + + LibereTenseur(); + }; + + // calcul des contraintes a t+dt et de ses variations +void Hypo_hooke1D::Calcul_DsigmaHH_tdt (TenseurHH& sigHH_t,TenseurBB& ,DdlElement & tab_ddl + ,BaseB& giB_t,TenseurBB & gijBB_t_,TenseurHH & gijHH_t_,BaseB& giB_tdt,Tableau & d_giB_tdt + ,BaseH& giH_tdt,Tableau & d_giH_tdt, + TenseurBB & epsBB_tdt,Tableau & d_epsBB,TenseurBB & delta_epsBB_,TenseurBB & gijBB_tdt, + TenseurHH & gijHH_tdt,Tableau & d_gijBB_tdt, + Tableau & d_gijHH_tdt,double& jacobien_0,double& jacobien, + Vecteur& d_jacobien_tdt,TenseurHH& sigHH_,Tableau & d_sigHH + ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement + ,const Met_abstraite::Impli& ex) +{ + #ifdef MISE_AU_POINT + if (delta_epsBB_.Dimension() != 1) + { cout << "\nErreur : la dimension devrait etre 1 !\n"; + cout << " Hypo_hooke1D::Calcul_DsigmaHH_tdt\n"; + Sortie(1); + }; + if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille()) + { cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt !\n"; + cout << " Hypo_hooke1D::Calcul_DsigmaHH_tdt\n"; + Sortie(1); + }; + #endif + const Tenseur1BB & epsBB = *((Tenseur1BB*) &epsBB_tdt); // passage en dim 1 + const Tenseur1BB & delta_epsBB = *((Tenseur1BB*) &delta_epsBB_); // passage en dim 1 + const Tenseur1HH & gijHH = *((Tenseur1HH*) &gijHH_tdt); // " " " " + const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_tdt); // " " " " + Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_); // " " " " + Tenseur1HH & sigHH_nn = *((Tenseur1HH*) &sigHH_t); // " " " " + + const Tenseur1HH & gijHH_t = *((Tenseur3HH*) &gijHH_t_); // " " " " + const Tenseur1BB & gijBB_t = *((Tenseur3BB*) &gijBB_t_); // " " " " + + SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveResul); + // tenseur intermédiaires utilisées selon les cas (par forcément dans tous les cas !!) + Tenseur1BH sigBH_n;Tenseur1HH sigHH_n; Tenseur1BB sigBB_n;Tenseur1BB sig_interBB_n; + switch (type_derive) + {case -1: // cas d'une dérivée de jauman: 1/2 deux fois covariant + deux fois contra + {sig_interBB_n = gijBB_t * sigHH_nn * gijBB_t; + sigBH_n = 0.5*( sig_interBB_n * gijHH + gijBB * sigHH_nn) ; + sigHH_n = gijHH * sigBH_n ; + sigBB_n = sigBH_n * gijBB; + // pour la déformation cumulée associée, + Tenseur1HH eps_interHH_n = gijHH_t * save_resul.eps_cumulBB_t * gijHH_t; + Tenseur1BB epsBB_n = 0.5*( gijBB * eps_interHH_n * gijBB + save_resul.eps_cumulBB_t); + save_resul.eps_cumulBB = delta_epsBB + epsBB_n; + break; + } + case 0: // cas d'une dérivée de Lie deux fois covariantes + {sigBB_n = gijBB_t * sigHH_nn * gijBB_t; + sigBH_n = sigBB_n * gijHH ; + sigHH_n = gijHH * sigBH_n ; + // pour la déformation cumulée associée, on sait que cela donne directement la def d'Almansi + // mais la def passée en paramètre pourrait ne pas être d'Almansi, donc on calcule + // quand même la def cumulée + save_resul.eps_cumulBB = save_resul.eps_cumulBB_t + delta_epsBB; + break; + } + case 1: // cas d'une dérivée de Lie deux fois contravariantes + {sigHH_n = sigHH_nn; + sigBH_n = gijBB * sigHH_n; + sigBB_n = sigBH_n * gijBB; + // pour la déformation cumulée associée, + Tenseur1HH eps_interHH_n = gijHH_t * save_resul.eps_cumulBB_t * gijHH_t; + Tenseur1BB epsBB_n = gijBB * eps_interHH_n * gijBB ; + save_resul.eps_cumulBB = delta_epsBB + epsBB_n; + break; + } + }; + + // opération de transmission de la métrique + 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; + + // cas de la thermo dépendance, on calcul les grandeurs + if (f_temperature != NULL) + { f = f_temperature->Valeur(*temperature);} + else if (f_nD != NULL) + // là il faut calculer la fonction nD + { // on utilise la méthode générique de loi abstraite + list list_save; // inter pour l'appel de la fonction + list_save.push_back(saveResul); + Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee + (f_nD,1 // une seule valeur attendue en retour + ,ex_impli,ex_expli_tdt,ex_expli + ,NULL + ,NULL + ,&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 + List_io & li_enu_scal = f_nD->Li_enu_etendu_scalaire(); + List_io & li_quelc = f_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 + Tableau val_ddl_enum(Valeur_multi_interpoler_ou_calculer + (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,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,ex_impli,ex_expli_tdt,ex_expli,NULL); + // calcul de la valeur et retour dans tab_ret + Tableau & tab_val = f_nD->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL); + #ifdef MISE_AU_POINT + if (tab_val.Taille() != 1) + { cout << "\nErreur : la fonction nD relative a la fonction f " + << " doit calculer un scalaire or le tableau de retour est de taille " + << tab_val.Taille() << " ce n'est pas normal !\n"; + cout << " Hypo_hooke1D::Calcul_DsigmaHH_tdt \n"; + Sortie(1); + }; + #endif + */ + // on récupère le premier élément du tableau uniquement + f = tab_val(1); + }; + + if (!compress_thermophysique) // sinon kc est mis à jour avec la méthode CalculGrandeurTravail(..) + {if (Kc_temperature != NULL) + {Kc = Kc_temperature->Valeur(*temperature);} + else if (Kc_nD != NULL) + // là il faut calculer la fonction nD + { // on utilise la méthode générique de loi abstraite + list list_save; // inter pour l'appel de la fonction + list_save.push_back(saveResul); + Tableau & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee + (Kc_nD,1 // une seule valeur attendue en retour + ,ex_impli,ex_expli_tdt,ex_expli + ,NULL + ,NULL + ,&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 + List_io & li_enu_scal = Kc_nD->Li_enu_etendu_scalaire(); + List_io & li_quelc = Kc_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 + Tableau val_ddl_enum(Valeur_multi_interpoler_ou_calculer + (absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,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,ex_impli,ex_expli_tdt,ex_expli,NULL); + // calcul de la valeur et retour dans tab_ret + Tableau & tab_val = Kc_nD->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL); + #ifdef MISE_AU_POINT + if (tab_val.Taille() != 1) + { cout << "\nErreur : la fonction nD relative a la fonction Kc " + << " doit calculer un scalaire or le tableau de retour est de taille " + << tab_val.Taille() << " ce n'est pas normal !\n"; + cout << " Hypo_hooke1D::Calcul_SigmaHH\n"; + Sortie(1); + }; + #endif + */ + // on récupère le premier élément du tableau uniquement + Kc = tab_val(1); + }; + }; + + // cas d'une non-linéarité en fonction de la déformation + Tenseur1BH epsBH; Courbe1D::ValDer valder_Kc,valder_f; + double Kc_base=Kc;double Kc_use=Kc;// pour ne pas changer les valeurs à chaque passage !! + double f_base=f;double f_use=f;// idem + if ((Kc_IIeps != NULL)||(f_IIeps != NULL)) + { epsBH = epsBB * gijHH;double II_eps=epsBH.II(); + if (Kc_IIeps != NULL) + { valder_Kc = Kc_IIeps->Valeur_Et_derivee(II_eps); + Kc_use *= valder_Kc.valeur; + }; + if (f_IIeps != NULL) + { valder_f = f_IIeps->Valeur_Et_derivee(II_eps); + f_use *= valder_f.valeur; + }; + }; + + + // sauvegarde des paramètres matériau + save_resul.Kc = Kc_use; + save_resul.f = f_use; + + // calcul de la contrainte unidirectionnelle + Tenseur1BH delta_epsBH = delta_epsBB * gijHH; + Tenseur1BH delta_sigmaBH(delta_epsBH*f_use); + Tenseur1BH sigBH = sigBH_n + delta_sigmaBH; + sigHH = gijHH * sigBH; + + + // on applique éventuellement une restriction sur la contrainte + // =0 -> pas de restriction + // = -1 : traction uniquement autorisée, la compression est mise à 0 + // = 1 : compression uniquement autorisée, la traction est mise à 0 + int restriction_effectuee=0; + switch (restriction_traction_compression) + { case 0: break; // cas sans restriction + case -1: // traction uniquement autorisée, la compression est mise à 0 + {if (sigHH(1,1)<0.) + sigHH.Coor(1,1) = 0.; + restriction_effectuee = -1; + break; + } + case 1: // compression uniquement autorisée, la traction est mise à 0 + {if (sigHH(1,1)>0.) + sigHH.Coor(1,1) = 0.; + restriction_effectuee = 1; + break; + } + default: + cout << "\n erreur Hypo_hooke1D::Calcul_DsigmaHH_tdt " + << " restriction_traction_compression= "<< restriction_traction_compression + << " ce qui n'est pas exploitable !! "; + Sortie(1); + }; + + // traitement des énergies + // on incrémente l'énergie élastique + energ.ChangeEnergieElastique(energ_t.EnergieElastique() + 0.5 * ((sigHH + sigHH_nn) && delta_epsBB)); + + // récup de la compressibilité (-p_point = compress * I_D, S_point = cisaille * D_barre) + module_compressibilite = Kc_use/3.; // en fait il s'agit ici de la compressibilité tangente + + // en faisant l'analogie : K == Kc/3. = E/3/(1-2nu) et f==E + // on obtient: nu = 1/2 (1-f/Kc) et G = Kc*f/(3*Kc - f) + module_cisaillement = Kc_use * f_use/(3. * Kc_use - f_use); + + // calcul des déformations transversales + // \varepsilon_i^{.i}~(t+\Delta t) = \varepsilon_i^{.i} (t) + + //\frac{1}{2} \left ( \frac{{\Delta_t^{t+\Delta t} \sigma_1^{.1}}}{3~K_c - \Delta_t^{t+\Delta t} \varepsilon_1^{.1} } \right ) + save_resul.eps22 = save_resul.eps22_t + 0.5 *(delta_sigmaBH(1,1)/(Kc_use) - delta_epsBH(1,1) ); + save_resul.eps33 = save_resul.eps22; + +////--- debug +//cout << "\n Hypo_hooke1D::Calcul_DsigmaHH_tdt"; +//cout << "\n save_resul.eps22= "<Valeur(*temperature); + if (!compress_thermophysique) // sinon kc est mis à jour avec la méthode CalculGrandeurTravail(..) + {if (Kc_temperature != NULL) Kc = Kc_temperature->Valeur(*temperature);}; + // cas d'une non-linéarité en fonction de la déformation + Tenseur3BH DepsBH; Tenseur3BB Deps_barre_BB;Tenseur3HH Deps_barre_HH; // init + Tenseur3BH Deps_barre_BH; + double IDeps=0.; // " + Tenseur3BH epsBH; + + Courbe1D::ValDer valder_Kc,valder_f; + double Kc_base=Kc;double Kc_use=Kc;// pour ne pas changer les valeurs à chaque passage !! + double f_base=f;double f_use=f;// idem + static const double untier=1./3.; + + + if (en_base_orthonormee) + {DepsBH = DepsBB.MonteDernierIndice(); + IDeps=DepsBH.Trace(); + Deps_barre_BH = DepsBH - (untier * IDeps) * IdBH3; + Deps_barre_BB = Deps_barre_BH * IdBB3; + Deps_barre_HH = IdHH3 * Deps_barre_BH ; + epsBH = epsBB.MonteDernierIndice(); // deformation en mixte + } + else + {DepsBH = DepsBB * gijHH; + IDeps=DepsBH.Trace(); + Deps_barre_BH = DepsBH - (untier * IDeps) * IdBH3; + Deps_barre_BB = Deps_barre_BH * gijBB; + Deps_barre_HH = gijHH * Deps_barre_BH ; + epsBH = epsBB * gijHH; // deformation en mixte + }; + + if ((Kc_IIeps != NULL)||(f_IIeps != NULL)) + { + double Ieps = epsBH.Trace(); + Tenseur3BH eps_barre_BH = epsBH - (untier * Ieps) * IdBH3; + + + double II_eps=epsBH.II(); + if (Kc_IIeps != NULL) + { valder_Kc = Kc_IIeps->Valeur_Et_derivee(II_eps); + Kc_use *= valder_Kc.valeur; + }; + if (f_IIeps != NULL) + { valder_f = f_IIeps->Valeur_Et_derivee(II_eps); + f_use *= valder_f.valeur; + }; + }; + + // recup de l'incrément de temps + double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); + // le calcul de la contrainte correspond à l'intégration d'une équation différencielle + // delta_eps_barre=1/f delta_sigma/dt pour la partie déviatoire + // et pour la partie sphérique + // I_D = 1/(Kc) I_Sig_point + double Isig_n = sigBH_n.Trace(); + Tenseur3BH SBH_n = sigBH_n - (untier * Isig_n) * IdBH3; + + 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 + { // non 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; + }; + // cas de la partie sphérique scalaire + double Isigma = Isig_n + deltat * Kc_use * IDeps; + + // cas de la partie déviatorique + +// if (en_base_orthonormee) +// { switch (cas_calcul) +// { case 0: // calcul normal (tous les termes) +// { // la partie sphérique est déjà calculé, cas de la partie déviatorique +// Tenseur3HH SHH = IdHH3 * (SBH_n + (deltat * f_use) * Deps_barre_BH); +// sigHH = SHH + (untier * Isigma) * IdHH3; +// break; +// } +// case 1: // calcul de la partie déviatorique seule +// { sigHH = IdHH3 * (SBH_n + (deltat * f_use) * Deps_barre_BH); +// break; +// } +// case 2: // calcul de la partie sphérique seule +// { sigHH = (untier * Isigma) * IdHH3; +// break; +// } +// default: +// { cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! " +// << "\n Hypo_hooke1D::Calcul_dsigma_deps (.... "; +// Sortie(1); +// } +// }; +// } +// else +// { switch (cas_calcul) +// { case 0: // calcul normal (tous les termes) +// { // la partie sphérique est déjà calculé, cas de la partie déviatorique +// Tenseur3HH SHH = gijHH * (SBH_n + deltat * f_use * Deps_barre_BH); +// sigHH = SHH + (untier * Isigma) * gijHH ; +// break; +// } +// case 1: // calcul de la partie déviatorique seule +// { sigHH = gijHH * (SBH_n + deltat * f_use * Deps_barre_BH); +// break; +// } +// case 2: // calcul de la partie sphérique seule +// { sigHH = (untier * Isigma) * gijHH ; // *IdBH3; +// break; +// } +// default: +// { cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! " +// << "\n Hypo_hooke1D::Calcul_dsigma_deps (.... "; +// Sortie(1); +// } +// }; +// }; + + +// +// +// +////--- debug +//cout << "\n Hypo_hooke1D::Calcul_dsigma_deps"; +//cout << "\n IDeps= "<Monte4Indices(); + +// switch (cas_calcul) +// { case 0: // calcul normal (tous les termes) +// { d_sigma_depsHHHH = f_use * d_deltat_Deps_barre_HHHH; +// if (f_IIeps != NULL) +// d_sigma_depsHHHH += ((Tenseur3BBBB*) &((f_base * valder_f.derivee * 2. * deltat) +// * Tenseur3BBBB::Prod_tensoriel((Deps_barre_BH * IdBB3),epsBB)))->Monte4Indices(); +// d_sigma_depsHHHH += (untier * Kc) * IdHHHH3; +// if (Kc_IIeps != NULL) +// d_sigma_depsHHHH += ((Tenseur3BBBB*) &((untier * Kc_base * valder_Kc.derivee * 2. ) +// * Tenseur3BBBB::Prod_tensoriel(IdBB3, epsBB)))->Monte4Indices(); +// break; +// } +// case 1: // calcul de la partie déviatorique seule +// { d_sigma_depsHHHH = f_use * d_deltat_Deps_barre_HHHH; +// if (f_IIeps != NULL) +// d_sigma_depsHHHH += ((Tenseur3BBBB*) &((f_base * valder_f.derivee * 2. * deltat) +// * Tenseur3BBBB::Prod_tensoriel((Deps_barre_BH * IdBB3),epsBB)))->Monte4Indices(); +// break; +// } +// case 2: // calcul de la partie sphérique seule +// { d_sigma_depsHHHH = (untier * Kc) * IdHHHH3; +// if (Kc_IIeps != NULL) +// d_sigma_depsHHHH += ((Tenseur3BBBB*) &((untier * Kc_base * valder_Kc.derivee * 2. ) +// * Tenseur3BBBB::Prod_tensoriel(IdBB3, epsBB)))->Monte4Indices(); +// break; +// } +// default: +// { cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! " +// << "\n Hypo_hooke1D::Calcul_DsigmaHH_tdt (.... "; +// Sortie(1); +// } +// }; + } + else // sinon cas où les bases sont curvilignes + { // calcul de variables intermédiaires + I_x_I_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH,gijHH); + I_xbarre_I_HHHH=Tenseur3HHHH::Prod_tensoriel_barre(gijHH,gijHH); + Tenseur3HH epsHH(gijHH * epsBH); + I_x_eps_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH,epsHH); + Tenseur3HH Deps_HH(gijHH * DepsBH); + I_x_D_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH,Deps_HH); + I_xbarre_D_HHHH=Tenseur3HHHH::Prod_tensoriel_barre(gijHH,Deps_HH); + + // variation de sigma_n et de la trace en fonction du transport + switch (type_derive) + { case -1: // // cas d'une dérivée de jauman: 1/2 deux fois covariant + deux fois contra + {// pour info sigBH_n = 0.5*(gijBB * sigHH_n + (gijBB_t * sigHH_n * gijBB_t) * gijHH) + d_sig_t_HHHH = Tenseur3HHHH::Prod_tensoriel_barre(gijHH,sigHH_n); + d_spherique_sig_t_HHHH = (1. /6.) * + (Tenseur3HHHH::Prod_tensoriel(gijHH,sigHH_n)+ Isig_n * I_xbarre_I_HHHH); + break; + } + case 0: // cas d'une dérivée de Lie deux fois covariantes + {// pour info sigBH_n = (gijBB_t * sigHH_n * gijBB_t) * gijHH + d_sig_t_HHHH = 2.*Tenseur3HHHH::Prod_tensoriel_barre(gijHH,sigHH_n); + d_spherique_sig_t_HHHH = untier * + (Tenseur3HHHH::Prod_tensoriel(gijHH,sigHH_n)+ Isig_n * I_xbarre_I_HHHH); + break; + } + case 1: // cas d'une dérivée de Lie deux fois contravariantes + {// pour info sigBH_n = gijBB * sigHH_n, ici il n'y a aucune conséquence + // sur sigma(t), par contre il y en a une sur la trace + d_sig_t_HHHH.Inita(0.); + d_spherique_sig_t_HHHH = untier * + (Tenseur3HHHH::Prod_tensoriel(gijHH,sigHH_n)+ Isig_n * I_xbarre_I_HHHH); + break; + } + }; + + // on choisit entre les différents cas + double e01,e02;e01=e02=0.; // les coefficients pour la partie transportée + double b01,b02,b03,b04,b05;b01=b02=b03=b04=b05=0.; // les coefficients pour le reste de l'équation + // l'expression générique est: + // d_sigma_depsHHHH = b01 * I_x_I_HHHH + b02 * I_x_D_HHHH + b03 * I_xbarre_I_HHHH + // + b04 * I_xbarre_D_HHHH + // + e01 * d_sig_t_HHHH + e02 * d_spherique_sig_t_HHHH; + // en fonction des coefficients nuls on simplifie + +// switch (cas_calcul) +// { case 0: // calcul normal (tous les termes) // cas complet +// { b01= (Kc_use - f_use) * untier; +// b02= -2.* (Kc_use - f_use) * untier * deltat; +// b03= -2. * (Kc_use - f_use) * untier * deltat * IDeps + f_use; +// b04= -4. * f_use * deltat; +// e01= f_use; +// e02= untier; +// d_sigma_depsHHHH = b01 * I_x_I_HHHH + b02 * I_x_D_HHHH +// + b03 * I_xbarre_I_HHHH + b04 * I_xbarre_D_HHHH +// + e01 * d_sig_t_HHHH + e02 * d_spherique_sig_t_HHHH; +// break; +// } +// case 1: // calcul de la partie déviatorique seule +// { b01= - f_use * untier; +// b02= 2.* f_use * untier * deltat; +// b03= 2. * f_use * untier * deltat * IDeps + f_use; +// b04= -4. * f_use * deltat; +// e01= f_use; +// e02= untier; +// d_sigma_depsHHHH = b01 * I_x_I_HHHH + b02 * I_x_D_HHHH +// + b03 * I_xbarre_I_HHHH + b04 * I_xbarre_D_HHHH +// + e01 * d_sig_t_HHHH + e02 * d_spherique_sig_t_HHHH; +// break; +// } +// case 2: // calcul de la partie sphérique seule +// { b01= Kc_use * untier; +// b02= -2.*Kc_use * untier * deltat; +// b03= -2.*Kc_use * untier * deltat * IDeps; +// b04= 0.; +// e01= 0.; +// e02= untier; +// d_sigma_depsHHHH = b01 * I_x_I_HHHH + b02 * I_x_D_HHHH +// + b03 * I_xbarre_I_HHHH +// + e02 * d_spherique_sig_t_HHHH; +// break; +// } +// default: +// { cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! " +// << "\n Hypo_hooke1D::Calcul_dsigma_deps (.... "; +// Sortie(1); +// } +// }; + }; + + // traitement des énergies + // on incrémente l'énergie élastique + energ.ChangeEnergieElastique(energ_t.EnergieElastique() + 0.5 * deltat * ((sigHH + sigHH_nn) && DepsBB)); + + // récup de la compressibilité (-p_point = compress * I_D, S_point = cisaille * D_barre) + module_compressibilite = Kc_use/3.; + module_cisaillement = 0.5 * f_use; + + LibereTenseur(); + LibereTenseurQ(); + }; + diff --git a/comportement/Hypo_elastique/Hypo_hooke1D.h b/comportement/Hypo_elastique/Hypo_hooke1D.h new file mode 100755 index 0000000..2699c42 --- /dev/null +++ b/comportement/Hypo_elastique/Hypo_hooke1D.h @@ -0,0 +1,367 @@ + // FICHIER : Hypo1D.h + // CLASSE : Hypo1D + + + // 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) . + // + // Herezh++ is distributed under GPL 3 license ou ultérieure. + // + // Copyright (C) 1997-2022 Université Bretagne Sud (France) + // AUTHOR : Gérard Rio + // E-MAIL : gerardrio56@free.fr + // + // This program is free software: you can redistribute it and/or modify + // it under the terms of the GNU General Public License as published by + // the Free Software Foundation, either version 3 of the License, + // or (at your option) any later version. + // + // This program is distributed in the hope that it will be useful, + // but WITHOUT ANY WARRANTY; without even the implied warranty + // of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + // See the GNU General Public License for more details. + // + // You should have received a copy of the GNU General Public License + // along with this program. If not, see . + // + // For more information, please consult: . + +/************************************************************************ + * DATE: 18/07/2020 * + * $ * + * AUTEUR: G RIO (mailto:gerardrio56@free.fr) * + * $ * + * PROJET: Herezh++ * + * $ * + ************************************************************************ + * BUT: La classe Hypo1D definit une loi 1D hypo-élastique * + * qui sous forme intégrée peut dans certain cas être * + * équivalente à hooke. * + * On a donc : * + * sigma_point = f(..) * D * + * $ * + * '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * * + * VERIFICATION: * + * * + * ! date ! auteur ! but ! * + * ------------------------------------------------------------ * + * ! ! ! ! * + * $ * + * '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * + * MODIFICATIONS: * + * ! date ! auteur ! but ! * + * ------------------------------------------------------------ * + * $ * + ************************************************************************/ + + +#ifndef HYPO_ELAS1D_H +#define HYPO_ELAS1D_H + + +#include "Loi_comp_abstraite.h" + +/** @defgroup Les_lois_hypoelastiques +* +* BUT: groupe des lois hypoélastiques +* +* +* \author Gérard Rio +* \version 1.0 +* \date 28/06/2004 +* \brief Définition des lois hypoélastiques +* +*/ + +/// @addtogroup Les_lois_hypoelastiques +/// @{ +/// + + +class Hypo_hooke1D : public Loi_comp_abstraite +{ + + + public : + + + // CONSTRUCTEURS : + + // Constructeur par defaut + Hypo_hooke1D (); + + + // Constructeur de copie + Hypo_hooke1D (const Hypo_hooke1D& loi) ; + + // DESTRUCTEUR : + + ~Hypo_hooke1D (); + + + // initialise les donnees particulieres a l'elements + // de matiere traite ( c-a-dire au pt calcule) + // Il y a creation d'une instance de SaveResul particuliere + // a la loi concernee + // la SaveResul classe est remplie par les instances heritantes + // le pointeur de SaveResul est sauvegarde au niveau de l'element + // c'a-d que les info particulieres au point considere sont stocke + // au niveau de l'element et non de la loi. + class SaveResulLoi_Hypo1D: public SaveResul + { public : + SaveResulLoi_Hypo1D(); // constructeur par défaut à ne pas utiliser + // le constructeur courant + SaveResulLoi_Hypo1D (SaveResul* ); + // de copie + SaveResulLoi_Hypo1D(const SaveResulLoi_Hypo1D& sav): // de copie + Kc(sav.Kc),Kc_t(sav.Kc_t),f(sav.f),f_t(sav.f_t) + ,eps22(sav.eps22),eps22_t(sav.eps22_t),eps33(sav.eps33),eps33_t(sav.eps33_t) + ,eps_cumulBB(sav.eps_cumulBB),eps_cumulBB_t(sav.eps_cumulBB_t) + {}; + + virtual ~SaveResulLoi_Hypo1D() {}; // destructeur + // définition d'une nouvelle instance identique + // appelle du constructeur via new + SaveResul * Nevez_SaveResul() const{return (new SaveResulLoi_Hypo1D(*this));}; + // affectation + virtual SaveResul & operator = ( const SaveResul & a) + { SaveResulLoi_Hypo1D& sav = *((SaveResulLoi_Hypo1D*) &a); + Kc=sav.Kc;Kc_t=sav.Kc_t;f=sav.f;f_t=sav.f_t; + eps22=sav.eps22;eps22_t=sav.eps22_t;eps33=sav.eps33;eps33_t=sav.eps33_t; + eps_cumulBB=sav.eps_cumulBB;eps_cumulBB_t=sav.eps_cumulBB_t; + return *this; + }; + //============= lecture écriture dans base info ========== + // 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 Lecture_base_info (ifstream& ent,const int cas); + + // cas donne le niveau de sauvegarde + // = 1 : on sauvegarde tout + // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) + void Ecriture_base_info(ofstream& sort,const int cas); + + // mise à jour des informations transitoires + void TdtversT() + {Kc_t = Kc; f_t=f; eps22_t=eps22;eps33_t=eps33; + eps_cumulBB_t = eps_cumulBB; + }; + void TversTdt() + {Kc = Kc_t; f=f_t;eps22=eps22_t;eps33=eps33_t; + eps_cumulBB = eps_cumulBB_t; + }; + + // affichage à l'écran des infos + void Affiche() const + { cout <<"\n Kc= "<< Kc << " f= " << f + << " eps22= "<< eps22 << " eps33= "<< eps33 + << " eps_cumulBB= " << eps_cumulBB + << " "; + }; + + //changement de base de toutes les grandeurs internes tensorielles stockées + // beta(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB + // gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne + // gpH(i) = gamma(i,j) * gH(j) + virtual void ChBase_des_grandeurs(const Mat_pleine& beta,const Mat_pleine& gamma) ; + + // procedure permettant de completer éventuellement les données particulières + // de la loi stockées + // au niveau du point d'intégration par exemple: exemple: un repère d'anisotropie + // completer est appelé apres sa creation avec les donnees du bloc transmis + // peut etre appeler plusieurs fois + virtual SaveResul* Complete_SaveResul(const BlocGen & bloc, const Tableau & tab_coor + ,const Loi_comp_abstraite* loi) {return NULL;}; + + + //------------------------------------------------------------------- + // données + //------------------------------------------------------------------- + double Kc,Kc_t; // les paramètres matériaux réellement utilisés + double f,f_t; + double eps33,eps22; // déformations transversale courantes + double eps33_t,eps22_t; // les dernières enregistrées + + Tenseur1BB eps_cumulBB,eps_cumulBB_t; // déformation cumulée associée à la loi + + }; + + SaveResul * New_et_Initialise(); + + friend class SaveResulLoi_Hypo1D; + + // Lecture des donnees de la classe sur fichier + void LectureDonneesParticulieres (UtilLecture * ,LesCourbes1D& lesCourbes1D + ,LesFonctions_nD& lesFonctionsnD); + // affichage de la loi + void Affiche() const ; + // test si la loi est complete + // = 1 tout est ok, =0 loi incomplete + int TestComplet(); + + //----- lecture écriture de restart ----- + // 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 Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D + ,LesFonctions_nD& lesFonctionsnD); + // cas donne le niveau de sauvegarde + // = 1 : on sauvegarde tout + // = 2 : on sauvegarde uniquement les données variables (supposées comme telles) + void Ecriture_base_info_loi(ofstream& sort,const int cas); + + // récupération des grandeurs particulière (hors ddl ) + // correspondant à liTQ + // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière + void Grandeur_particuliere + (bool absolue,List_io& ,Loi_comp_abstraite::SaveResul * ,list& decal) const; + // récupération de la liste de tous les grandeurs particulières + // ces grandeurs sont ajoutées à la liste passées en paramètres + // absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière + void ListeGrandeurs_particulieres(bool absolue,List_io& ) const; + + // calcul d'un module d'young équivalent à la loi + double Module_young_equivalent(Enum_dure temps,const Deformation & ,SaveResul * saveResul); + + // récupération d'un module de compressibilité équivalent à la loi, ceci pour un chargement nul + // il s'agit ici de la relation -pression = sigma_trace/3. = module de compressibilité * I_eps + // >>> en fait ici il s'agit du dernier module tangent calculé !! + double Module_compressibilite_equivalent(Enum_dure temps,const Deformation & def,SaveResul * saveResul); + + // 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 Eps33BH(SaveResul * saveDon) const + { SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + return save_resul.eps33; + }; + + // 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 Eps22BH(SaveResul * saveDon) const + { SaveResulLoi_Hypo1D & save_resul = *((SaveResulLoi_Hypo1D*) saveDon); + return save_resul.eps22; + }; + + // récupération de la variation relative d'épaisseur calculée: h/h0 + // 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 HsurH0(SaveResul * saveResul) const; + + // 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 +// pour l'instant en attente *** virtual double d_HsurH0(SaveResul * saveResul,Vecteur & d_hsurh0) const ; + + + // 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 BsurB0(SaveResul * saveResul) const ; + + // création d'une loi à l'identique et ramène un pointeur sur la loi créée + Loi_comp_abstraite* Nouvelle_loi_identique() const { return (new Hypo_hooke1D(*this)); }; + + // affichage et definition interactive des commandes particulières à chaques lois + void Info_commande_LoisDeComp(UtilLecture& lec); + // calcul de grandeurs de travail aux points d'intégration via la def et autres + // ici permet de récupérer la compressibilité + // fonction surchargée dans les classes dérivée si besoin est + virtual void CalculGrandeurTravail + (const PtIntegMecaInterne& ptintmeca + ,const Deformation & def,Enum_dure temps,const ThermoDonnee& dTP + ,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* exclure_dd_etend + ,const List_io* exclure_Q + ) + {if (compress_thermophysique) Kc = 3./dTP.Compressibilite(); }; + + protected : + // donnée de la loi + double f; // coef de proportionalité entre sig_1^{.1} et D_1^{.1} + Courbe1D* f_temperature; // courbe éventuelle d'évolution de f en fonction de la température + Courbe1D* f_IIeps; // courbe éventuelle d'évolution de f en fonction du deuxième invariant d'epsilon + Fonction_nD* f_nD; // fonction nD éventuelle pour f + + double Kc; // 3 * coefficient de compressibilité tangent + Courbe1D* Kc_temperature; // courbe éventuelle d'évolution de Kc en fonction de la température + Courbe1D* Kc_IIeps; // courbe éventuelle d'évolution de Kc en fonction du deuxième invariant d'epsilon + Fonction_nD * Kc_nD; // fonction nD éventuelle pour Kc + + bool compress_thermophysique; // indique si oui ou non la compressibilité est calculée par une loi + // thermophysique et donc + // récupéré par la fonction "CalculGrandeurTravail" + int type_derive; // type de dérivée objective utilisée pour sigma + + int restriction_traction_compression; // =0 -> pas de restriction + // = -1 : traction uniquement autorisée, la compression est mise à 0 + // = 1 : compression uniquement autorisée, la traction est mise à 0 + + // on introduit un certain nombre de tenseur du quatrième ordre, qui vont nous servir pour + // Calcul_dsigma_deps, dans le cas où on n'est pas en orthonormee + Tenseur3HHHH I_x_I_HHHH,I_xbarre_I_HHHH,I_x_eps_HHHH,I_x_D_HHHH,I_xbarre_D_HHHH,d_sig_t_HHHH; + Tenseur3HHHH d_spherique_sig_t_HHHH; + + // codage des METHODES VIRTUELLES protegees: + // calcul des contraintes a t+dt + // calcul des contraintes + void Calcul_SigmaHH (TenseurHH & sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl + ,TenseurBB & gijBB_t,TenseurHH & gijHH_t,BaseB& giB,BaseH& gi_H, TenseurBB & epsBB_ + ,TenseurBB & delta_epsBB_ + ,TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau & d_gijBB_ + ,double& jacobien_0,double& jacobien,TenseurHH & sigHH + ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement + ,const Met_abstraite::Expli_t_tdt& ex); + + // calcul des contraintes et de ses variations a t+dt + void Calcul_DsigmaHH_tdt (TenseurHH & sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl + ,BaseB& giB_t,TenseurBB & gijBB_t,TenseurHH & gijHH_t + ,BaseB& giB_tdt,Tableau & d_giB_tdt,BaseH& giH_tdt,Tableau & d_giH_tdt + ,TenseurBB & epsBB_tdt,Tableau & d_epsBB + ,TenseurBB & delta_epsBB,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt + ,Tableau & d_gijBB_tdt + ,Tableau & d_gijHH_tdt,double& jacobien_0,double& jacobien + ,Vecteur& d_jacobien_tdt,TenseurHH& sigHH,Tableau & d_sigHH + ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement + ,const Met_abstraite::Impli& ex); + + // 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 orthonormee + // 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 Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB& DepsBB + ,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien + ,TenseurHH& sigHH,TenseurHHHH& d_sigma_deps + ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement + ,const Met_abstraite::Umat_cont& ex) ; //= 0; + + +}; +/// @} // end of group + + +#endif + + + diff --git a/contact/ElContact.cc b/contact/ElContact.cc index d72c84a..3ef697d 100644 --- a/contact/ElContact.cc +++ b/contact/ElContact.cc @@ -705,48 +705,28 @@ bool ElContact::Contact( bool init) return false; }; +// **** il y a des choses à faire au niveau de la droite tangente et du point d'impact +// lorsque l'on a une ligne d'angle mort, contrairement au cas du point, le point d'intersection +// doit-être calculé, et c'est important car il va servir pour RecupInfo et ensuite pour +// le calcul de la raideur et du second membre... +// bref il faut revoir la chose: c'est un peu moins simple que dans le cas d'un point_G + // 1) === le cas des angles morts est particulier, on commence par le traiter if (elfront->Angle_mort()) { ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite - if (elfro.TypeFrontiere() == "FrontPointF") - {// s'il s'agit d'un point il n'y a pas d'existance de normale - // on remplace la notion de normale par le vecteur qui joint le noeud à la position - // du point frontière - // Coordonnee V = elfro.Ref()- noeud->Coord2(); - const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); - const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); - N = gr.ConteneurCoordonnee_const(); - - // on utilise la normale au noeud -// oui, mais ici ne sert donc on ne récupère pas -// const TypeQuelconque& tiq = noeud->Grandeur_quelconque(N_FRONT_t); -// const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); -// Coordonnee V = gr.ConteneurCoordonnee_const(); -// -// // l'intersection avec la frontière est systématiquement le point frontière -// Coordonnee M1 = elfro.Ref(); - // la seule chose que l'on peut faire pour valider le fait qu'il y a contact ou pas - // c'est de tester si le noeud est à l'intérieur de l'élément d'angle mort ou pas - // c'est plus long qu'un test de surface mais il n'y a pas vraiment d'autre solution - int cas = elfront->PtEI()->Interne_tdt(noeud->Coord2()); - if (Permet_affichage() > 5) - {cout << "\n -- ElContact::Contact: "; - cout << " Interne_tdt ? cas= " << cas <<", "; - this->Affiche(1); - }; - - // cout << "\n debug ElContact::Actualisation, angle mort: cas = " << cas ; - // elfront->Affiche(); - // cout << "\n----------------------------------"; - - return cas; // si != 0 le point est interne - } - else - { cout << "\n *** cas element Angle_mort actuellement non implante: , " - << elfro.TypeFrontiere() - << "\n ElContact::Contact(..."; - Sortie(1); - }; + int dim = noeud->Dimension(); + Plan pl(dim); Droite dr(dim); int indic; // def des variables de tangence + N = Calcul_Normale(dim,pl,dr,indic); + // la seule chose que l'on peut faire pour valider le fait qu'il y a contact ou pas + // c'est de tester si le noeud est à l'intérieur de l'élément d'angle mort ou pas + // c'est plus long qu'un test de surface mais il n'y a pas vraiment d'autre solution + int cas = elfront->PtEI()->Interne_tdt(noeud->Coord2()); + if (Permet_affichage() > 5) + {cout << "\n -- ElContact::Contact: "; + cout << " Interne_tdt ? cas= " << cas <<", "; + this->Affiche(1); + }; + return cas; // si != 0 le point est interne }; // 2) === suite pour le cas différent des angles morts @@ -1213,29 +1193,9 @@ Coordonnee ElContact::Trajectoire(int & test) int dim = noeud->Dimension(); // dimension des points // ----- on traite les cas particuliers des éléments d'angle mort - if (elfront->Angle_mort()) - { if (elfro.TypeFrontiere() == "FrontPointF") - {// s'il s'agit d'un point il n'y a pas d'existance de normale -// // on remplace la notion de normale par le vecteur qui joint le noeud à la position -// // du point frontière -// V = elfro.Ref()- noeud->Coord2(); - // on utilise la normale au noeud - const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); - const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); - N = gr.ConteneurCoordonnee_const(); - V=N; - test = type_trajectoire_tdt = 5; - } - else - { cout << "\n *** cas element Angle_mort actuellement non implante: , " - << elfro.TypeFrontiere() - << "\n ElContact::Trajectoire(int & test)"; - Sortie(1); - }; - - } - // ---- maintenant on regarde si on veut une normale lissée - else if (normale_lisser) + // ---- on regarde si on veut une normale lissée + // ou s'il s'agit d'angle mort + if ((normale_lisser) || (elfront->Angle_mort())) {// dans le cas d'une normale lissée on va utiliser une interpolation des normales aux noeuds ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // recup et actualisation du dernier plan tangent, coordonnées locale etc. diff --git a/contact/ElContact_2.cc b/contact/ElContact_2.cc index 3b8daf3..16d7549 100755 --- a/contact/ElContact_2.cc +++ b/contact/ElContact_2.cc @@ -1444,28 +1444,23 @@ Element::ResRaid ElContact::SM_K_charge_contact() Coordonnee ElContact::Calcul_Normale(int dim, Plan & pl, Droite & dr,int indic) { Coordonnee N(dim); // le vecteur normal initialisé à 0 + ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // --- le cas des angles morts est particulier, on commence par le traiter - if (elfront->Angle_mort()) - { ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite - if (elfro.TypeFrontiere() == "FrontPointF") - {// s'il s'agit d'un point il n'y a pas d'existance de normale -// // on remplace la notion de normale par le vecteur qui joint le noeud à la position -// // du point frontière -// N = elfro.Ref()- noeud->Coord2(); - // on utilise la normale au noeud - const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); - const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); - N = gr.ConteneurCoordonnee_const(); - } - else - { cout << "\n *** cas element Angle_mort actuellement non implante: , " - << elfro.TypeFrontiere() - << "\n ElContact::Calcul_Normale(..."; - Sortie(1); - }; + if ((elfront->Angle_mort()) && (elfro.Type_geom_front() == POINT_G)) + { // s'il s'agit d'un point il n'y a pas d'existance de normale + // on utilise la normale au noeud + const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); + const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); + N = gr.ConteneurCoordonnee_const(); } - // ---- maintenant on regarde si on veut une normale lissée - else if (normale_lisser) + // ---- maintenant on regarde 1) si on veut une normale lissée + // 2) qui représente aussi le cas d'angle mort avec une ligne normalement qu'en 3D, + // pour laquelle il n'y a pas d'existance d'une normale classique + + // on utilise ici la normale interpolée aux noeuds + else if ( (normale_lisser) + || ((elfront->Angle_mort()) && (elfro.Type_geom_front() == LIGNE)) + ) {// dans le cas d'une normale lissée on va utiliser une interpolation des normales aux noeuds ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // recup et actualisation du dernier plan tangent, coordonnées locale etc. @@ -1509,7 +1504,6 @@ Coordonnee ElContact::Calcul_Normale(int dim, Plan & pl, Droite & dr,int indic) else if ( (dim == 3) && (indic == 2)) // cas 3d avec une surface { N = pl.Vecplan(); - return N; } else if ( ((dim == 2) && (indic == 1)) @@ -1587,8 +1581,9 @@ Coordonnee ElContact::Calcul_Normale(int dim, Plan & pl, Droite & dr,int indic) cout << " une droite "; else cout << " un plan "; - cout << " n\'est pas actuellement traite " - << "ElContact::Normal(int dim, Plan & pl, Droite & dr, indic)" << endl; + cout << " n\'est pas actuellement traite \n" ; + Affiche(1); + cout << "ElContact::Normal(.." << endl; Sortie(1); }; return N; @@ -1610,30 +1605,26 @@ Tableau * ElContact::RecupInfo(Coordonnee& N,Coordonnee& M_impact, // different cas N.Change_dim(dim); // le vecteur normal // on commence par traiter le cas des éléments particuliers qui décrivent les angles morts - if (elfront->Angle_mort()) - { if (elfro.TypeFrontiere() == "FrontPointF") - {// s'il s'agit d'un point il n'y a pas d'existance de normale -// // on remplace la notion de normale par le vecteur qui joint le noeud à la position -// // du point frontière -// N = elfro.Ref()- noeud->Coord2(); - // on utilise la normale au noeud - elfro.DernierTangent(dr,pl,indic,false); - // on ne considère pas la variation de la normale - d_T = NULL; - const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); - const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); - N = gr.ConteneurCoordonnee_const(); - M_impact = elfro.Ref(); - } - else - { cout << "\n *** cas element Angle_mort actuellement non implante: , " - << elfro.TypeFrontiere() - << "\n ElContact::RecupInfo(.."; - Sortie(1); - }; + if ((elfront->Angle_mort()) && (elfro.Type_geom_front() == POINT_G)) + { // s'il s'agit d'un point il n'y a pas d'existance de normale + // on utilise la normale au noeud + elfro.DernierTangent(dr,pl,indic,false); + // on ne considère pas la variation de la normale + d_T = NULL; + const TypeQuelconque& tiq = tabNoeud(2)->Grandeur_quelconque(N_FRONT_t); + const Grandeur_coordonnee& gr= *((Grandeur_coordonnee*) (tiq.Const_Grandeur_pointee())); + N = gr.ConteneurCoordonnee_const(); + M_impact = elfro.Ref(); } // ---- maintenant on regarde si on veut une normale lissée - else if (normale_lisser) + // ---- maintenant on regarde 1) si on veut une normale lissée + // 2) qui représente aussi le cas d'angle mort avec une ligne normalement qu'en 3D, + // pour laquelle il n'y a pas d'existance d'une normale classique + + // on utilise ici la normale interpolée aux noeuds + else if ( (normale_lisser) + || ((elfront->Angle_mort()) && (elfro.Type_geom_front() == LIGNE)) + ) {// dans le cas d'une normale lissée on va utiliser une interpolation des normales aux noeuds ElFrontiere & elfro = *(elfront->Eleme()); // pour commodite // recup et actualisation du dernier plan tangent, coordonnées locale etc. @@ -1657,17 +1648,17 @@ Tableau * ElContact::RecupInfo(Coordonnee& N,Coordonnee& M_impact, N += phii(inoe)*Nnoe; }; N.Normer(); -//debug -if (Permet_affichage() > 4) - {cout << "\n debug ElContact::RecupInfo N= "< 4) +// {cout << "\n debug ElContact::RecupInfo N= "<Coord2();; break; + case 0: M_impact = tabNoeud(2)->Coord2(); break; default: {cout << "\n *** erreur indic =" << indic <<" c-a-d diff de 2 ou 1 ou 0 " << " on ne peut pas continuer !"; @@ -1784,7 +1775,9 @@ if (Permet_affichage() > 4) cout << " une droite "; else cout << " un plan "; - cout << " n\'est pas actuellement traite " << "ElContact::RecupInfo(.."<< endl; + cout << " n\'est pas actuellement traite \n" ; + Affiche(1); + cout << "ElContact::RecupInfo(.."<< endl; Sortie(1); }; }; @@ -1983,76 +1976,92 @@ double ElContact::CalFactPenal(const double& gap,double & d_beta_gap,int contact //******** partie en construction case LIGNE: - if (ParaGlob::AxiSymetrie()) - // cas où l'on est en axisymétrique, - { // on regarde la dimension de l'élément associé - if (elem->ElementGeometrie(X1).Dimension()==2) - // élément 2D en axi donc 3D en réalité - { // récup de la longueur de la frontière - double longueur = elfrontiere->LongueurApprox(); - // on calcul la surface associée - // M_impact(1) = x donc = le rayon car on est axi autour de y - //// double surf = longueur * M_impact(1) * ConstMath::Pi * 2.; - // 9 nov 2016: je crois qu'ici il y a une incompréhension. en fait il s'agit de la surface du maître - // donc il ne faut pas tenir compte du point d'impact ?? sinon par exemple si le point est au centre - // M_impact(1) = 0 et on n'a pas de surf !! - // mais plus généralement, le pt d'impact ne doit pas interférer à mon avis - double surf = longueur * ConstMath::Pi * 2.; - // on récupère la surface de l'élément qui a créé la frontière - double volume = elem->Volume(); - if (Dabs(volume) <= ConstMath::pasmalpetit) - { if (Permet_affichage() > 1) - cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " - << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; - if (Permet_affichage() > 4) - cout << "\n ElContact::CalFactPenal()"; - facteur = elfrontiere->MaxDiagonale_tdt(); - } - else // sinon c'est bon - {facteur = surf*surf/volume;} - } - else // sinon ce n'est pas pris en charge - {cout << "\n *** erreur: cas non pris en charge pour l'instant: " - << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) - << " contact sur une ligne axisyme trique" - << " \n ElContact::CalFactPenal() "; - Sortie(1); - }; + // --- on traite les cas particuliers des éléments d'angle mort + if (elfront->Angle_mort()) + { // il n'y a pas de notion de surface de contact et l'idée est de récupérer + // une grandeur du type module de compressibilité * une longueur caractéristique de l'élément + // on retient une expression simple + double longueur = elem->LongueurGeometrique_mini(1); + facteur = longueur; + //%%%% essai %%%% + longueur = elem->LongueurGeometrique_mini(1); + double volume = elem->Volume(); + facteur = volume / (longueur * longueur); + //%%%% fin essai %%%% } - else if (elem->ElementGeometrie(X1).Dimension()==3) - // cas d'une ligne frontière d'un élément coque ou plaque - { - // il faut finir la fonction ElFrontiere::LongueurApprox() - // cf le laïus qui est écrit dans la méthode - - - double surf = elfrontiere->SurfaceApprox(); - double volume = elem->Volume(); + else + {// --- cas d'un élément de contact autre que d'angle mort + if (ParaGlob::AxiSymetrie()) + // cas où l'on est en axisymétrique, + { // on regarde la dimension de l'élément associé + if (elem->ElementGeometrie(X1).Dimension()==2) + // élément 2D en axi donc 3D en réalité + { // récup de la longueur de la frontière + double longueur = elfrontiere->LongueurApprox(); + // on calcul la surface associée + // M_impact(1) = x donc = le rayon car on est axi autour de y + //// double surf = longueur * M_impact(1) * ConstMath::Pi * 2.; + // 9 nov 2016: je crois qu'ici il y a une incompréhension. en fait il s'agit de la surface du maître + // donc il ne faut pas tenir compte du point d'impact ?? sinon par exemple si le point est au centre + // M_impact(1) = 0 et on n'a pas de surf !! + // mais plus généralement, le pt d'impact ne doit pas interférer à mon avis + double surf = longueur * ConstMath::Pi * 2.; + // on récupère la surface de l'élément qui a créé la frontière + double volume = elem->Volume(); + if (Dabs(volume) <= ConstMath::pasmalpetit) + { if (Permet_affichage() > 1) + cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " + << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; + if (Permet_affichage() > 4) + cout << "\n ElContact::CalFactPenal()"; + facteur = elfrontiere->MaxDiagonale_tdt(); + } + else // sinon c'est bon + {facteur = surf*surf/volume;} + } + else // sinon ce n'est pas pris en charge + {cout << "\n *** erreur: cas non pris en charge pour l'instant: " + << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) + << " contact sur une ligne axisyme trique" + << " \n ElContact::CalFactPenal() "; + Sortie(1); + }; + } + else if (elem->ElementGeometrie(X1).Dimension()==3) + // cas d'une ligne frontière d'un élément coque ou plaque + { + // il faut finir la fonction ElFrontiere::LongueurApprox() + // cf le laïus qui est écrit dans la méthode - // virtual double EpaisseurMoyenne(Enum_dure ) - if (Dabs(volume) <= ConstMath::pasmalpetit) - { if (Permet_affichage() > 1) - cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " - << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; - if (Permet_affichage() > 4) - cout << "\n ElContact::CalFactPenal()"; - facteur = elfrontiere->MaxDiagonale_tdt(); - } - else // sinon c'est bon - {facteur = surf*surf/volume;} - } - else if (elem->ElementGeometrie(X1).Dimension()==2) - // cas d'une ligne frontière d'une plaque - { double longueur = elfrontiere->LongueurApprox(); - double epais = 0.; // init - if (!ParaGlob::AxiSymetrie()) // si on n'est pas en axi - {epais = ((ElemMeca*) elem)->EpaisseurMoyenne(TEMPS_tdt );} - else // si on est en axi - {epais = 1.;}; - double surf = elem->Volume() / epais; - facteur = surf/longueur; - }; + double surf = elfrontiere->SurfaceApprox(); + double volume = elem->Volume(); + + // virtual double EpaisseurMoyenne(Enum_dure ) + + if (Dabs(volume) <= ConstMath::pasmalpetit) + { if (Permet_affichage() > 1) + cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " + << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; + if (Permet_affichage() > 4) + cout << "\n ElContact::CalFactPenal()"; + facteur = elfrontiere->MaxDiagonale_tdt(); + } + else // sinon c'est bon + {facteur = surf*surf/volume;} + } + else if (elem->ElementGeometrie(X1).Dimension()==2) + // cas d'une ligne frontière d'une plaque + { double longueur = elfrontiere->LongueurApprox(); + double epais = 0.; // init + if (!ParaGlob::AxiSymetrie()) // si on n'est pas en axi + {epais = ((ElemMeca*) elem)->EpaisseurMoyenne(TEMPS_tdt );} + else // si on est en axi + {epais = 1.;}; + double surf = elem->Volume() / epais; + facteur = surf/longueur; + }; + }; break; //******** fin partie en construction @@ -2397,70 +2406,86 @@ double ElContact::CalFactPenalTangentiel(const double& dep_T,double & d_beta_dep //******** partie en construction case LIGNE: - if (ParaGlob::AxiSymetrie()) - // cas où l'on est en axisymétrique, - { // on regarde la dimension de l'élément associé - if (elem->ElementGeometrie(X1).Dimension()==2) - // élément 2D en axi donc 3D en réalité - { // récup de la longueur de la frontière - double longueur = elfrontiere->LongueurApprox(); - // on calcul la surface associée - double surf = longueur * ConstMath::Pi * 2.; - // on récupère la surface de l'élément qui a créé la frontière - double volume = elem->Volume(); - if (Dabs(volume) <= ConstMath::pasmalpetit) - { if (Permet_affichage() > 1) - cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " - << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; - if (Permet_affichage() > 4) - cout << "\n ElContact::CalFactPenalTangentiel()"; - facteur = elfrontiere->MaxDiagonale_tdt(); - } - else // sinon c'est bon - {facteur = surf*surf/volume;} - } - else // sinon ce n'est pas pris en charge - {cout << "\n *** erreur: cas non pris en charge pour l'instant: " - << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) - << " contact sur une ligne axisyme trique" - << " \n ElContact::CalFactPenalTangentiel() "; - Sortie(1); - }; + // --- on traite les cas particuliers des éléments d'angle mort + if (elfront->Angle_mort()) + { // il n'y a pas de notion de surface de contact et l'idée est de récupérer + // une grandeur du type module de compressibilité * une longueur caractéristique de l'élément + // on retient une expression simple + double longueur = elem->LongueurGeometrique_mini(1); + facteur = longueur; + //%%%% essai %%%% + longueur = elem->LongueurGeometrique_mini(1); + double volume = elem->Volume(); + facteur = volume / (longueur * longueur); + //%%%% fin essai %%%% } - else if (elem->ElementGeometrie(X1).Dimension()==3) - // cas d'une ligne frontière d'un élément coque ou plaque - { - // il faut finir la fonction ElFrontiere::LongueurApprox() - // cf le laïus qui est écrit dans la méthode - - - double surf = elfrontiere->SurfaceApprox(); - double volume = elem->Volume(); + else + {// --- cas d'un élément de contact autre que d'angle mort + if (ParaGlob::AxiSymetrie()) + // cas où l'on est en axisymétrique, + { // on regarde la dimension de l'élément associé + if (elem->ElementGeometrie(X1).Dimension()==2) + // élément 2D en axi donc 3D en réalité + { // récup de la longueur de la frontière + double longueur = elfrontiere->LongueurApprox(); + // on calcul la surface associée + double surf = longueur * ConstMath::Pi * 2.; + // on récupère la surface de l'élément qui a créé la frontière + double volume = elem->Volume(); + if (Dabs(volume) <= ConstMath::pasmalpetit) + { if (Permet_affichage() > 1) + cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " + << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; + if (Permet_affichage() > 4) + cout << "\n ElContact::CalFactPenalTangentiel()"; + facteur = elfrontiere->MaxDiagonale_tdt(); + } + else // sinon c'est bon + {facteur = surf*surf/volume;} + } + else // sinon ce n'est pas pris en charge + {cout << "\n *** erreur: cas non pris en charge pour l'instant: " + << " type de frontiere : " << Nom_type_geom(elfrontiere->Type_geom_front()) + << " contact sur une ligne axisyme trique" + << " \n ElContact::CalFactPenalTangentiel() "; + Sortie(1); + }; + } + else if (elem->ElementGeometrie(X1).Dimension()==3) + // cas d'une ligne frontière d'un élément coque ou plaque + { + // il faut finir la fonction ElFrontiere::LongueurApprox() + // cf le laïus qui est écrit dans la méthode - // virtual double EpaisseurMoyenne(Enum_dure ) - if (Dabs(volume) <= ConstMath::pasmalpetit) - { if (Permet_affichage() > 1) - cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " - << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; - if (Permet_affichage() > 4) - cout << "\n ElContact::CalFactPenalTangentiel()"; - facteur = elfrontiere->MaxDiagonale_tdt(); - } - else // sinon c'est bon - {facteur = surf*surf/volume;} + double surf = elfrontiere->SurfaceApprox(); + double volume = elem->Volume(); + + // virtual double EpaisseurMoyenne(Enum_dure ) + + if (Dabs(volume) <= ConstMath::pasmalpetit) + { if (Permet_affichage() > 1) + cout << "\n *** attention, le volume de l'element " << elem->Num_elt_const() << " du maillage " + << elem->Num_maillage() << " est nul, on utilise le max de diagonal "; + if (Permet_affichage() > 4) + cout << "\n ElContact::CalFactPenalTangentiel()"; + facteur = elfrontiere->MaxDiagonale_tdt(); + } + else // sinon c'est bon + {facteur = surf*surf/volume;} + } + else if (elem->ElementGeometrie(X1).Dimension()==2) + // cas d'une ligne frontière d'une plaque + { double longueur = elfrontiere->LongueurApprox(); + double epais = 0.; // init + if (!ParaGlob::AxiSymetrie()) // si on n'est pas en axi + {epais = ((ElemMeca*) elem)->EpaisseurMoyenne(TEMPS_tdt );} + else // si on est en axi + {epais = 1.;}; + double surf = elem->Volume() / epais; + facteur = surf/longueur; + }; } - else if (elem->ElementGeometrie(X1).Dimension()==2) - // cas d'une ligne frontière d'une plaque - { double longueur = elfrontiere->LongueurApprox(); - double epais = 0.; // init - if (!ParaGlob::AxiSymetrie()) // si on n'est pas en axi - {epais = ((ElemMeca*) elem)->EpaisseurMoyenne(TEMPS_tdt );} - else // si on est en axi - {epais = 1.;}; - double surf = elem->Volume() / epais; - facteur = surf/longueur; - }; break; //******** fin partie en construction diff --git a/contact/LesContacts_3.cc b/contact/LesContacts_3.cc index 4c46001..12f078c 100755 --- a/contact/LesContacts_3.cc +++ b/contact/LesContacts_3.cc @@ -1420,13 +1420,23 @@ void LesContacts::ElementAngleMort(LesMaillages& lesMail) // on regarde si cet élément ne fait pas partie de la liste d'élément // qui contient déjà une frontière if (find(li_Elem_avec_front.begin(),li_Elem_avec_front.end(),el) == li_Elem_avec_front.end()) - // on n'a pas trouvé d'élément, on ajoute l'élément front - // on doit indiquer qu'il s'agit d'un élément d'angle mort - // on se sert d'un élément intermédiaire - {Front inter((*jl)); - inter.Change_angle_mort(1); - la_list_des_front.push_back(inter); - }; + // on n'a pas trouvé d'élément, + //on ajoute l'élément front sauf cas particulier + {// on fait quelques tests pour écarter des cas non recevable + bool continuer = true; // par défaut + // si on est en 3D et que l'élément est lui même un élément ligne, ou un + // élément de surface, il ne constitue pas un angle mort puisqu'il n'y à pas de volume mort associé + Enum_type_geom enu_type_geom_ele = Type_geom_generique(el->Id_geometrie()); + if ((enu_type_geom_ele == LIGNE) && (enu_type_geom_ele == SURFACE)) + {continuer = false;} + if (continuer) + // on doit indiquer qu'il s'agit d'un élément d'angle mort + // on se sert d'un élément intermédiaire + {Front inter((*jl)); + inter.Change_angle_mort(1); + la_list_des_front.push_back(inter); + }; + } }; int nb_front_fin = la_list_des_front.size(); // la taille finale if (niveau_commentaire_lescontacts >2 ) @@ -1465,13 +1475,24 @@ void LesContacts::ElementAngleMort(LesMaillages& lesMail) // on regarde si cet élément ne fait pas partie de la liste d'élément // qui contient déjà une frontière if (find(li_Elem_avec_front.begin(),li_Elem_avec_front.end(),el) == li_Elem_avec_front.end()) - // on n'a pas trouvé d'élément, on ajoute l'élément front - // on doit indiquer qu'il s'agit d'un élément d'angle mort - // on se sert d'un élément intermédiaire - {Front inter((*jl)); - inter.Change_angle_mort(1); - la_list_des_front.push_back(inter); - }; + // on n'a pas trouvé d'élément, + //on ajoute l'élément front sauf cas particulier + {// on fait quelques tests pour écarter des cas non recevable + bool continuer = true; // par défaut + // si on est en 3D et que l'élément est lui même un élément ligne, ou un + // élément de surface, il ne constitue pas un angle mort puisqu'il n'y à pas de volume mort associé + Enum_type_geom enu_type_geom_ele = Type_geom_generique(el->Id_geometrie()); + if ((enu_type_geom_ele == LIGNE) && (enu_type_geom_ele == SURFACE)) + {continuer = false;} + if (continuer) + // on n'a pas trouvé d'élément, on ajoute l'élément front + // on doit indiquer qu'il s'agit d'un élément d'angle mort + // on se sert d'un élément intermédiaire + {Front inter((*jl)); + inter.Change_angle_mort(1); + la_list_des_front.push_back(inter); + }; + }; }; int nb_front_fin = la_list_des_front.size(); // la taille finale if (niveau_commentaire_lescontacts >2 ) @@ -1564,36 +1585,52 @@ void LesContacts::ElementAngleMort(LesMaillages& lesMail) // on passe par l'intermédiaire de lesMail pour pouvoir modifier l'élément // car c'est l'élément qui stocke et gère la frontière Element& elem_non_const = lesMail.Element_LesMaille(el->Num_maillage(), el->Num_elt_const()); - // 3. on crée l'élément de frontière ad hoc - ElFrontiere* const elfpoint = elem_non_const.Frontiere_points(num_noe_local,true); - // ramene le numero de la frontière passée en argument si elle existe actuellement au niveau de l'élément - // sinon ramène 0 - // ramene également type_front: qui indique le type de frontière: POINT_G, LIGNE ou SURFACE - // c'est une méthode très longue, a utiliser avec précaution - Enum_type_geom type_front; - int num_front = el->Num_frontiere(*elfpoint,type_front); - // 4. on crée l'élément front ad hoc et on l'ajoute à la liste - const DdlElement & ddlElem = el->TableauDdl(); // récup des ddl - int angle_mort = 1; - la_list_des_front.push_front(Front(*elfpoint,&elem_non_const,num_front,angle_mort)); - LaLIST_io ::iterator iik = la_list_des_front.begin(); // récup de l'iterator - // on met à jour les frontières mitoyennes - int num_noeud = ttn(1)->Num_noeud(); // le num global du noeud - // récup de la liste des frontières initiales contenant le noeud - list & t_liFront_noeud_ijk = t_liFront_noeud(ilistfront)(1)(num_noeud); - list ::iterator km,kmfin=t_liFront_noeud_ijk.end(); - for (km=t_liFront_noeud_ijk.begin();km != kmfin;km++) - {(*km)->AjoutMitoyen(&(*iik));// ajout frontière pour la frontière initiale - (*iik).AjoutMitoyen(*km); // ajout frontière pour la nouvelle frontière d'angle mort - }; - // on sauvegarde le nouvelle front dans une liste spécifique - front_angle_mort.push_front(&(*iik)); + // on fait quelques tests pour écarter des cas non recevable + bool continuer = true; // par défaut + // a) si on est en 2D axi et que l'élément est un élément noeud ou un élément ligne, cela veut dire + // qu'il génére par rotation en 3D : une ligne ou une surface, or lignes ou surfaces n'intègrent pas + // d'élément de contact d'angle mort puisqu'il n'y à pas de volume mort associé + Enum_type_geom enu_type_geom_ele = Type_geom_generique(elem_non_const.Id_geometrie()); + if (ParaGlob::AxiSymetrie() && ((enu_type_geom_ele == POINT_G) || (enu_type_geom_ele == LIGNE))) + {continuer = false;} + // b) si on est en 3D et que l'élément est un élément noeud, ou ligne, ou surface, même sentence + else if ((ParaGlob::Dimension() == 3) + && ((enu_type_geom_ele == POINT_G) || (enu_type_geom_ele == LIGNE) || (enu_type_geom_ele == SURFACE)) + ) + {continuer = false;} + // on continue sauf cas particulier + if (continuer) + {// 3. on crée l'élément de frontière ad hoc + ElFrontiere* const elfpoint = elem_non_const.Frontiere_points(num_noe_local,true); + // ramene le numero de la frontière passée en argument si elle existe actuellement au niveau de l'élément + // sinon ramène 0 + // ramene également type_front: qui indique le type de frontière: POINT_G, LIGNE ou SURFACE + // c'est une méthode très longue, a utiliser avec précaution + Enum_type_geom type_front; + int num_front = el->Num_frontiere(*elfpoint,type_front); + // 4. on crée l'élément front ad hoc et on l'ajoute à la liste + const DdlElement & ddlElem = el->TableauDdl(); // récup des ddl + int angle_mort = 1; + la_list_des_front.push_front(Front(*elfpoint,&elem_non_const,num_front,angle_mort)); + LaLIST_io ::iterator iik = la_list_des_front.begin(); // récup de l'iterator + // on met à jour les frontières mitoyennes + int num_noeud = ttn(1)->Num_noeud(); // le num global du noeud + // récup de la liste des frontières initiales contenant le noeud + list & t_liFront_noeud_ijk = t_liFront_noeud(ilistfront)(1)(num_noeud); + list ::iterator km,kmfin=t_liFront_noeud_ijk.end(); + for (km=t_liFront_noeud_ijk.begin();km != kmfin;km++) + {(*km)->AjoutMitoyen(&(*iik));// ajout frontière pour la frontière initiale + (*iik).AjoutMitoyen(*km); // ajout frontière pour la nouvelle frontière d'angle mort + }; + // on sauvegarde le nouvelle front dans une liste spécifique + front_angle_mort.push_front(&(*iik)); - if (niveau_commentaire_lescontacts >4 ) - { cout << "\n ajout front point: noeud "<Num_noeud(num_noe_local) - << " mail: " << jf //<< " zone " - << " element: "<Num_elt_const() ; - }; + if (niveau_commentaire_lescontacts >4 ) + { cout << "\n ajout front point: noeud "<Num_noeud(num_noe_local) + << " mail: " << jf //<< " zone " + << " element: "<Num_elt_const() ; + }; + }; }; }; }; @@ -1715,38 +1752,54 @@ void LesContacts::ElementAngleMort(LesMaillages& lesMail) // on passe par l'intermédiaire de lesMail pour pouvoir modifier l'élément // car c'est l'élément qui stocke et gère la frontière Element& elem_non_const = lesMail.Element_LesMaille(el->Num_maillage(), el->Num_elt_const()); - // 3. on crée l'élément de frontière ad hoc - ElFrontiere* const elfpoint = elem_non_const.Frontiere_points(num_noe_local,true); - // ramene le numero de la frontière passée en argument si elle existe actuellement au niveau de l'élément - // sinon ramène 0 - // ramene également type_front: qui indique le type de frontière: POINT_G, LIGNE ou SURFACE - // c'est une méthode très longue, a utiliser avec précaution - Enum_type_geom type_front; - int num_front = el->Num_frontiere(*elfpoint,type_front); - // 4. on crée l'élément front ad hoc et on l'ajoute à la liste - const DdlElement & ddlElem = el->TableauDdl(); // récup des ddl - int angle_mort = 1; - la_list_des_front.push_front(Front(*elfpoint,&elem_non_const,num_front,angle_mort)); - LaLIST_io ::iterator iik = la_list_des_front.begin(); // récup de l'iterator - // on met à jour les frontières mitoyennes - // // t_liFront_noeud(i)(j)(k) : la liste des frontières initiales qui contiennent le noeud k de la zone j du maillage i - // Tableau < Tableau < Tableau < list > > > t_liFront_noeud; - int num_noeud = ttn(1)->Num_noeud(); // le num global du noeud - // récup de la liste des frontières initiales contenant le noeud - list & t_liFront_noeud_ijk = t_liFront_noeud(ilistfront)(i_zone)(num_noeud); - list ::iterator km,kmfin=t_liFront_noeud_ijk.end(); - for (km=t_liFront_noeud_ijk.begin();km != kmfin;km++) - {(*km)->AjoutMitoyen(&(*iik));// ajout frontière pour la frontière initiale - (*iik).AjoutMitoyen(*km); // ajout frontière pour la nouvelle frontière d'angle mort - }; - // on sauvegarde le nouvelle front dans une liste spécifique - front_angle_mort.push_front(&(*iik)); + // on fait quelques tests pour écarter des cas non recevable + bool continuer = true; // par défaut + // a) si on est en 2D axi et que l'élément est un élément noeud ou un élément ligne, cela veut dire + // qu'il génére par rotation en 3D : une ligne ou une surface, or lignes ou surfaces n'intègrent pas + // d'élément de contact d'angle mort puisqu'il n'y à pas de volume mort associé + Enum_type_geom enu_type_geom_ele = Type_geom_generique(elem_non_const.Id_geometrie()); + if (ParaGlob::AxiSymetrie() && ((enu_type_geom_ele == POINT_G) || (enu_type_geom_ele == LIGNE))) + {continuer = false;} + // b) si on est en 3D et que l'élément est un élément noeud, ou ligne, ou surface, même sentence + else if ((ParaGlob::Dimension() == 3) + && ((enu_type_geom_ele == POINT_G) || (enu_type_geom_ele == LIGNE) || (enu_type_geom_ele == SURFACE)) + ) + {continuer = false;} + // on continue sauf cas particulier + if (continuer) + {// 3. on crée l'élément de frontière ad hoc + ElFrontiere* const elfpoint = elem_non_const.Frontiere_points(num_noe_local,true); + // ramene le numero de la frontière passée en argument si elle existe actuellement au niveau de l'élément + // sinon ramène 0 + // ramene également type_front: qui indique le type de frontière: POINT_G, LIGNE ou SURFACE + // c'est une méthode très longue, a utiliser avec précaution + Enum_type_geom type_front; + int num_front = el->Num_frontiere(*elfpoint,type_front); + // 4. on crée l'élément front ad hoc et on l'ajoute à la liste + const DdlElement & ddlElem = el->TableauDdl(); // récup des ddl + int angle_mort = 1; + la_list_des_front.push_front(Front(*elfpoint,&elem_non_const,num_front,angle_mort)); + LaLIST_io ::iterator iik = la_list_des_front.begin(); // récup de l'iterator + // on met à jour les frontières mitoyennes + // // t_liFront_noeud(i)(j)(k) : la liste des frontières initiales qui contiennent le noeud k de la zone j du maillage i + // Tableau < Tableau < Tableau < list > > > t_liFront_noeud; + int num_noeud = ttn(1)->Num_noeud(); // le num global du noeud + // récup de la liste des frontières initiales contenant le noeud + list & t_liFront_noeud_ijk = t_liFront_noeud(ilistfront)(i_zone)(num_noeud); + list ::iterator km,kmfin=t_liFront_noeud_ijk.end(); + for (km=t_liFront_noeud_ijk.begin();km != kmfin;km++) + {(*km)->AjoutMitoyen(&(*iik));// ajout frontière pour la frontière initiale + (*iik).AjoutMitoyen(*km); // ajout frontière pour la nouvelle frontière d'angle mort + }; + // on sauvegarde le nouvelle front dans une liste spécifique + front_angle_mort.push_front(&(*iik)); - if (niveau_commentaire_lescontacts >4 ) - { cout << "\n ajout front point: noeud "<Num_noeud(num_noe_local) - << " mail: " << jf << " zone "<< i_zone - << " element: "<Num_elt_const() ; - }; + if (niveau_commentaire_lescontacts >4 ) + { cout << "\n ajout front point: noeud "<Num_noeud(num_noe_local) + << " mail: " << jf << " zone "<< i_zone + << " element: "<Num_elt_const() ; + }; + }; }; }; };