rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/comportement/Hyper_elastique/Hyper3D.cc
Gérard Rio 3ccffa8899 V 7.017 
- intro para limite_inf_bIIb_ et limite_inf_Qeps_ dans les potentiels de grenoble
- cor bug diff finie pour Qeps petit sur loi grenoble
- modif prise en compte de la regularisation loi grenoble (passage systématique en différence finie pour Qeps très faible)
2023-07-06 08:42:57 +02:00

1386 lines
68 KiB
C++
Raw Permalink Blame History

// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
#include "Hyper3D.h"
//#include "ComLoi_comp_abstraite.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
//================== initialisation des variables static ======================
// ---------- classe Hyper3D ---------
Hyper3D::Hyper3D () : // Constructeur par defaut
HyperD(RIEN_COMP,RIEN_CATEGORIE_LOI_COMP,3,false)
,indic_Verif_Potentiel_et_var(0),limite_inf_Qeps(sqrt(2)*6.e-5)
,limite_inf_bIIb(36.e-10)
{};
// Constructeur utile si l'identificateur du nom de la loi
// de comportement et le paramètre phase sont connus
Hyper3D::Hyper3D (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp,bool avec_ph) :
HyperD(id_compor,categorie_comp,3,avec_ph)
,indic_Verif_Potentiel_et_var(0),limite_inf_Qeps(sqrt(2)*6.e-5)
,limite_inf_bIIb(36.e-10)
{} ;
// Constructeur utile si l'identificateur du nom de la loi
// de comportement et le paramètre phase sont connus
Hyper3D::Hyper3D (char* nom,Enum_categorie_loi_comp categorie_comp,bool avec_ph) :
HyperD(nom,categorie_comp,3,avec_ph)
,indic_Verif_Potentiel_et_var(0),limite_inf_Qeps(sqrt(2)*6.e-5)
,limite_inf_bIIb(36.e-10)
{} ;
// Constructeur de copie
Hyper3D::Hyper3D (const Hyper3D& loi) :
HyperD (loi)
,indic_Verif_Potentiel_et_var(loi.indic_Verif_Potentiel_et_var)
,limite_inf_Qeps(loi.limite_inf_Qeps)
,limite_inf_bIIb(loi.limite_inf_bIIb)
{};
// =========== METHODES Protégées dérivant de virtuelles : ==============
// Calcul des invariants, de epsBH,
// retour de IdGBH qui pointe sur le bon tenseur
//++ TensBH * Invariants (TenseurBB& epsBB_t,TenseurBB& gijBB_t,
TenseurBH * Hyper3D::Invariants (const TenseurBB& epsBB_t,const TenseurBB& ,
const TenseurHH & gijHH_t,const double& jacobien_0,const double& jacobien_t,
//++ Invariant & invariant,TensBH & epsBH);
Invariant & inv,TenseurBH & epsBH_t)
{
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_t); // passage en dim 3
Tenseur3BH & epsBH = *((Tenseur3BH*) &epsBH_t); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_t); // passage en dim 3
epsBH = epsBB * gijHH; // deformation en mixte
inv.Ieps = epsBH.Trace(); // trace de la déformation
Tenseur3BH eps_barreBH = epsBH - (inv.Ieps/3.) * IdBH3;
inv.V = (jacobien_t/jacobien_0);
// dans le cas de la prise en compte de la dilatation il faut changer V
if (dilatation)
// on considère que la variation de volume thermique = la trace du tenseur de def thermique
{double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace();
inv.V -= var_vol_thermique;
};
double inv_V = inv.V; // init pour être éventuellement modifié
// s'il n'y a pas de régularisation : fact_regularisation = 0.
if (Dabs(inv_V) < Dabs(fact_regularisation))
inv_V += fact_regularisation;
inv.bIIb = eps_barreBH.II() / 2.; // second invariant
// la première méthode de calcul est la méthode traditionnelle, mais elle nécessite beaucoup
// plus de calcul que la seconde (qui donne le même résultat à 10^{-18} près environ
// une fois que V est calculé
// inv.bIIIb = eps_barreBH.III() / 3.; // troisième invariant
inv.bIIIb = -0.125/(inv_V * inv_V ) +0.125 - 0.25*inv.Ieps
+ 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
- 0.5*inv.bIIb +1./3. * inv.bIIb * inv.Ieps; // troisième invariant
// if (avec_regularisation)
// {// inv.V += fact_regularisation;
// inv.bIIb += fact_regularisation;
// inv.bIIIb += fact_regularisation;
// };
return & IdBH3;
};
// calcul des invariants et de leurs variations, de epsBH,
// et de sa variation, puis retour de IdGBH qui pointe sur le bon tenseur
//++ TensBH * Invariants_et_var (TenseurBB& epsBB_tdt,
TenseurBH * Hyper3D::Invariants_et_var (const TenseurBB& epsBB_tdt,
const TenseurBB& ,const Tableau <TenseurBB *>& d_gijBB_tdt,
const TenseurHH & gijHH_tdt,const Tableau <TenseurHH *>& d_gijHH_tdt,
const double& jacobien_0,const double& jacobien_tdt,const Vecteur& d_jacobien_tdt,
InvariantVarDdl& inv,
//++ TensBH & epsBH_tdt,Tableau<TensBH> & depsBH_tdt);
TenseurBH & epsBH_tdt,Tableau<TenseurBH*> & depsBH_tdt)
{
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
Tenseur3BH & epsBH = *((Tenseur3BH*) &epsBH_tdt); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // passage en dim 3
epsBH = epsBB * gijHH; // deformation en mixte
inv.Ieps = epsBH.Trace(); // trace de la déformation
Tenseur3BH eps_barreBH = epsBH - (inv.Ieps/3.) * IdBH3;
inv.V = (jacobien_tdt/jacobien_0);
// dans le cas de la prise en compte de la dilatation il faut changer V
if (dilatation)
// on considère que la variation de volume thermique = la trace du tenseur de def thermique
{double var_vol_thermique = (gijHH_tdt * (*epsBB_therm)).Trace();
inv.V -= var_vol_thermique;
};
double inv_V = inv.V; // init pour être éventuellement modifié
// s'il n'y a pas de régularisation : fact_regularisation = 0.
if (Dabs(inv_V) < Dabs(fact_regularisation))
inv_V += fact_regularisation;
inv.bIIb = eps_barreBH.II() / 2.; // second invariant
// la première méthode de calcul est la méthode traditionnelle, mais elle nécessite beaucoup
// plus de calcul que la seconde (qui donne le même résultat à 10^{-18} près environ
// une fois que V est calculé
// inv.bIIIb = eps_barreBH.III() / 3.; // troisième invariant
// cout << "\n inv.bIIIb= "<< inv.bIIIb;
inv.bIIIb = -0.125/(inv_V * inv_V ) +0.125 - 0.25*inv.Ieps
+ 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
- 0.5*inv.bIIb +1./3. * inv.bIIb * inv.Ieps; // troisième invariant
// cout << " , "<< inv.bIIIb << " debug Hyper3D::Invariants_et_var " << endl;
// if (avec_regularisation)
// { //inv.V += fact_regularisation;
// inv.bIIb += fact_regularisation;
// inv.bIIIb += fact_regularisation;
// };
// --- debug
/// pour info la bonne expression de bIIIb est donnée par tutu !!
// mais le résultat est identique (au erreur de calcul près) avec titi (c'est d'après la formule de denis
// thèse chapitre, II.2.3, formule 2.19
// et idem avec toto1 , en fait il n'y a que toto2 qui est complètement faux
//
// double toto2 = inv.bIIIb - (-1./(8.*inv.V * inv.V ) +1./8. - 0.25*inv.Ieps
// + 1./3. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps*11/54.
// -0.5*inv.bIIb + inv.bIIb * inv.Ieps);
//
// double je = 1.-2./3.* inv.Ieps;
// double je3=je*je*je;
// double titi = 1./(inv.V*inv.V) - (je3 - 4.*inv.bIIb*je-8.*inv.bIIIb);
// double tutu = 1./(8.*inv.V*inv.V) - 1./8.+1./4.*inv.Ieps-1./6.*inv.Ieps*inv.Ieps
// +1./27.*inv.Ieps*inv.Ieps*inv.Ieps
// + 0.5*inv.bIIb*1.-1./3.*inv.bIIb* inv.Ieps+inv.bIIIb;
// double toto1 = inv.bIIIb -( -1./(8.*inv.V * inv.V ) +1./8. - 0.25*inv.Ieps
// + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
// - 0.5*inv.bIIb +1./3. * inv.bIIb * inv.Ieps);
//
// if ((Abs(toto1)>ConstMath::petit) || (Abs(toto2)>ConstMath::petit)|| (Abs(titi)>ConstMath::petit))
// cout << "\n ** toto1= " << toto1 << " toto2= " << toto2 << " titi= " << titi<< " tutu= " << tutu;
//
//// -- fin debug
//
// cas le la variation
int nbddl = d_gijBB_tdt.Taille();
Tenseur3BH deps_barreBH; // grandeurs intermédiaires
for (int i = 1; i<= nbddl; i++)
{const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
const Tenseur3BB & dgijBB = *((Tenseur3BB*)(d_gijBB_tdt(i))); // passage en dim 3
Tenseur3BH & depsBH = *((Tenseur3BH*)(depsBH_tdt(i))); // passage en dim 3
depsBH = 0.5 * dgijBB * gijHH_tdt + epsBB * dgijHH;
double dIeps = (depsBH).Trace();inv.dIeps(i) = dIeps;
deps_barreBH = depsBH - (dIeps/3.) * IdBH3; // en mixte les coordonnées de IdBH3 sont fixes
double dbIIb = inv.dbIIb(i) = deps_barreBH && eps_barreBH;
double dV = inv.dV(i) = (d_jacobien_tdt(i)/jacobien_0) ;
// Pour minimiser le temps de calcul on calcul la dérivée à l'aide de la formule scalaire
// inv.dbIIIb(i) = (deps_barreBH * eps_barreBH) && eps_barreBH;
// -- formule scalaire de base
// inv.bIIIb = -0.125/(inv.V * inv.V ) +0.125 - 0.25*inv.Ieps
// + 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
// - 0.5*inv.bIIb +1./3. * inv.bIIb * inv.Ieps; // troisième invariant
// d'où la variation
inv.dbIIIb(i) = 0.25*dV/(inv_V * inv_V * inv_V) - 0.25*dIeps
+ 1./3. * dIeps * inv.Ieps - dIeps*inv.Ieps*inv.Ieps/9.
- 0.5*dbIIb +1./3. * ( dbIIb * inv.Ieps + inv.bIIb * dIeps ) ; // troisième invariant
};
// --- debug
// affichage des invariants pour comparer
// cout << "\n Hyper3D::Invariants_et_var ( V:" << inv.V <<" Ieps:" << inv.Ieps
// << " bIIb:" << inv.bIIb << " bIIIb:" << inv.bIIIb
// << "\n dV:" << inv.dV <<"\n dIeps:" << inv.dIeps
// << "\n dbIIb:" << inv.dbIIb << "\n dbIIIb:" << inv.dbIIIb
// << endl;
//
//// -- fin debug
return & IdBH3;
};
// calcul des invariants et de leurs variations par rapport aux déformations
// et de sa variation, puis retour de IdGBH qui pointe sur le bon tenseur
TenseurBH * Hyper3D::Invariants_et_varEps (const TenseurBB& epsBB_tdt,
const TenseurBB& ,const TenseurHH & gijHH_tdt,
const double& jacobien_0,const double& jacobien_tdt,
InvariantVarEps& inv,TenseurBH & epsBH_tdt)
{
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
Tenseur3BH & epsBH = *((Tenseur3BH*) &epsBH_tdt); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // passage en dim 3
epsBH = epsBB * gijHH; // deformation en mixte
inv.Ieps = epsBH.Trace(); // trace de la déformation
Tenseur3BH eps_barreBH = epsBH - (inv.Ieps/3.) * IdBH3;
inv.V = (jacobien_tdt/jacobien_0);
// dans le cas de la prise en compte de la dilatation il faut changer V
if (dilatation)
// on considère que la variation de volume thermique = la trace du tenseur de def thermique
{double var_vol_thermique = (gijHH * (*epsBB_therm)).Trace();
inv.V -= var_vol_thermique;
};
double inv_V = inv.V; // init pour être éventuellement modifié
// s'il n'y a pas de régularisation : fact_regularisation = 0.
if (Dabs(inv_V) < Dabs(fact_regularisation))
inv_V += fact_regularisation;
inv.bIIb = eps_barreBH.II() / 2.; // second invariant
// la première méthode de calcul est la méthode traditionnelle, mais elle nécessite beaucoup
// plus de calcul que la seconde (qui donne le même résultat à 10^{-18} près environ
// une fois que V est calculé
inv.bIIIb = eps_barreBH.III() / 3.; // troisième invariant
// cout << "\n inv.bIIIb= "<< inv.bIIIb;
inv.bIIIb = -0.125/(inv_V * inv_V ) +0.125 - 0.25*inv.Ieps
+ 1./6. * inv.Ieps * inv.Ieps - inv.Ieps*inv.Ieps*inv.Ieps/27.
- 0.5*inv.bIIb +1./3. * inv.bIIb * inv.Ieps; // troisième invariant
// cout << " , "<< inv.bIIIb << " debug Hyper3D::Invariants_et_varEps " << endl;
// if (avec_regularisation)
// { //inv.V += fact_regularisation;
// inv.bIIb += fact_regularisation;
// inv.bIIIb += fact_regularisation;
// };
// cas le la variation
Tenseur3HH epsHH = gijHH_tdt * epsBH;
inv.dIeps_deps_HH = gijHH_tdt - 2. * epsHH;
// on a : bIIb = bII - 1/6 Ieps^2 ; et dbII = epsHH - 2. * epsHH * epsBH
Tenseur3HH d_bII_HH = epsHH - 2. * epsHH * epsBH ;
inv.dbIIb_deps_HH = d_bII_HH - (inv.Ieps/3.) * inv.dbIIb_deps_HH;
// on a = bIIIb = bIII - 2/3 (I * bII - I^3/9) : thèse denis, page II.2.1
Tenseur3HH d_bIII = epsHH * epsBH - 2. * epsHH * epsBH * epsBH;
inv.dbIIIb_deps_HH = d_bIII + ((inv.Ieps * inv.Ieps / 3.)- (2.*inv.bIIb/3.) ) // mise en facteur du coeff
* inv.dIeps_deps_HH // scalaire
- (2.*inv.Ieps/3.) * inv.dbIIb_deps_HH;
// formule identique au cas de Mooney Rivlin , voir document de synthèse sur les lois hyper-élastiques chap 9
inv.dV_deps_HH = inv.V * gijHH_tdt;
// on libère les tenseurs
LibereTenseur();LibereTenseurQ();
return & IdBH3;
};
// calcul du potentiel et de ses dérivées non compris la phase
Hyper3D::PotenSansPhaseSansVar Hyper3D::Potentiel(const Invariant & inv,const double& jacobien0)
{ // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
// mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
// directement une technique de différence finie pour obtenir une dérivée numérique
// def du retour
PotenSansPhaseSansVar potret;
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
// if ((Dabs(inv.bIIb) > limite_inf_bIIb)||(avec_regularisation))
if (Dabs(inv.bIIb) > limite_inf_bIIb)
{ // calcul du nouvel invariant sous la nouvelle forme
InvariantQeps inv_spes = InvariantSpecifSansPhase(inv);
// calcul du potentiel sous la forme polaire
PoGrenobleSansPhaseSansVar pot = PoGrenoble(inv_spes,inv);
// calcul de potret
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
potret.EbIIb = jacobien0 * inv.V * pot.EQ * inv_spes.dQepsdbIIb;
// stockage éventuel des invariants et du potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= inv_spes.Qeps;
save_resulHyperD.invP->cos3phi= 0.;
save_resulHyperD.invP->potentiel= pot.E;
};
}
else
// cas d'un calcul par différence finie
{ // calcul des invariants sous la nouvelle forme
double Qeps = Invariant0Qeps(inv);
// calcul du potentiel sous la forme polaire au point actuel
// ainsi que la variation suivant V, la phase ne peut pas être prise en compte
// car Qeps est nulle
PoGrenoble_V pot = PoGrenoble_et_V(Qeps,inv);
// calcul du potentiel pour une petite variation suivant bIIb
Invariant inv_n = inv; inv_n.bIIb += limite_inf_bIIb;
double Qeps_n = Invariant0Qeps(inv_n);
double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
// stockage éventuel des invariants et potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv_n.V;
save_resulHyperD.invP->Qeps= Qeps_n;
save_resulHyperD.invP->cos3phi= 0.;
save_resulHyperD.invP->potentiel= pot.E;
};
double E_n = PoGrenoble(Qeps_n,inv_n);
// calcul de potret
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
{double inter = (E_n - pot.E );
// on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
if (inter < ConstMath::trespetit)
{ potret.EbIIb = 0.; }
else
{potret.EbIIb = jacobien0 * inv.V * inter * unsurdelta;};
};
}
return potret;
};
// calcul du potentiel et de ses dérivées avec la phase
Hyper3D::PotenAvecPhaseSansVar Hyper3D::PotentielPhase(const Invariant& inv,const double& jacobien0)
{ // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
// mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
// directement une technique de différence finie pour obtenir une dérivée numérique
// def du retour
PotenAvecPhaseSansVar potret;
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
// if ((Dabs(inv.bIIb) > limite_inf_bIIb)||(avec_regularisation))
if (Dabs(inv.bIIb) > limite_inf_bIIb)
{ // calcul des invariants sous la nouvelle forme
InvariantQepsCosphi inv_spes = InvariantSpecif(inv);
// calcul du potentiel sous la forme polaire
PoGrenobleAvecPhaseSansVar pot = PoGrenoblePhase(inv_spes,inv);
// stockage éventuel des invariants et potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= inv_spes.Qeps;
save_resulHyperD.invP->cos3phi= inv_spes.cos3phi;
save_resulHyperD.invP->potentiel= pot.E;
};
// calcul de potret
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * ( pot.EV + pot.Ecos3phi * inv_spes.dcos3phidV) + pot.E);
potret.EIeps = jacobien0 * inv.V * (pot.Ecos3phi * inv_spes.dcos3phidIeps);
potret.EbIIb = jacobien0 * inv.V * (pot.EQ * inv_spes.dQepsdbIIb
+ pot.Ecos3phi * inv_spes.dcos3phidbIIb);
}
else
// cas d'un calcul par différence finie, on considère que toute les dérivées suivant la phase
// sont nulles dans ce cas
{ // calcul des invariants sous la nouvelle forme
double Qeps = Invariant0Qeps(inv);
// calcul du potentiel sous la forme polaire au point actuel
// ainsi que la variation suivant V
PoGrenoble_V pot = PoGrenoble_et_V(Qeps,inv);
// calcul du potentiel pour une petite variation suivant bIIb
Invariant inv_n = inv; inv_n.bIIb += limite_inf_bIIb;
Invariant0QepsCosphi inv_spes_n = Invariant0Specif(inv_n);
double E_n = PoGrenoble(inv_spes_n,inv_n);
double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
// stockage éventuel des invariants et du potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= inv_spes_n.Qeps;
save_resulHyperD.invP->cos3phi= inv_spes_n.cos3phi;
save_resulHyperD.invP->potentiel= pot.E;
};
// calcul de potret
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
potret.EIeps = 0.;
// potret.EbIIb = jacobien0 * inv.V * (-pot.E + E_n)/limite_inf_bIIb;
{double inter = (E_n-pot.E );
// on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
if (inter < ConstMath::trespetit)
{ potret.EbIIb = 0.; }
else
{potret.EbIIb = jacobien0 * inv.V * inter *unsurdelta;};
};
}
return potret;
};
// calcul du potentiel sans phase et dérivées avec ses variations par rapport aux invariants
Hyper3D::PotenSansPhaseAvecVar Hyper3D::Potentiel_et_var(const Invariant& inv,const double& jacobien0 )
{ // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
// mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
// directement une technique de différence finie pour obtenir une dérivée numérique
// def du retour
PotenSansPhaseAvecVar potret;
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
// if ((Dabs(inv.bIIb) > ConstMath::petit)||(avec_regularisation))
if (Dabs(inv.bIIb) > limite_inf_bIIb)
{ // calcul des invariants sous la nouvelle forme
Invariant2Qeps inv_spes = Invariant2SpecifSansPhase(inv);
// calcul du potentiel sous la forme polaire
PoGrenobleSansPhaseAvecVar pot = PoGrenoble_et_var(inv_spes,inv);
//debug
//cout << "\n debug: Hyper3D::"
// << " E="<<pot.E << " pot.EV=" << pot.EV << " pot.EVV="<<pot.EVV<<endl;
//if (pot.E < 0)
// cout << "\n pb ";
//fin debug
// calcul de potret
//dérivées premières
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
potret.EbIIb = jacobien0 * inv.V * pot.EQ * inv_spes.dQepsdbIIb;
// dérivées secondes
potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
potret.EbIIb2 = jacobien0 * inv.V * ( pot.EQQ * inv_spes.dQepsdbIIb * inv_spes.dQepsdbIIb
+ pot.EQ * inv_spes.dQepsdbIIb2 );
potret.EbIIbV = jacobien0 * ((pot.EQ + inv.V * pot.EQV) * inv_spes.dQepsdbIIb);
// stockage éventuel des invariants et du potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= inv_spes.Qeps;
save_resulHyperD.invP->cos3phi= 0.;
save_resulHyperD.invP->potentiel= pot.E;
};
// appel d'une vérification éventuelle
//+++verif*** Verif_Potentiel_et_var(inv,jacobien0,potret);
}
else
// cas d'un calcul par différence finie, on considère que toute les dérivées suivant la phase
// sont nulles dans ce cas
{ // calcul des invariants sous la nouvelle forme
Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb);
double Qeps = Invariant0Qeps(inv_n0);
// calcul du potentiel sous la forme polaire au point actuel
// ainsi que la variation première et seconde suivant V
PoGrenoble_VV pot = PoGrenoble_et_VV(Qeps,inv_n0);
// calcul du potentiel pour une petite variation suivant bIIb
Invariant inv_n = inv_n0;
inv_n.bIIb += limite_inf_bIIb;
double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
double Qeps_n = Invariant0Qeps(inv_n);
// stockage éventuel des invariants et du potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= Qeps_n;
save_resulHyperD.invP->cos3phi= 0.;
save_resulHyperD.invP->potentiel= pot.E;
};
PoGrenoble_V pot_V = PoGrenoble_et_V(Qeps_n,inv_n);
// et pour la dérivée seconde de nouveau une petite variation
inv_n.bIIb += limite_inf_bIIb;
Qeps_n = Invariant0Qeps(inv_n);
double E_n = PoGrenoble(Qeps_n,inv_n);
// calcul de potret
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
{double inter = (pot_V.E - pot.E);
// on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
if (inter < ConstMath::trespetit)
{ potret.EbIIb = 0.; }
else
{potret.EbIIb = jacobien0 * inv.V * inter *unsurdelta;};
};
// potret.EbIIb = jacobien0 * inv.V * (pot_V.E - pot.E)/limite_inf_bIIb;
// dérivées secondes
potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
{double inter = ((E_n + pot.E) - 2.*pot_V.E );
// on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// la dérivée seconde du potentiel doit-être elle aussi > 0 (au moins à l'origine)
if (inter < ConstMath::trespetit)
{ potret.EbIIb2 = 0.; }
else
{potret.EbIIb2 = jacobien0 * inv.V * inter
* unsurdelta * unsurdelta;};
//cout << "\n debug ((E_n + pot.E) - 2.*pot_V.E )= " << ((E_n + pot.E) - 2.*pot_V.E ) << endl;
};
// pour la variation de la dérivée / V, a priori on peut avoir les deux signes
// donc on ne fait pas de filtrage
// {double inter = (pot_V.EV - pot.EV);
// // on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// if (inter < ConstMath::trespetit)
// { potret.EbIIbV = potret.EbIIb / inv.V; }
// else
// {potret.EbIIbV = potret.EbIIb / inv.V
// + jacobien0 * inv.V * inter * unsurdelta;};
// };
potret.EbIIbV = potret.EbIIb / inv.V
+ jacobien0 * inv.V * (pot_V.EV - pot.EV) * unsurdelta;
// potret.EbIIbV = potret.EbIIb / inv.V
// + jacobien0 * inv.V * (pot_V.EV - pot.EV) / limite_inf_bIIb;
// appel d'une vérification éventuelle
//+++verif*** Verif_Potentiel_et_var(inv,jacobien0,potret);
}
// retour
return potret;
};
// calcul du potentiel avec phase et dérivées avec ses variations par rapport aux invariants
Hyper3D::PotenAvecPhaseAvecVar Hyper3D::PotentielPhase_et_var(const Invariant& inv,const double& jacobien0)
{ // en fait il s'agit du transfert du calcul du potentiel de la forme grenobloise à la forme standard
// mais ce transfert n'est possible que si bIIb est non nulle, dans le cas contraire on utilise `
// directement une technique de différence finie pour obtenir une dérivée numérique
// def du retour
PotenAvecPhaseAvecVar potret; // mise à 0 de tous les termes
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyperD & save_resulHyperD = *((SaveResulHyperD*) saveResul);
// cout << "\n " << potret.E; // pour avoir a cces en debug ??
// if ((Dabs(inv.bIIb) > limite_inf_bIIb)||avec_regularisation)
if (Dabs(inv.bIIb) > limite_inf_bIIb)
// if (Dabs(inv.bIIb) > limite_inf_Qeps*limite_inf_Qeps)
{ // calcul des invariants sous la nouvelle forme
Invariant2QepsCosphi inv_spes = Invariant2Specif(inv);
// calcul du potentiel sous la forme polaire
PoGrenobleAvecPhaseAvecVar pot = PoGrenoblePhase_et_var(inv_spes,inv);
// stockage éventuel des invariants et du potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= inv_spes.Qeps;
save_resulHyperD.invP->cos3phi= inv_spes.cos3phi;
save_resulHyperD.invP->potentiel= pot.E;
};
//cout << "\n calcul normal";
//debug
//cout << "\n debug: Hyper3D::"
// << " E="<<pot.E << " pot.EV=" << pot.EV << " pot.EVV="<<pot.EVV<<endl;
//if (pot.E < 0)
// cout << "\n pb ";
//fin debug
// calcul de potret
//dérivées premières
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
potret.EV = jacobien0 * (inv.V * ( pot.EV + pot.Ecos3phi * inv_spes.dcos3phidV) + pot.E);
potret.EIeps = jacobien0 * inv.V * (pot.Ecos3phi * inv_spes.dcos3phidIeps);
potret.EbIIb = jacobien0 * inv.V * (pot.EQ * inv_spes.dQepsdbIIb
+ pot.Ecos3phi * inv_spes.dcos3phidbIIb);
// dérivées secondes
double dcos3phidV_2 = inv_spes.dcos3phidV * inv_spes.dcos3phidV; // simplification
double dcos3phidIeps_2 = inv_spes.dcos3phidIeps * inv_spes.dcos3phidIeps; // "
double dcos3phidbIIb_2 = inv_spes.dcos3phidbIIb * inv_spes.dcos3phidbIIb; // "
potret.EVV = jacobien0 * (2. * ( pot.Ecos3phi + inv.V * pot.EVcos3phi) * inv_spes.dcos3phidV
+ 2. * pot.EV + inv.V * pot.EVV + inv.V * pot.Ecos3phi2 * dcos3phidV_2
+ inv.V * pot.Ecos3phi * inv_spes.dcos3phidV2);
potret.EIeps2 = jacobien0 * inv.V * ( pot.Ecos3phi2 * dcos3phidIeps_2 + pot.Ecos3phi * inv_spes.dcos3phidIeps2);
potret.EbIIb2 = jacobien0 * inv.V * ( pot.EQQ * inv_spes.dQepsdbIIb * inv_spes.dQepsdbIIb
+ 2. * pot.EQcos3phi * inv_spes.dQepsdbIIb * inv_spes.dcos3phidbIIb
+ pot.Ecos3phi2 * dcos3phidbIIb_2
+ pot.EQ * inv_spes.dQepsdbIIb2 + pot.Ecos3phi * inv_spes.dcos3phidbIIb2);
potret.EIepsV = jacobien0 * (( pot.Ecos3phi + inv.V * pot.EVcos3phi) * inv_spes.dcos3phidIeps
+ inv.V * pot.Ecos3phi2 * inv_spes.dcos3phidV * inv_spes.dcos3phidIeps);
potret.EbIIbV = jacobien0 * ((pot.EQ + inv.V * pot.EQV) * inv_spes.dQepsdbIIb
+ (pot.Ecos3phi + inv.V * pot.EVcos3phi ) * inv_spes.dcos3phidbIIb
+ inv.V * pot.EQcos3phi * inv_spes.dcos3phidV * inv_spes.dQepsdbIIb
+ inv.V * pot.Ecos3phi2 * inv_spes.dcos3phidV * inv_spes.dcos3phidbIIb
+ inv.V * pot.Ecos3phi * inv_spes.dcos3phidVdbIIb);
potret.EbIIbIeps = jacobien0 * inv.V * ( pot.EQcos3phi * inv_spes.dcos3phidIeps * inv_spes.dQepsdbIIb
+ pot.Ecos3phi2 * inv_spes.dcos3phidIeps * inv_spes.dcos3phidbIIb
+ pot.Ecos3phi * inv_spes.dcos3phidIepsdbIIb);
// PotenAvecPhaseAvecVar toto = potret; // pour le debug
// cout << "\n " << toto.E << " " << toto.EV;
// cout << "\n cas analytique ";
// --- pour le débug
// switch(fpclassify(toto.EV))
// { case FP_NAN: case FP_INFINITE:
// { // cas d'un nombre infini
// cout << "\n erreur valeur nan pour EV! ";
// cout << " poret= " << potret.E; // pour que potret existe pour le débug
// };
// break;
// };
// --- fin débug
// appel d'une vérification éventuelle
// Verif_Potentiel_et_var(inv,jacobien0,potret);
}
// ---------------------------------------------------------------
// non, le calcul par différence finies vraiment ne marche pas !!! c'est la merde
// car le calcul normale est extrèmement sensible aux limites initiales limite_inf_Qeps et limite_inf_bIIb
// la ça marche pour le cu mais peut-être pour un autre matériau ça ne marchera pas !!!
// peut-être mettre en paramètre, mais c'est vraiment difficile à utiliser
// bref à suivre ...
/*
else if (Dabs(inv.bIIb) > limite_inf_bIIb /10.)
{ // calcul des dérivées par différences finies
cout << "\n cas difference fini avec phase ";
Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb); // recopie des invariants sans les varddl
Invariant0QepsCosphi is_sp0 = Invariant0Specif(inv_n0);
double E_n = PoGrenoble(is_sp0,inv_n0); // la valeur du potentiel de référence
// ici on est avec phase donc trois invariants indépendant V et Qeps et cos3phi,
// 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/100.;
Invariant inv_et_dbIIb = inv_n0;inv_et_dbIIb.bIIb += delta;
Invariant0QepsCosphi is_et_dbIIb = InvariantDe_V_bIIb_Ieps(inv_et_dbIIb);
double E_et_dbIIb = PoGrenoble(is_et_dbIIb,inv_et_dbIIb);
Invariant inv_et_dbIIb2 = inv_n0;inv_et_dbIIb2.bIIb += 2. * delta;
Invariant0QepsCosphi is_et_dbIIb2 = InvariantDe_V_bIIb_Ieps(inv_et_dbIIb2);
double E_et_dbIIb2 = PoGrenoble(is_et_dbIIb2,inv_et_dbIIb2);
Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
Invariant0QepsCosphi is_et_dV = InvariantDe_V_bIIb_Ieps(inv_et_dV);
double E_et_dV = PoGrenoble(is_et_dV,inv_et_dV);
Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
Invariant0QepsCosphi is_et_dV2 = InvariantDe_V_bIIb_Ieps(inv_et_dV2);
double E_et_dV2 = PoGrenoble(is_et_dV2,inv_et_dV2);
Invariant inv_et_dVdbIIb = inv_n0;
inv_et_dVdbIIb.V += delta;
inv_et_dVdbIIb.bIIb += delta;
Invariant0QepsCosphi is_et_dVdbIIb = InvariantDe_V_bIIb_Ieps(inv_et_dVdbIIb);
double E_et_dVdbIIb = PoGrenoble(is_et_dVdbIIb,inv_et_dVdbIIb);
// cas avec Ieps
Invariant inv_et_dIeps = inv_n0;inv_et_dIeps.Ieps += delta;
Invariant0QepsCosphi is_et_dIeps = InvariantDe_V_bIIb_Ieps(inv_et_dIeps);
double E_et_dIeps = PoGrenoble(is_et_dIeps,inv_et_dIeps);
Invariant inv_et_dIeps2 = inv_n0;inv_et_dIeps2.Ieps += 2. * delta;
Invariant0QepsCosphi is_et_dIeps2 = InvariantDe_V_bIIb_Ieps(inv_et_dIeps2);
double E_et_dIeps2 = PoGrenoble(is_et_dIeps2,inv_et_dIeps2);
Invariant inv_et_dIepsdbIIb = inv_n0;
inv_et_dIepsdbIIb.Ieps += delta;
inv_et_dIepsdbIIb.bIIb += delta;
Invariant0QepsCosphi is_et_dIepsdbIIb = InvariantDe_V_bIIb_Ieps(inv_et_dIepsdbIIb);
double E_et_dIepsdbIIb = PoGrenoble(is_et_dIepsdbIIb,inv_et_dIepsdbIIb);
Invariant inv_et_dIepsdV = inv_n0;
inv_et_dIepsdV.Ieps += delta;
inv_et_dIepsdV.V += delta;
Invariant0QepsCosphi is_et_dIepsdV = InvariantDe_V_bIIb_Ieps(inv_et_dIepsdV);
double E_et_dIepsdV = PoGrenoble(is_et_dIepsdV,inv_et_dIepsdV);
// calcul des dérivées premières
potret.EV = jacobien0 *(E_et_dV * inv_et_dV.V - E_n * inv.V)/delta;
potret.EbIIb = (jacobien0 * inv.V)*(E_et_dbIIb - E_n)/delta;
potret.EIeps = (jacobien0 * inv.V)*(E_et_dIeps - E_n )/delta;
// calcul des dérivées secondes hors phases
potret.EVV = jacobien0 *(E_et_dV2 * inv_et_dV2.V - 2.*E_et_dV * inv_et_dV.V + E_n * inv.V)
/(delta * delta);
potret.EbIIb2 = (jacobien0 * inv.V)*(E_et_dbIIb2 - 2.*E_et_dbIIb + E_n) /(delta * delta);
potret.EIeps2 = (jacobien0 * inv.V) *(E_et_dIeps2 - 2.*E_et_dIeps + E_n )/(delta * delta);
double E_V_a_dbIIb = jacobien0 * (E_et_dVdbIIb * inv_et_dVdbIIb.V - E_et_dbIIb * inv_et_dbIIb.V) /delta;
potret.EbIIbV = ( E_V_a_dbIIb - potret.EV)/delta;
double E_Ieps_a_dbIIb = (jacobien0 * inv.V) * (E_et_dIepsdbIIb - E_et_dbIIb) /delta;
potret.EbIIbIeps = ( E_Ieps_a_dbIIb - potret.EIeps)/delta;
double E_Ieps_a_dV = jacobien0 * (E_et_dIepsdV * inv_et_dIepsdV.V - E_et_dV * inv.V) /delta;
potret.EIepsV = ( E_Ieps_a_dV - potret.EIeps)/delta;
}*/
/* else if (Dabs(inv.bIIb) > ConstMath::petit)
{ // ici on calcul comme si on n'avait pas de phase
// calcul des invariants sous la nouvelle forme
Invariant2Qeps inv_spes = Invariant2SpecifSansPhase(inv);
// calcul du potentiel sous la forme polaire
PoGrenobleSansPhaseAvecVar pot = PoGrenoble_et_var(inv_spes,inv);
// cout << "\n cas analytique sans phase ";
// calcul de potret
//dérivées premières
potret.E = pot.E * inv.V * jacobien0;
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
potret.EbIIb = jacobien0 * inv.V * pot.EQ * inv_spes.dQepsdbIIb;
// dérivées secondes
potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
potret.EbIIb2 = jacobien0 * inv.V * ( pot.EQQ * inv_spes.dQepsdbIIb * inv_spes.dQepsdbIIb
+ pot.EQ * inv_spes.dQepsdbIIb2 );
potret.EbIIbV = jacobien0 * ((pot.EQ + inv.V * pot.EQV) * inv_spes.dQepsdbIIb);
// appel d'une vérification éventuelle
//+++verif*** Verif_Potentiel_et_var(inv,jacobien0,potret);
}*/
else // dernier cas avec un inv.bIIb vraiment petit
// cas d'un calcul par différence finie, on considère que toute les dérivées suivant la phase
// sont nulles dans ce cas, car de toutes manière on ne peut pas calculer la phase
{ // calcul des invariants sous la nouvelle forme
// cout << "\n cas difference fini sans phase ";
Invariant inv_n0(inv.Ieps,inv.V,inv.bIIb,inv.bIIIb);
double Qeps = Invariant0Qeps(inv_n0);
// calcul du potentiel sous la forme polaire au point actuel
// ainsi que la variation première et seconde suivant V
PoGrenoble_VV pot = PoGrenoble_et_VV(Qeps,inv_n0);
// calcul du potentiel pour une petite variation suivant bIIb
Invariant inv_n = inv_n0;
double delta = limite_inf_bIIb; double unsurdelta = 1./delta;
inv_n.bIIb += delta;
Invariant0QepsCosphi inv_spes_n = Invariant0Specif(inv_n);
// stockage éventuel des invariants et du potentiel
if (sortie_post)
{ save_resulHyperD.invP->V= inv.V;
save_resulHyperD.invP->Qeps= inv_spes_n.Qeps;
save_resulHyperD.invP->cos3phi= inv_spes_n.cos3phi;
save_resulHyperD.invP->potentiel= pot.E;
};
PoGrenoble_V pot_V = PoGrenoble_et_V(inv_spes_n,inv_n);
// et pour la dérivée seconde de nouveau une petite variation
inv_n.bIIb += delta;
inv_spes_n = Invariant0Specif(inv_n);
double E_n = PoGrenoble(inv_spes_n,inv_n);
// calcul de potret
potret.E = pot.E * inv.V * jacobien0;
potret.Ks = pot.Ks;
// dérivées premières
potret.EV = jacobien0 * (inv.V * pot.EV + pot.E);
potret.EbIIb = jacobien0 * inv.V * ( pot_V.E - pot.E)*unsurdelta;
{double inter = (pot_V.E - pot.E);
// on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// le potentiel doit-être croissant, on doit donc avoir toujours inter > 0
if (inter < ConstMath::trespetit)
{ potret.EbIIb = 0.; }
else
{potret.EbIIb = jacobien0 * inv.V * inter * unsurdelta;};
};
// dérivées secondes
potret.EVV = jacobien0 * (2. * pot.EV + inv.V * pot.EVV );
{double inter = ((E_n + pot.E) - 2.*pot_V.E ); //(E_n - 2.*pot_V.E + pot.E)
// on est obligé de faire ainsi car sinon on a des erreurs d'arrondi
// la dérivée seconde du potentiel doit-être elle aussi > 0 (au moins à l'origine)
if (inter < ConstMath::trespetit)
{ potret.EbIIb2 = 0.; }
else
{potret.EbIIb2 = jacobien0 * inv.V * inter
* unsurdelta * unsurdelta;};
//cout << "\n debug ((E_n + pot.E) - 2.*pot_V.E )= " << ((E_n + pot.E) - 2.*pot_V.E ) << endl;
}
// pour la variation de la dérivée / V, a priori on peut avoir les deux signes
// donc on ne fait pas de filtrage
potret.EbIIbV = potret.EbIIb / inv.V
+ jacobien0 * inv.V * (pot_V.EV - pot.EV) *unsurdelta;
// pour les dérivées selon l'angle tous les termes sont considéré nul
potret.EIeps = 0.;
potret.EIeps2 = 0.; potret.EIepsV = 0.; potret.EbIIbIeps = 0.;
// appel d'une vérification éventuelle
// Verif_Potentiel_et_var(inv,jacobien0,potret);
};
// retour
// --- debug --------
//cout << "\n E= " << potret.E << " EV= " << potret.EV << " EbIIb= " << potret.EbIIb << " EIeps= " << potret.EIeps
// << "\n EVV= " << potret.EVV << " EbIIb2= " << potret.EbIIb2 << " EIeps2= " << potret.EIeps2
// << "\n EbIIbIeps = " << potret.EbIIbIeps << " EIepsV= " << potret.EIepsV << " EbIIbV= " << potret.EbIIbV << endl;
// ---- fin débug -------
return potret;
};
// calcul de l'invariant Qeps seul valable même si Qeps est nul
double Hyper3D::Invariant0Qeps(const Invariant & inv)
{ // les nouveaux invariants
double Qeps = sqrt(2. * inv.bIIb);
// if (avec_regularisation)
// {Qeps += fact_regularisation;}
return Qeps;
};
// calcul des invariants Qeps,cos3phi seuls
Hyper3D::Invariant0QepsCosphi Hyper3D::Invariant0Specif(const Invariant & inv)
{ Hyper3D::Invariant0QepsCosphi inv_ret; // def du retour
// les nouveaux invariants
inv_ret.Qeps = sqrt(2. * inv.bIIb);
double inv_ret_Qeps = inv_ret.Qeps; // init pour être éventuellement modifié
// ensuite avec la régularisation
if (avec_regularisation)
{inv_ret_Qeps += fact_regularisation;}
double unsurQeps = 1./inv_ret_Qeps; // pour simplifier
// --- angle de phase
double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
double troisracine6 = 3. * sqrt(6.);
inv_ret.cos3phi = troisracine6 * inv.bIIIb * unsurQeps3;
if (inv_ret.cos3phi > 1.) inv_ret.cos3phi=1.;
if (inv_ret.cos3phi < -1.) inv_ret.cos3phi=-1.;
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::Invariant0Specif(" << endl;
//fin debug
// retour
return inv_ret;
};
// calcul des invariants Qeps,cos3phi en fonction de V, bIIb Ieps
// on utilise pas l'invariant bIIIb,
Hyper3D::Invariant0QepsCosphi Hyper3D::InvariantDe_V_bIIb_Ieps(const Invariant & inv)
{ Hyper3D::Invariant0QepsCosphi inv_ret; // def du retour
///******** en fait ça merdouille, le bIIIb calculé donne n'importe quoi quand c'est petit donc
/// a éviter !!!
/// a défaut d'autre chose pour l'instant on utilise Invariant0Specif
///
// les nouveaux invariants
inv_ret.Qeps = sqrt(2. * inv.bIIb);
double inv_ret_Qeps = inv_ret.Qeps; // init pour être éventuellement modifié
double inv_V = inv.V; // init pour être éventuellement modifié
// ensuite avec la régularisation
if (avec_regularisation)
{inv_ret_Qeps += fact_regularisation;
inv_V += fact_regularisation;
};
double unsurQeps = 1./inv_ret_Qeps; // pour simplifier
// on recalcul bIIIb
double bIIIb=-1./(8.*inv_V*inv_V) + 1./8.-1./4.*inv.Ieps+1./6.*inv.Ieps*inv.Ieps
-1./27.*inv.Ieps*inv.Ieps*inv.Ieps
- 0.5*inv.bIIb+1./3.*inv.bIIb* inv.Ieps;
// --- angle de phase avec la valeur locale de bIIIb
double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
double troisracine6 = 3. * sqrt(6.);
inv_ret.cos3phi = troisracine6 * bIIIb * unsurQeps3;
// on limite les valeurs au normal
if (inv_ret.cos3phi > 1.) { inv_ret.cos3phi = 1.;}
else if (inv_ret.cos3phi < -1.) { inv_ret.cos3phi = -1.;};
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::InvariantDe_V_bIIb_Ieps(" << endl;
//fin debug
// retour
return inv_ret;
};
// calcul de l'invariants Qeps ainsi que sa variation par rapport
// à bIIb (Qeps doit être non nul)
Hyper3D::InvariantQeps Hyper3D::InvariantSpecifSansPhase(const Invariant & inv)
{ Hyper3D::InvariantQeps inv_ret; // def du retour
// tout d'abord le nouvel invariant
inv_ret.Qeps = sqrt(2. * inv.bIIb);
// maintenant la variation
if (avec_regularisation)
{inv_ret.dQepsdbIIb = 1./(inv_ret.Qeps+fact_regularisation);}
else
{inv_ret.dQepsdbIIb = 1./inv_ret.Qeps;};
// retour
return inv_ret;
};
// calcul des invariants Qeps,cos3phi, ainsi que leur variation par rapport
// aux trois invariants Ieps, V et bIIb
Hyper3D::InvariantQepsCosphi Hyper3D::InvariantSpecif(const Invariant & inv)
{ Hyper3D::InvariantQepsCosphi inv_ret; // def du retour, init à 0
// tout d'abord les nouveaux invariants
inv_ret.Qeps = sqrt(2. * inv.bIIb);
double inv_V = inv.V; // init pour être éventuellement modifié
double unsurQeps; // init = 1./inv_ret.Qeps; // pour simplifier
if (avec_regularisation)
{unsurQeps = 1./(inv_ret.Qeps+fact_regularisation);
inv_V += fact_regularisation;
}
else {unsurQeps = 1./inv_ret.Qeps;};
// maintenant les variations / pour Qeps
double untiers = 1./3.;
inv_ret.dQepsdbIIb = unsurQeps;
// --- angle de phase
// dans le cas ou Qeps est trop petit on laisse les 0
// if ((inv_ret.Qeps >= limite_inf_Qeps)||(avec_regularisation))
if (inv_ret.Qeps >= limite_inf_Qeps)
{ double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
double unsurQeps5 = unsurQeps3 * unsurQeps * unsurQeps;
double troisracine6 = 3. * sqrt(6.);
double neufracine6 = 3. * troisracine6;
inv_ret.cos3phi = troisracine6 * inv.bIIIb * unsurQeps3;
if (inv_ret.cos3phi > 1.) inv_ret.cos3phi=1.;
if (inv_ret.cos3phi < -1.) inv_ret.cos3phi=-1.;
// puis les variations
inv_ret.dcos3phidV = troisracine6 / 4. * unsurQeps3 /(inv_V * inv_V * inv_V);
inv_ret.dcos3phidIeps = troisracine6 * unsurQeps3 * (-0.25 + inv.bIIb*untiers
+ inv.Ieps*untiers - inv.Ieps * inv.Ieps/9.);
inv_ret.dcos3phidbIIb = - neufracine6 * unsurQeps5 * inv.bIIIb
+ troisracine6 * unsurQeps3 * (inv.Ieps*untiers - 0.5);
};
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::InvariantSpecif(" << endl;
//fin debug
// retour
return inv_ret;
};
// calcul de l'invariants Qeps ainsi que ses variations première et seconde
// par rapport à bIIb (Qeps doit être non nul)
Hyper3D::Invariant2Qeps Hyper3D::Invariant2SpecifSansPhase(const Invariant & inv)
{ Hyper3D::Invariant2Qeps inv_ret; // def du retour et init à 0
// tout d'abord les nouveaux invariants
inv_ret.Qeps = sqrt(2. * inv.bIIb);
double unsurQeps; // init = 1./inv_ret.Qeps; // pour simplifier
if (avec_regularisation)
{unsurQeps = 1./(inv_ret.Qeps+fact_regularisation);}
else {unsurQeps = 1./inv_ret.Qeps;};
// double unsurQeps = 1./inv_ret.Qeps; // pour simplifier
// maintenant les variations
inv_ret.dQepsdbIIb = unsurQeps;
// et les variations secondes
inv_ret.dQepsdbIIb2 = -1. * unsurQeps * unsurQeps * unsurQeps;
// retour
return inv_ret;
};
// calcul des invariants Qeps,cos3phi, ainsi que leur variation première et seconde
// par rapport aux trois invariants Ieps, V et bIIb
Hyper3D::Invariant2QepsCosphi Hyper3D::Invariant2Specif(const Invariant & inv)
{ Hyper3D::Invariant2QepsCosphi inv_ret; // def du retour, init à 0
// tout d'abord les nouveaux invariants
inv_ret.Qeps = sqrt(2. * inv.bIIb);
double unsurQeps; // init = 1./inv_ret.Qeps; // pour simplifier
double inv_V = inv.V; // init pour être éventuellement modifié
if (avec_regularisation)
{unsurQeps = 1./(inv_ret.Qeps+fact_regularisation);
inv_V += fact_regularisation;
}
else {unsurQeps = 1./inv_ret.Qeps;};
// double unsurQeps = 1./inv_ret.Qeps; // pour simplifier
// et les variations / Qeps
double untiers = 1./3.;
inv_ret.dQepsdbIIb = unsurQeps;
double unsurQeps3 = unsurQeps * unsurQeps * unsurQeps;
inv_ret.dQepsdbIIb2 = -1. * unsurQeps3;
// --- angle de phase
// dans le cas ou Qeps est trop petit on laisse les 0
// if ((inv_ret.Qeps >= limite_inf_Qeps)||(avec_regularisation))
if (inv_ret.Qeps >= limite_inf_Qeps)
{ double unsurQeps5 = unsurQeps3 * unsurQeps * unsurQeps;
double unsurQeps7 = unsurQeps5 * unsurQeps * unsurQeps;
double racine6=sqrt(6.);
double troisracine6 = 3. * racine6;
double neufracine6 = 3. * troisracine6;
inv_ret.cos3phi = troisracine6 * inv.bIIIb * unsurQeps3;
if (inv_ret.cos3phi > 1.) inv_ret.cos3phi=1.;
if (inv_ret.cos3phi < -1.) inv_ret.cos3phi=-1.;
double unsurV3 = 1. / (inv_V * inv_V * inv_V);
// maintenant les variations
inv_ret.dcos3phidV = troisracine6 * unsurQeps3 * unsurV3 / 4. ;
inv_ret.dcos3phidIeps = troisracine6 * unsurQeps3 * (-0.25 + inv.bIIb*untiers
+ inv.Ieps*untiers - inv.Ieps * inv.Ieps/9.);
inv_ret.dcos3phidbIIb = - neufracine6 * unsurQeps5 * inv.bIIIb
+ troisracine6 * unsurQeps3 * (inv.Ieps*untiers - 0.5);
// les variations secondes
inv_ret.dcos3phidV2 = - neufracine6 * unsurQeps3 * unsurV3 / (4. * inv_V);
inv_ret.dcos3phidIeps2 = troisracine6 * unsurQeps3 * (1./3. - 2./9. * inv.Ieps);
inv_ret.dcos3phidbIIb2 = 15. * troisracine6 * inv.bIIIb * unsurQeps7
- 2.*neufracine6 * unsurQeps5 * (untiers * inv.Ieps - 0.5);
inv_ret.dcos3phidVdbIIb = - neufracine6 * unsurQeps5 * unsurV3 * 0.25;
inv_ret.dcos3phidIepsdbIIb = racine6 * unsurQeps3
- neufracine6 * unsurQeps5
* (-0.25 + inv.bIIb*untiers + untiers * inv.Ieps
- 1./9. * inv.Ieps * inv.Ieps );
};
//debug
//cout << "\n inv_ret.cos3phi= "<<inv_ret.cos3phi << "Hyper3D::Invariant2Specif(" << endl;
//fin debug
// retour
return inv_ret;
};
// Lecture des paramètre specifique sur fichier
void Hyper3D::LectParaSpecifiques(UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D,LesFonctions_nD& lesFonctionsnD)
{ string nom_class_methode("Hyper3D::LectParaSpecifiques");
// cas limite_inf_Qeps (variable stockée dans Hyper3D.h)
string nom1;
if (strstr(entreePrinc->tablcar,"limite_inf_Qeps_")!=NULL)
{ *(entreePrinc->entree) >> nom1 ;
if (nom1 != "limite_inf_Qeps_")
{ cout << "\n erreur en lecture du drapeau de limite_inf_Qeps, on attendait le mot cle "
<< " limite_inf_Qeps_ on a lue: " << nom1 ;
entreePrinc->MessageBuffer("**Hyper3D::LectParaSpecifiques((.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> limite_inf_Qeps;
};
// cas limite_inf_Qeps (variable stockée dans Hyper3D.h)
if (strstr(entreePrinc->tablcar,"limite_inf_bIIb_")!=NULL)
{ *(entreePrinc->entree) >> nom1 ;
if (nom1 != "limite_inf_bIIb_")
{ cout << "\n erreur en lecture du drapeau de limite_inf_bIIb, on attendait le mot cle "
<< " limite_inf_bIIb_ on a lue: " << nom1 ;
entreePrinc->MessageBuffer("**Hyper3D::LectParaSpecifiques((.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> limite_inf_bIIb;
};
// cas avec régularisation (variable stockée dans Hyper3D.h)
avec_regularisation = false; // par défaut
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 cle "
<< " avec_regularisation_ on a lue: " << nom1 ;
entreePrinc->MessageBuffer("**Hyper3D::LectParaSpecifiques((.....**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> fact_regularisation;
avec_regularisation=true;
};
// lecture de l'indication éventuelle du post traitement
string le_mot_cle = "sortie_post_";
entreePrinc->Lecture_un_parametre_int(0,nom_class_methode,0,1,le_mot_cle,sortie_post);
};
// I/O des paramètres spécifiques: ex: régularisation, limite inf ...
void Hyper3D::Ecriture_para_specifiques(ofstream& sort,const int cas)
{ if (cas == 1)
{// la régularisation éventuelle
if (avec_regularisation)
sort << " avec_regularisation " << fact_regularisation;
else
sort << " sans_regularisation ";
// limites
sort << " limite_inf_Qeps= " << limite_inf_Qeps
<< " limite_inf_bIIb= " << limite_inf_bIIb << " " ;
};
};
void Hyper3D::Lecture_para_specifiques(ifstream& ent,const int cas)
{ if (cas == 1)
{// la régularisation éventuelle
string nom;
ent >> nom;
if (nom == "avec_regularisation")
{ent >> fact_regularisation;avec_regularisation=true;}
else if (nom == "sans_regularisation")
avec_regularisation=false;
else
{cout << "\n erreur en lecture du facteur de regularisation (loi hyperelastique) "
<< " on a lue: " << nom << " au lieu du mot cle: avec_regularisation ou sans_regularisation";
Sortie(1);
};
// limites
ent >> nom;
if (nom == "limite_inf_Qeps=")
ent >> limite_inf_Qeps;
else
{cout << "\n erreur en lecture de limite_inf_Qeps (loi hyperelastique) "
<< " on a lue: " << nom << " au lieu du mot cle: limite_inf_Qeps= ";
Sortie(1);
};
ent >> nom;
if (nom == "limite_inf_bIIb=")
ent >> limite_inf_bIIb;
else
{cout << "\n erreur en lecture de limite_inf_bIIb (loi hyperelastique) "
<< " on a lue: " << nom << " au lieu du mot cle: limite_inf_bIIb= ";
Sortie(1);
};
};
};
///=============== 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 Hyper3D::Verif_Potentiel_et_var(const InvariantVarDdl& inv,const double& jacobien0
,const PotenSansPhaseAvecVar& potret )
{ // dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_Potentiel_et_var == 0)
{ cout << "\n ****verification des derivees du potentiels par rapport aux invariants****";
cout << "\n Hyper3D::Verif_Potentiel_et_var \n";
}
indic_Verif_Potentiel_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 Qeps = Invariant0Qeps(inv_n0);
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 = limite_inf_bIIb;
double mini_val = ConstMath::pasmalpetit;
Invariant inv_et_dbIIb = inv_n0;inv_et_dbIIb.bIIb += delta;
double Qeps_et_dbIIb = Invariant0Qeps(inv_et_dbIIb);
double E_et_dbIIb = PoGrenoble(Qeps_et_dbIIb,inv_et_dbIIb);
Invariant inv_et_dbIIb2 = inv_n0;inv_et_dbIIb2.bIIb += 2. * delta;
double Qeps_et_dbIIb2 = Invariant0Qeps(inv_et_dbIIb2);
double E_et_dbIIb2 = PoGrenoble(Qeps_et_dbIIb2,inv_et_dbIIb2);
Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
double Qeps_et_dV = Invariant0Qeps(inv_et_dV);
double E_et_dV = PoGrenoble(Qeps_et_dV,inv_et_dV);
Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
double Qeps_et_dV2 = Invariant0Qeps(inv_et_dV2);
double E_et_dV2 = PoGrenoble(Qeps_et_dV2,inv_et_dV2);
Invariant inv_et_dVdbIIb = inv_n0;
inv_et_dVdbIIb.V += delta;
inv_et_dVdbIIb.bIIb += delta;
double Qeps_et_dVdbIIb = Invariant0Qeps(inv_et_dVdbIIb);
double E_et_dVdbIIb = PoGrenoble(Qeps_et_dVdbIIb,inv_et_dVdbIIb);
// calcul des dérivées premières
double E_V = jacobien0 *(E_et_dV * inv_et_dV.V - E_n * inv.V)/delta;
double E_bIIb = (jacobien0 * inv.V)*(E_et_dbIIb - E_n)/delta;
// calcul des dérivées secondes
double E_V2 = jacobien0 *(E_et_dV2 * inv_et_dV2.V - 2.*E_et_dV * inv_et_dV.V + E_n * inv.V)
/(delta * delta);
double E_bIIb2 = (jacobien0 * inv.V)*(E_et_dbIIb2 - 2.*E_et_dbIIb + E_n)
/(delta * delta);
double E_V_a_dbIIb = jacobien0 * (E_et_dVdbIIb * inv_et_dVdbIIb.V - E_et_dbIIb * inv.V)
/delta;
double E_VbIIb = ( E_V_a_dbIIb - E_V)/delta;
// comparaison avec les valeurs de dérivées analytiques
int erreur = 0;
if (diffpourcent(potret.EV,E_V,E_V,0.05))
{ if (MiN(Dabs(potret.EV),Dabs(E_V)) <= mini_val)
{if( MaX(Dabs(potret.EV),Dabs(E_V)) > 50.*delta)
if (Dabs(potret.EV) == 0.)
// on regarde si la dérivée numérique tend vers zéro lorque l'incrément tend vers 0
{
Invariant inv_et_ddV = inv_n0;inv_et_ddV.V += delta/10.;
double Qeps_et_ddV = Invariant0Qeps(inv_et_ddV);
double E_et_ddV = PoGrenoble(Qeps_et_ddV,inv_et_ddV);
double E_ddV = jacobien0 *(E_et_ddV * inv_et_ddV.V - E_n * inv.V)/(delta*10.);
if (10.* Dabs(E_ddV) > Dabs(E_V)) erreur = 5;
}
else {erreur = 5;cout << " \n err5: potret.EV= " << potret.EV << " E_V= " << E_V;}
else {erreur = 52;cout << " \n err52: potret.EV= " << potret.EV << " E_V= " << E_V;};
}
else {erreur = 53;cout << " \n err53: potret.EV= " << potret.EV << " E_V= " << E_V;};
};
// cout << " \n info: potret.EV= " << potret.EV << " E_V= " << E_V;
if (diffpourcent(potret.EbIIb,E_bIIb,E_bIIb,0.05))
{if (MiN(Dabs(potret.EbIIb),Dabs(E_bIIb)) <= mini_val)
{if( MaX(Dabs(potret.EbIIb),Dabs(E_bIIb)) > 50.*delta)
{erreur = 1;cout << " \n err 1: potret.EbIIb= " << potret.EbIIb << " E_bIIb= " << E_bIIb;}
}
else {erreur = 11;cout << " \n err11: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
};
if (diffpourcent(potret.EbIIb2,E_bIIb2,E_bIIb2,0.05))
{if (MiN(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) <= mini_val)
{if( MaX(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) > 50.*delta)
{erreur = 2;cout << " \n err 2: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;}
}
else {erreur = 21;cout << " \n err21: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;};
};
if (diffpourcent(potret.EVV,E_V2,E_V2,0.05))
{if (MiN(Dabs(potret.EVV),Dabs(E_V2)) <= mini_val)
{if( MaX(Dabs(potret.EVV),Dabs(E_V2)) > 50.*delta)
{erreur = 3;cout << " \n err 3: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;}
}
else {erreur = 31;cout << " \n err31: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;};
};
if (diffpourcent(potret.EbIIbV,E_VbIIb,E_VbIIb,0.05))
{if (MiN(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) <= mini_val)
{if( MaX(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) > 50.*delta)
{erreur = 4;cout << " \n err 4: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;}
}
else {erreur = 41;cout << " \n err41: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
};
// if (erreur)
// { cout << "\n erreur dans le calcul analytique des dérivées du potentiel";
// cout << "\n Hyper3D::Verif_Potentiel_et_var, numéro d'appel = "
// << indic_Verif_Potentiel_et_var;
// Sortie(1);
// }
};
// vérif des dérivées du potentiels par rapport aux invariants, ceci par différences finies
// cas avec la phase
void Hyper3D::Verif_Potentiel_et_var(const InvariantVarDdl& inv,const double& jacobien0
,const PotenAvecPhaseAvecVar& potret )
{ // dans le cas du premier passage on indique qu'il y a vérification
if (indic_Verif_Potentiel_et_var == 0)
{ cout << "\n ****verification des derivees du potentiels par rapport aux invariants****";
cout << "\n Hyper3D::Verif_Potentiel_et_var \n";
}
indic_Verif_Potentiel_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 is_sp0 = Invariant0Specif(inv_n0);
double E_n = PoGrenoble(is_sp0,inv_n0); // la valeur du potentiel de référence
// ici on est avec phase donc trois invariants indépendant V et Qeps et cos3phi,
// 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 = limite_inf_bIIb;
double mini_val = ConstMath::pasmalpetit;
Invariant inv_et_dbIIb = inv_n0;inv_et_dbIIb.bIIb += delta;
Invariant0QepsCosphi is_et_dbIIb = Invariant0Specif(inv_et_dbIIb);
double E_et_dbIIb = PoGrenoble(is_et_dbIIb,inv_et_dbIIb);
Invariant inv_et_dbIIb2 = inv_n0;inv_et_dbIIb2.bIIb += 2. * delta;
Invariant0QepsCosphi is_et_dbIIb2 = Invariant0Specif(inv_et_dbIIb2);
double E_et_dbIIb2 = PoGrenoble(is_et_dbIIb2,inv_et_dbIIb2);
Invariant inv_et_dV = inv_n0;inv_et_dV.V += delta;
Invariant0QepsCosphi is_et_dV = Invariant0Specif(inv_et_dV);
double E_et_dV = PoGrenoble(is_et_dV,inv_et_dV);
Invariant inv_et_dV2 = inv_n0;inv_et_dV2.V += 2. * delta;
Invariant0QepsCosphi is_et_dV2 = Invariant0Specif(inv_et_dV2);
double E_et_dV2 = PoGrenoble(is_et_dV2,inv_et_dV2);
Invariant inv_et_dVdbIIb = inv_n0;
inv_et_dVdbIIb.V += delta;
inv_et_dVdbIIb.bIIb += delta;
Invariant0QepsCosphi is_et_dVdbIIb = Invariant0Specif(inv_et_dVdbIIb);
double E_et_dVdbIIb = PoGrenoble(is_et_dVdbIIb,inv_et_dVdbIIb);
// cas avec Ieps
Invariant inv_et_dIeps = inv_n0;inv_et_dIeps.Ieps += delta;
Invariant0QepsCosphi is_et_dIeps = Invariant0Specif(inv_et_dIeps);
double E_et_dIeps = PoGrenoble(is_et_dIeps,inv_et_dIeps);
Invariant inv_et_dIeps2 = inv_n0;inv_et_dIeps2.Ieps += 2. * delta;
Invariant0QepsCosphi is_et_dIeps2 = Invariant0Specif(inv_et_dIeps2);
double E_et_dIeps2 = PoGrenoble(is_et_dIeps2,inv_et_dIeps2);
Invariant inv_et_dIepsdbIIb = inv_n0;
inv_et_dIepsdbIIb.Ieps += delta;
inv_et_dIepsdbIIb.bIIb += delta;
Invariant0QepsCosphi is_et_dIepsdbIIb = Invariant0Specif(inv_et_dIepsdbIIb);
double E_et_dIepsdbIIb = PoGrenoble(is_et_dIepsdbIIb,inv_et_dIepsdbIIb);
Invariant inv_et_dIepsdV = inv_n0;
inv_et_dIepsdV.Ieps += delta;
inv_et_dIepsdV.V += delta;
Invariant0QepsCosphi is_et_dIepsdV = Invariant0Specif(inv_et_dIepsdV);
double E_et_dIepsdV = PoGrenoble(is_et_dIepsdV,inv_et_dIepsdV);
// calcul des dérivées premières
// calcul des dérivées premières
double E_V = jacobien0 *(E_et_dV * inv_et_dV.V - E_n * inv.V)/delta;
double E_bIIb = (jacobien0 * inv.V)*(E_et_dbIIb - E_n)/delta;
double E_Ieps = (jacobien0 * inv.V)*(E_et_dIeps - E_n )/delta;
// calcul des dérivées secondes hors phases
double E_V2 = jacobien0 *(E_et_dV2 * inv_et_dV2.V - 2.*E_et_dV * inv_et_dV.V + E_n * inv.V)
/(delta * delta);
double E_bIIb2 = (jacobien0 * inv.V)*(E_et_dbIIb2 - 2.*E_et_dbIIb + E_n) /(delta * delta);
double E_Ieps2 = (jacobien0 * inv.V) *(E_et_dIeps2 - 2.*E_et_dIeps + E_n )/(delta * delta);
double E_V_a_dbIIb = jacobien0 * (E_et_dVdbIIb * inv_et_dVdbIIb.V - E_et_dbIIb * inv_et_dbIIb.V) /delta;
double E_VbIIb = ( E_V_a_dbIIb - potret.EV)/delta;
double E_Ieps_a_dbIIb = (jacobien0 * inv.V) * (E_et_dIepsdbIIb - E_et_dbIIb) /delta;
double E_bIIbIeps = ( E_Ieps_a_dbIIb - potret.EIeps)/delta;
double E_Ieps_a_dV = jacobien0 * (E_et_dIepsdV * inv_et_dIepsdV.V - E_et_dV * inv.V) /delta;
double E_IepsV = ( E_Ieps_a_dV - potret.EIeps)/delta;
// comparaison avec les valeurs de dérivées analytiques
int erreur = 0;
if (diffpourcent(potret.EV,E_V,E_V,0.05))
{if (MiN(Dabs(potret.EV),Dabs(E_V)) <= mini_val)
{if( MaX(Dabs(potret.EV),Dabs(E_V)) > 50.*delta)
if (Dabs(potret.EV) == 0.)
// on regarde si la dérivée numérique tend vers zéro lorque l'incrément tend vers 0
{
Invariant inv_et_ddV = inv_n0;inv_et_ddV.V += delta/10.;
double Qeps_et_ddV = Invariant0Qeps(inv_et_ddV);
double E_et_ddV = PoGrenoble(Qeps_et_ddV,inv_et_ddV);
double E_ddV = jacobien0 *(E_et_ddV * inv_et_ddV.V - E_n * inv.V)/(delta*10.);
if (10.* Dabs(E_ddV) > Dabs(E_V)) erreur = 5;
}
else {erreur = 5;cout << " \n err5: potret.EV= " << potret.EV << " E_V= " << E_V;}
else {erreur = 52;cout << " \n err52: potret.EV= " << potret.EV << " E_V= " << E_V;};
}
else
{erreur = 53;cout << " \n err53: potret.EV= " << potret.EV << " E_V= " << E_V;};
};
if (diffpourcent(potret.EbIIb,E_bIIb,E_bIIb,0.05))
{if (MiN(Dabs(potret.EbIIb),Dabs(E_bIIb)) <= mini_val)
{if( MaX(Dabs(potret.EbIIb),Dabs(E_bIIb)) > 50.*delta)
{erreur = 1;cout << " \n err 1: potret.EbIIb= " << potret.EbIIb << " E_bIIb= " << E_bIIb;}
}
else {erreur = 11;cout << " \n err 11: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
};
if (diffpourcent(potret.EbIIb2,E_bIIb2,E_bIIb2,0.05))
{if (MiN(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) <= mini_val)
{if( MaX(Dabs(potret.EbIIb2),Dabs(E_bIIb2)) > 50.*delta)
{erreur = 2;cout << " \n err 2: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;}
}
else {erreur = 21;cout << " \n err 21: potret.EbIIb2= " << potret.EbIIb2 << " E_bIIb2= " << E_bIIb2;};
};
if (diffpourcent(potret.EVV,E_V2,E_V2,0.05))
{if (MiN(Dabs(potret.EVV),Dabs(E_V2)) <= mini_val)
{if( MaX(Dabs(potret.EVV),Dabs(E_V2)) > 50.*delta)
{erreur = 3;cout << " \n err 3: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;}
}
else {erreur = 31;cout << " \n err 31: potret.EVV= " << potret.EVV << " E_V2= " << E_V2;};
};
if (diffpourcent(potret.EbIIbV,E_VbIIb,E_VbIIb,0.05))
{if (MiN(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) <= mini_val)
{if( MaX(Dabs(potret.EbIIbV),Dabs(E_VbIIb)) > 50.*delta)
{erreur = 4;cout << " \n err 4: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;}
}
else {erreur = 41;cout << " \n err 41: potret.EbIIbV= " << potret.EbIIbV << " E_VbIIb= " << E_VbIIb;};
};
if (diffpourcent(potret.EIeps,E_Ieps,E_Ieps,0.05))
{if (MiN(Dabs(potret.EIeps),Dabs(E_Ieps)) <= mini_val)
{if( MaX(Dabs(potret.EIeps),Dabs(E_Ieps)) > 50.*delta)
{erreur = 6;cout << " \n err 6: potret.EIeps= " << potret.EIeps << " E_Ieps= " << E_Ieps;}
}
else {erreur = 61;cout << " \n err 61: potret.EIeps= " << potret.EIeps << " E_Ieps= " << E_Ieps;};
};
if (diffpourcent(potret.EIeps2,E_Ieps2,E_Ieps2,0.05))
{if (MiN(Dabs(potret.EIeps2),Dabs(E_Ieps2)) <= mini_val)
{if( MaX(Dabs(potret.EIeps2),Dabs(E_Ieps2)) > 50.*delta)
{erreur = 7;cout << " \n err 7: potret.EIeps2= " << potret.EIeps2 << " E_Ieps2= " << E_Ieps2;}
}
else {erreur = 71;cout << " \n err 71: potret.EIeps2= " << potret.EIeps2 << " E_Ieps2= " << E_Ieps2;};
};
if (diffpourcent(potret.EbIIbIeps,E_bIIbIeps,E_bIIbIeps,0.05))
{if (MiN(Dabs(potret.EbIIbIeps),Dabs(E_bIIbIeps)) <= mini_val)
{if( MaX(Dabs(potret.EbIIbIeps),Dabs(E_bIIbIeps)) > 50.*delta)
{erreur = 8;cout << " \n err 8: potret.EbIIbIeps= " << potret.EbIIbIeps << " E_bIIbIeps= " << E_bIIbIeps;}
}
else {erreur = 81;cout << " \n err 81: potret.EbIIbIeps= " << potret.EbIIbIeps << " E_bIIbIeps= " << E_bIIbIeps;};
};
if (diffpourcent(potret.EIepsV,E_IepsV,E_IepsV,0.05))
{if (MiN(Dabs(potret.EIepsV),Dabs(E_IepsV)) <= mini_val)
{if( MaX(Dabs(potret.EIepsV),Dabs(E_IepsV)) > 50.*delta)
{erreur = 9;cout << " \n err 9: potret.EIepsV= " << potret.EIepsV << " E_IepsV= " << E_IepsV;}
}
else {erreur = 91;cout << " \n err 91: potret.EIepsV= " << potret.EIepsV << " E_IepsV= " << E_IepsV;};
};
/* if (erreur)
{ cout << "\n erreur dans le calcul analytique des derivees du potentiel";
cout << "\n Hyper3D::Verif_Potentiel_et_var, numero d'appel = "
<< indic_Verif_Potentiel_et_var << " erreur = " << erreur << endl;
Sortie(1);
}*/
};
Reference in a new issue View git blame Copy permalink