// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
#include "IsoHyper3DFavier3.h"
//#include "ComLoi_comp_abstraite.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "ConstMath.h"
#include "CharUtil.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"
//================== initialisation des variables static ======================
// indicateur utilisé par Verif_Potentiel_et_var
int IsoHyper3DFavier3::indic_Verif_PoGrenoble_et_var = 0;
double IsoHyper3DFavier3::limite_co2=700.; // limite à partir de laquelle on considère que cosh(co2) = infini
// ici cosh(700.) = 5.0712e+303 !!
//================== fin d'initialisation des variables static ================
//================== stockage particulier SaveResul ================
IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::SaveResulIsoHyper3DFavier3(bool avec_para): // constructeur par défaut :
SaveResulHyperD(),para_loi(NULL),para_loi_t(NULL)
{if (avec_para)
{para_loi = new Vecteur(4,-ConstMath::trespetit);
para_loi_t = new Vecteur(4,-ConstMath::trespetit);
};
};
IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::SaveResulIsoHyper3DFavier3(const SaveResulIsoHyper3DFavier3& sav): // de copie
SaveResulHyperD(sav),para_loi(NULL),para_loi_t(NULL)
{if (sav.para_loi != NULL)
{para_loi = new Vecteur(*sav.para_loi);
para_loi_t = new Vecteur(*sav.para_loi_t);
};
};
IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::~SaveResulIsoHyper3DFavier3() // destructeur
{ if (para_loi != NULL) delete para_loi;
if (para_loi_t != NULL) delete para_loi_t;
};
// affectation
Loi_comp_abstraite::SaveResul & IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3
::operator = ( const Loi_comp_abstraite::SaveResul & a)
{ SaveResulHyperD::operator=(a);
SaveResulIsoHyper3DFavier3& sav = *((SaveResulIsoHyper3DFavier3*) &a);
if (sav.para_loi != NULL)
{if (para_loi == NULL)
{para_loi = new Vecteur(4,-ConstMath::trespetit);
para_loi_t = new Vecteur(4,-ConstMath::trespetit);
};
*para_loi = *(sav.para_loi);
*para_loi_t = *(sav.para_loi_t);
};
map_type_quelconque = sav.map_type_quelconque;
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 IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::Lecture_base_info (istream& ent,const int cas)
{SaveResulHyperD::Lecture_base_info (ent,cas);
string toto;
ent >> toto;
if (toto == "sans_para")
{if (para_loi != NULL)
{delete para_loi;para_loi=NULL;
delete para_loi_t;para_loi_t=NULL;
};
}
else
{if (para_loi == NULL)
{para_loi = new Vecteur(4);
para_loi_t = new Vecteur(4);
};
ent >> *para_loi;
*para_loi_t=*para_loi;
};
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::Ecriture_base_info(ostream& sort,const int cas)
{SaveResulHyperD::Ecriture_base_info (sort,cas);
if (para_loi == NULL)
{sort << " sans_para ";}
else
{sort << " avec_para ";
sort << *para_loi;
};
};
// mise à jour des informations transitoires
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::TdtversT()
{SaveResulHyperD::TdtversT();
if (para_loi != NULL)
{*para_loi_t = *para_loi;};
// mise à jour de la liste des grandeurs quelconques internes
Mise_a_jour_map_type_quelconque();
};
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::TversTdt()
{SaveResulHyperD::TversTdt();
if (para_loi != NULL)
{*para_loi = *para_loi_t;};
// mise à jour de la liste des grandeurs quelconques internes
Mise_a_jour_map_type_quelconque();
};
// affichage à l'écran des infos
void IsoHyper3DFavier3::SaveResulIsoHyper3DFavier3::Affiche() const
{SaveResulHyperD::Affiche();
if (para_loi == NULL)
{cout << " sans_para ";}
else
{cout << " avec_para ";
cout << *para_loi;
};
};
//================== fin stockage particulier SaveResul ================
IsoHyper3DFavier3::IsoHyper3DFavier3 () : // Constructeur par defaut
Hyper3D(ISOHYPER3DFAVIER3,CAT_MECANIQUE,false),K(ConstMath::trespetit),Qor(ConstMath::trespetit)
,mur(ConstMath::trespetit),mu_inf(ConstMath::trespetit)
,nQor(ConstMath::trespetit),gammaQor(ConstMath::trespetit)
,n_mu_inf(ConstMath::trespetit),gamma_mu_inf(ConstMath::trespetit)
// dépendance nD éventuelle
,K_nD(NULL),Qor_nD(NULL),mur_nD(NULL),mu_inf_nD(NULL)
{};
// Constructeur de copie
IsoHyper3DFavier3::IsoHyper3DFavier3 (const IsoHyper3DFavier3& loi) :
Hyper3D (loi),K(loi.K),Qor(loi.Qor),mur(loi.mur),mu_inf(loi.mu_inf)
,nQor(loi.nQor),gammaQor(loi.gammaQor)
,n_mu_inf(loi.n_mu_inf),gamma_mu_inf(loi.gamma_mu_inf)
// dépendance nD éventuelle
,K_nD(loi.K_nD),Qor_nD(loi.Qor_nD),mur_nD(loi.mur_nD),mu_inf_nD(loi.mu_inf_nD)
{// on regarde s'il s'agit d'une fonction locale ou d'une fonction globale
if (K_nD != NULL) // 1)
if (K_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
K_nD = Fonction_nD::New_Fonction_nD(*loi.K_nD);
};
if (Qor_nD != NULL) // 2)
if (Qor_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
Qor_nD = Fonction_nD::New_Fonction_nD(*loi.Qor_nD);
};
if (mur_nD != NULL) // 3)
if (mur_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
mur_nD = Fonction_nD::New_Fonction_nD(*loi.mur_nD);
};
if (mu_inf_nD != NULL) // 4)
if (mu_inf_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
mu_inf_nD = Fonction_nD::New_Fonction_nD(*loi.mu_inf_nD);
};
// --- on remplit avec les grandeurs succeptible d'être utilisées
// // acces à listdeTouslesQuelc_dispo_localement
// list <EnumTypeQuelconque >& list_tousQuelc = ListdeTouslesQuelc_dispo_localement();
// list_tousQuelc.push_back(E_YOUNG);
// list_tousQuelc.push_back(NU_YOUNG);
// // on supprime les doublons localement
// list_tousQuelc.sort(); // on ordonne la liste
// list_tousQuelc.unique(); // suppression des doublons
};
// Lecture des donnees de la classe sur fichier
void IsoHyper3DFavier3::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ // --- lecture des quatres coefficients de la loi
string nom_class_methode("IsoHyper3DFavier3::LectureDonneesParticulieres");
// on regarde si K dépend d'une fonction nD
if(strstr(entreePrinc->tablcar,"K_fonction_nD:")!=0)
{ string nom;
string mot_cle2="K_fonction_nD:";
// on lit le nom de la fonction
string nom_fonct;
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer("**erreur01 en lecture** "+mot_cle2);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{K_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
K_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
K_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (K_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << K_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur02** \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 = K_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
}
// sinon si le module K n'a pas été lue, on le récupère en version scalaire
else
{ // lecture du module
*(entreePrinc->entree) >> K ;
};
// idem si Qor dépend d'une fonction nD
if(strstr(entreePrinc->tablcar,"Qor_fonction_nD:")!=0)
{ string nom;
string mot_cle2="Qor_fonction_nD:";
// on lit le nom de la fonction
string nom_fonct;
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer("**erreur03 en lecture** "+mot_cle2);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{Qor_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
Qor_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
Qor_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (Qor_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << Qor_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur04** \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 = Qor_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
}
// sinon si le module Qor n'a pas été lue, on le récupère en version scalaire
else
{ // lecture du module
*(entreePrinc->entree) >> Qor ;
};
// idem si mur dépend d'une fonction nD
if(strstr(entreePrinc->tablcar,"mur_fonction_nD:")!=0)
{ string nom;
string mot_cle2="mur_fonction_nD:";
// on lit le nom de la fonction
string nom_fonct;
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer("**erreur05 en lecture** "+mot_cle2);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{mur_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
mur_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
mur_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (mur_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << mur_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur06** \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 = mur_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
}
// sinon si le module mur n'a pas été lue, on le récupère en version scalaire
else
{ // lecture du module
*(entreePrinc->entree) >> mur ;
};
// idem si mu_inf dépend d'une fonction nD
if(strstr(entreePrinc->tablcar,"mu_inf_fonction_nD:")!=0)
{ string nom;
string mot_cle2="mu_inf_fonction_nD:";
// on lit le nom de la fonction
string nom_fonct;
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer("**erreur07 en lecture** "+mot_cle2);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{mu_inf_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
mu_inf_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
mu_inf_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (mu_inf_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << mu_inf_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur08** \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 = mu_inf_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
}
// sinon si le module mu_inf n'a pas été lue, on le récupère en version scalaire
else
{ // lecture du module
*(entreePrinc->entree) >> mu_inf ;
};
// *(entreePrinc->entree) >> K >> Qor >> mur >> mu_inf ;
// on regarde ensuite s'il y a la phase
if (strstr(entreePrinc->tablcar,"avec_phase")!=NULL)
{ // lecture des 4 paramètres de phase
entreePrinc->NouvelleDonnee(); // passage d'une ligne
*(entreePrinc->entree) >> nQor >> gammaQor >> n_mu_inf >> gamma_mu_inf ;
avec_phase=true;
};
// Lecture des paramètre specifique sur fichier (exe: regularisation sortie_post...
Hyper3D::LectParaSpecifiques(entreePrinc,lesCourbes1D,lesFonctionsnD);
// appel au niveau de la classe mère
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD);
};
// affichage de la loi
void IsoHyper3DFavier3::Affiche() const
{ cout << " \n loi de comportement 3D hyperelastique isotrope favier3 : " << Nom_comp(id_comp)
<< " paramètres : \n";
if (K_nD != NULL)
{cout << "\n K fonction nD: ";
cout << "K_fonction_nD:" << " ";
if (K_nD->NomFonction() != "_")
cout << K_nD->NomFonction();
else
K_nD->Affiche();
}
else
{ cout << " K= " << K ;};
if (Qor_nD != NULL)
{cout << "\n Qor fonction nD: ";
cout << "Qor_fonction_nD:" << " ";
if (Qor_nD->NomFonction() != "_")
cout << Qor_nD->NomFonction();
else
Qor_nD->Affiche();
}
else
{ cout << " Qor= " << Qor ;};
if (mur_nD != NULL)
{cout << "\n mur fonction nD: ";
cout << "mur_fonction_nD:" << " ";
if (mur_nD->NomFonction() != "_")
cout << mur_nD->NomFonction();
else
mur_nD->Affiche();
}
else
{ cout << " mur_nD= " << mur_nD ;};
if (mu_inf_nD != NULL)
{cout << "\n mu_inf fonction nD: ";
cout << "mu_inf_fonction_nD:" << " ";
if (mu_inf_nD->NomFonction() != "_")
cout << mu_inf_nD->NomFonction();
else
mu_inf_nD->Affiche();
}
else
{ cout << " mu_inf= " << mu_inf ;};
// cout << " K= " << K << " ; Qor = " << Qor << " ; mur = " << mur << " ; mu_inf = " << mu_inf ;
if (avec_phase)
{ cout << "\n cas de la loi avec phase, parametre de la phase: "
<< " nQor= " << nQor << " gammaQor= " << gammaQor << " n_mu_inf= "
<< n_mu_inf << " gamma_mu_inf= " << gamma_mu_inf;
}
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 IsoHyper3DFavier3::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# ....... loi de comportement 3D hyperelastique isotrope favier3 ........";
if ((rep != "o") && (rep != "O" ) && (rep != "0") )
{ sort
<< "\n# **** exemple de cas AVEC la phase :"
<< "\n#------------------------------------------------------------------------------------"
<< "\n# K | Qor | mur | mu_inf | avec_phase(falculatif) |"
<< "\n#------------------------------------------------------------------------------------"
<< "\n 160000 150 17000 220 avec_phase"
<< "\n#-------------------------------------------------"
<< "\n# n_Q | gamma_Q | n_mu | gamma_mu |"
<< "\n#-------------------------------------------------"
<< "\n 0.25 0.4 0.25 0.4 "
<< "\n# **** exemple de cas AVEC utilisation de fonction nD :"
<< "\n#------------------------------------------------------------------------------------"
<< "\n# K | Qor | mur | mu_inf | avec_phase(falculatif) |"
<< "\n#------------------------------------------------------------------------------------"
<< "\n# 16000 Qor_fonction_nD: fct1 17000 mu_inf_fonction_nD: fct2"
<< "\n#"
<< "\n# une fonction nD est possible pour chacun des 4 coefficients, via les mots cle "
<< "\n# K_fonction_nD: "
<< "\n# Qor_fonction_nD: "
<< "\n# mur_fonction_nD: "
<< "\n# mu_inf_fonction_nD: "
<< "\n# on peut utiliser 1 ou 2 ou 3 ou 4 fct nD, au choix "
<< "\n# ( bien noter que la loi obtenue peut-etre quelconque, en particulier plus du tout "
<< "\n# hyperelastique, tout depend des choix de fct nD !) "
<< "\n#"
<< "\n#------------------------------------------------------------------------------------"
<< "\n#------------------------------------------------------------------------------------"
<< "\n# il est possible d'indiquer 2 parametres specifiques de limite inf "
<< "\n# limite_inf_Qeps_ : si Qeps < limite_inf_Qeps on considere que l'angle de phase et ses derivees sont nulles "
<< "\n# limite_inf_bIIb_ : si Dabs(inv.bIIb) < limite_inf_bIIb "
<< "\n# les derivees du potentiel sont calculees pas difference fini, sauf ceux relatif a la phase qui sont mises a 0 "
<< "\n# les mots cle limite_inf_Qeps_ suivi du facteur voulu et limite_inf_bIIb_ suivi du facteur "
<< "\n# doivent etre dans cet ordre "
<< "\n# ex: "
<< "\n# limite_inf_Qeps_ 8.5e-5 limite_inf_bIIb_ 36.e-10 "
<< "\n# ces mots cle doivent se situer avant le mot cle avec_regularisation_ "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# il est possible d'indiquer un facteur de regularisation qui permet d'eviter "
<< "\n# de potentiels problemes de NaN, de type division par 0 par exemple "
<< "\n# 1/a est remplace par 1/(a+fact_regularisation), par defaut fact_regularisation = 1.e-12 "
<< "\n# pour indiquer un facteur de regulation non nul on indique en dernier parametre "
<< "\n# le mot cle avec_regularisation_ suivi du facteur voulu "
<< "\n# ex: "
<< "\n# avec_regularisation_ 1.e-12 "
<< "\n# ce mot cle doit se situer avant le mot cle sortie_post_ "
<< "\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#------------------------------------------------------------------------------------";
};
sort << "\n# **** exemple de cas SANS la phase :"
<< "\n#------------------------------------------------------------------------------------"
<< "\n# K | Qor | mur | mu_inf | avec_phase(falculatif) |"
<< "\n#------------------------------------------------------------------------------------"
<< "\n 160000 150 17000 220 "
<< endl;
// appel de la classe Hyper3D
Hyper3D::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 IsoHyper3DFavier3::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
if ((K == ConstMath::trespetit) && (K_nD == NULL))
{ cout << " \n Le parametre K n'est pas defini pour la loi " << Nom_comp(id_comp)
<< '\n';
Affiche();
ret = 0;
}
if ((Qor == ConstMath::trespetit) && (Qor_nD == NULL))
{ cout << " \n Le parametre Qor n'est pas defini pour la loi " << Nom_comp(id_comp)
<< '\n';
Affiche();
ret = 0;
}
if ((mur == ConstMath::trespetit) && (mur_nD == NULL))
{ cout << " \n Le parametre mur n'est pas defini pour la loi " << Nom_comp(id_comp)
<< '\n';
Affiche();
ret = 0;
}
if ((mu_inf == ConstMath::trespetit) && (mu_inf_nD == NULL))
{ cout << " \n Le parametre mu_inf n'est pas defini pour la loi " << Nom_comp(id_comp)
<< '\n';
Affiche();
ret = 0;
}
if (avec_phase)
if ((nQor == ConstMath::trespetit) || (gammaQor == ConstMath::trespetit)
|| (n_mu_inf == ConstMath::trespetit) || (gamma_mu_inf == ConstMath::trespetit))
{ cout << " \n Les paramètres de la phase ne sont pas défini pour la loi " << Nom_comp(id_comp)
<< '\n';
Affiche();
ret = 0;
}
return ret;
};
// 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 IsoHyper3DFavier3::Grandeur_particuliere
(bool absolue,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list<int>& decal) const
{ // tout d'abord on appelle la class mère
HyperD::Grandeur_particuliere(absolue,liTQ,saveDon,decal);
// --- puis la loi elle-même
// 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 PARA_LOI_FAVIER:
// ----- cas des paramètres de la loi
{ SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveDon);
Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
if (save_resul.para_loi != NULL)
{tyTQ(1+(*idecal))= *(save_resul.para_loi);}
else {tyTQ(1+(*idecal)).Zero();};
(*idecal)++; break;
}
default: ;// on ne fait rien
};
};
};
// 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 IsoHyper3DFavier3::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
{ // tout d'abord on appelle la class mère
HyperD::ListeGrandeurs_particulieres(absolue,liTQ);
// --- puis la loi elle-même
// $$$ cas des paramètres de la loi de comportement
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == PARA_LOI_FAVIER)
{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)
{Vecteur tens(4); // les 4 paramètres
Tab_Grandeur_Vecteur gr_vec_para(tens,1);
// def d'un type quelconque représentatif
TypeQuelconque typQ(PARA_LOI_FAVIER,SIG11,gr_vec_para);
liTQ.push_back(typQ);
};
};
};
//----- 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 IsoHyper3DFavier3::Lecture_base_info_loi(istream& ent,const int cas,LesReferences& lesRef
,LesCourbes1D& lesCourbe1D,LesFonctions_nD& lesFonctionsnD)
{ string toto,nom;
if (cas == 1)
{ ent >> toto >> thermo_dependant >> toto;
// tout d'abord la lecture de K
int type =0;
ent >> type;
switch (type)
{ case 1 :
{ent >> nom;
if (nom != " K_fonction_nD: ")
{ cout << "\n erreur en lecture de la fonction nD, on attendait "
<< " K_fonction_nD: et on a lue " << nom
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
};
K_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,K_nD);
// on regarde si la fonction nD intègre la température
// car l'utilisateur peut éventuellement lors d'un restart changer
// de paramètres pour la fonction nD
const Tableau <Ddl_enum_etendu>& tab_enu = K_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
break;
};
case 2 :
{ent >> K;
break;
};
default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
<< " reconnus sont uniquement: 1 c-a-d K fonction nD, "
<< " 2 K une valeur fixe "
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
break;
};
// puis lecture de Qor
type =0;
ent >> toto >> type;
switch (type)
{ case 1 :
{ent >> nom;
if (nom != " Qor_fonction_nD: ")
{ cout << "\n erreur en lecture de la fonction nD, on attendait "
<< " Qor_fonction_nD: et on a lue " << nom
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
};
Qor_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,Qor_nD);
// on regarde si la fonction nD intègre la température
// car l'utilisateur peut éventuellement lors d'un restart changer
// de paramètres pour la fonction nD
const Tableau <Ddl_enum_etendu>& tab_enu = Qor_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
break;
};
case 2 :
{ent >> Qor;
break;
};
default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
<< " reconnus sont uniquement: 1 c-a-d Qor fonction nD, "
<< " 2 Qor une valeur fixe "
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
break;
};
// puis lecture de mur
type =0;
ent >> toto >> type;
switch (type)
{ case 1 :
{ent >> nom;
if (nom != " mur_fonction_nD: ")
{ cout << "\n erreur en lecture de la fonction nD, on attendait "
<< " mur_fonction_nD: et on a lue " << nom
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
};
mur_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mur_nD);
// on regarde si la fonction nD intègre la température
// car l'utilisateur peut éventuellement lors d'un restart changer
// de paramètres pour la fonction nD
const Tableau <Ddl_enum_etendu>& tab_enu = mur_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
break;
};
case 2 :
{ent >> mur;
break;
};
default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
<< " reconnus sont uniquement: 1 c-a-d mur fonction nD, "
<< " 2 mur une valeur fixe "
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
break;
};
// puis lecture de mu_inf
type =0;
ent >> toto >> type;
switch (type)
{ case 1 :
{ent >> nom;
if (nom != " mu_inf_fonction_nD: ")
{ cout << "\n erreur en lecture de la fonction nD, on attendait "
<< " mu_inf_fonction_nD: et on a lue " << nom
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
};
mu_inf_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,mu_inf_nD);
// on regarde si la fonction nD intègre la température
// car l'utilisateur peut éventuellement lors d'un restart changer
// de paramètres pour la fonction nD
const Tableau <Ddl_enum_etendu>& tab_enu = mu_inf_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
break;
};
case 2 :
{ent >> mu_inf;
break;
};
default: cout << "\n erreur type " << type << " non prevu, pour l'instant, les types "
<< " reconnus sont uniquement: 1 c-a-d mu_inf fonction nD, "
<< " 2 mu_inf une valeur fixe "
<< "\n IsoHyper3DFavier3::Lecture_base_info_loi(...";
Sortie(1);
break;
};
// ent >> toto >> K >> toto >> Qor >> toto >> mur >> toto >> mu_inf >> avec_phase;
// avec phase éventuelle
ent >> avec_phase;
if (avec_phase)
{ ent >> toto >> nQor >> toto >> gammaQor >> toto
>> n_mu_inf >> toto >> gamma_mu_inf;
};
// prise en compte éventuelle de paramètres particulier: ex: régularisation
// géré par HyperD et Hyper3D
Hyper3D::Lecture_para_specifiques(ent,cas);
// gestion du post-traitement
ent >> toto >> sortie_post ;
};
// appel class mère
Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbe1D,lesFonctionsnD);
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void IsoHyper3DFavier3::Ecriture_base_info_loi(ostream& sort,const int cas)
{ if (cas == 1)
{ sort << " ISOHYPER3DFAVIER3,thermodependance= " << thermo_dependant ;
sort << " module_K ";
if (K_nD != NULL)
{sort << " 1 K_fonction_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,K_nD);
}
else
{ sort << " 2 " << K ; };
sort << " seuil ";
if (Qor_nD != NULL)
{sort << " 1 Qor_fonction_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,Qor_nD);
}
else
{ sort << " 2 " << Qor ; };
sort << " pente_origine ";
if (mur_nD != NULL)
{sort << " 1 mur_fonction_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,mur_nD);
}
else
{ sort << " 2 " << mur ; };
sort << " pente_infini ";
if (mu_inf_nD != NULL)
{sort << " 1 mu_inf_fonction_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,mu_inf_nD);
}
else
{ sort << " 2 " << mu_inf ; };
// sort << " module_dilatation " << K << " seuil " << Qor
// << " pente_origine " << mur << " pente_infini " << mu_inf << " avec_phase " << avec_phase;
sort << " avec_phase " << avec_phase;
if (avec_phase)
{ sort << "\n nQor= " << nQor << " gammaQor= " << gammaQor << " n_mu_inf= "
<< n_mu_inf << " gamma_mu_inf= " << gamma_mu_inf;
};
// prise en compte éventuelle de paramètres particulier: ex: régularisation
// géré par HyperD et Hyper3D
Hyper3D::Ecriture_para_specifiques(sort,cas);
// gestion du post-traitement
sort << " sortie_post= "<< sortie_post << " ";
};
// appel de la classe mère
Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);
};
// calcul d'un module d'young équivalent à la loi, ceci pour un
// chargement nul si rien n'est calculé, sinon c'est
// une valeur correspondante au dernier calcul
double IsoHyper3DFavier3::Module_young_equivalent(Enum_dure temps,const Deformation & def,SaveResul * saveDon)
{ // 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
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
double K_sortie=K; // init
double mur_sortie=mur; // init
double mu_inf_sortie=mu_inf; // init
// traitement du cas particulier de l'existence de fct nD
if (save_resul.para_loi != NULL)
{bool recup = false;
switch (temps)
{
case TEMPS_t : // on utilise les grandeurs stockées à t
if ((*save_resul.para_loi_t)(1) != (-ConstMath::trespetit))
{K_sortie= (*save_resul.para_loi_t)(1);recup = true;}
break;
case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
if ((*save_resul.para_loi)(1) != (-ConstMath::trespetit))
{K_sortie= (*save_resul.para_loi)(1);recup = true;}
break;
case TEMPS_0 : // rien n'a été calculé
recup = false;
};
// dans tous les autres cas on utilise soit les fonctions
// si elles existent sinon on concerve les valeurs par défaut
if (!recup)
{//on essaie un calcul
// cas des courbes d'évolution
if (K_nD != NULL)
// là il faut calcul 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 = K_nD->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = K_nD->Li_equi_Quel_evolue();
// on initialise les grandeurs du tableau pour les valeurs numériques
Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
// init par défaut des types quelconques
List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
for (il = li_quelc.begin();il != ilfin; il++)
(*il).Grandeur_pointee()->InitParDefaut();
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val
= K_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 au module d'Young "
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
K_sortie = tab_val(1);
};
};
// -- idem pour mur
recup = false;
switch (temps)
{
case TEMPS_t : // on utilise les grandeurs stockées à t
if ((*save_resul.para_loi_t)(3) != (-ConstMath::trespetit))
{mur_sortie= (*save_resul.para_loi_t)(3);recup = true;}
break;
case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
if ((*save_resul.para_loi)(3) != (-ConstMath::trespetit))
{mur_sortie= (*save_resul.para_loi)(3);recup = true;}
break;
case TEMPS_0 : // rien n'a été calculé
recup = false;
};
// dans tous les autres cas on utilise soit les fonctions
// si elles existent sinon on concerve les valeurs par défaut
if (!recup)
{//on essaie un calcul
// cas des courbes d'évolution
if (mur_nD != NULL)
// là il faut calcul 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 = mur_nD->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = mur_nD->Li_equi_Quel_evolue();
// on initialise les grandeurs du tableau pour les valeurs numériques
Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
// init par défaut des types quelconques
List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
for (il = li_quelc.begin();il != ilfin; il++)
(*il).Grandeur_pointee()->InitParDefaut();
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val
= mur_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 au module d'Young "
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
mur_sortie = tab_val(1);
};
};
// -- idem pour mu_inf
recup = false;
switch (temps)
{
case TEMPS_t : // on utilise les grandeurs stockées à t
if ((*save_resul.para_loi_t)(4) != (-ConstMath::trespetit))
{mu_inf_sortie= (*save_resul.para_loi_t)(4);recup = true;}
break;
case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
if ((*save_resul.para_loi)(4) != (-ConstMath::trespetit))
{mu_inf_sortie= (*save_resul.para_loi)(4);recup = true;}
break;
case TEMPS_0 : // rien n'a été calculé
recup = false;
};
// dans tous les autres cas on utilise soit les fonctions
// si elles existent sinon on concerve les valeurs par défaut
if (!recup)
{//on essaie un calcul
// cas des courbes d'évolution
if (mu_inf_nD != NULL)
// là il faut calcul 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 = mu_inf_nD->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = mu_inf_nD->Li_equi_Quel_evolue();
// on initialise les grandeurs du tableau pour les valeurs numériques
Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
// init par défaut des types quelconques
List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
for (il = li_quelc.begin();il != ilfin; il++)
(*il).Grandeur_pointee()->InitParDefaut();
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val
= mu_inf_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 au module d'Young "
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
mu_inf_sortie = tab_val(1);
};
};
};
// retour du module
double module = (9.*K_sortie*(mur_sortie+mu_inf_sortie)/((mur_sortie+mu_inf_sortie)+3.*K_sortie));
return module;
};
// récupération d'un module de compressibilité équivalent à la loi, ceci pour un chargement nul
// si rien n'est calculé, sinon c'est
// une valeur correspondante au dernier calcul
double IsoHyper3DFavier3::Module_compressibilite_equivalent(Enum_dure temps,const Deformation & def,SaveResul * saveDon)
{ // 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
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
double K_sortie=K; // init
// traitement du cas particulier de l'existence de fct nD
if (save_resul.para_loi != NULL)
{bool recup = false;
switch (temps)
{
case TEMPS_t : // on utilise les grandeurs stockées à t
if ((*save_resul.para_loi_t)(1) != (-ConstMath::trespetit))
{K_sortie= (*save_resul.para_loi_t)(1);recup = true;}
break;
case TEMPS_tdt : // on utilise les grandeurs stockées à tdt
if ((*save_resul.para_loi)(1) != (-ConstMath::trespetit))
{K_sortie= (*save_resul.para_loi)(1);recup = true;}
break;
case TEMPS_0 : // rien n'a été calculé
recup = false;
};
// dans tous les autres cas on utilise soit les fonctions
// si elles existent sinon on concerve les valeurs par défaut
if (!recup)
{//on essaie un calcul
// cas des courbes d'évolution
if (K_nD != NULL)
// là il faut calcul 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 = K_nD->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = K_nD->Li_equi_Quel_evolue();
// on initialise les grandeurs du tableau pour les valeurs numériques
Tableau <double> val_ddl_enum(li_enu_scal.size(),0.);
// init par défaut des types quelconques
List_io <TypeQuelconque >::iterator il,ilfin=li_quelc.end();
for (il = li_quelc.begin();il != ilfin; il++)
(*il).Grandeur_pointee()->InitParDefaut();
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val
= K_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 au module d'Young "
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " IsoHyper3DFavier3::Module_compressibilite_equivalent(..\n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
K_sortie = tab_val(1);
};
};
};
// retour du module
double K= K_sortie/3.;
return K;
};
// =========== METHODES Protégées dérivant de virtuelles : ==============
// METHODES internes spécifiques à l'hyperélasticité isotrope découlant de
// méthodes virtuelles de Hyper3D
// calcul du potentiel tout seul sans la phase car Qeps est nul
// ou très proche de 0
double IsoHyper3DFavier3::PoGrenoble
(const double & Qeps,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
// des variables intermédiaires
double a = Qor/2./mur;
double co1 = Qor*a;
double co2 = Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées
double E = K/6. * (logV)*(logV)
+ co1*log_A + mu_inf * Qeps * Qeps;
// retour
return E;
};
// calcul du potentiel tout seul avec la phase donc dans le cas où Qeps est non nul
double IsoHyper3DFavier3::PoGrenoble
(const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
// dans le cas de l'existence de la phase, on modifie Qr et muinf
double bmuinf,bQr,muinfF,Qr;
if (avec_phase)
{bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
bQr = (1.+gammaQor*inv_spec.cos3phi);
muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
Qr = Qor/pow(bQr,nQor);
}
else
{bmuinf = 1.;
bQr = 1.;
muinfF = mu_inf;
Qr = Qor;
}
// des variables intermédiaires
double a = Qr/2./mur;
double co1 = Qr*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées
double E = K/6. * (logV)*(logV)
+ co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
// retour
return E;
};
// calcul du potentiel tout seul sans la phase car Qeps est nul
// ou très proche de 0, et de sa variation suivant V uniquement
Hyper3D::PoGrenoble_V IsoHyper3DFavier3::PoGrenoble_et_V
(const double & Qeps,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
PoGrenoble_V ret;
// des variables intermédiaires
double a = Qor/2./mur;
double co1 = Qor*a;
double co2 = Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + mu_inf * Qeps * Qeps;
ret.EV = K/3.*logV/inv.V;
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel et de sa variation suivant V uniquement
Hyper3D::PoGrenoble_V IsoHyper3DFavier3::PoGrenoble_et_V
(const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
Hyper3D::PoGrenoble_V ret;
// dans le cas de l'existence de la phase, on modifie Qr et muinf
double bmuinf,bQr,muinfF,Qr;
if (avec_phase)
{bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
bQr = (1.+gammaQor*inv_spec.cos3phi);
muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
Qr = Qor/pow(bQr,nQor);
}
else
{bmuinf = 1.;
bQr = 1.;
muinfF = mu_inf;
Qr = Qor;
}
// des variables intermédiaires
double a = Qr/2./mur;
double co1 = Qr*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
ret.EV = K/3.*logV/inv.V;
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel tout seul sans la phase car Qeps est nul
// ou très proche de 0, et de ses variations première et seconde suivant V uniquement
Hyper3D::PoGrenoble_VV IsoHyper3DFavier3::PoGrenoble_et_VV
(const double & Qeps,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
PoGrenoble_VV ret;
// des variables intermédiaires
double a = Qor/2./mur;
double co1 = Qor*a;
double co2 = Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + mu_inf * Qeps * Qeps;
ret.EV = K/3.*logV/inv.V;
// dérivées secondes
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel et de sa variation première et seconde suivant V uniquement
Hyper3D::PoGrenoble_VV IsoHyper3DFavier3::PoGrenoble_et_VV
(const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
Hyper3D::PoGrenoble_VV ret;
// dans le cas de l'existence de la phase, on modifie Qr et muinf
double bmuinf,bQr,muinfF,Qr;
if (avec_phase)
{bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
bQr = (1.+gammaQor*inv_spec.cos3phi);
muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
Qr = Qor/pow(bQr,nQor);
}
else
{bmuinf = 1.;
bQr = 1.;
muinfF = mu_inf;
Qr = Qor;
}
// des variables intermédiaires
double a = Qr/2./mur;
double co1 = Qr*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
ret.EV = K/3.*logV/inv.V;
// dérivées secondes
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel et de ses dérivées non compris la phase
Hyper3D::PoGrenobleSansPhaseSansVar IsoHyper3DFavier3::PoGrenoble
(const InvariantQeps & inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
PoGrenobleSansPhaseSansVar ret;
// des variables intermédiaires
double a = Qor/2./mur;
double co1 = Qor*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + mu_inf * inv_spec.Qeps * inv_spec.Qeps;
ret.EV = K/3.*logV/inv.V;
// // dans le cas où l'on est très près de l'origine (voir nul)
// // il faut un traitement particulier
// if (co2 > ConstMath::unpeupetit)
// cas normal
ret.EQ = Qor * tanh(co2) + 2.* mu_inf * inv_spec.Qeps;
/* else
// cas ou Qeps est très proche de zéro, on utilise un développement
// limité
ret.EQ = (Qor/a + mu_inf)* inv_spec.Qeps
- 1./3. * Qor/(a*a*a)
* inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps; */
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PoGrenobleAvecPhaseSansVar IsoHyper3DFavier3::PoGrenoblePhase
(const InvariantQepsCosphi& inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
PoGrenobleAvecPhaseSansVar ret;
// dans le cas de l'existence de la phase, on modifie Qr et muinf
double bmuinf,bQr,muinfF,Qr;
if (avec_phase)
{bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
bQr = (1.+gammaQor*inv_spec.cos3phi);
muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
Qr = Qor/pow(bQr,nQor);
}
else
{bmuinf = 1.;
bQr = 1.;
muinfF = mu_inf;
Qr = Qor;
}
// des variables intermédiaires
double a = Qr/2./mur;
double co1 = Qr*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
}
else
{ log_A = log(A);};
double logV = log(inv.V);
// le potentiel et ses dérivées
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + muinfF * inv_spec.Qeps * inv_spec.Qeps;
ret.EV = K/3.*logV/inv.V;
// // dans le cas où l'on est très près de l'origine (voir nul)
// // il faut un traitement particulier
// if (co2 > ConstMath::unpeupetit)
// cas normal
ret.EQ = Qr * tanh(co2) + 2.* muinfF * inv_spec.Qeps;
/* else
// cas ou Qeps est très proche de zéro, on utilise un développement
// limité
ret.EQ = (Qor/a + mu_inf)* inv_spec.Qeps
- 1./3. * Qor/(a*a*a)
* inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps; */
// cas de la variation faisant intervenir la phase
double dEdmuinf= inv_spec.Qeps * inv_spec.Qeps;
double dmuinf_dcos3phi = (- n_mu_inf * gamma_mu_inf * muinfF/bmuinf);
double dQdcos3phi = (- nQor * gammaQor * Qr / bQr);
ret.Ecos3phi = dEdmuinf * dmuinf_dcos3phi + ret.EQ * dQdcos3phi;
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel sans phase et dérivées avec ses variations par rapport aux invariants
Hyper3D::PoGrenobleSansPhaseAvecVar IsoHyper3DFavier3::PoGrenoble_et_var
(const Invariant2Qeps& inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
PoGrenobleSansPhaseAvecVar ret;
// des variables intermédiaires
double a = Qor/2./mur;
double co1 = Qor*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
double tanhco2 = 1.; // init par défaut
if ((Abs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est // if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
tanhco2 = 1.; // et dans ce cas on est à la limite : th(A) =(environ)= exp(co2)/exp(co2) = 1
}
else
{ log_A = log(A);
tanhco2 = tanh(co2);
};
double logV = log(inv.V);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV)
+ co1*log_A + mu_inf * inv_spec.Qeps * inv_spec.Qeps;
ret.EV = K/3.*logV/inv.V;
// dérivées secondes
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
ret.EQV = 0.;
// dans le cas où l'on est très près de l'origine (voir nul)
// il faut un traitement particulier
// if (co2 > ConstMath::unpeupetit)
// cas normal
{ ret.EQ = Qor * tanhco2 + 2.* mu_inf * inv_spec.Qeps;
// dérivée seconde
ret.EQQ = 2.*mur*(1. - tanhco2 * tanhco2) + 2.* mu_inf;
}
/* else
// cas ou Qeps est très proche de zéro, on utilise un développement
// limité
{ ret.EQ = (Qor/a + mu_inf)* inv_spec.Qeps
- 1./3. * Qor/(a*a*a)
* inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;
// dérivée seconde
ret.EQQ = (Qor/a + mu_inf)
- Qor/(a*a*a)
* inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;
}*/
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
// --- Verif_PoGrenoble_et_var(inv_spec.Qeps,inv,ret );
//--- débug: autre vérification à décommenter si on veut vérifier un nan événtuel !! --------------------
// if ((Dabs(ret.E) > ConstMath::tresgrand) || (Dabs(ret.EV) > ConstMath::tresgrand)
// || (Dabs(ret.EVV) > ConstMath::tresgrand) || (Dabs(ret.EQV) > ConstMath::tresgrand)
// || (Dabs(ret.EQQ) > ConstMath::tresgrand) || (Dabs(ret.EQ) > ConstMath::tresgrand) )
// { cout << "\n attention *** on a detecter un comportement bizarre sur le potentiel de la loi Favier3D "
// << " potentiel= " << ret.E << " var/V= " << ret.EV << " var/Q= " << ret.EQ
// << " var/VV= " << ret.EVV << " var/QQ= " << ret.EQQ << " var/QV= " << ret.EQV << endl;
// };
//---- fin débug -----------------------------------------------------------------------------------------
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel avec phase et dérivées avec ses variations par rapport aux invariants
Hyper3D::PoGrenobleAvecPhaseAvecVar IsoHyper3DFavier3::PoGrenoblePhase_et_var
(const Invariant2QepsCosphi& inv_spec,const Invariant & inv)
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
Hyper3D::PoGrenobleAvecPhaseAvecVar ret;
// dans le cas de l'existence de la phase, on modifie Qr et muinf
double bmuinf,bQr,muinfF,Qr;
if (avec_phase)
{bmuinf = (1.+gamma_mu_inf*inv_spec.cos3phi);
bQr = (1.+gammaQor*inv_spec.cos3phi);
muinfF = mu_inf/ pow(bmuinf,n_mu_inf);
Qr = Qor/pow(bQr,nQor);
}
else
{bmuinf = 1.;
bQr = 1.;
muinfF = mu_inf;
Qr = Qor;
}
//debug
//cout << "\n debug: IsoHyper3DFavier3::PoGrenoblePhase_et_var "
// << " Qr="<<Qr << " Qor=" << Qor << " muinfF="<<muinfF <<" mu_inf="<<mu_inf<<endl;
//fin debug
// des variables intermédiaires
double a = Qr/2./mur;
double co1 = Qr*a;
double co2 = inv_spec.Qeps/a ;
double A = cosh(co2);
double log_A=0.;
double tanhco2 = 1.; // init par défaut
if ((Dabs(co2) > limite_co2) || (!isfinite(A))) // en fait A est "toujours" très grand, donc c'est normal, et il vaut mieux regarder s'il est infini
// if ((Dabs(A)>ConstMath::pasmalgrand)||(!isfinite(A)))
// if (Dabs(A)>ConstMath::tresgrand)
{// si co2 est grand (par exemple > 800 uniquement !!) ou très petit, comme cosh(co2) = (exp(co2) + exp(-co2))/2, on a exp(co2) = inf
// ou exp(-co2) = inf
// mais en fait on n'utilise que le log(A) donc on peut faire une approximation
// soit co2 est grand
// on fait l'approximation que expo(-co2) est pratiquement nulle d'où cosh(co2) =(environ)= exp(co2)/2
// d'où log(A) =(environ)= co2-log(2)
log_A = MaX(0.,Dabs(co2)-log(2.));
tanhco2 = 1.; // et dans ce cas on est à la limite : th(A) =(environ)= exp(co2)/exp(co2) = 1
}
else
{ log_A = log(A);
tanhco2 = tanh(co2);
};
//debug
//cout << "\n debug: 3DFavier3::Po..ePhase_et_var "
// << " A="<< A << " co2="<< co2<< " log_A="<<log_A << " log(A)=" << log(A) << " tanhco2="<<tanhco2 <<" tanh(co2)="<<tanh(co2) <<endl;
//fin debug
double logV = log(inv.V);
double logCosh = log_A;
double unsurmur=1./mur;
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV)
+ co1*logCosh + muinfF * inv_spec.Qeps * inv_spec.Qeps;
ret.EV = K/3.*logV/inv.V;
// dérivées secondes
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
ret.EQV = 0.;
// dans le cas où l'on est très près de l'origine (voir nul)
// il faut un traitement particulier
// if (co2 > ConstMath::unpeupetit)
// cas normal
{ ret.EQ = Qr * tanhco2 + 2.* muinfF * inv_spec.Qeps;
// dérivée seconde
ret.EQQ = 2.*mur*(1. - tanhco2 * tanhco2) + 2.* muinfF;
}
/* else
// cas ou Qeps est très proche de zéro, on utilise un développement
// limité
{ ret.EQ = (Qr/a + muinfF)* inv_spec.Qeps
- 1./3. * Qr/(a*a*a)
* inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;
// dérivée seconde
ret.EQQ = (Qr/a + muinfF)
- Qr/(a*a*a)
* inv_spec.Qeps * inv_spec.Qeps * inv_spec.Qeps;
}*/
// cas de la variation première par rapport à la phase
double dEdQr= Qr*unsurmur * logCosh - inv_spec.Qeps * tanhco2;
double QepsSurQr = inv_spec.Qeps/Qr;
double tanhco2_2=tanhco2*tanhco2;
double dEdmuinf= inv_spec.Qeps * inv_spec.Qeps;
double dmuinf_dcos3phi = (- n_mu_inf * gamma_mu_inf * muinfF/bmuinf);
double dQdcos3phi = (- nQor * gammaQor * Qr / bQr);
ret.Ecos3phi = dEdmuinf * dmuinf_dcos3phi + dEdQr * dQdcos3phi;
// cas des variations seconde faisant intervenir la phase la phase
double d2EdQr2=unsurmur*logCosh-2.* QepsSurQr * tanhco2 + co2*QepsSurQr*(1.-tanhco2_2);
double d2muinf_dcos3phi_2 = n_mu_inf*(1.+n_mu_inf)* muinfF*gamma_mu_inf*gamma_mu_inf/bmuinf/bmuinf;
double d2Qr_dcos3phi_2 = nQor*(1.+nQor)*Qr*gammaQor*gammaQor/bQr/bQr;
double d2E_dQ_dmuinf = 2. * inv_spec.Qeps ;
double d2E_dQ_dQr = tanhco2 - co2 * (1.-tanhco2_2);
ret.EVcos3phi= 0.;
ret.EQcos3phi= d2E_dQ_dmuinf * dmuinf_dcos3phi + d2E_dQ_dQr * dQdcos3phi;
ret.Ecos3phi2= d2EdQr2 * ret.Ecos3phi * ret.Ecos3phi + dEdmuinf * d2muinf_dcos3phi_2
+ dEdQr * d2Qr_dcos3phi_2;
/* ret.EVcos3phi= 0.;
ret.EQcos3phi= 0.;
ret.Ecos3phi2= 0.;
ret.Ecos3phi = 0.;*/
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
// --- Verif_PoGrenoble_et_var(inv_spec.Qeps,inv,inv_spec.cos3phi,ret );
// --- Invariant2Qeps toto(inv_spec.Qeps,inv_spec.dQepsdbIIb,inv_spec.dQepsdbIIb2);
// --- Hyper3D::PoGrenobleSansPhaseAvecVar truc = PoGrenoble_et_var(toto,inv);
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
///=============== fonctions pour la vérification et la mise au point ===============
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DFavier3::Verif_PoGrenoble_et_var
(const double & Qeps,const Invariant & inv
,const PoGrenobleSansPhaseAvecVar& potret )
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
// dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_PoGrenoble_et_var == 0)
{ cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var \n";
}
indic_Verif_PoGrenoble_et_var++;
// potret contiend le potentiel et ses variations première et seconde
// on va calculer ces mêmes variations par différences finies et comparer les deux résultats
// calcul des invariants sous la nouvelle forme
double toto = potret.E; // pour que le débugger affiche potret
Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl
double E_n = PoGrenoble(Qeps,inv_n0); // la valeur du potentiel de référence
// ici on est sans phase donc deux invariants indépendant V et Qeps,
// pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
// pour les dérivées premières et secondes
double delta = ConstMath::unpeupetit;
double Qeps_et_dQeps = Qeps+delta;
double E_et_dQeps = PoGrenoble(Qeps_et_dQeps,inv_n0);
double Qeps_et_dQeps2 = Qeps + 2.*delta;
double E_et_dQeps2 = PoGrenoble(Qeps_et_dQeps2,inv_n0);
Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
double E_et_dV = PoGrenoble(Qeps,inv_et_dV);
Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
double E_et_dV2 = PoGrenoble(Qeps,inv_et_dV2);
Invariant inv_et_dVdQeps = inv_n0;
inv_et_dVdQeps.V += delta;
double Qeps_et_dVdQeps = Qeps+delta;
double E_et_dVdQeps = PoGrenoble(Qeps_et_dVdQeps,inv_et_dVdQeps);
// calcul des dérivées premières
double E_V = (E_et_dV - E_n)/delta;
double E_Qeps = (E_et_dQeps - E_n)/delta;
// calcul des dérivées secondes
double E_V2 = (E_et_dV2 - 2.*E_et_dV + E_n ) /(delta * delta);
double E_Qeps2 = (E_et_dQeps2 - 2.*E_et_dQeps + E_n)/(delta * delta);
double E_V_a_dQeps = (E_et_dVdQeps - E_et_dQeps )/delta;
double E_VQeps = ( E_V_a_dQeps - E_V)/delta;
// comparaison avec les valeurs de dérivées analytiques
bool erreur = false;
if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
{if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
{if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = true;}
else
erreur = true;}
if (diffpourcent(potret.EQ,E_Qeps,MaX(Dabs(E_Qeps),Dabs(potret.EQ)),0.05))
{if (MiN(Dabs(E_Qeps),Dabs(potret.EQ)) == 0.)
{if ( MaX(Dabs(E_Qeps),Dabs(potret.EQ)) > 50.*delta) erreur = true;}
else
erreur = true;}
if (diffpourcent(potret.EQQ,E_Qeps2,MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)),0.05))
{if (MiN(Dabs(E_Qeps2),Dabs(potret.EQQ)) == 0.)
{if ( MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)) > 50.*delta) erreur = true;}
else
erreur = true;}
if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05))
{if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
{if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = true;}
else
erreur = true;}
if (diffpourcent(potret.EQV,E_VQeps,MaX(Dabs(E_VQeps),Dabs(potret.EQV)),0.05))
{if (MiN(Dabs(E_VQeps),Dabs(potret.EQV)) == 0.)
{if (Dabs(potret.EQV) == 0.)
// si la valeur de la dérivée numérique est nulle cela signifie peut-être
// que l'on est à un extréma, la seule chose que le l'on peut faire est de
// vérifier que l'on tend numériquement vers cette extréma ici un minima
{
Invariant inv_et_ddVddQeps = inv_n0;
inv_et_ddVddQeps.V += delta/10.;
double Qeps_et_ddVddQeps = Qeps+delta/10.;
double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
double Qeps_et_ddQeps = Qeps+delta/10.;
double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
double E_V_a_ddQeps = (E_et_ddVddQeps - E_et_ddQeps )/(delta*10.);
if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps)) erreur = true;
}
else
if ( MaX(Dabs(E_VQeps),Dabs(potret.EQV)) > 150.*delta) erreur = true;
}
else
erreur = true;
}
if (erreur)
{ cout << "\n erreur dans le calcul analytique des derivees du potentiel";
cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var(.."
<< " , numero d'increment = " << indic_Verif_PoGrenoble_et_var;
Sortie(1);
}
};
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DFavier3::Verif_PoGrenoble_et_var
(const double & Qeps,const Invariant & inv,const double& cos3phi
,const PoGrenobleAvecPhaseAvecVar& potret )
{ // calcul éventuelle des paramètres de la loi
if (K_nD != NULL)
{K = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (Qor_nD != NULL)
{Qor = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(Qor_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mur_nD != NULL)
{mur = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mur_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
if (mu_inf_nD != NULL)
{mu_inf = (Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(mu_inf_nD,1,ex_impli_hyper,ex_expli_tdt_hyper,ex_umat_hyper))(1);
};
// sauvegarde éventuelle des paramètres matériau
SaveResulIsoHyper3DFavier3 & save_resul = *((SaveResulIsoHyper3DFavier3*) saveResul);
if (save_resul.para_loi != NULL)
{(*save_resul.para_loi)(1) = K;
(*save_resul.para_loi)(2) = Qor;
(*save_resul.para_loi)(3) = mur;
(*save_resul.para_loi)(4) = mu_inf;
};
// dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_PoGrenoble_et_var == 0)
{ cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var \n";
}
indic_Verif_PoGrenoble_et_var++;
// potret contiend le potentiel et ses variations première et seconde
// on va calculer ces mêmes variations par différences finies et comparer les deux résultats
// calcul des invariants sous la nouvelle forme
double toto = potret.E; // pour que le débugger affiche potret
Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl
Invariant0QepsCosphi invcos3phi(Qeps,cos3phi);
double E_n = PoGrenoble(invcos3phi,inv_n0); // la valeur du potentiel de référence
// ici on est sans phase donc deux invariants indépendant V et Qeps,
// pour calculer les variations on définit des points distants d'un incrément puis de deux incréments
// pour les dérivées premières et secondes
double delta = 10.*ConstMath::unpeupetit;
double Qeps_et_dQeps = Qeps+delta;
Invariant0QepsCosphi invcos3phi1(Qeps_et_dQeps,cos3phi);
double E_et_dQeps = PoGrenoble(invcos3phi1,inv_n0);
double Qeps_moins_dQeps = Qeps - delta;
Invariant0QepsCosphi invcos3phi2(Qeps_moins_dQeps,cos3phi);
double E_moins_dQeps = PoGrenoble(invcos3phi2,inv_n0);
Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
double E_et_dV = PoGrenoble(invcos3phi,inv_et_dV);
Invariant inv_moins_dV = inv_n0;inv_moins_dV.V -= delta;
double E_moins_dV = PoGrenoble(invcos3phi,inv_moins_dV);
Invariant inv_et_dVdQeps = inv_n0;
inv_et_dVdQeps.V += delta;
double Qeps_et_dVdQeps = Qeps+delta;
Invariant0QepsCosphi invcos3phi3(Qeps_et_dVdQeps,cos3phi);
double E_et_dVdQeps = PoGrenoble(invcos3phi3,inv_et_dVdQeps);
// cas des variations avec cos3phi
double cos3phi_et_dcos3phi = cos3phi+delta;
Invariant0QepsCosphi invcos3phi_et_dcos3phi(Qeps,cos3phi_et_dcos3phi);
double E_et_dcos3phi = PoGrenoble(invcos3phi_et_dcos3phi,inv_n0);
double cos3phi_moins_dcos3phi = cos3phi-delta;
Invariant0QepsCosphi invcos3phi_moins_dcos3phi(Qeps,cos3phi_moins_dcos3phi);
double E_moins_dcos3phi = PoGrenoble(invcos3phi_moins_dcos3phi,inv_n0);
Invariant0QepsCosphi invcos3phi_et_dcos3phi_a_dQeps(Qeps_et_dQeps,cos3phi_et_dcos3phi);
double E_et_dcos3phi_a_dQeps = PoGrenoble(invcos3phi_et_dcos3phi_a_dQeps,inv_n0);
double E_et_dcos3phi_a_dV = PoGrenoble(invcos3phi_et_dcos3phi,inv_et_dV);
// calcul des dérivées premières
double E_V = (E_et_dV - E_n)/delta;
double E_Qeps = (E_et_dQeps - E_n)/delta;
double E_cos3phi = (E_et_dcos3phi -E_n)/delta;
// calcul des dérivées secondes
double E_V2 = (E_et_dV - 2.*E_n + E_moins_dV ) /(delta * delta);
double E_Qeps2 = (E_et_dQeps - 2.*E_n + E_moins_dQeps)/(delta * delta);
double E_V_a_dQeps = (E_et_dVdQeps - E_et_dQeps )/delta;
double E_VQeps = ( E_V_a_dQeps - E_V)/delta;
double E_cos3phi_2 = (E_et_dcos3phi - 2.*E_n + E_moins_dcos3phi ) /(delta * delta);
double E_dQeps_a_dcos3phi = (E_et_dcos3phi_a_dQeps - E_et_dcos3phi )/delta;
double E_dcos3phi_dQeps = (E_dQeps_a_dcos3phi - E_Qeps )/delta;
double E_dV_a_dcos3phi = (E_et_dcos3phi_a_dV - E_et_dcos3phi )/delta;
double E_dcos3phi_dV = (E_dV_a_dcos3phi - E_V )/delta;
// dans les dérivées secondes par rapport à cos3phi on a des erreurs, sans doute à cause des petites valeurs ??
// comparaison avec les valeurs de dérivées analytiques
int erreur = 0;
if (diffpourcent(potret.EV,E_V,MaX(Dabs(E_V),Dabs(potret.EV)),0.05))
{if (MiN(Dabs(E_V),Dabs(potret.EV)) == 0.)
{if ( MaX(Dabs(E_V),Dabs(potret.EV)) > 50.*delta) erreur = 1;}
else
erreur = 2; }
if (diffpourcent(potret.EQ,E_Qeps,MaX(Dabs(E_Qeps),Dabs(potret.EQ)),0.05))
{if (MiN(Dabs(E_Qeps),Dabs(potret.EQ)) == 0.)
{if ( MaX(Dabs(E_Qeps),Dabs(potret.EQ)) > 50.*delta) erreur = 3;}
else
erreur = 4;}
if (diffpourcent(potret.Ecos3phi,E_cos3phi,MaX(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)),0.05))
{if (MiN(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)) == 0.)
{if ( MaX(Dabs(E_cos3phi),Dabs(potret.Ecos3phi)) > 50.*delta) erreur = 31;}
else
erreur = 41;}
if (diffpourcent(potret.EQQ,E_Qeps2,MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)),0.05))
{if (MiN(Dabs(E_Qeps2),Dabs(potret.EQQ)) == 0.)
{if ( MaX(Dabs(E_Qeps2),Dabs(potret.EQQ)) > 50.*delta) erreur = 5;}
else
erreur = 6;}
if (diffpourcent(potret.EVV,E_V2,MaX(Dabs(E_V2),Dabs(potret.EVV)),0.05))
{if (MiN(Dabs(E_V2),Dabs(potret.EVV)) == 0.)
{if ( MaX(Dabs(E_V2),Dabs(potret.EVV)) > 50.*delta) erreur = 7;}
else
erreur = 8;}
if (diffpourcent(potret.Ecos3phi2,E_cos3phi_2,MaX(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)),0.05))
{if (MiN(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)) == 0.)
{if ( MaX(Dabs(E_cos3phi_2),Dabs(potret.Ecos3phi2)) > 50.*delta) erreur = 71;}
else
erreur = 81;}
if (diffpourcent(potret.EQV,E_VQeps,MaX(Dabs(E_VQeps),Dabs(potret.EQV)),0.05))
{if (MiN(Dabs(E_VQeps),Dabs(potret.EQV)) == 0.)
{if (Dabs(potret.EQV) == 0.)
// si la valeur de la dérivée numérique est nulle cela signifie peut-être
// que l'on est à un extréma, la seule chose que le l'on peut faire est de
// vérifier que l'on tend numériquement vers cette extréma ici un minima
{
Invariant inv_et_ddVddQeps = inv_n0;
inv_et_ddVddQeps.V += delta/10.;
double Qeps_et_ddVddQeps = Qeps+delta/10.;
double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
double Qeps_et_ddQeps = Qeps+delta/10.;
double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
double E_V_a_ddQeps = (E_et_ddVddQeps - E_et_ddQeps )/(delta*10.);
if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps)) erreur = 9;
}
else
if ( MaX(Dabs(E_VQeps),Dabs(potret.EQV)) > 150.*delta) erreur = 10;
}
else
erreur = 11;
}
if (diffpourcent(potret.EVcos3phi,E_dcos3phi_dV,MaX(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)),0.05))
{if (MiN(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)) == 0.)
{if (Dabs(potret.EVcos3phi) == 0.)
// si la valeur de la dérivée numérique est nulle cela signifie peut-être
// que l'on est à un extréma, la seule chose que le l'on peut faire est de
// vérifier que l'on tend numériquement vers cette extréma ici un minima
{
/*Invariant inv_et_ddVddQeps = inv_n0;
inv_et_ddVddQeps.V += delta/10.;
double Qeps_et_ddVddQeps = Qeps+delta/10.;
double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
double Qeps_et_ddQeps = Qeps+delta/10.;
double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
double E_V_a_ddQeps = (E_et_ddVddQeps - E_et_ddQeps )/(delta*10.);
if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps))*/ erreur = 91;
}
else
if ( MaX(Dabs(E_dcos3phi_dV),Dabs(potret.EVcos3phi)) > 150.*delta) erreur = 101;
}
else
erreur = 111;
}
if (diffpourcent(potret.EQcos3phi,E_dcos3phi_dQeps,MaX(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)),0.05))
{if (MiN(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)) == 0.)
{if (Dabs(potret.EQcos3phi) == 0.)
// si la valeur de la dérivée numérique est nulle cela signifie peut-être
// que l'on est à un extréma, la seule chose que le l'on peut faire est de
// vérifier que l'on tend numériquement vers cette extréma ici un minima
{
/*Invariant inv_et_ddVddQeps = inv_n0;
inv_et_ddVddQeps.V += delta/10.;
double Qeps_et_ddVddQeps = Qeps+delta/10.;
double E_et_ddVddQeps = PoGrenoble(Qeps_et_ddVddQeps,inv_et_ddVddQeps);
double Qeps_et_ddQeps = Qeps+delta/10.;
double E_et_ddQeps = PoGrenoble(Qeps_et_ddQeps,inv_n0);
double E_V_a_ddQeps = (E_et_ddVddQeps - E_et_ddQeps )/(delta*10.);
if (10.* Dabs(E_V_a_ddQeps) > Dabs(E_V_a_dQeps))*/ erreur = 92;
}
else
if ( MaX(Dabs(E_dcos3phi_dQeps),Dabs(potret.EQcos3phi)) > 150.*delta) erreur = 102;
}
else
erreur = 112;
}
if (erreur > 0)
{ cout << "\n erreur dans le calcul analytique des derivees du potentiel";
cout << "\n IsoHyper3DFavier3::Verif_PoGrenoble_et_var(.."
<< " , numero d'increment = " << indic_Verif_PoGrenoble_et_var
<< " , numero d'erreur : " << erreur ;
// Sortie(1);
}
};