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

1363 lines
59 KiB
C++
Raw Permalink Normal View History

ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// FICHIER : Hyper_externe_W.cc
// CLASSE : Hyper_externe_W
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
//#include "Debug.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "CharUtil.h"
#include "Hyper_externe_W.h"
#include "MathUtil.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"
// CONSTRUCTEURS :
// ---------- classe de stockage des grandeurs spécifiques pour la loi ---------
Hyper_externe_W::SaveResulHyper_externe_W::SaveResulHyper_externe_W(): // constructeur par défaut
SaveResulHyper_W_gene_3D()
,module_compressibilite(0.),module_compressibilite_t(0.)
,W_poten(NULL),W_poten_t(NULL)
,inv_B(NULL),inv_B_t(NULL)
,inv_J(NULL),inv_J_t(NULL)
{};
Hyper_externe_W::SaveResulHyper_externe_W::SaveResulHyper_externe_W(int sortie_post): // avec init ou pas
SaveResulHyper_W_gene_3D(sortie_post)
,module_compressibilite(0.),module_compressibilite_t(0.)
,W_poten(NULL),W_poten_t(NULL)
,inv_B(NULL),inv_B_t(NULL)
,inv_J(NULL),inv_J_t(NULL)
{ if (sortie_post)
{W_poten = new Vecteur(9);
W_poten_t = new Vecteur(9);
inv_B = new Vecteur(3);
inv_B_t = new Vecteur(3);
inv_J = new Vecteur(3);
inv_J_t = new Vecteur(3);
};
};
Hyper_externe_W::SaveResulHyper_externe_W::SaveResulHyper_externe_W(const SaveResulHyper_externe_W& sav): // de copie
SaveResulHyper_W_gene_3D(sav) // par défaut on sauvegarde les infos
,module_compressibilite(sav.module_compressibilite),module_compressibilite_t(sav.module_compressibilite)
,W_poten(NULL),W_poten_t(NULL)
,inv_B(NULL),inv_B_t(NULL)
,inv_J(NULL),inv_J_t(NULL)
{if (sav.W_poten != NULL)
{W_poten = new Vecteur(*sav.W_poten);
W_poten_t = new Vecteur(*sav.W_poten_t);
inv_B = new Vecteur(*sav.inv_B);
inv_B_t = new Vecteur(*sav.inv_B_t);
inv_J = new Vecteur(*sav.inv_J);
inv_J_t = new Vecteur(*sav.inv_J_t);
};
};
Hyper_externe_W::SaveResulHyper_externe_W::~SaveResulHyper_externe_W() // destructeur
{if (W_poten != NULL)
{delete W_poten;
delete W_poten_t;
delete inv_B;
delete inv_B_t;
delete inv_J;
delete inv_J_t;
};
};
// affectation
Loi_comp_abstraite::SaveResul & Hyper_externe_W::SaveResulHyper_externe_W::operator = ( const SaveResul & a)
{ // on s'occupe tout d'abord de la classe mère
SaveResulHyper_W_gene_3D::operator = (a);
// puis les nouvelles variables
SaveResulHyper_externe_W& sav = *((SaveResulHyper_externe_W*) &a);
module_compressibilite = sav.module_compressibilite;
module_compressibilite_t = sav.module_compressibilite_t;
if (sav.W_poten != NULL)
{if (W_poten != NULL)
{*W_poten = (*sav.W_poten);
*W_poten_t = (*sav.W_poten_t);
*inv_B = (*sav.inv_B);
*inv_B_t = (*sav.inv_B_t);
*inv_J = (*sav.inv_J);
*inv_J_t = (*sav.inv_J_t);
}
else // sinon il faut définir les conteneurs
{W_poten = new Vecteur(*sav.W_poten);
W_poten_t = new Vecteur(*sav.W_poten_t);
inv_B = new Vecteur(*sav.inv_B);
inv_B_t = new Vecteur(*sav.inv_B_t);
inv_J = new Vecteur(*sav.inv_J);
inv_J_t = new Vecteur(*sav.inv_J_t);
};
}
else
{if (W_poten != NULL)
// il faut détruire les conteneurs
{delete W_poten;W_poten=NULL;
delete W_poten_t;W_poten_t=NULL;
delete inv_B;inv_B=NULL;
delete inv_B_t;inv_B_t=NULL;
delete inv_J;inv_J=NULL;
delete inv_J_t;inv_J_t=NULL;
};
// sinon rien n'a faire c'est ok
};
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 Hyper_externe_W::SaveResulHyper_externe_W::Lecture_base_info (ifstream& ent,const int cas)
{ // entête via la classe mère
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Lecture_base_info(ent,cas);
// puis les données spécifiques
string toto;int test;
ent >> toto >> test;
if (test)
{if (W_poten == NULL)
{ W_poten = new Vecteur(9);
W_poten_t = new Vecteur(9);
inv_B = new Vecteur(3);
inv_B_t = new Vecteur(3);
inv_J = new Vecteur(3);
inv_J_t = new Vecteur(3);
};
ent >> toto >>(*W_poten)
>> toto >>(*inv_B)
>> toto >>(*inv_J);
*W_poten_t = (*W_poten);
*inv_B_t = (*inv_B);
*inv_J_t = (*inv_J);
}
else
{if (W_poten != NULL)
// il faut détruire les conteneurs
{delete W_poten;W_poten=NULL;
delete W_poten_t;W_poten_t=NULL;
delete inv_B;inv_B=NULL;
delete inv_B_t;inv_B_t=NULL;
delete inv_J;inv_J=NULL;
delete inv_J_t;inv_J_t=NULL;
};
// sinon rien n'a faire c'est ok
};
//module_compressibilite
ent >> toto >> module_compressibilite;
module_compressibilite_t = module_compressibilite;
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables
//(supposées comme telles)
void Hyper_externe_W::SaveResulHyper_externe_W::Ecriture_base_info(ofstream& sort,const int cas)
{ // entête via la classe mère
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Ecriture_base_info(sort,cas);
// puis les données spécifiques
sort << "\n dat_Hyper_externe_W: ";
if (W_poten != NULL)
{sort << " 1 \n fonction_potentiel " << (*W_poten);
sort << "\n invar_B "<< (*inv_B) << " invar_J "<< (*inv_J);
}
else
{ sort << " 0 ";};
// module_compressibilite
sort << " module_compressibilite "<<module_compressibilite;
};
// mise à jour des informations transitoires
void Hyper_externe_W::SaveResulHyper_externe_W::TdtversT()
{ SaveResulHyper_W_gene_3D::TdtversT();
if (W_poten != NULL)
{*W_poten_t = (*W_poten);
*inv_B_t = (*inv_B);
*inv_J_t = (*inv_J);
};
module_compressibilite_t = module_compressibilite;
};
void Hyper_externe_W::SaveResulHyper_externe_W::TversTdt()
{ SaveResulHyper_W_gene_3D::TversTdt();
if (W_poten != NULL)
{*W_poten = (*W_poten_t);
*inv_B = (*inv_B_t);
*inv_J = (*inv_J_t);
};
module_compressibilite = module_compressibilite_t;
};
// affichage à l'écran des infos
void Hyper_externe_W::SaveResulHyper_externe_W::Affiche() const
{ // entête via la classe mère
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Affiche();
// puis les données spécifiques
cout << "\n dat_Hyper_externe_W: ";
if (W_poten != NULL)
{cout << "\n fonction_potentiel " << (*W_poten);
cout << "\n invar_B "<< (*inv_B) << " invar_J "<< (*inv_J);
};
// module_compressibilite
cout << " module_compressibilite "<<module_compressibilite << " ";
};
// ---------- classe Hyper_externe_W ---------
Hyper_externe_W::Hyper_externe_W () : // Constructeur par defaut
Hyper_W_gene_3D(HYPER_EXTERNE_W,CAT_MECANIQUE,3)
,W_potentiel(NULL)
,W_d(0.),W_v(0.)
,W_d_J1(0.),W_d_J2(0.),W_v_J3(0.),W_v_J3J3(0.)
,W_d_J1_2(0.),W_d_J1_J2(0.),W_d_J2_2(0.)
// variables de travail
,invar(),exclure_ddenum()
{ // on remplit invar, pour l'appel de la fonction potentiel nD
Grandeur_scalaire_double grand_courant(0.);
{TypeQuelconque typQ(INVAR_B1,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_B2,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_B3,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_J1,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_J2,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_J3,EPS11,grand_courant);
invar.push_back(typQ);
}
// on rempli exclure_ddenum
exclure_ddenum.push_back(Ddl_enum_etendu("I_B"));
exclure_ddenum.push_back(Ddl_enum_etendu("II_B"));
exclure_ddenum.push_back(Ddl_enum_etendu("III_B"));
exclure_ddenum.push_back(Ddl_enum_etendu("J1"));
exclure_ddenum.push_back(Ddl_enum_etendu("J2"));
exclure_ddenum.push_back(Ddl_enum_etendu("J3"));
};
// Constructeur de copie
Hyper_externe_W::Hyper_externe_W (const Hyper_externe_W& loi) :
Hyper_W_gene_3D(loi)
,W_potentiel(loi.W_potentiel)
,W_d(loi.W_d),W_v(loi.W_v)
,W_d_J1(loi.W_d_J1),W_d_J2(loi.W_d_J2),W_v_J3(loi.W_v_J3),W_v_J3J3(loi.W_v_J3J3)
,W_d_J1_2(loi.W_d_J1_2),W_d_J1_J2(loi.W_d_J1_J2),W_d_J2_2(loi.W_d_J2_2)
// variables de travail
,invar(),exclure_ddenum()
{// on regarde s'il s'agit d'une fonction nD locale ou globale
if (W_potentiel != NULL)
if (W_potentiel->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
W_potentiel = Fonction_nD::New_Fonction_nD(*loi.W_potentiel);
};
// on remplit invar, pour l'appel de la fonction potentiel nD
Grandeur_scalaire_double grand_courant(0.);
{TypeQuelconque typQ(INVAR_B1,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_B2,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_B3,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_J1,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_J2,EPS11,grand_courant);
invar.push_back(typQ);
}
{TypeQuelconque typQ(INVAR_J3,EPS11,grand_courant);
invar.push_back(typQ);
}
// on rempli exclure_ddenum
exclure_ddenum.push_back(Ddl_enum_etendu("I_B"));
exclure_ddenum.push_back(Ddl_enum_etendu("II_B"));
exclure_ddenum.push_back(Ddl_enum_etendu("III_B"));
exclure_ddenum.push_back(Ddl_enum_etendu("J1"));
exclure_ddenum.push_back(Ddl_enum_etendu("J2"));
exclure_ddenum.push_back(Ddl_enum_etendu("J3"));
};
Hyper_externe_W::~Hyper_externe_W ()
// Destructeur
{ if (W_potentiel != NULL)
if (W_potentiel->NomFonction() == "_") delete W_potentiel;
};
// Lecture des donnees de la classe sur fichier
void Hyper_externe_W::LectureDonneesParticulieres
(UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string nom_class_methode("Hyper_externe_W::LectureDonneesParticulieres");
string nom;
string mot_cle1="W_potentiel_fonction_nD:";
// on lit le nom de la fonction
string nom_fonct;
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle1,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer("**erreur02 en lecture** "+mot_cle1);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{W_potentiel = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
W_potentiel = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
W_potentiel->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (W_potentiel->NbComposante() != 9 )
{ cout << "\n erreur en lecture de W_potentiel_fonction_nD:, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << W_potentiel->NbComposante()
<< " composantes alors qu'elle devrait en avoir 9 ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur03** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// on regarde si la fonction nD intègre la température
const Tableau <Ddl_enum_etendu>& tab_enu = W_potentiel->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
// lecture de l'indication éventuelle du post traitement
string le_mot_cle = "sortie_post_";
entreePrinc->Lecture_un_parametre_int(0,nom_class_methode,0,1,le_mot_cle,sortie_post);
// --- appel au niveau de la classe mère
// ici il n'y a pas de type de déformation associé
// mais on prend la def standart d'almansi, pour les fonctions associées éventuelles
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD,true);
};
// affichage de la loi
void Hyper_externe_W::Affiche() const
{ cout << " \n loi de comportement hyper elastique 3D Hyper_externe_W \n";
cout << "\n la fonction potentiel calculee avec une fonction nD: ";
cout << "W_potentiel_fonction_nD:" << " ";
if (W_potentiel->NomFonction() != "_")
cout << W_potentiel->NomFonction();
else
W_potentiel->Affiche();
// appel de la classe mère
Loi_comp_abstraite::Affiche_don_classe_abstraite();
};
// affichage et definition interactive des commandes particulières à chaques lois
void Hyper_externe_W::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");
sort << "\n# ---------------------------------------------------------------------"
<< "\n# |... loi hyper elastique 3D Hyper_externe_W .. |"
<< "\n# ---------------------------------------------------------------------"
<< "\n\n# exemple de definition de loi"
<< "\n W_potentiel_fonction_nD: ma_fct_nD_W "
<< "\n# .. fin de la definition de la loi Hyper_externe_W \n" << endl;
if ((rep != "o") && (rep != "O" ) && (rep != "0") )
{ sort << "\n# "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# la loi se construit via l'utilisation d'une fonction vectorielle nD "
<< "\n# celle-ci doit ramener un tableau tab_val de 9 valeurs tels que: "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# (1) -> W_d c-a-d la valeur de la partie deviatorique du potentiel "
<< "\n# (2) -> W_v c-a-d la valeur de la partie spherique du potentiel "
<< "\n# (3) -> W_d_J1 c-a-d variation de W_d / a J1 "
<< "\n# (4) -> W_d_J2 c-a-d variation de W_d / a J2 "
<< "\n# (5) -> W_d_J1_2 c-a-d variation seconde de W_d / a J1 "
<< "\n# (6) -> W_d_J1_J2 c-a-d variation de W_d / a J1 et J2 "
<< "\n# (7) -> W_d_J2_2 c-a-d variation seconde de W_d / a J2 "
<< "\n# (8) -> W_v_J3 c-a-d variation de W_v / J3 "
<< "\n# (9) -> W_v_J3J3 c-a-d variation seconde de W_v / J3 "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# les variables de passage a la fonction sont celles definies lors de la definition "
<< "\n# de la fonction nD, avec en particulier: "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# INVAR_B1 -> I_B "
<< "\n# INVAR_B2 -> II_B "
<< "\n# INVAR_B3 -> III_B "
<< "\n# INVAR_J1 -> J1 "
<< "\n# INVAR_J2 -> J2 "
<< "\n# INVAR_J3 -> J3 "
<< "\n# on peut aussi y ajouter toutes les grandeurs locales ou globales generiques "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# il est possible de recuperer differentes grandeurs de travail par exemple "
<< "\n# l'intensite du potentiel, comme ces grandeurs sont calculees au moment de la resolution "
<< "\n# et ne sont pas stockees, il faut indiquer le mot cle sortie_post_ suivi de 1 (par defaut = 0) "
<< "\n# ensuite au moment de la constitution du .CVisu on aura acces aux grandeurs de travail "
<< "\n# ex: "
<< "\n# sortie_post_ 1 "
<< "\n# ce mot cle est le dernier des parametres specifiques de la loi il doit se situe "
<< "\n# a la fin de la derniere ligne de donnees "
<< "\n#"
<< "\n#------------------------------------------------------------------------------------"
<< "\n# NB: pour la definition de la fct nD il est aussi possible"
<< "\n# d'introduire directement la fonction a la place de son nom, comme pour "
<< "\n# toutes les autres lois "
<< endl;
};
// appel de la classe Hyper_W_gene_3D
Hyper_W_gene_3D::Info_commande_LoisDeComp_hyper3D(entreePrinc);
// appel de la classe mère
Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);
};
// test si la loi est complete
int Hyper_externe_W::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
if (W_potentiel == NULL)
{ret = 0;};
//
if (ret == 0)
{this-> Affiche();
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 Hyper_externe_W::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string nom; bool test;
if (cas == 1)
{ ent >> nom;
if (nom != "HYPER_EXTERNE_W")
{ cout << "\n erreur en lecture de la loi : Hyper_externe_W, on attendait le mot cle : HYPER_EXTERNE_W "
<< "\n HYPER_EXTERNE_W::Lecture_base_info_loi(...";
Sortie(1);
};
ent >> nom >> sortie_post;
int test; // sert pour le test de l'existence de la fonction
// compressibilité
ent >> test ;
if (test == 1)
{W_potentiel = lesFonctionsnD.Lecture_pour_base_info(ent,cas,W_potentiel);
}
else
{// cas où on ne peut pas lire la fonction
if (W_potentiel != NULL)
if (W_potentiel->NomFonction() == "_")
{delete W_potentiel;W_potentiel=NULL;
cout << "\n *** attention, on supprime la fonction potentiel, car"
<< " on ne peut pas la lire, l'utilisation de la loi est impossible par la suite ";
this->Affiche();
};
};
};
// appel de la classe mère
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 Hyper_externe_W::Ecriture_base_info_loi(ofstream& sort,const int cas)
{ if (cas == 1)
{ sort << " HYPER_EXTERNE_W sortie_post_ "<< sortie_post ;
if (W_potentiel != NULL)
{sort << " 1 W_potentiel: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,W_potentiel);
}
else
{ sort << " 0 W_potentiel: " << " ";};
}
// appel de la classe mère
Loi_comp_abstraite::Ecriture_don_base_info(sort,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 Hyper_externe_W::Grandeur_particuliere
(bool absolue,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon ,list<int>& decal ) const
{// tout d'abord on récupère le conteneur
SaveResulHyper_externe_W & sav = *((SaveResulHyper_externe_W*) saveDon);
// on passe en revue la liste
List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end();
list<int>::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
{switch (tipParticu.EnuTypeQuelconque().EnumTQ())
{
case MODULE_COMPRESSIBILITE:
{ Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
tyTQ(1+(*idecal))=sav.module_compressibilite;(*idecal)++;
break;
}
case INVAR_B1 :
{Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) ((*itq).Grandeur_pointee()));
if (sortie_post) {tyTQ(1+(*idecal))=(*sav.inv_B)(1);}
else {tyTQ(1+(*idecal))=0.;}
(*idecal)++;
break;
}
case INVAR_B2 :
{Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) ((*itq).Grandeur_pointee()));
if (sortie_post) {tyTQ(1+(*idecal))=(*sav.inv_B)(2);}
else {tyTQ(1+(*idecal))=0.;}
(*idecal)++;
break;
}
case INVAR_B3 :
{Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) ((*itq).Grandeur_pointee()));
if (sortie_post) {tyTQ(1+(*idecal))=(*sav.inv_B)(3);}
else {tyTQ(1+(*idecal))=0.;}
(*idecal)++;
break;
}
case INVAR_J1 :
{Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) ((*itq).Grandeur_pointee()));
if (sortie_post) {tyTQ(1+(*idecal))=(*sav.inv_J)(1);}
else {tyTQ(1+(*idecal))=0.;}
(*idecal)++;
break;
}
case INVAR_J2 :
{Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) ((*itq).Grandeur_pointee()));
if (sortie_post) {tyTQ(1+(*idecal))=(*sav.inv_J)(2);}
else {tyTQ(1+(*idecal))=0.;}
(*idecal)++;
break;
}
case INVAR_J3 :
{Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) ((*itq).Grandeur_pointee()));
if (sortie_post) {tyTQ(1+(*idecal))=(*sav.inv_J)(3);}
else {tyTQ(1+(*idecal))=0.;}
(*idecal)++;
break;
}
case FCT_POTENTIEL_ND :
{Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) ((*itq).Grandeur_pointee()));
if (sortie_post)
{tyTQ(1+(*idecal))=(*sav.W_poten);}
else
{tyTQ(1+(*idecal)).Zero();}
(*idecal)++;
break;
}
case MODULE_CISAILLEMENT:
{ Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
tyTQ(1+(*idecal))=0.;(*idecal)++;
break;
}
default: ;// on ne fait rien
};
}; // 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 Hyper_externe_W::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
{//on commence par définir une grandeur_scalaire_double
Tableau <double> 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: intéressant quand il dépend de la température
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == MODULE_COMPRESSIBILITE)
{ 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(MODULE_COMPRESSIBILITE,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ1);
};
};
// $$$ cas de MODULE_CISAILLEMENT: intéressant quand il dépend de la température
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == MODULE_CISAILLEMENT)
{ 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_CISAILLEMENT,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
// $$$ cas des invariants de B
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == INVAR_B1)
{ 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(INVAR_B1,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == INVAR_B2)
{ 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(INVAR_B2,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == INVAR_B3)
{ 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(INVAR_B3,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
// $$$ cas des invariants de J
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == INVAR_J1)
{ 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(INVAR_J1,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == INVAR_J2)
{ 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(INVAR_J2,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == INVAR_J3)
{ 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(INVAR_J3,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ2);
};
};
// $$$ cas de la fonction potentiel et de ses dérivées
{Vecteur v(9);
List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == FCT_POTENTIEL_ND)
{ Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
int taille = tyTQ.Taille()+1;
tyTQ.Change_taille(taille); nexistePas = false;
};
if (nexistePas)
{Tab_Grandeur_Vecteur tab(v,1);
TypeQuelconque typQ2(FCT_POTENTIEL_ND,enu_ddl_type_pt,tab);
liTQ.push_back(typQ2);
};
};
};
// calcul d'un module d'young équivalent à la loi: pour l'instant à faire !!
double Hyper_externe_W::Module_young_equivalent(Enum_dure temps,const Deformation & ,SaveResul * )
{// on met un message d'erreur
cout << "\n *** erreur, pour l'instant le module d'Young n'est pas calcule !! "
<< "\n Hyper_externe_W::Module_young_equivalent(... "<<flush;
Sortie(1);
return 0.;
};
// 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 Hyper_externe_W::Module_compressibilite_equivalent(Enum_dure temps,const Deformation & def,SaveResul * saveResul)
{ SaveResulHyper_externe_W& sav = *((SaveResulHyper_externe_W*) saveResul);
return sav.module_compressibilite;
};
// ========== codage des METHODES VIRTUELLES protegees:================
// calcul des contraintes a t+dt
void Hyper_externe_W::Calcul_SigmaHH (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl,
TenseurBB & gijBB_t,TenseurHH & gijHH_t,BaseB& giB,BaseH& gi_H,TenseurBB & epsBB_,
TenseurBB & ,
TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,
double& jacobien_0,double& jacobien,TenseurHH & sigHH_
,EnergieMeca & energ,const EnergieMeca & ,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Expli_t_tdt& ex)
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 3)
{ cout << "\n Hyper_externe_W::Calcul_SigmaHH: avant critere ";
cout << "\n ele= "<< (Ptintmeca_en_cours()->Nb_ele()) << ", nbpti= "<< (Ptintmeca_en_cours()->Nb_pti());
};
if (epsBB_.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hyper_externe_W::Calcul_SigmaHH\n";
Sortie(1);
};
if (tab_ddl.NbDdl() != d_gijBB_.Taille())
{ cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_ !\n";
cout << " Hyper_externe_W::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage explicite en tenseur dim 3
// calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D)
Invariants_et_var1(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_,gijHH_,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
// le potentiel: il faut calculer la fonction nD
// 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_umat = NULL;
Potentiel_et_var2(ex_impli,ex_expli_tdt,ex_umat);
// calcul du tenseur des contraintes
sigHH = ((W_d_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3)/V) * d_J_r_epsBB_HH(3);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{TenseurHH* ptHH = NevezTenseurHH(sigHH);
sigHH.BaseAbsolue(*ptHH,giB);
cout << "\n base giB(1) "<< giB(1)
<< "\n base giB(2) "<< giB(2);
cout << "\n sigHH= en giB "; sigHH.Ecriture(cout);
cout << "\n sigma en absolu: ";
ptHH->Ecriture(cout);
delete ptHH;
};
#endif
// calcul du module de compressibilité
module_compressibilite = 2. * W_v_J3;
SaveResulHyper_externe_W& sav = *((SaveResulHyper_externe_W*) saveResul);
sav.module_compressibilite = module_compressibilite;
// pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
module_cisaillement = 0.;
// traitement des énergies
energ.Inita(0.);
energ.ChangeEnergieElastique((W_d)/V);
// energ.ChangeEnergieElastique((W_d+W_v)/V);
LibereTenseur();
};
// calcul des contraintes a t+dt et de ses variations
void Hyper_externe_W::Calcul_DsigmaHH_tdt (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl
,BaseB& giB_t,TenseurBB & gijBB_t,TenseurHH & gijHH_t
,BaseB& giB_tdt,Tableau <BaseB> & d_giB_tdt,BaseH& giH_tdt,Tableau <BaseH> & d_giH_tdt
,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,TenseurBB &
,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt_
,Tableau <TenseurBB *>& d_gijBB_tdt
,Tableau <TenseurHH *>& ,double& jacobien_0,double& jacobien
,Vecteur& ,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
,EnergieMeca & energ,const EnergieMeca & ,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Impli& ex )
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 3)
{ cout << "\n Hyper_externe_W::Calcul_DsigmaHH_tdt: ";
cout << "\n ele= "<< (Ptintmeca_en_cours()->Nb_ele()) << ", nbpti= "<< (Ptintmeca_en_cours()->Nb_pti());
};
if (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hyper_externe_W::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 << " Hyper_externe_W::Calcul_SDsigmaHH_tdt\n";
Sortie(1);
};
#endif
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // passage en dim 3 explicite
Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // passage en dim 3 explicite
// calcul des invariants et de leurs variations premières et secondes
Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
// le potentiel: il faut calculer la fonction nD
// 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_umat = NULL;
Potentiel_et_var2(ex_impli,ex_expli_tdt,ex_umat);
// calcul du tenseur des contraintes
double unSurV=1./V;
sigHH = ((W_d_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3)/V) * d_J_r_epsBB_HH(3);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{TenseurHH* ptHH = NevezTenseurHH(sigHH);
sigHH.BaseAbsolue(*ptHH,giB_tdt);
cout << "\n base giB(1) "<< giB_tdt(1)
<< "\n base giB(2) "<< giB_tdt(2);
cout << "\n sigHH= en giB "; sigHH.Ecriture(cout);
cout << "\n sigma en absolu: ";
ptHH->Ecriture(cout);
delete ptHH;
if (Permet_affichage() > 5)
cout << "\n d_J_r_epsBB_HH1 " << d_J_r_epsBB_HH(1)
<< "\n d_J_r_epsBB_HH2 " << d_J_r_epsBB_HH(2)
<< "\n d_J_r_epsBB_HH3 " << d_J_r_epsBB_HH(3)
<< "\n (W_d_J1*unSurV) " << (W_d_J1*unSurV)
<< " (W_d_J2*unSurV) " << (W_d_J2*unSurV)
<< " (W_d_J3*unSurV) " << (W_v_J3*unSurV);
};
#endif
// calcul de la variation seconde du potentiel par rapport à epsij epskl
Tenseur3HHHH d2W_d2epsHHHH
= // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
(W_d_J1_2 ) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
+ W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
+ W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
// puis les dérivées premières du potentiel déviatoire + courbure éventuellement
+ (W_d_J1 ) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
// enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
+ (W_v_J3 ) * d_J_3_eps2BB_HHHH
+ (W_v_J3J3 ) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV);
// Tenseur3HHHH dSigdepsHHHH(1,interHH,d_J_r_epsBB_HH(3));
Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt);
// dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH;
Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH); // cas des tenseurs généraux
dSigdepsHHHH += interHHHH; // cas des tenseurs généraux
//---------------------------------------------------------------------------------
// vérif numérique de l'opérateur tangent
// Cal_dsigma_deps_num (*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt
// ,gijHH_tdt,jacobien_0,jacobien,dSigdepsHHHH
// ,const Met_abstraite::Impli& ex );
//---------------------------------------------------------------------------------
// calcul des variations / aux ddl
int nbddl = d_gijBB_tdt.Taille();
for (int i = 1; i<= nbddl; i++)
{ // on fait uniquement une égalité d'adresse et de ne pas utiliser
// le constructeur d'ou la profusion d'* et de ()
Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i))); // passage en dim 3
const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); // "
dsigHH = dSigdepsHHHH && depsBB;
};
// calcul du module de compressibilité
module_compressibilite = 2. * W_v_J3;
SaveResulHyper_externe_W& sav = *((SaveResulHyper_externe_W*) saveResul);
sav.module_compressibilite = module_compressibilite;
// pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
module_cisaillement = 0.;
// traitement des énergies
energ.Inita(0.);
energ.ChangeEnergieElastique((W_d+W_v)/V);
LibereTenseur();
LibereTenseurQ();
};
// 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 Hyper_externe_W::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & ,TenseurBB&
,TenseurBB & epsBB_tdt,TenseurBB &, double& jacobien_0,double& jacobien
,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps_
,EnergieMeca & energ,const EnergieMeca & ,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Umat_cont& ex )
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 3)
{ cout << "\n Hyper_externe_W::Calcul_dsigma_deps: ";
cout << "\n ele= "<< (Ptintmeca_en_cours()->Nb_ele()) << ", nbpti= "<< (Ptintmeca_en_cours()->Nb_pti());
};
if (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hyper_externe_W::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // // passage en dim 3 explicite
Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) ex.gijHH_tdt); // passage en dim 3 explicite
// calcul des invariants et de leurs variations premières et secondes
Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),*(ex.gijBB_tdt),gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
// le potentiel: il faut calculer la fonction nD
// opération de transmission de la métrique
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_umat = &ex;
Potentiel_et_var2(ex_impli,ex_expli_tdt,ex_umat);
// calcul du tenseur des contraintes, on travaille ici dans le repère matériel finale correspondant
// aux coordonnées initiales X_(0)^a, on obtient donc un tenseur dans la base naturelle finale
double unSurV=1./V;
Tenseur3HH sig_localeHH = ((W_d_J1)*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2)
+ ((W_v_J3)*unSurV) * d_J_r_epsBB_HH(3);
// passage éventuelle dans la base I_a
if (en_base_orthonormee)
{sig_localeHH.BaseAbsolue(sigHH,*(ex.giB_tdt));}
else {sigHH = sig_localeHH;}; // sinon la base locale est la bonne
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{cout << "\n sigHH en local" << sigHH;
Tenseur3HH sig_3D_inter;
sigHH.BaseAbsolue(sig_3D_inter,*(ex.giB_tdt));
cout << "\n sigHH en absolue 3D : "<<sig_3D_inter;
if (Permet_affichage() > 5)
cout << "\n d_J_r_epsBB_HH1 " << d_J_r_epsBB_HH(1)
<< "\n d_J_r_epsBB_HH2 " << d_J_r_epsBB_HH(2)
<< "\n d_J_r_epsBB_HH3 " << d_J_r_epsBB_HH(3)
<< "\n (W_d_J1*unSurV) " << (W_d_J1*unSurV)
<< " (W_d_J2*unSurV) " << (W_d_J2*unSurV)
<< " (W_d_J3*unSurV) " << (W_v_J3*unSurV);
};
#endif
// calcul de la variation seconde du potentiel par rapport à epsij epskl
// calcul de la variation seconde du potentiel par rapport à epsij epskl
Tenseur3HHHH d2W_d2epsHHHH
= // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
(W_d_J1_2 ) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
+ W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
+ W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
// puis les dérivées premières du potentiel déviatoire + courbure éventuellement
+ (W_d_J1 ) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
// enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
+ (W_v_J3 ) * d_J_3_eps2BB_HHHH
+ (W_v_J3J3 ) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HH interHH = -sig_localeHH ; //* (-0.5*unSurV*unSurV);
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt);
Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH); // cas des tenseurs généraux
dSigdepsHHHH += interHHHH; // cas des tenseurs généraux
// transfert des informations: on pas d'un tenseur de 81 composantes à 36
// avec des symétries par rapport aux deux premiers indices et par rapport aux deux derniers
/// Tenseur3HHHH d_sigma_depsHHHH; d_sigma_depsHHHH.TransfertDunTenseurGeneral(dSigdepsHHHH.Symetrise1et2_3et4());
// calcul de la première partie de l'opérateur tangent (correspond au changement de repère
// gi_tdt -> Ia de l'opérateur calculer précédemment
Tenseur3HHHH & d_sigma_depsFinHHHH = *((Tenseur3HHHH*) &d_sigma_deps_); // pour accés directe
// passage éventuelle dans la base I_a
if (en_base_orthonormee)
{dSigdepsHHHH.ChangeBase(d_sigma_depsFinHHHH,*(ex.giB_tdt));}
else
{d_sigma_depsFinHHHH = dSigdepsHHHH;};
// calcul du module de compressibilité
module_compressibilite = 2. * W_v_J3;
SaveResulHyper_externe_W& sav = *((SaveResulHyper_externe_W*) saveResul);
sav.module_compressibilite = module_compressibilite;
// pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
module_cisaillement = 0.;
// traitement des énergies
energ.Inita(0.);
energ.ChangeEnergieElastique((W_d)/V);
// energ.ChangeEnergieElastique((W_d+W_v)/V);
LibereTenseur();
LibereTenseurQ();
};
//---------------------- méthodes privées -------------------------------
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
void Hyper_externe_W::Potentiel_et_var2
(const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat)
{ // calcul de grandeurs intermédiaires
// le potentiel: il faut calculer la fonction nD
// ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = W_potentiel->Li_enu_etendu_scalaire();
// on regarde s'il faut renseigner les invariants de B et/ou de J
List_io<Ddl_etendu> deja_calculer_etend;
{List_io<Ddl_enum_etendu>::iterator ipq,ipqfin=li_enu_scal.end();
for (ipq=li_enu_scal.begin();ipq!=ipqfin;ipq++)
{int posi = (*ipq).Position()-NbEnum_ddl();
switch (posi)
{
case 131 :
{deja_calculer_etend.push_back(Ddl_etendu((*ipq),I_B));
break;
}
case 132 :
{deja_calculer_etend.push_back(Ddl_etendu((*ipq),II_B));
break;
}
case 133 :
{deja_calculer_etend.push_back(Ddl_etendu((*ipq),III_B));
break;
}
case 134 :
{deja_calculer_etend.push_back(Ddl_etendu((*ipq),J_r(1)));
break;
}
case 135 :
{deja_calculer_etend.push_back(Ddl_etendu((*ipq),J_r(2)));
break;
}
case 136 :
{deja_calculer_etend.push_back(Ddl_etendu((*ipq),J_r(3)));
break;
}
default: // on ne fait rien
break;
};
};
}
// on utilise la méthode générique de loi abstraite
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(W_potentiel,9 // 9 valeurs attendue en retour
,ex_impli,ex_expli_tdt,ex_umat
,&deja_calculer_etend
,NULL
,NULL
);
/*
// Tableau <double> & Loi_comp_Valeur_FnD_Evoluee
// (Fonction_nD* fct,int nb_retour
// ,const Met_abstraite::Impli* ex_impli
// ,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
// ,const Met_abstraite::Umat_cont* ex_umat
// ,const List_io<Ddl_enum_etendu>* exclure_dd_etend = NULL
// ,const List_io<EnumTypeQuelconque>* exclure_Q = NULL
// ,list <SaveResul*>* list_save = 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_umat,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = W_potentiel->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 9)
{ cout << "\nErreur : la fonction nD relative a la fonction potentielle "
<< " doit calculer un vecteur de dimention 9 or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " Hyper_externe_W::Potentiel_et_var\n";
Sortie(1);
};
#endif
*/
// on récupère les valeurs du tableau
W_d= tab_val(1); W_v= tab_val(2); // le potentiel: partie déviatorique, partie sphérique
W_d_J1= tab_val(3);W_d_J2= tab_val(4); // dérivées premières du potentiel déviatoire par rapport aux J_1 et J_2
W_d_J1_2= tab_val(5);W_d_J1_J2= tab_val(6);W_d_J2_2= tab_val(7);; // dérivées secondes
W_v_J3= tab_val(8);W_v_J3J3= tab_val(9); // dérivées premières et seconde du potentiel volumique / J3
// stockage éventuel pour du post-traitement
if (sortie_post)
{ // récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
save_resulHyper_W.invP->potentiel= W_d + W_v;
// puis les nouvelles variables
SaveResulHyper_externe_W& sav = *((SaveResulHyper_externe_W*) saveResul);
Vecteur& inv__B = *sav.inv_B; // pour simplifier
inv__B(1) = I_B;
inv__B(2) = II_B;
inv__B(3) = III_B;
Vecteur& inv__J = *sav.inv_J; // pour simplifier
inv__J = J_r;
// cout << "\n inv__B= "<<inv__B << " inv__J= "<<inv__J<<flush;
Vecteur& W__poten = *sav.W_poten; // pour simplifier
W__poten(1) = W_d;
W__poten(2) = W_v;
W__poten(3) = W_d_J1;
W__poten(4) = W_d_J2;
W__poten(5) = W_d_J1_2;
W__poten(6) = W_d_J1_J2;
W__poten(7) = W_d_J2_2;
W__poten(8) = W_v_J3;
W__poten(9) = W_v_J3J3;
};
};
// calcul de la dérivée numérique de la contrainte
void Hyper_externe_W::Cal_dsigma_deps_num
(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
,const double& jacobien_0,const double& jacobien
,Tenseur3HHHH& dSigdepsHHHH
,const Met_abstraite::Impli& ex )
{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); // "
const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); // "
const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // "
Tenseur3BB gijBBtdt_N; // tenseur modifié
Tenseur3HH gijHHtdt_N; // idem_0
double delta = ConstMath::unpeupetit*10.;
double unSurDelta = 1./delta;
// cas des contraintes et de ses variations analytiques
// Tenseur3HHHH dSigdepsHHHH; // le tenseur contenant les dérivées analytiques
Tenseur3HH SigmaHH_deb;
Cal_sigmaEtDer_pour_num(gijBB_0_,gijHH_0_,gijBB_tdt_,gijHH_tdt_
,jacobien_0,jacobien,SigmaHH_deb,dSigdepsHHHH,ex);
// dimensionnement pour la matrice numérique
Tenseur3HHHH dSigdepsHHHH_num;
// on va boucler sur les composantes de gijBB
for (int i=1;i<=3;i++)
for (int j=1;j<=3;j++)
{ gijBBtdt_N = gijBB_tdt;
gijBBtdt_N.Coor(i,j) += delta;
// en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i)
// car le tenseur est symétrique
// on a donc en variation numérique la somme des deux dérivées
// on définit un coeff multiplicatif qui vaut 1 ou 0.5
double coef=1.; if (i != j) coef = 0.5;
gijHHtdt_N = gijBBtdt_N.Inverse();
double jacobien_N=sqrt(gijBBtdt_N.Det());
// cas des contraintes
Tenseur3HH SigmaHH_N;
Cal_sigma_pour_num(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N
,(const TenseurHH &) gijHHtdt_N
,jacobien_0,jacobien_N,SigmaHH_N,ex);
// calcul des dérivées numériques et comparaisons
for (int k=1;k<=3;k++)
for (int l=1;l<=3;l++)
{ //
double derSigNum = coef * 2.*(SigmaHH_N(k,l) - SigmaHH_deb(k,l) )*unSurDelta;
dSigdepsHHHH_num.Change(k,l,i,j,derSigNum);
double derSigAna = dSigdepsHHHH(k,l,i,j);//0.5*(dSigdepsHHHH(k,l,i,j) + dSigdepsHHHH(k,l,j,i));
bool erreur = false;
if (diffpourcent(derSigNum,derSigAna,MaX(Dabs(derSigNum),Dabs(derSigAna)),0.1))
if (MaX(Dabs(derSigNum),Dabs(derSigAna)) > 200.)
{if (MiN(Dabs(derSigNum),Dabs(derSigAna)) == 0.)
{if ( MaX(Dabs(derSigNum),Dabs(derSigAna)) > 50.*delta) erreur = true;}
else erreur = true;
};
// erreur = false; // a virer
if (erreur)
{
// calcul des dérivées d'éléments intermédiaires pour voir
//
cout << "\n erreur dans le calcul analytique de l'operateur tangent "
<< "\n derSigNum= " << derSigNum << " derSigAna= " << derSigAna
<< " klij= " << k << " " << l << " " << i << " " << j
<< " SigmaHH_N(k,l)= " << SigmaHH_N(k,l);
cout << "\n Hyper_externe_W::Calcul_derivee_numerique(..";
cout << "\n un caractere ";
// --- pour le débug ----
// calcul des invariants et de leurs variations premières en numérique
Invariants_et_var1_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
,jacobien_0,jacobien_N);
// calcul des invariants et de leurs variations premières et secondes
Invariants_et_var2_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
,jacobien_0,jacobien_N);
string toto;
toto= lect_chaine();
};
};
};
// passage des dérivées numériques aux dérivées finales
dSigdepsHHHH= dSigdepsHHHH_num;
};
// calcul de la contrainte avec le minimum de variable de passage, utilisé pour le numérique
void Hyper_externe_W::Cal_sigma_pour_num
(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0
,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt
,const double& jacobien_0,const double& jacobien,TenseurHH & sigHH_
,const Met_abstraite::Impli& ex)
{
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage en dim 3 explicite
// calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D)
// Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// pour vérif on appelle var2, mais c'est à virer
Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
// le potentiel: il faut calculer la fonction nD
// 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_umat = NULL;
Potentiel_et_var2(ex_impli,ex_expli_tdt,ex_umat);
// calcul du tenseur des contraintes
double unSurV=1./V;
sigHH = ((W_d_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3)/V) * d_J_r_epsBB_HH(3);
};
// idem avec la variation
void Hyper_externe_W::Cal_sigmaEtDer_pour_num
(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0
,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt
,const double& jacobien_0,const double& jacobien
,TenseurHH & sigHH_,Tenseur3HHHH& dSigdepsHHHH
,const Met_abstraite::Impli& ex)
{
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage en dim 3 explicite
// calcul des invariants et de leurs variations premières et seconde (méthode de Hyper_W_gene_3D)
Invariants_et_var2(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
// le potentiel: il faut calculer la fonction nD
// 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_umat = NULL;
Potentiel_et_var2(ex_impli,ex_expli_tdt,ex_umat);
// calcul du tenseur des contraintes
double unSurV=1./V;
sigHH = ((W_d_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3)/V) * d_J_r_epsBB_HH(3);
// calcul de la variation seconde du potentiel par rapport à epsij epskl
Tenseur3HHHH d2W_d2epsHHHH
= // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
(W_d_J1_2 ) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
+ W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
+ W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
// puis les dérivées premières du potentiel déviatoire + courbure éventuellement
+ (W_d_J1 ) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
// enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
+ (W_v_J3 ) * d_J_3_eps2BB_HHHH
+ (W_v_J3J3 ) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV);
Tenseur3HHHH d_igdepsHHHH(1,interHH,gijHH_tdt);
d_igdepsHHHH += (unSurV) * d2W_d2epsHHHH;
dSigdepsHHHH = d_igdepsHHHH;
};
Reference in a new issue Copy permalink