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

2253 lines
111 KiB
C++
Raw Permalink Normal View History

ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
#include "IsoHyper3DOrgeas1.h"
#include "ComLoi_comp_abstraite.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "CharUtil.h"
//================== initialisation des variables static ======================
// indicateur utilisé par Verif_Potentiel_et_var
int IsoHyper3DOrgeas1::indic_Verif_PoGrenobleorgeas1_et_var = 0;
// indicateur utilisé par Verif_Potentiel_et_var_VV
int IsoHyper3DOrgeas1::indic_Verif_PoGrenobleorgeas1_et_var_VV = 0;
// indicateur utilisé par CompaFormel
int IsoHyper3DOrgeas1::indic_Verif_CompaFormel = 0;
//================== fin d'initialisation des variables static ================
IsoHyper3DOrgeas1::IsoHyper3DOrgeas1 () : // Constructeur par defaut
Hyper3D(ISOHYPER3DORGEAS1,CAT_MECANIQUE,false),K(ConstMath::tresgrand),Q0s(ConstMath::tresgrand)
,mu01(ConstMath::tresgrand),mu02(ConstMath::tresgrand),mu03(ConstMath::tresgrand)
,alpha1(ConstMath::tresgrand),alpha3(ConstMath::tresgrand),Q0e(ConstMath::tresgrand)
,alpha1_2(0.),alpha3_2(0.) //,mu1_2(0.)
,nQs(1.),gammaQs(0.),nQe(1.),gammaQe(0.),nMu1(1.),gammaMu1(0.),nMu2(1.),gammaMu2(0.),nMu3(1.),gammaMu3(0.)
,cas_Q0s_thermo_dependant(0),cas_Q0e_thermo_dependant(0)
,T0rQe(ConstMath::tresgrand),Qe_T0rQe(ConstMath::tresgrand),h1(ConstMath::tresgrand)
,T0r(ConstMath::tresgrand),Gr(ConstMath::tresgrand)
,Qrmax(ConstMath::tresgrand),ar(ConstMath::tresgrand),Tc(0.),Ts(0.),Qs_T0rQe(0.)
,Q0s_temperature(NULL),Q0e_temperature(NULL)
{};
// Constructeur de copie
IsoHyper3DOrgeas1::IsoHyper3DOrgeas1 (const IsoHyper3DOrgeas1& loi) :
Hyper3D (loi),K(loi.K),Q0s(loi.Q0s),mu01(loi.mu01),mu02(loi.mu02)
,mu03(loi.mu03),alpha1(loi.alpha1),alpha3(loi.alpha3),Q0e(loi.Q0e)
,alpha1_2(loi.alpha1_2)
,alpha3_2(loi.alpha3_2) //,mu1_2(loi.mu1_2)
,nQs(loi.nQs),gammaQs(loi.gammaQs),nQe(loi.nQe),gammaQe(loi.gammaQe)
,nMu1(loi.nMu2),gammaMu1(loi.gammaMu2),nMu2(loi.nMu2),gammaMu2(loi.gammaMu2),nMu3(loi.nMu3),gammaMu3(loi.gammaMu3)
,cas_Q0s_thermo_dependant(loi.cas_Q0s_thermo_dependant),cas_Q0e_thermo_dependant(loi.cas_Q0e_thermo_dependant)
,T0rQe(loi.T0rQe),Qe_T0rQe(loi.Qe_T0rQe),h1(loi.h1)
,T0r(loi.T0r),Gr(loi.Gr),Qrmax(loi.Qrmax),ar(loi.ar)
,Tc(loi.Tc),Ts(loi.Ts),Qs_T0rQe(loi.Qs_T0rQe)
,Q0s_temperature(loi.Q0s_temperature),Q0e_temperature(loi.Q0e_temperature)
{// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale pour Qs(T)
if (Q0s_temperature != NULL)
if (Q0s_temperature->NomCourbe() == "_")
{// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
Q0s_temperature = Courbe1D::New_Courbe1D(*loi.Q0s_temperature);
};
// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale pour Qe(T)
if (Q0e_temperature != NULL)
if (Q0e_temperature->NomCourbe() == "_")
{// comme il s'agit d'une courbe locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
Q0e_temperature = Courbe1D::New_Courbe1D(*loi.Q0e_temperature);
};
};
IsoHyper3DOrgeas1::~IsoHyper3DOrgeas1 ()
// Destructeur
{ if (Q0s_temperature != NULL)
if (Q0s_temperature->NomCourbe() == "_") delete Q0s_temperature;
if (Q0e_temperature != NULL)
if (Q0e_temperature->NomCourbe() == "_") delete Q0e_temperature;
};
// Lecture des donnees des paramètres matériau sur fichier
void IsoHyper3DOrgeas1::LectureDonneesParticulieres
(UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D,LesFonctions_nD& lesFonctionsnD)
{ // lecture des quatres coefficients de la loi
*(entreePrinc->entree) >> K ;
// on regarde s'il y a une thermo-dépendance pour Q0s
if (strstr(entreePrinc->tablcar,"Q0s_thermo_dependant_")!=NULL)
{ string nom;
*(entreePrinc->entree) >> nom;
if (nom != "Q0s_thermo_dependant_") // vérification
{ cout << "\n erreur en lecture du parametre Q0s, on attendait le mot cle "
<< " Q0s_thermo_dependant_ et on a lue: " << nom;
entreePrinc->MessageBuffer("**1--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
Q0s=0.; // init néanmoins de la valeur
cas_Q0s_thermo_dependant =1; // on signale qu'il y a thermo dépendance (= 1 n'est pas définitif)
}
else // sinon cas sans thermo-dépendance
{ *(entreePrinc->entree) >> Q0s; };
// puis lecture des autres paramètres
*(entreePrinc->entree) >> mu01 >> mu02 >> mu03
>> alpha1 >> alpha3;
if ((alpha1 <= 0.) || (alpha3 <= 0.) )
{ cout << "\n erreur en lecture des parametres alpha1= " << alpha1 << " et alpha3 = " << alpha3
<< " , ces deux parametres doivent etre strictement positifs pour eviter des / par 0 ";
entreePrinc->MessageBuffer("**2--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
if (mu01 < ConstMath::petit)
{ cout << "\n erreur en lecture du parametre mu1 " << mu01 << " ce potentiel n'est pas prevu pour "
<< " fonctionner avec un mu1 quasi-nul ";
entreePrinc->MessageBuffer("**3--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
if (mu03 < ConstMath::petit)
{ cout << "\n erreur en lecture du parametre mu3 " << mu03 << " ce potentiel n'est pas prevu pour "
<< " fonctionner avec un mu3 quasi-nul ";
entreePrinc->MessageBuffer("**4--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde s'il y a une thermo-dépendance pour Q0e
if (strstr(entreePrinc->tablcar,"Q0e_thermo_dependant_")!=NULL)
{ string nom;
*(entreePrinc->entree) >> nom;
if (nom != "Q0e_thermo_dependant_") // vérification
{ cout << "\n erreur en lecture du parametre Q0e, on attendait le mot cle "
<< " Q0e_thermo_dependant_ et on a lue: " << nom;
entreePrinc->MessageBuffer("**1_1--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
Q0e=0.; // init néanmoins de la valeur
cas_Q0e_thermo_dependant =1; // on signale qu'il y a thermo dépendance (= 1 n'est pas définitif)
}
else // sinon cas sans thermo-dépendance
{ *(entreePrinc->entree) >> Q0e; };
// on regarde s'il y a dépendance à la phase
if (strstr(entreePrinc->tablcar,"avec_phase")!=NULL)
{ string nom1,nom2;
// lecture des paramètres de phase
entreePrinc->NouvelleDonnee(); // passage d'une ligne
// cas de Qsigrev
if (strstr(entreePrinc->tablcar,"nQs=")!=NULL)
{ *(entreePrinc->entree) >> nom1 >> nQs >> nom2 >> gammaQs;
if ((nom1 != "nQs=") || (nom2 != "gammaQs="))
{ cout << "\n erreur en lecture des parametres de phase pour Q0sig_rev, on attendait les mots cles "
<< " nQs= et gammaQs= et on a lue: " << nom1 << " " << nom2;
entreePrinc->MessageBuffer("**5--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// cas de Q0epsrev
if (strstr(entreePrinc->tablcar,"nQe=")!=NULL)
{ *(entreePrinc->entree) >> nom1 >> nQe >> nom2 >> gammaQe;
if ((nom1 != "nQe=") || (nom2 != "gammaQe="))
{ cout << "\n erreur en lecture des parametres de phase pour Qeps_rev, on attendait les mots cles "
<< " nQe= et gammaQe= et on a lue: " << nom1 << " " << nom2;
entreePrinc->MessageBuffer("**6--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// cas de mu1rev
if (strstr(entreePrinc->tablcar,"nMu1=")!=NULL)
{ *(entreePrinc->entree) >> nom1 >> nMu1 >> nom2 >> gammaMu1;
if ((nom1 != "nMu1=") || (nom2 != "gammaMu1="))
{ cout << "\n erreur en lecture des parametres de phase pour mu1_rev, on attendait les mots cles "
<< " nMu1= et gammaMu1= et on a lue: " << nom1 << " " << nom2;
entreePrinc->MessageBuffer("**71--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// cas de mu2rev
if (strstr(entreePrinc->tablcar,"nMu2=")!=NULL)
{ *(entreePrinc->entree) >> nom1 >> nMu2 >> nom2 >> gammaMu2;
if ((nom1 != "nMu2=") || (nom2 != "gammaMu2="))
{ cout << "\n erreur en lecture des parametres de phase pour mu2_rev, on attendait les mots cles "
<< " nMu2= et gammaMu2= et on a lue: " << nom1 << " " << nom2;
entreePrinc->MessageBuffer("**7--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// cas de mu3rev
if (strstr(entreePrinc->tablcar,"nMu3=")!=NULL)
{ *(entreePrinc->entree) >> nom1 >> nMu3 >> nom2 >> gammaMu3;
if ((nom1 != "nMu3=") || (nom2 != "gammaMu3="))
{ cout << "\n erreur en lecture des parametres de phase pour mu3_rev, on attendait les mots cles "
<< " nMu3= et gammaMu3= et on a lue: " << nom1 << " " << nom2;
entreePrinc->MessageBuffer("**8--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
avec_phase=true;
};
// ---- maintenant on traite le cas des thermo-dépendances
// .... cas de Q0s .......
if (cas_Q0s_thermo_dependant != 0)
{ thermo_dependant=true; // on indique que la loi est thermo-dépendant
string nom;
entreePrinc->NouvelleDonnee(); // passage d'une ligne
*(entreePrinc->entree) >> nom ;
if (nom != "cas_Q0s_thermo_dependant_") // vérification
{ cout << "\n erreur en lecture du parametre cas_Q0s_thermo_dependant, on attendait le mot cle "
<< " cas_Q0s_thermo_dependant_ et on a lue: " << nom;
entreePrinc->MessageBuffer("**9--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> cas_Q0s_thermo_dependant;
if (cas_Q0s_thermo_dependant == 1)
{// cas de T0r
*(entreePrinc->entree) >> nom ;
if (nom != "T0r=") // vérification
{ cout << "\n erreur en lecture du parametre T0r, on attendait le mot cle "
<< " T0r= et on a lue: " << nom;
entreePrinc->MessageBuffer("**10--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> T0r;
// cas de Gr
*(entreePrinc->entree) >> nom ;
if (nom != "Gr=") // vérification
{ cout << "\n erreur en lecture du parametre Gr, on attendait le mot cle "
<< " Gr= et on a lue: " << nom;
entreePrinc->MessageBuffer("**11--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> Gr;
// cas de Qrmax
*(entreePrinc->entree) >> nom ;
if (nom != "Qrmax=") // vérification
{ cout << "\n erreur en lecture du parametre Qrmax, on attendait le mot cle "
<< " Qrmax= et on a lue: " << nom;
entreePrinc->MessageBuffer("**12--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> Qrmax;
// cas de ar
*(entreePrinc->entree) >> nom ;
if (nom != "ar=") // vérification
{ cout << "\n erreur en lecture du parametre ar, on attendait le mot cle "
<< " ar= et on a lue: " << nom;
entreePrinc->MessageBuffer("**13--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> ar;
// on calcul maintenant les variables intermédiaires fixes
Tc=T0r+(Qrmax/Gr/2.);
Ts=(ar*ar- Sqr(Qrmax/Gr) )/2./(Qrmax/Gr);
}
else if (cas_Q0s_thermo_dependant == 2)
{// lecture de la loi d'évolution de Q0s en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ Q0s_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
Q0s_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
Q0s_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
}
else
{ cout << "\n erreur en lecture du type de dependance de Q0s de la temperature "
<< " actuellement seules les cas 1 et 2 sont possible et on a lue: " << cas_Q0s_thermo_dependant;
entreePrinc->MessageBuffer("**14--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// .... cas de Q0e .......
if (cas_Q0e_thermo_dependant != 0)
{ thermo_dependant=true; // on indique que la loi est thermo-dépendant
string nom;
entreePrinc->NouvelleDonnee(); // passage d'une ligne
*(entreePrinc->entree) >> nom ;
if (nom != "cas_Q0e_thermo_dependant_") // vérification
{ cout << "\n erreur en lecture du parametre cas_Q0e_thermo_dependant, on attendait le mot cle "
<< " cas_Q0e_thermo_dependant_ et on a lue: " << nom;
entreePrinc->MessageBuffer("**15--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> cas_Q0e_thermo_dependant;
if (cas_Q0e_thermo_dependant == 1)
{// cas de T0rQe
*(entreePrinc->entree) >> nom ;
if (nom != "T0rQe=") // vérification
{ cout << "\n erreur en lecture du parametre T0rQe, on attendait le mot cle "
<< " T0rQe= et on a lue: " << nom;
entreePrinc->MessageBuffer("**16--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> T0rQe;
// cas de Qe_T0rQe
*(entreePrinc->entree) >> nom ;
if (nom != "Qe_T0rQe=") // vérification
{ cout << "\n erreur en lecture du parametre Qe_T0rQe, on attendait le mot cle "
<< " Qe_T0rQe= et on a lue: " << nom;
entreePrinc->MessageBuffer("**17--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> Qe_T0rQe;
// cas de h1
*(entreePrinc->entree) >> nom ;
if (nom != "h1=") // vérification
{ cout << "\n erreur en lecture du parametre h1, on attendait le mot cle "
<< " h1= et on a lue: " << nom;
entreePrinc->MessageBuffer("**17--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> h1;
}
else if (cas_Q0e_thermo_dependant == 2)
{// lecture de la loi d'évolution de Q0e en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ Q0e_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
Q0e_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
Q0e_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
}
else
{ cout << "\n erreur en lecture du type de dependance de Q0e de la temperature "
<< " actuellement seules les cas 1 et 2 sont possible et on a lue: " << cas_Q0e_thermo_dependant;
entreePrinc->MessageBuffer("**18--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// maintenant on calcule la valeur de Qs_T0rQr
{Qs_T0rQe = Q0s; // au cas où il n'y a pas de dépendance de Qs à la température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if (T0rQe >= Tc)
{Qs_T0rQe=0.5*Gr*( (T0rQe-Tc+Ts)-sqrt(Sqr(T0rQe-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs_T0rQe=0.5*Gr*( (T0rQe-Tc-Ts)+sqrt(Sqr(T0rQe-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs_T0rQe = Q0s_temperature->Valeur(T0rQe);
break;
}
};
};
};
// cas avec régularisation (variable stockée dans Hyper3D.h)
string nom1;
if (strstr(entreePrinc->tablcar,"avec_regularisation_")!=NULL)
{ *(entreePrinc->entree) >> nom1 ;
if (nom1 != "avec_regularisation_")
{ cout << "\n erreur en lecture du drapeau de regularisation, on attendait le mot cles "
<< " avec_regularisation_ on a lue: " << nom1 ;
entreePrinc->MessageBuffer("**9--IsoHyper3DOrgeas1::LectureDonneesParticulieres(.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> fact_regularisation;
avec_regularisation=true;
};
// lecture de l'indication éventuelle du post traitement
string nom_class_methode = "IsoHyper3DOrgeas1";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
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD);
// mise à jour et même calcul des variables intermédiaires
// les grandeurs de base ici sans phases pour une première vérification
alpha1_2 = Sqr(alpha1);
alpha3_2 = Sqr(alpha3);
// mu1_2 = Sqr(mu1);
};
// affichage de la loi
void IsoHyper3DOrgeas1::Affiche() const
{ cout << " \n loi de comportement 3D hyperelastique isotrope Orgeas 1 : " << Nom_comp(id_comp)
<< " parametres : \n";
cout << " K= " << K ;
if (cas_Q0s_thermo_dependant != 0)
{ cout << " Q0s thermo dependant, cas= " << cas_Q0s_thermo_dependant;
switch (cas_Q0s_thermo_dependant)
{ case 1: cout << " T0r= " << T0r << " Gr= " << Gr << " Qrmax= " << Qrmax << " ar= " << ar << " ";
break;
case 2: { cout << " courbe Q0s=f(T): " << Q0s_temperature->NomCourbe() <<" ";
if ( Q0s_temperature->NomCourbe() == "_") Q0s_temperature->Affiche();
}
break;
default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
<< " cas non pris en compte !!!! dans Affiche() "
<< "\n IsoHyper3DOrgeas1::Affiche()";
Sortie(1);
};
}
else
{cout << " ; Q0s = " << Q0s ;};
// puis les autres para
cout << " ; mu01 = " << mu01 << " ; mu02 = " << mu02 << " ; mu03 = " << mu03
<< " ; alpha1 = " << alpha1 << " ; alpha3 = " << alpha3 ;
// Qe thermo dépendant ou pas
if (cas_Q0e_thermo_dependant != 0)
{ cout << " Q0e thermo dependant, cas= " << cas_Q0e_thermo_dependant;
switch (cas_Q0e_thermo_dependant)
{ case 1: cout << " T0rQe= " << T0rQe << " Qe_T0rQe= " << Qe_T0rQe << " h1= " << h1 << " ";
break;
case 2: { cout << " courbe Q0e=f(T): " << Q0e_temperature->NomCourbe() <<" ";
if ( Q0e_temperature->NomCourbe() == "_") Q0e_temperature->Affiche();
}
break;
default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
<< " cas non pris en compte !!!! dans Affiche() "
<< "\n IsoHyper3DOrgeas1::Affiche()";
Sortie(1);
};
}
else
{cout << " ; Q0e = " << Q0e ;};
// les dépendances à la phase
if (avec_phase)
{ cout << "\n nQs= " << nQs << " gammaQs= " << gammaQs
<< "\n nQe= " << nQe << " gammaQe= " << gammaQe
<< "\n nMu1= " << nMu1 << " gammaMu1= " << gammaMu1
<< "\n nMu2= " << nMu2 << " gammaMu2= " << gammaMu2
<< "\n nMu3= " << nMu3 << " gammaMu3= " << gammaMu3;
};
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 IsoHyper3DOrgeas1::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 Orgeas 1 ........"
<< "\n#-----------------------------------------------------------------------------------------"
<< "\n# K | Q0s | mu01 | mu02 | mu03 | alpha1 | alpha2 | Q0e |"
<< "\n#-----------------------------------------------------------------------------------------"
<< "\n 160000 150 17000 220 17000 0.002 0.004 0.01 ";
if ((rep != "o") && (rep != "O" ) && (rep != "0") )
{ sort << "\n# il est egalement possible de definir une dependance a la phase des parametres"
<< "\n# dans ce cas, a la suite du parametre Q0e on met le mot cle: avec_phase "
<< "\n# puis sur la ligne qui suit on met les 10 coefficients de controle de la dependance "
<< "\n# suit un exemple de declaration (il est possible d'ommettre un groupe de 2 "
<< "\n# mais il faut respecter l'ordre ) "
<< "\n#-----------------------------------------------------------------------"
<< "\n# K | Q0s | mu01 | mu02 | mu03 | alpha1| alpha2 | Q0e |"
<< "\n#-----------------------------------------------------------------------"
<< "\n# 160000 150 17000 220 17000 0.002 0.004 0.01 avec_phase "
<< "\n # nQs= 0.2 gammaQs= 0.5 nQe= 0.2 gammaQe= 0.5 nMu1= 0.2 gammaMu1= 0.5 nMu2= 0.2 gammaMu2= 0.5 nMu3= 0.2 gammaMu3= 0.5 "
<< "\n # "
<< "\n # noter qu'il ne faut pas que mu3 devienne nulle"
<< "\n # "
<< "\n # il est egalement possible de faire dependre certains coefficients a la temperature"
<< "\n # pour l'instant seule les parametres Q0s et Q0e peuvent-etre thermo-dependant, un exemple:"
<< "\n#-----------------------------------------------------------------------"
<< "\n# K | Q0s | mu01 | mu02 | mu03 | alpha1| alpha2 | Q0e |"
<< "\n#-----------------------------------------------------------------------"
<< "\n# 160000 Q0s_thermo_dependant_ 17000 220 17000 0.002 0.004 Q0e_thermo_dependant_ avec_phase "
<< "\n # nQs= 1.1 gammaQs= 1. nQe= 1.1 gammaQe= 1. nMu2= 1.1 gammaMu2= 1. nMu3= 1.1 gammaMu3= 1. "
<< "\n # cas_Q0s_thermo_dependant_ 1 T0r= 308 Gr= 7.5 Qrmax= 570 ar= 9 "
<< "\n # cas_Q0e_thermo_dependant_ 1 T0rQe= 308 Qe_T0rQe= 0.01 h1= 1.2"
<< "\n # "
<< "\n # la valeur de Q0s est remplassee par le mot cle: Q0s_thermo_dependant_,"
<< "\n # idem (eventuellement) pour Q0e qui est remplace par Q0e_thermo_dependant_, ensuite apres la definition"
<< "\n # de tous les autres parametres, sur la ligne qui suit, on definit la thermodependance:"
<< "\n # 1) tout d'abord pour Q0s"
<< "\n # cas_Q0s_thermo_dependant_ : indique le cas de la thermodependance, "
<< "\n # = 1: indique que l'on a une fonction pre-programmee (voir doc) donnee par Laurent orgeas, "
<< "\n # dependant des parametre qui suivent "
<< "\n # = 2: indique Q0s depend d'une fonction Q0s(T) que l'on indique ensuite ex: "
<< "\n # cas_Q0s_thermo_dependant_ 2 courbe_evol_Q0s "
<< "\n # ici courbe_evol_Q0s est le nom d'une courbe qui doit avoir ete definit auparavant "
<< "\n # "
<< "\n # 2) puis eventuellement pour Q0e, dans ce cas on passe une ligne puis: "
<< "\n # cas_Q0e_thermo_dependant_ : indique le cas de la thermodependance, "
<< "\n # = 1: indique que l'on a la fonction pre-programmee suivante: "
<< "\n # la thermo-dépendance est fonction de la pente mu3+mu2 -> Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2)) "
<< "\n # dependant des parametres T0rQe Qe_T0rQe et h1 "
<< "\n # ex: cas_Q0e_thermo_dependant_ 1 T0rQe= 273 Qe_T0rQe= 0.07 h1= 1.2 "
<< "\n # second cas: = 2: indique Q0e depend d'une fonction Q0e(T) que l'on indique ensuite ex: "
<< "\n # cas_Q0e_thermo_dependant_ 2 courbe_evol_Q0e "
<< "\n # ici courbe_evol_Q0e est le nom d'une courbe qui doit avoir ete definit auparavant "
<< "\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#------------------------------------------------------------------------------------"
<< "\n # ";
};
sort << 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 IsoHyper3DOrgeas1::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
if ((K == ConstMath::tresgrand) || (Q0s == ConstMath::tresgrand) ||
(mu01 == ConstMath::tresgrand) || (mu02 == ConstMath::tresgrand)
|| (mu03 == ConstMath::tresgrand) || (alpha1 == ConstMath::tresgrand)
|| (alpha3 == ConstMath::tresgrand) || (Q0e == ConstMath::tresgrand))
{ cout << " \n Les parametres ne sont pas defini pour la loi " << Nom_comp(id_comp)
<< '\n';
Affiche();
ret = 0;
}
if (cas_Q0s_thermo_dependant != 0)
{ switch (cas_Q0s_thermo_dependant)
{ case 1: if ((T0r == ConstMath::tresgrand) || (Gr == ConstMath::tresgrand) ||
(Qrmax == ConstMath::tresgrand) || (ar == ConstMath::tresgrand))
{ cout << " \n un des parametres T0r, Gr, Qrmax ou ar n'est pas defini pour la loi "
<< Nom_comp(id_comp) << '\n';
Affiche();ret = 0;
}
break;
case 2: if (Q0s_temperature == NULL)
{ cout << " \n la loi Q0s(T) n'est pas defini pour la loi "
<< Nom_comp(id_comp) << '\n';
Affiche();ret = 0;
}
break;
default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
<< " cas non pris en compte !!!! dans Affiche() "
<< "\n IsoHyper3DOrgeas1::TestComplet()";
ret = 0.; Sortie(1);
};
};
if (cas_Q0e_thermo_dependant != 0)
{ switch (cas_Q0e_thermo_dependant)
{ case 1: if ((T0rQe == ConstMath::tresgrand) || (Qe_T0rQe == ConstMath::tresgrand)
|| (h1 == ConstMath::tresgrand) )
{ cout << " \n un des parametres T0rQe, Qe_T0rQe ou h1 n'est pas defini pour la loi "
<< Nom_comp(id_comp) << '\n';
Affiche();ret = 0;
}
break;
case 2: if (Q0e_temperature == NULL)
{ cout << " \n la loi Q0e(T) n'est pas defini pour la loi "
<< Nom_comp(id_comp) << '\n';
Affiche();ret = 0;
}
break;
default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
<< " cas non pris en compte !!!! dans Affiche() "
<< "\n IsoHyper3DOrgeas1::TestComplet()";
ret = 0.; Sortie(1);
};
};
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 IsoHyper3DOrgeas1::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string toto;
if (cas == 1)
{ ent >> toto >> K >> toto >> cas_Q0s_thermo_dependant;
if (cas_Q0s_thermo_dependant != 0)
{ thermo_dependant=true; // on indique que la loi est thermo-dépendant
switch (cas_Q0s_thermo_dependant)
{ case 1:{ent >> toto >> T0r >> toto >> Gr >> toto >> Qrmax >> toto >> ar ;
// on calcul maintenant les variables intermédiaires fixes
Tc=T0r+(Qrmax/Gr/2.);
Ts=(ar*ar- Sqr(Qrmax/Gr) )/2./(Qrmax/Gr);
break;
}
case 2:{ ent >> toto;
if (toto != " Q0s(T) ")
{ cout << "\n erreur en lecture de la fonction temperature, on attendait Q0s(T) et on a lue "
<< toto
<< "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(...";
Sortie(1);
};
Q0s_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Q0s_temperature);
break;
};
default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
<< " cas non pris en compte !!!! dans Lecture_base_info_loi(.. "
<< "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(..";
Sortie(1);
};
}
else
{ent >> toto >> Q0s ;};
// puis les autres para sauf Qe
ent >> toto >> Q0s >> toto >> mu01 >> toto >> mu02 >> toto >> mu03
>> toto >> alpha1 >> toto >> alpha3 ;
// Qe dépendant éventuellement de la température
ent >> toto >> cas_Q0e_thermo_dependant;
if (cas_Q0e_thermo_dependant != 0)
{ thermo_dependant=true; // on indique que la loi est thermo-dépendant
switch (cas_Q0e_thermo_dependant)
{ case 1:{ent >> toto >> T0rQe >> toto >> Qe_T0rQe >> toto >> h1 ;
break;
}
case 2:{ ent >> toto;
if (toto != " Q0e(T) ")
{ cout << "\n erreur en lecture de la fonction temperature, on attendait Q0e(T) et on a lue "
<< toto
<< "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(...";
Sortie(1);
};
Q0e_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,Q0e_temperature);
break;
};
default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
<< " cas non pris en compte !!!! dans Lecture_base_info_loi(.. "
<< "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(..";
Sortie(1);
};
// maintenant on calcule la valeur de Qs_T0rQr
{Qs_T0rQe = Q0s ; // au cas où il n'y a pas de dépendance de Qs à la température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if (T0rQe >= Tc)
{Qs_T0rQe=0.5*Gr*( (T0rQe-Tc+Ts)-sqrt(Sqr(T0rQe-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs_T0rQe=0.5*Gr*( (T0rQe-Tc-Ts)+sqrt(Sqr(T0rQe-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs_T0rQe = Q0s_temperature->Valeur(T0rQe);
break;
}
default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
<< " cas non pris en compte !!!! dans Lecture_base_info_loi(.. "
<< "\n IsoHyper3DOrgeas1::Lecture_base_info_loi(..";
Sortie(1);
};
};
}
else
{ent >> toto >> Q0e ;};
// puis les dépendances éventuelles à la phase
ent >> toto >> avec_phase;
if (avec_phase)
{ ent >> toto >> nQs >> toto >> gammaQs
>> toto >> nQe >> toto >> gammaQe
>> toto >> nMu1 >> toto >> gammaMu1
>> toto >> nMu2 >> toto >> gammaMu2
>> toto >> nMu3 >> toto >> gammaMu3;
};
// calcul des constantes particulières
alpha1_2 = Sqr(alpha1);
alpha3_2 = Sqr(alpha3);
// mu1_2 = Sqr(mu1);
// gestion du post-traitement
ent >> toto >> sortie_post ;
};
// appel class 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 IsoHyper3DOrgeas1::Ecriture_base_info_loi(ofstream& sort,const int cas)
{ if (cas == 1)
{ sort << " module_dilatation " << K << " cas_Q0s_thermo_dependant " << cas_Q0s_thermo_dependant;
if (cas_Q0s_thermo_dependant != 0)
{ sort << " cas_Q0s_thermo_dependant " << cas_Q0s_thermo_dependant;
switch (cas_Q0s_thermo_dependant)
{ case 1: sort << " T0r= " << T0r << " Gr= " << Gr << " Qrmax= " << Qrmax << " ar= " << ar << " ";
break;
case 2: {sort << " Q0s(T) ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,Q0s_temperature);
};
break;
default: cout << "\n erreur cas_Q0s_thermo_dependant= " << cas_Q0s_thermo_dependant
<< " cas non pris en compte !!!! dans Ecriture_base_info_loi(.. "
<< "\n IsoHyper3DOrgeas1::Ecriture_base_info_loi(..";
Sortie(1);
};
}
else
{sort << " Q0s " << Q0s ;};
// puis les autres para sauf Qe
sort << " mu01 " << mu01 << " mu02 " << mu02 << " mu03 " << mu03
<< " alpha1 " << alpha1 << " alpha3 " << alpha3 << " ";
// Qe eventuellement thermo dépendant
sort << " cas_Q0e_thermo_dependant " << cas_Q0e_thermo_dependant;
if (cas_Q0e_thermo_dependant != 0)
{ sort << " cas_Q0e_thermo_dependant " << cas_Q0e_thermo_dependant;
switch (cas_Q0e_thermo_dependant)
{ case 1: sort << " T0rQe= " << T0rQe << " Qe_T0rQe= " << Qe_T0rQe << " h1= " << h1 << " ";
break;
case 2: {sort << " Q0e(T) ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,Q0e_temperature);
};
break;
default: cout << "\n erreur cas_Q0e_thermo_dependant= " << cas_Q0e_thermo_dependant
<< " cas non pris en compte !!!! dans Ecriture_base_info_loi(.. "
<< "\n IsoHyper3DOrgeas1::Ecriture_base_info_loi(..";
Sortie(1);
};
}
else
{sort << " Q0e " << Q0e ;};
// puis les dépendances éventuelles à la phase
sort << " avec_phase " << avec_phase;
if (avec_phase)
{ sort << "\n nQs= " << nQs << " gammaQs= " << gammaQs
<< "\n nQe= " << nQe << " gammaQe= " << gammaQe
<< "\n nMu1= " << nMu1 << " gammaMu1= " << gammaMu1
<< "\n nMu2= " << nMu2 << " gammaMu2= " << gammaMu2
<< "\n nMu3= " << nMu3 << " gammaMu3= " << gammaMu3;
};
// 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
double IsoHyper3DOrgeas1::Module_young_equivalent(Enum_dure ,const Deformation &,SaveResul * )
{ return (9.*K*(mu01+mu02)/((mu01+mu02)+3.*K));
};
// 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
double IsoHyper3DOrgeas1::Module_compressibilite_equivalent(Enum_dure ,const Deformation & def)
{ return K/3.;};
// =========== 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 IsoHyper3DOrgeas1::PoGrenoble
(const double & Qeps,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01;double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = Qeps * Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
double E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(Qeps*(Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
// retour
return E;
};
// calcul du potentiel tout seul avec la phase donc dans le cas où Qeps est non nul
double IsoHyper3DOrgeas1::PoGrenoble
(const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
if (avec_phase)
{Qs = Q0sT / pow((1.+gammaQs * inv_spec.cos3phi),nQs);
Qe = Q0eT / pow((1.+gammaQe * inv_spec.cos3phi),nQe);
mu1 = mu01 / pow((1.+gammaMu1 * inv_spec.cos3phi),nMu1);
mu2 = mu02 / pow((1.+gammaMu2 * inv_spec.cos3phi),nMu2);
mu3 = mu03 / pow((1.+gammaMu3 * inv_spec.cos3phi),nMu3);
};
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = inv_spec.Qeps * inv_spec.Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (inv_spec.Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = inv_spec.Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel
double E = K/6. * (logV)*(logV) + Qs * inv_spec.Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* inv_spec.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(inv_spec.Qeps*(inv_spec.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
// 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 IsoHyper3DOrgeas1::PoGrenoble_et_V
(const double & Qeps,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
PoGrenoble_V ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = Qeps * Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(Qeps*(Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
// dérivées premières
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 IsoHyper3DOrgeas1::PoGrenoble_et_V
(const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
if (avec_phase)
{Qs = Q0sT / pow((1.+gammaQs * inv_spec.cos3phi),nQs);
Qe = Q0eT / pow((1.+gammaQe * inv_spec.cos3phi),nQe);
mu1 = mu01 / pow((1.+gammaMu1 * inv_spec.cos3phi),nMu1);
mu2 = mu02 / pow((1.+gammaMu2 * inv_spec.cos3phi),nMu2);
mu3 = mu03 / pow((1.+gammaMu3 * inv_spec.cos3phi),nMu3);
};
PoGrenoble_V ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps = inv_spec.Qeps; // pour simplifier la notation
double Qeps2 = Qeps * Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(Qeps*(Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
// dérivées premières
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 IsoHyper3DOrgeas1::PoGrenoble_et_VV
(const double & Qeps,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
PoGrenoble_VV ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = Qeps * Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(Qeps*(Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
/* double toto = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2;
toto = + mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(QplusQec/alpha1_2 * S2) - log(Qec/alpha1_2 * S1)));
cout << "toto = " << toto;
toto = + mu3/2.*(Qeps*(Qeps-2. * S4));
cout << "toto = " << toto;
toto= + alpha3_2*(log(QmoinsQe/alpha3_2 * S3) - log(-Qe/alpha3_2 * S4));
cout << "toto = " << toto;
toto = + QmoinsQe * S3 + Qe * S4;
cout << "toto = " << toto; */
// dérivées premières
ret.EV = K/3.*logV/inv.V;
// dérivées secondes
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
//++++verif***
//Verif_PoGrenoble_et_var_VV(Qeps,inv,ret );
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel et de sa variation première et seconde suivant V uniquement
// on prend en compte une dépendance éventuelle à la phase
Hyper3D::PoGrenoble_VV IsoHyper3DOrgeas1::PoGrenoble_et_VV
(const Invariant0QepsCosphi & inv_spec,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
if (avec_phase)
{Qs = Q0sT / pow((1.+gammaQs * inv_spec.cos3phi),nQs);
Qe = Q0eT / pow((1.+gammaQe * inv_spec.cos3phi),nQe);
mu1 = mu01 / pow((1.+gammaMu1 * inv_spec.cos3phi),nMu1);
mu2 = mu02 / pow((1.+gammaMu2 * inv_spec.cos3phi),nMu2);
mu3 = mu03 / pow((1.+gammaMu3 * inv_spec.cos3phi),nMu3);
};
PoGrenoble_VV ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps = inv_spec.Qeps; // pour simplifier la notation
double Qeps2 = Qeps * Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(Qeps*(Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
// dérivées premières
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 premières non compris la phase
Hyper3D::PoGrenobleSansPhaseSansVar IsoHyper3DOrgeas1::PoGrenoble
(const InvariantQeps & i_s,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
PoGrenobleSansPhaseSansVar ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = i_s.Qeps * i_s.Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (i_s.Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = i_s.Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4); // ok
// dérivées premières
ret.EV = K/3.*logV/inv.V;
double dS2_Q = QplusQec / S2;
double dS3_Q = QmoinsQe / S3;
ret.EQ = Qs + 2.*mu2*i_s.Qeps
+ 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2
-alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
+ 0.5*mu3*(2.*(i_s.Qeps - S4)
+alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
+S3+QmoinsQe*dS3_Q);
ret.Ks = K / 3. *(logV/2. + 1.);
// retour
return ret;
};
// calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PoGrenobleAvecPhaseSansVar IsoHyper3DOrgeas1::PoGrenoblePhase
(const InvariantQepsCosphi& i_s,const Invariant & inv)
{ // les grandeurs de base,
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
if (avec_phase)
{Qs = Q0sT / pow((1.+gammaQs * i_s.cos3phi),nQs);
Qe = Q0eT / pow((1.+gammaQe * i_s.cos3phi),nQe);
mu1 = mu01 / pow((1.+gammaMu1 * i_s.cos3phi),nMu1);
mu2 = mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2);
mu3 = mu03 / pow((1.+gammaMu3 * i_s.cos3phi),nMu3);
};
double mu1_2 = Sqr(mu1);
PoGrenobleAvecPhaseSansVar ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = i_s.Qeps * i_s.Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (i_s.Qeps + Qec);
double QplusQec_2 = Sqr(QplusQec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = i_s.Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4); // ok
// dérivées premières
ret.EV = K/3.*logV/inv.V;
double dS2_Q = QplusQec / S2;
double dS3_Q = QmoinsQe / S3;
ret.EQ = Qs + 2.*mu2*i_s.Qeps
+ 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2
-alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
+ 0.5*mu3*(2.*(i_s.Qeps - S4)
+alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
+S3+QmoinsQe*dS3_Q);
// + mu1 *(QplusQec - S2) + mu3 * (i_s.Qeps - S4 + S3);
// cas de la variation faisant intervenir la phase
// 1-- calcul des variations des grandeurs de bases
double dQs_dcos3phi = - nQs * gammaQs * Qs / (1.+gammaQs * i_s.cos3phi);
double dQe_dcos3phi = - nQe * gammaQe * Qe / (1.+gammaQe * i_s.cos3phi);
double dMu1_dcos3phi = - nMu1 * gammaMu1 * mu1 / (1.+gammaMu1 * i_s.cos3phi);
double dMu2_dcos3phi = - nMu2 * gammaMu2 * mu2 / (1.+gammaMu2 * i_s.cos3phi);
double dMu3_dcos3phi = - nMu3 * gammaMu3 * mu3 / (1.+gammaMu3 * i_s.cos3phi);
// 2-- variations des grandeurs intermédiaires par rapport aux grandeurs de bases
double Qs2 = Qs * Qs;
double dQec_Qs = - (Qs2 + mu1_2 * alpha1_2) / (2. * mu1 * Qs2); // OK
double dS2_Qs = QplusQec * dQec_Qs / S2; // OK
double dS1_Qs = Qec * dQec_Qs / S1; // OK
double dS4_Qe = Qe / S4; // ok
double dS3_Qe = -(QmoinsQe) / S3; // OK
double dQec_Mu1 = (Qs2 + mu1_2 * alpha1_2) / (2. * mu1_2 * Qs);
double dQec2_Mu1 = 2. * Qec * dQec_Mu1;
double dQplusQec_Mu1 = dQec_Mu1; // pourra être supprimé par la suite mais pour la mise au point ... c'est plus cohérent
double dS1_Mu1 = dQec_Mu1 * Qec/ S1;
double dS2_Mu1 = dQec_Mu1 * QplusQec/ S2;
// 3-- variations du potentiel par rapport aux grandeurs de bases
// A_x,B_x,C_x sont des variables intermédiaires
double A_x = (dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
-alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
);
double dE_Qs = i_s.Qeps + 0.5*mu1*A_x;
// + 0.5*mu1*(dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
// -alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
// ); // OK
double dE_Qe = 0.5*mu3*(-2.*i_s.Qeps*dS4_Qe
+alpha3_2*((-1.+dS3_Qe)/(QmoinsQe + S3)-(-1.+dS4_Qe)/(-Qe+S4))
- S3 + S4 + Qe*dS4_Qe + QmoinsQe * dS3_Qe); // OK
double B_x = 2.* i_s.Qeps * dQec_Mu1 - (dQplusQec_Mu1 * S2 + QplusQec * dS2_Mu1) + dQec_Mu1* S1 + Qec* dS1_Mu1
- alpha1_2*( ((dQplusQec_Mu1+dS2_Mu1 )) / (QplusQec+S2)
- ( (dQec_Mu1+dS1_Mu1)) / (Qec+S1)
);
double C_x = Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Dabs(QplusQec + S2)) - log(Dabs(Qec + S1)));
double dE_Mu1 = 0.5 * ( C_x + mu1 * B_x ); // ok
double dE_Mu2 = i_s.Qeps * i_s.Qeps; // OK
double dE_Mu3 = 0.5 * (i_s.Qeps*(i_s.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4); // OK
// 4-- variations du potentiel par rapport à cos3phi
ret.Ecos3phi = dE_Qs * dQs_dcos3phi + dE_Qe * dQe_dcos3phi
+ dE_Mu1 * dMu1_dcos3phi + dE_Mu2 * dMu2_dcos3phi + dE_Mu3 * dMu3_dcos3phi; // OK
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 IsoHyper3DOrgeas1::PoGrenoble_et_var
(const Invariant2Qeps& i_s,const Invariant & inv)
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
PoGrenobleSansPhaseAvecVar ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = i_s.Qeps * i_s.Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (i_s.Qeps + Qec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double QmoinsQe = i_s.Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4);
// dérivées premières
ret.EV = K/3.*logV/inv.V;
double dS2_Q = QplusQec / S2;
double dS3_Q = QmoinsQe / S3;
ret.EQ = Qs + 2.*mu2*i_s.Qeps
+ 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2
-alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
+ 0.5*mu3*(2.*(i_s.Qeps - S4)
+alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
+S3+QmoinsQe*dS3_Q);
// ret.EQ = Qs + 2.*mu2*i_s.Qeps
// + mu1 *(QplusQec - S2) + mu3 * (i_s.Qeps - S4 + S3);
// dérivées secondes
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V);
ret.EQV = 0.;
// double dS2_QQ = (S2-dS2_Q)/(S2*S2);
// double dS3_QQ = (S3-dS3_Q)/(S3*S3);
double dS2_QQ = (1. - dS2_Q*dS2_Q)/S2; // OK
double dS3_QQ = (1. - dS3_Q*dS3_Q)/S3; // OK
ret.EQQ = 2. * mu2
+ 0.5*mu1*( 2.*(1-dS2_Q) - QplusQec * dS2_QQ
-alpha1_2*(dS2_QQ/(QplusQec + S2)-Sqr((1.+dS2_Q)/(QplusQec + S2))))
+ 0.5*mu3*(2.
+alpha3_2*(dS3_QQ/(QmoinsQe + S3)-Sqr((1.+dS3_Q)/(QmoinsQe + S3)))
+2.* dS3_Q + QmoinsQe * dS3_QQ);
// ret.EQQ = 2. * mu2 + mu1*(1.-QplusQec*S2) + mu3*(1. + QmoinsQe*S3);
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
//++++verif***
// Verif_PoGrenoble_et_var(i_s.Qeps,inv,ret );
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 IsoHyper3DOrgeas1::PoGrenoblePhase_et_var
(const Invariant2QepsCosphi& i_s,const Invariant & inv)
{ // les grandeurs de base,
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
double Q0sT=Q0s; double Q0eT=Q0e; // pour la dépendance éventuelle à la phase et température
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs = Q0sT=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0sT = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Q0eT = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0eT = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas de l'existence de la phase, on modifie les variables qui dépendent de la phase
if (avec_phase)
{Qs = Q0sT / pow((1.+gammaQs * i_s.cos3phi),nQs);
Qe = Q0eT / pow((1.+gammaQe * i_s.cos3phi),nQe);
mu1 = mu01 / pow((1.+gammaMu1 * i_s.cos3phi),nMu1);
mu2 = mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2);
mu3 = mu03 / pow((1.+gammaMu3 * i_s.cos3phi),nMu3);
//cout << "\n cos(3phi)= " << i_s.cos3phi << " Qs= " << Qs << " mu2= " << mu2;
};
double mu1_2 = Sqr(mu1);
double Qs_2 = Sqr(Qs);
PoGrenobleAvecPhaseAvecVar ret;
// des variables intermédiaires
double logV = log(inv.V);
double Qeps2 = i_s.Qeps * i_s.Qeps;
double Qe2 = Qe * Qe;
double Qec = (mu1 * alpha1 * mu1 * alpha1 - Qs * Qs) / (2. * mu1 * Qs);
double Qec2 = Sqr(Qec);
double QplusQec = (i_s.Qeps + Qec);
double QplusQec_2 = Sqr(QplusQec);
double S1 = sqrt(alpha1_2 + Qec2);
double S2 = sqrt(alpha1_2 + Sqr(QplusQec));
double S2_2 = Sqr(S2);
double QmoinsQe = i_s.Qeps - Qe;
double S3 = sqrt(alpha3_2 + Sqr(QmoinsQe));
double S4 = sqrt(alpha3_2+Qe2);
// le potentiel et ses dérivées premières
ret.E = K/6. * (logV)*(logV) + Qs * i_s.Qeps + mu2 * Qeps2
+ mu1/2. * (Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Abs(QplusQec + S2)) - log(Abs(Qec + S1))))
+ mu3/2.*(i_s.Qeps*(i_s.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4); // ok
// +++++++++++++++++++ dérivées premières +++++++++++++++++++
// a) la partie sans phase
ret.EV = K/3.*logV/inv.V; // ok
// --- pour le débug
// switch(fpclassify(ret.EV))
// { case FP_NAN: case FP_INFINITE:
// { // cas d'un nombre infini
// cout << "\n erreur valeur nan pour EV! ";
// };
// break;
// };
// --- fin débug
double dS2_Q = QplusQec / S2; // OK
double dS3_Q = QmoinsQe / S3; // OK
// D_x est une variable intermédiaire, utilisée par la suite
double D_x = QplusQec*(2.-dS2_Q) -S2 - alpha1_2*((1.+dS2_Q)/(QplusQec + S2));
ret.EQ = Qs + 2.*mu2*i_s.Qeps
+ 0.5 * mu1 * D_x
// + 0.5*mu1*(QplusQec*(2.-dS2_Q) -S2
// -alpha1_2*((1.+dS2_Q)/(QplusQec + S2)))
+ 0.5*mu3*(2.*(i_s.Qeps - S4)
+alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
+S3+QmoinsQe*dS3_Q); // OK
// b) cas de la variation faisant intervenir la phase
// 1-- calcul des variations des grandeurs de bases
double dQs_dcos3phi = - nQs * gammaQs * Qs / (1.+gammaQs * i_s.cos3phi); // 0K
double dQe_dcos3phi = - nQe * gammaQe * Qe / (1.+gammaQe * i_s.cos3phi); // OK
double dMu1_dcos3phi = - nMu1 * gammaMu1 * mu1 / (1.+gammaMu1 * i_s.cos3phi); // OK
double dMu2_dcos3phi = - nMu2 * gammaMu2 * mu2 / (1.+gammaMu2 * i_s.cos3phi); // OK
// double toto = - nMu2 * gammaMu2 * mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2+1);
// if (Abs(toto -dMu2_dcos3phi) > ConstMath::unpeupetit)
// cout << "\n pb dMu2_dcos3phi" << toto << " " << dMu2_dcos3phi << endl;
double dMu3_dcos3phi = - nMu3 * gammaMu3 * mu3 / (1.+gammaMu3 * i_s.cos3phi); // OK
// 2-- variations des grandeurs intermédiaires par rapport aux grandeurs de bases
double Qs2 = Qs * Qs;
double dQec_Qs = - (Qs2 + mu1_2 * alpha1_2) / (2. * mu1 * Qs2); // OK
double dS2_Qs = QplusQec * dQec_Qs / S2; // OK
double dS1_Qs = Qec * dQec_Qs / S1; // OK
double dS4_Qe = Qe / S4; // ok
double dS3_Qe = -(QmoinsQe) / S3; // OK
double dQec_Mu1 = (Qs2 + mu1_2 * alpha1_2) / (2. * mu1_2 * Qs);
double dQec2_Mu1 = 2. * Qec * dQec_Mu1;
double dQplusQec_Mu1 = dQec_Mu1; // pourra être supprimé par la suite mais pour la mise au point ... c'est plus cohérent
double dS1_Mu1 = dQec_Mu1 * Qec/ S1;
double dS2_Mu1 = dQec_Mu1 * QplusQec/ S2;
// 3-- variations du potentiel par rapport aux grandeurs de bases
// A_x,B_x,C_x sont des variables intermédiaires
double A_x = (dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
-alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
);
double dE_Qs = i_s.Qeps + 0.5*mu1*A_x;
// + 0.5*mu1*(dQec_Qs*(2.*i_s.Qeps-S2+S1)-QplusQec*dS2_Qs+Qec*dS1_Qs
// -alpha1_2*((dQec_Qs+dS2_Qs)/(QplusQec + S2)-(dQec_Qs+dS1_Qs)/(Qec+S1) )
// ); // OK
double dE_Qe = 0.5*mu3*(-2.*i_s.Qeps*dS4_Qe
+alpha3_2*((-1.+dS3_Qe)/(QmoinsQe + S3)-(-1.+dS4_Qe)/(-Qe+S4))
- S3 + S4 + Qe*dS4_Qe + QmoinsQe * dS3_Qe); // OK
double B_x = 2.* i_s.Qeps * dQec_Mu1 - (dQplusQec_Mu1 * S2 + QplusQec * dS2_Mu1) + dQec_Mu1* S1 + Qec* dS1_Mu1
- alpha1_2*( ((dQplusQec_Mu1+dS2_Mu1 )) / (QplusQec+S2)
- ( (dQec_Mu1+dS1_Mu1)) / (Qec+S1)
);
double C_x = Qeps2 + 2.* i_s.Qeps * Qec - QplusQec * S2 + Qec* S1
- alpha1_2*(log(Dabs(QplusQec + S2)) - log(Dabs(Qec + S1)));
double dE_Mu1 = 0.5 * ( C_x + mu1 * B_x ); // ok
double dE_Mu2 = i_s.Qeps * i_s.Qeps; // OK
double dE_Mu3 = 0.5 * (i_s.Qeps*(i_s.Qeps-2. * S4)
+ alpha3_2*(log(Abs(QmoinsQe + S3)) - log(Abs(-Qe + S4)))
+ QmoinsQe * S3 + Qe * S4); // OK
// 4-- variations du potentiel par rapport à cos3phi
ret.Ecos3phi = dE_Qs * dQs_dcos3phi + dE_Qe * dQe_dcos3phi
+ dE_Mu1 * dMu1_dcos3phi + dE_Mu2 * dMu2_dcos3phi + dE_Mu3 * dMu3_dcos3phi; // OK
// --- pour le débug
// switch(fpclassify(ret.Ecos3phi))
// { case FP_NAN: case FP_INFINITE:
// { // cas d'un nombre infini
// cout << "\n erreur valeur nan pour EV! ";
// };
// break;
// };
// --- fin débug
// +++++++++++++++++++ dérivées secondes +++++++++++++++++++
// a) la partie sans phase
ret.EVV = K/3. * (1.-logV) / (inv.V * inv.V); // ok
ret.EQV = 0.;
double dS2_QQ = (1. - dS2_Q*dS2_Q)/S2; // OK
double dS3_QQ = (1. - dS3_Q*dS3_Q)/S3; // OK
ret.EQQ = 2. * mu2
+ 0.5*mu1*( 2.*(1-dS2_Q) - QplusQec * dS2_QQ
-alpha1_2*(dS2_QQ/(QplusQec + S2)-Sqr((1.+dS2_Q)/(QplusQec + S2))))
+ 0.5*mu3*(2.
+alpha3_2*(dS3_QQ/(QmoinsQe + S3)-Sqr((1.+dS3_Q)/(QmoinsQe + S3)))
+2.* dS3_Q + QmoinsQe * dS3_QQ); // OK
// ret.EQQ = 2. * mu2 + mu1*(1.-QplusQec*S2) + mu3*(1. + QmoinsQe*S3);
// b) la partie dépendant de la phase
// 1-- calcul des variations seconde des grandeurs de bases
double dQs_dcos3phi2 = - (nQs+1) * gammaQs * dQs_dcos3phi / (1.+gammaQs * i_s.cos3phi); // OK
double dQe_dcos3phi2 = - (nQe+1) * gammaQe * dQe_dcos3phi / (1.+gammaQe * i_s.cos3phi); // OK
double dMu1_dcos3phi2 = - (nMu1+1) * gammaMu1 * dMu1_dcos3phi / (1.+gammaMu1 * i_s.cos3phi); // OK
double dMu2_dcos3phi2 = - (nMu2+1) * gammaMu2 * dMu2_dcos3phi / (1.+gammaMu2 * i_s.cos3phi); // OK
// toto = nMu2 *(nMu2+1)* gammaMu2*gammaMu2 * mu02 / pow((1.+gammaMu2 * i_s.cos3phi),nMu2+2);
// if (Abs(toto -dMu2_dcos3phi2) > ConstMath::unpeupetit)
// cout << "\n pb dMu2_dcos3phi2" << toto << " " << dMu2_dcos3phi2 << endl;
double dMu3_dcos3phi2 = - (nMu3+1) * gammaMu3 * dMu3_dcos3phi / (1.+gammaMu3 * i_s.cos3phi); // OK
// 2-- variations seconde des grandeurs intermédiaires par rapport aux grandeurs de bases
double dQec_QsQs = mu1 * alpha1_2 / (Qs2 * Qs); // OK
double dS2_QsQs = ((dQec_Qs*dQec_Qs + QplusQec * dQec_QsQs) - Sqr(QplusQec*dQec_Qs/S2))/S2; // OK
double dS1_QsQs = ((dQec_Qs*dQec_Qs + Qec * dQec_QsQs) - Sqr(Qec*dQec_Qs/S1))/S1; // OK
double dS4_QeQe = (1.-Sqr(Qe/S4))/S4; // OK
double dS3_QeQe = (1.-Sqr(QmoinsQe/S3))/S3; // OK
double dS2_QQs = (dQec_Qs - (QplusQec) * dS2_Qs / S2) / S2; // OK
double dS3_QQe = (-1. - QmoinsQe * dS3_Qe / S3) / S3 ; // OK
double dQec_Mu1Mu1 = -Qs / (mu1 * mu1_2); // ok
double dQec_QsMu1 = -alpha1_2/Qs_2 + (Qs_2+mu1_2*alpha1_2) / (2.*mu1_2*Qs_2); // OK
double dQec_Mu1Q = 0.;
double dS1_Mu1Mu1 = ((dQec_Mu1*dQec_Mu1 + Qec*dQec_Mu1Mu1) -Sqr(Qec*dQec_Mu1)/Sqr(S1)) / S1; // ok
double dS1_QsMu1 = ((dQec_Mu1*dQec_Qs+Qec*dQec_QsMu1)-dQec_Mu1*dQec_Qs*Sqr(Qec/S1)) / S1; // ok
double dS2_Mu1Mu1 = dQec_Mu1 * dQec_Mu1 * (1.-QplusQec_2/S2_2)/S2 + QplusQec/S2 * dQec_Mu1Mu1; // ok
double dS2_QsMu1 = dQec_Mu1 * dQec_Qs * (1.-QplusQec_2/S2_2)/S2 + QplusQec/S2 * dQec_QsMu1; // ok
double dS2_QMu1 = dQec_Mu1 * (1.-QplusQec_2/S2_2)/S2; // ok
// 3-- variations seconde du potentiel par rapport aux grandeurs de bases
double dE_QsQs = 0.5*mu1*(dQec_QsQs*(2.*i_s.Qeps-S2+S1) - 2.*dQec_Qs*dS2_Qs -QplusQec*dS2_QsQs
+2.*dQec_Qs*dS1_Qs + Qec * dS1_QsQs
-alpha1_2*((dQec_QsQs+dS2_QsQs)/(QplusQec+S2)-Sqr((dQec_Qs+dS2_Qs)/(QplusQec+S2))
-(dQec_QsQs+dS1_QsQs)/(Qec+S1)+Sqr((dQec_Qs+dS1_Qs)/(Qec+S1))
)
); // OK
double dE_QeQe = 0.5*mu3*(-2.*i_s.Qeps*dS4_QeQe
+alpha3_2*((dS3_QeQe)/(QmoinsQe + S3)-Sqr((-1.+dS3_Qe)/(QmoinsQe + S3))
-(dS4_QeQe)/(-Qe + S4)+Sqr((-1.+dS4_Qe)/(-Qe + S4))
)
-2.*dS3_Qe + 2.*dS4_Qe + Qe*dS4_QeQe + QmoinsQe*dS3_QeQe
); // OK
double dE_Mu3Qe = dE_Qe/mu3; // OK
// puis pour la dérivée croisée
double dE_QQs = 1. + 0.5*mu1*(dQec_Qs*(2.-dS2_Q)-dS2_Qs-QplusQec*dS2_QQs
-alpha1_2*( (dS2_QQs-(1.+dS2_Q)*(dQec_Qs+dS2_Qs)/(QplusQec+S2))/(QplusQec+S2))
); // OK
double dE_QQe = 0.5*mu3*(-2.*dS4_Qe
+alpha3_2*((dS3_QQe-(1.+dS3_Q)*(-1.+dS3_Qe)/(QmoinsQe+S3))/(QmoinsQe+S3))
+dS3_Qe - dS3_Q + QmoinsQe*dS3_QQe
); // OK
double dE_QMu2 = 2.* i_s.Qeps; // OK
double dE_QMu3 = 0.5 * (2.*i_s.Qeps - 2. * S4
+alpha3_2*((1.+dS3_Q)/(QmoinsQe + S3))
+S3+QmoinsQe*dS3_Q); // OK
double dE_QsMu1 = 0.5 * A_x
+ 0.5*mu1 * (dQec_QsMu1*(2.*i_s.Qeps-S2+S1)+dQec_Qs*(-dS2_Mu1+dS1_Mu1)-dQec_Mu1*dS2_Qs-QplusQec*dS2_QsMu1
+dQec_Mu1*dS1_Qs+Qec*dS1_QsMu1
-alpha1_2*( (dQec_QsMu1+dS2_QsMu1)/(QplusQec + S2)-(dQec_Qs+dS2_Qs)*(dQec_Mu1+dS2_Mu1)/Sqr(QplusQec + S2)
- (dQec_QsMu1+dS1_QsMu1)/(Qec + S1) + (dQec_Qs+dS1_Qs)*(dQec_Mu1+dS1_Mu1)/Sqr(Qec+S1))
); // ok
double dE_QsQe = 0.;
double dE_Mu1Mu1 = B_x // en fait ici il y a la somme de deux termes 0.5 * B
+ 0.5 * mu1 * ( dQec_Mu1Mu1 * (2.* i_s.Qeps - S2 + S1) - 2.*dQec_Mu1*(dS2_Mu1 - dS1_Mu1)
- QplusQec * dS2_Mu1Mu1 + Qec * dS1_Mu1Mu1
- alpha1_2*( (dQec_Mu1Mu1+dS2_Mu1Mu1 ) / (QplusQec+S2) - Sqr( (dQec_Mu1+dS2_Mu1) / (QplusQec+S2))
- (dQec_Mu1Mu1+dS1_Mu1Mu1)/(Qec + S1) + Sqr((dQec_Mu1+dS1_Mu1)/(Qec+S1)))
); // ok
double dE_QMu1 = 0.5 * D_x
+ 0.5 * mu1 * (dQec_Mu1 * (2.-dS2_Q) - dS2_Mu1 - QplusQec * dS2_QMu1
-alpha1_2*( (dS2_QMu1)/(QplusQec+S2)
-(dQec_Mu1+dS2_Mu1)*(1.+dS2_Q)/Sqr(QplusQec+S2)
)
); // ok
// 4-- variations seconde du potentiel par rapport à cos3phi
// on tiens compte que dE_QsMu2=0, dE_QsMu3=0, dE_QeMu1=, dE_QeMu2=0, dE_Mu1Mu2=0, dE_Mu1Mu3=0,
// et également que dE_Mu2Mu2=0, dE_Mu3Mu3=0,
ret.Ecos3phi2 = dE_QsQs * Sqr(dQs_dcos3phi) + dE_Qs * dQs_dcos3phi2
+ dE_QeQe * Sqr(dQe_dcos3phi) + dE_Qe * dQe_dcos3phi2
+ dE_Mu1Mu1 * Sqr(dMu1_dcos3phi) + dE_Mu1 * dMu1_dcos3phi2
// + dE_Mu2Mu2 * Sqr(dMu2_dcos3phi) + dE_Mu2 * dMu2_dcos3phi2
+ dE_Mu2 * dMu2_dcos3phi2
// + dE_Mu3Mu3 * Sqr(dMu3_dcos3phi) + dE_Mu3 * dMu3_dcos3phi2
+ dE_Mu3 * dMu3_dcos3phi2
+ 2. * dE_Mu3Qe * dQe_dcos3phi * dMu3_dcos3phi
+ 2. * dE_QsQe * dQs_dcos3phi * dQe_dcos3phi
+ 2. * dE_QsMu1 * dQs_dcos3phi * dMu1_dcos3phi; // ok
// 5-- variations coisée seconde du potentiel par rapport à cos3phi et Qepsilon
ret.EVcos3phi= 0.;
ret.EQcos3phi= dE_QQs*dQs_dcos3phi + dE_QQe*dQe_dcos3phi + dE_QMu1*dMu1_dcos3phi+ dE_QMu2*dMu2_dcos3phi + dE_QMu3*dMu3_dcos3phi; // ok
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
//++++verif*** ++++++
// la vérification a été effectuée selon trois étapes:
// 1) développement des formules avec maxima, et comparaison des formules
// 2) vérification par rapport aux dérivées numériques (pas satisfaisant: ne permet pas de voir où est l'erreur)
// 3) calcul de la dérivée totale du potentiel par rapport aux invariants V,Q,cos(3phi) (en calcul formel), et évaluation numérique
// du résultat pour différentes valeurs de (V,Q,cos(3phi)) et comparaison avec le calcul donné par herezh
// --> extrémement puissant, à permis de trouver toutes les erreurs, et on peut être a peut-près certain
// qu'actuellement il n'y a plus d'erreur (jamais complètement certain ...)
// la méthode CompaFormel() est utilisée pour cela
// Verif_PoGrenoble_et_var(i_s.Qeps,inv,ret );
// Verif_PoGrenoble_et_var(i_s.Qeps,inv,i_s.cos3phi,ret );
// --- Invariant2Qeps toto(inv_spec.Qeps,inv_spec.dQepsdbIIb,inv_spec.dQepsdbIIb2);
// --- Hyper3D::PoGrenobleSansPhaseAvecVar truc = PoGrenoble_et_var(toto,inv);
// utilitaire de vérif pour comparaison avec des calculs effectués avec un programme
// de calcul formel: par exemple maxima
// pour la vérification des dérivées / mu1, on utilise l'outil yacas: les informations sont dans les fichiers remarques_verif_yacas
// associées à plusieurs autres fichiers... ça marche super !!!
//------ débug -------
/*
if (indic_Verif_CompaFormel==1)
{ cout << "\n\n----------------------------------------------------------------";
int prec = ParaGlob::NbdigdoCA();
cout << "\n V= " << setprecision(prec) << inv.V << " Qeps= "
<< setprecision(prec) << i_s.Qeps << " cos3phi= "
<< setprecision(prec) << i_s.cos3phi;
cout << "\n Qs= " << setprecision(prec) << Qs
<< " dQs_dcos3phi= " << setprecision(prec) << dQs_dcos3phi
<< " dQs_dcos3phi2= " << setprecision(prec) << dQs_dcos3phi2
<< "\n Qe= " << setprecision(prec) << Qe
<< " dQe_dcos3phi= " << setprecision(prec) << dQe_dcos3phi
<< " dQe_dcos3phi2= " << setprecision(prec) << dQe_dcos3phi2
<< "\n mu1= " << setprecision(prec) << mu1
<< " dMu1_dcos3phi= " << setprecision(prec) << dMu1_dcos3phi
<< " dMu1_dcos3phi2= " << setprecision(prec) << dMu1_dcos3phi2
<< "\n mu2= " << setprecision(prec) << mu2
<< " dMu2_dcos3phi= " << setprecision(prec) << dMu2_dcos3phi
<< " dMu2_dcos3phi2= " << setprecision(prec) << dMu2_dcos3phi2
<< "\n mu3= " << setprecision(prec) << mu3
<< " dMu3_dcos3phi= " << setprecision(prec) << dMu3_dcos3phi
<< " dMu3_dcos3phi2= " << setprecision(prec) << dMu3_dcos3phi2;
cout << "\n\n Qec= " << setprecision(prec) << Qec << " dQec_Qs= " << setprecision(prec) << dQec_Qs
<< " dQec_QsQs= " << setprecision(prec) << dQec_QsQs
<< "\n dQec_Mu1= " << setprecision(prec) << dQec_Mu1 << " dQec_Mu1Mu1= " << setprecision(prec) << dQec_Mu1Mu1
<< " dQec_QsMu1= " << setprecision(prec) << dQec_QsMu1
<< "\n S1= " << setprecision(prec) << S1 << " dS1_Qs= " << setprecision(prec) << dS1_Qs
<< " dS1_QsQs= " << setprecision(prec) << dS1_QsQs
<< "\n dS1_Mu1= " << setprecision(prec) << dS1_Mu1 << " dS1_Mu1Mu1= " << setprecision(prec) << dS1_Mu1Mu1
<< " dS1_QsMu1= " << setprecision(prec) << dS1_QsMu1
<< "\n S2= " << setprecision(prec) << S2 << " dS2_Qs= " << setprecision(prec) << dS2_Qs
<< " dS2_QsQs= " << setprecision(prec) << dS2_QsQs
<< "\n dS2_Q= " << setprecision(prec) << dS2_Q << " dS2_QQ= " << setprecision(prec) << dS2_QQ
<< " dS2_QQs= " << setprecision(prec) << dS2_QQs
<< "\n dS2_Mu1= " << setprecision(prec) << dS2_Mu1 << " dS2_Mu1Mu1= " << setprecision(prec) << dS2_Mu1Mu1
<< " dS2_QsMu1= " << setprecision(prec) << dS2_QsMu1 << " dS2_QMu1= " << setprecision(prec) << dS2_QMu1
<< "\n S3= " << setprecision(prec) << S3 << " dS3_Qe= " << setprecision(prec) << dS3_Qe
<< " dS3_QeQe= " << setprecision(prec) << dS3_QeQe
<< "\n dS3_Q= " << setprecision(prec) << dS3_Q << " dS3_QQ= " << setprecision(prec) << dS3_QQ
<< " dS3_QQe= " << setprecision(prec) << dS3_QQe
<< "\n S4= " << S4 << " dS4_Qe= " << dS4_Qe << " dS4_QeQe= " << dS4_QeQe;
cout << "\n\n dE_Qs= " << setprecision(prec) << dE_Qs << " dE_Qe= " << setprecision(prec) << dE_Qe
<< "\n dE_Mu1= " << setprecision(prec) << dE_Mu1 << "\n dE_Mu2= " << setprecision(prec) << dE_Mu2
<< " dE_mu3= " << setprecision(prec) << dE_Mu3
<< "\n dE_QsQs= " << setprecision(prec) << dE_QsQs << " dE_QeQe= " << setprecision(prec) << dE_QeQe
<< " dE_Mu3Qe= " << setprecision(prec) << dE_Mu3Qe
<< "\n dE_QsMu1= " << setprecision(prec) << dE_QsMu1 << " dE_Mu1Mu1= " << setprecision(prec) << dE_Mu1Mu1;
cout << "\n dE_QQs= " << setprecision(prec) << dE_QQs << " dE_QQe= " << setprecision(prec) << dE_QQe
<< " dE_QMu2= " << setprecision(prec) << dE_QMu2 << " dE_QMu3= " << setprecision(prec) << dE_QMu3
<< "\n dE_QMu1= " << setprecision(prec) << dE_QMu1;
indic_Verif_CompaFormel++; // pour éviter une boucle infinie
};
if (indic_Verif_CompaFormel==0)
CompaFormel();
*/
// ------- fin débug --------
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 IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var
(const double & Qeps,const Invariant & inv
,const PoGrenobleSansPhaseAvecVar& potret )
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_PoGrenobleorgeas1_et_var == 0)
{ cout << "\n ****verification des derivées du potentiels par rapport aux invariants****";
cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var \n";
};
indic_Verif_PoGrenobleorgeas1_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 dérivees du potentiel";
cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var(.."
<< " , numero d'increment = " << indic_Verif_PoGrenobleorgeas1_et_var;
Sortie(1);
};
};
// vérif des dérivées du potentiels par rapport à l'invariant V, ceci par différences finies
void IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var_VV
(const double & Qeps,const Invariant & inv
,const PoGrenoble_VV& potret )
{ // les grandeurs de base, qui ici ne dépendent pas de la phase
double Qs = Q0s; double Qe = Q0e; double mu1 = mu01; double mu2 = mu02; double mu3 = mu03;
switch (cas_Q0s_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon fct laurent
{ if ((*temperature) >= Tc)
{Qs=0.5*Gr*( ((*temperature)-Tc+Ts)-sqrt(Sqr((*temperature)-Tc+Ts) + ar*ar)) + Qrmax;}
else
{Qs=0.5*Gr*( ((*temperature)-Tc-Ts)+sqrt(Sqr((*temperature)-Tc-Ts) + ar*ar));};
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qs = Q0s_temperature->Valeur(*temperature);
break;
}
};
switch (cas_Q0e_thermo_dependant)
{case 1:// cas d'une dépendance à la température selon Qe(T) = Qe(T0rQe) + (Qs(T)-Qs(T0rQe))/(3.* h1 * (mu3+mu2))
{ Qe = Qe_T0rQe + (Qs-Qs_T0rQe)/(3.* h1 * (mu3+mu2));
break;
}
case 2:// cas d'une dépendance à la température selon une fonction d'évolution
{ Qe = Q0e_temperature->Valeur(*temperature);
break;
}
};
// dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_PoGrenobleorgeas1_et_var_VV == 0)
{ cout << "\n ****verification des derivées du potentiels par rapport aux invariants****";
cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var \n";
};
indic_Verif_PoGrenobleorgeas1_et_var_VV++;
// 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); // recopie des invariants
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;
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);
// calcul des dérivées premières
double E_V = (E_et_dV - E_n)/delta;
// calcul des dérivées secondes
double E_V2 = (E_et_dV2 - 2.*E_et_dV + E_n ) /(delta * 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.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 (erreur)
{ cout << "\n erreur dans le calcul analytique des derivees du potentiel selon V et VV";
cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var_VV(.."
<< " , numero d'increment = " << indic_Verif_PoGrenobleorgeas1_et_var_VV;
Sortie(1);
};
};
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
void IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var
(const double & Qeps,const Invariant & inv,const double& cos3phi
,const PoGrenobleAvecPhaseAvecVar& potret )
{ // dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_PoGrenobleorgeas1_et_var == 0)
{ cout << "\n ****vérification des dérivées du potentiels par rapport aux invariants****";
cout << "\n IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var \n";
};
indic_Verif_PoGrenobleorgeas1_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 = 100.*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;cout << " " << potret.Ecos3phi2 << " " <<E_cos3phi_2; }
}
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 IsoHyper3DOrgeas1::Verif_PoGrenoble_et_var(.."
<< " , numero d'increment = " << indic_Verif_PoGrenobleorgeas1_et_var
<< " , numero d'erreur : " << erreur ;
// Sortie(1);
}
};
// utilitaire de vérif pour comparaison avec des calculs effectués avec un programme
// de calcul formel: par exemple maxima
void IsoHyper3DOrgeas1::CompaFormel()
{ indic_Verif_CompaFormel++; // pour éviter une boucle infinie
// l'objectif ici est d'obtenir des résultats numériques à partir de données entrées interactivement
bool fin_interactif = false;
// on définit des invariants initialisées à 0 par défaut
Invariant2QepsCosphi i_s;
Invariant inv;
while (!fin_interactif)
{ cout << "\n variation relative de volume (rep V ) ? "
<< "\n Qeps (rep Qeps) ? "
<< "\n cos(3phi) (rep cos) ? "
<< "\n affichage (rep a ) ? "
<< "\n arret verif (rep f ) ? ";
string rep;
// procédure de lecture avec prise en charge d'un retour chariot
rep = lect_return_defaut(false,"o");
if (rep == "V")
{ cout << "\n valeur de V "; inv.V=lect_double();}
else if (rep == "Qeps")
{ cout << "\n valeur de Qeps "; i_s.Qeps=lect_double();}
else if (rep == "cos")
{ cout << "\n valeur de cos3phi "; i_s.cos3phi=lect_double();}
else if (rep=="f")
{fin_interactif=true;}
else if (rep=="a")
{
// appel du calcul des différentes valeurs
PoGrenobleAvecPhaseAvecVar pot = PoGrenoblePhase_et_var(i_s,inv);
// affichage
cout << "\n\n----------------------------------------------------------------"
<< "\n V= " << inv.V << " Qeps= " << i_s.Qeps << " cos3phi= " << i_s.cos3phi << "\n\n";
cout << "\n E= "<< pot.E; // valeur du potentiel
cout << "\n EV= "<< pot.EV; // variation du potentiel par rapport à V
cout << "\n EQ= "<< pot.EQ; // variation du potentiel par rapport à Q
cout << "\n Ecos3phi= "<< pot.Ecos3phi; // variation du potentiel par rapport à cos3phi
cout << "\n EVV= "<< pot.EVV; // variation seconde par rapport à V
cout << "\n EQQ= "<< pot.EQQ; // variation seconde par rapport à Q
cout << "\n Ecos3phi2= "<< pot.Ecos3phi2; // variation seconde par rapport à cos3phi
cout << "\n EQV= "<< pot.EQV; // variation seconde par rapport à Q et V
cout << "\n EVcos3phi= "<< pot.EVcos3phi; // variation seconde par rapport à V cos3phi
cout << "\n EQcos3phi= "<< pot.EQcos3phi; // variation seconde par rapport à Q cos3phi
cout << "\n\n\n" << endl;
}
};
};
Reference in a new issue Copy permalink