rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/herezh_pp/comportement/hysteresis/Hysteresis3D_2.cc

2818 lines
161 KiB
C++
Executable file
Raw Blame History

// FICHIER : Hysteresis3D.cc
// CLASSE : Hysteresis3D
// 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-2021 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 "Debug.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "ConstMath.h"
#include "MatLapack.h"
#include "ExceptionsLoiComp.h"
#include "TenseurQ3gene.h"
#include "Hysteresis3D.h"
// ========== codage des METHODES VIRTUELLES protegees:================
// calcul des contraintes
void Hysteresis3D::Calcul_SigmaHH (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl,
TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB& epsBB_,
TenseurBB& delta_epsBB_,TenseurBB & gijBB_,
TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,double& ,double& ,
TenseurHH & sigHH_
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Expli_t_tdt& )
{
#ifdef MISE_AU_POINT
if (epsBB_.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hysteresis3D::Calcul_SigmaHH\n";
Sortie(1);
};
if (tab_ddl.NbDdl() != d_gijBB_.Taille())
{ cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_ !\n";
cout << " Hysteresis3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_); // passage en dim 3
const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_); // " " " "
const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_); // " " " "
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // " " " "
// Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t); // " " " "
SaveResulHysteresis3D & save_resul = *((SaveResulHysteresis3D*) saveResul);
// initialisation des variables de travail
save_resul.Init_debut_calcul();
// initialisation éventuelle des variables thermo-dépendantes
Init_thermo_dependance();
// cas d'une dépendance à la phase de la déformation totale
if (avec_phase_paradefEqui || avec_phaseEps)
{ xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
ParaMateriauxAvecPhaseEpsilon(gijHH,epsBB,xnp,xmu,Qzero);
};
// === 1 == calcul de l'avancement temporel sur 1 pas,
int cas = 1; // cas normal avec dérivée de Jauman
// try { // débug
Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
// débug---------
/* }
catch (ErrNonConvergence_loiDeComportement excep)
{ // initialisation des variables de travail
save_resul.Init_debut_calcul();
// initialisation éventuelle des variables thermo-dépendantes
Init_thermo_dependance();
Avancement_temporel(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ);
}
catch ( ... )
{ // initialisation des variables de travail
save_resul.Init_debut_calcul();
// initialisation éventuelle des variables thermo-dépendantes
Init_thermo_dependance();
Avancement_temporel(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ);
}
*/
// Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// fin débug -------------
// on libère les tenseurs intermédiaires
LibereTenseur();LibereTenseurQ();
};
// calcul des contraintes a t+dt et de ses variations
void Hysteresis3D::Calcul_DsigmaHH_tdt (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl
,BaseB& ,TenseurBB & ,TenseurHH &
,BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> &
,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,TenseurBB & delta_epsBB_
,TenseurBB & gijBB_tdt ,TenseurHH & gijHH_tdt
,Tableau <TenseurBB *>& d_gijBB_tdt
,Tableau <TenseurHH *>& d_gijHH_tdt,double& ,double&
,Vecteur& ,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite
,double& module_cisaillement,const Met_abstraite::Impli& )
{
#ifdef MISE_AU_POINT
if (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hysteresis3D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille())
{ cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt !\n";
cout << " Hysteresis3D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
#endif
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // " " " "
const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_tdt); // " " " "
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " "
// Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t); // " " " "
SaveResulHysteresis3D & save_resul = *((SaveResulHysteresis3D*) saveResul);
// initialisation des variables de travail
save_resul.Init_debut_calcul();
// initialisation éventuelle des variables thermo-dépendantes
Init_thermo_dependance();
// cas d'une dépendance à la phase de la déformation totale
if (avec_phase_paradefEqui || avec_phaseEps)
{ xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
ParaMateriauxAvecPhaseEpsilon(gijHH,epsBB,xnp,xmu,Qzero);
};
// === 1 == calcul de l'avancement temporel sur 1 pas,
int cas = 1; // cas normal avec dérivée de Jauman
Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
// ==== 2 === calcul de l'opérateur tangent =============
Hysteresis3D::Dsig_d_ddl(d_sigHH,epsBB_tdt,gijHH_tdt,d_gijBB_tdt,d_gijHH_tdt,d_epsBB );
// débug---------
// Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// fin débug -------------
// on libère les tenseurs intermédiaires
LibereTenseur();LibereTenseurQ();
};
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormee
// le tenseur de déformation et son incrémentsont également en orthonormees
// si = false: les bases transmises sont utilisées
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void Hysteresis3D::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB&
,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB_,double& ,double&
,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Umat_cont& ex)
{
#ifdef MISE_AU_POINT
if (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hysteresis3D::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
const Tenseur3HH & gijHH = *((Tenseur3HH*) ex.gijHH_tdt); // " " " "
const Tenseur3BB & gijBB = *((Tenseur3BB*) ex.gijBB_tdt); // " " " "
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
const Tenseur3BB & delta_epsBB = *((Tenseur3BB*) &delta_epsBB_); // passage en dim 3
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " "
// Tenseur3HH & sigHH_i = *((Tenseur3HH*) &sigHH_t); // " " " "
SaveResulHysteresis3D & save_resul = *((SaveResulHysteresis3D*) saveResul);
// initialisation des variables de travail
save_resul.Init_debut_calcul();
// initialisation éventuelle des variables thermo-dépendantes
Init_thermo_dependance();
// cas d'une dépendance à la phase de la déformation totale
if (avec_phase_paradefEqui || avec_phaseEps)
{ xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
if (en_base_orthonormee)
ParaMateriauxAvecPhaseEpsilon(IdHH3,epsBB,xnp,xmu,Qzero);
else
ParaMateriauxAvecPhaseEpsilon(gijHH,epsBB,xnp,xmu,Qzero);
};
// === 1 == calcul de l'avancement temporel sur 1 pas,
// ici il n'y a pas de différence entre orthonormeee ou pas au niveau du type de calcul
int cas = 1; // donc cas avec dérivée de Jauman en mixte
if (en_base_orthonormee)
Avancement_temporel(delta_epsBB,IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ);
else
Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
// ==== 2 === calcul de l'opérateur tangent =============
// Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps);
// try { // débug
Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps,sigHH);
// débug---------
// }
// catch (ErrNonConvergence_loiDeComportement excep)
// { // initialisation des variables de travail
// save_resul.Init_debut_calcul();
// // initialisation éventuelle des variables thermo-dépendantes
// Init_thermo_dependance();
// if (en_base_orthonormee)
// Avancement_temporel(delta_epsBB,IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ);
// else
// Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
// Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps);
// }
// catch ( ... )
// { // initialisation des variables de travail
// racine.Change_taille(9);
// save_resul.Init_debut_calcul();
// // initialisation éventuelle des variables thermo-dépendantes
// Init_thermo_dependance();
// if (en_base_orthonormee)
// Avancement_temporel(delta_epsBB,IdHH3,IdBB3,cas,save_resul,sigHH,energ_t,energ);
// else
// Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
// Hysteresis3D::Dsig_depsilon(en_base_orthonormee,gijHH,d_sigma_deps);
// };
// Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// fin débug -------------
// on libère les tenseurs intermédiaires
LibereTenseur();LibereTenseurQ();
};
// =============== fonction protegee ============
// calcul de la fonction résidu de la résolution de l'équation constitutive
// l'argument test ramène
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
// fatal, qui invalide le calcul du résidu
Vecteur& Hysteresis3D::Residu_constitutif (const double & alpha,const Vecteur & x, int& test)
{ // récupération du vecteur delta_sigma
for (int i=1;i<=9;i++) // accroissement de t à tdt
rdelta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = x(i);
test = 1; // init par défaut
// variation de sigma de R à tdt
rdelta_sigma_barre_BH_Ratdt = sigma_i_barre_BH + rdelta_sigma_barre_tdt_BH - sigma_barre_BH_R;
delta_barre_alpha_epsBH = alpha * delta_barre_epsBH; // limitation de la charge
// -- calcul de QdeltaSigma
QdeltaSigmaRatdt = sqrt(rdelta_sigma_barre_BH_Ratdt.II());
// -- calculs relatifs aux angles de phases
Cal_wprime(QdeltaSigmaRatdt,rdelta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR
,wBase,rdelta_sigma_barre_tdt_BH,wprime,wprime_point,trajet_neutre);
// --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
if (avec_phaseDeltaSig)
{xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,rdelta_sigma_barre_BH_Ratdt,xmu,Qzero);
};
// -- modif 7 mai 2012 ---
// dans le cas où le phidt est négatif, l'équation constitutive n'est plus valide, dans ce cas on considère
// que phi est nulle, c-a-d que l'on évolue suivant un trajet neutre ce qui n'est sans doute pas l'idéal, mais permet d'éviter l'explosion
// -- calcul de phidt
double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH) ; // init partie radial // -- modif 16 juin 2012 --
if (!trajet_neutre) // -- modif 16 juin 2012 --
phidt -= QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime); // -- modif 16 juin 2012 --
// if ((phidt < 0) && force_phi_zero_si_negatif) // -- modif 16 juin 2012 --
if (phidt < 0) // -- modif 9 fev 2015 --
{if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 8)) || (Permet_affichage() > 8) )
cout << "\n phidt= " << phidt << "\n Hysteresis3D::Residu_constitutif(..." << endl;
if ( Dabs(phidt) < force_phi_zero_si_negatif) // -- modif 9 fev 2015 --
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5) )
{ cout << "\n attention on utilise une valeur de phi negative dans pour le calcul de l'equation constitutive "
<< " phidt= " << phidt;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::Residu_constitutif(..." << endl;
};
phidt = 0.;trajet_neutre = true;
test = -1; // c'est un pb mais pas fatal
};
};
// fin modif modif 7 mai 2012 ---
if (trajet_neutre)
{ // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
// -- calcul du résidu
residuBH = rdelta_sigma_barre_tdt_BH - 2. * xmu * delta_barre_alpha_epsBH;
// --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
double phitau = 0.;
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
}
else
{ // cas d'un trajet non neutre
// modif 7 mai 2012 --- on ne recalcule plus le phi qui a été calculé plus haut pour savoir s'il était négatif
// // -- calcul de phidt
// double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH)
// - QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime);
// --- si dépendance à la déformation équivalente
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phidt,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
// -- calcul de Beta
double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
double beta=0.; // init dans le cas où ne le calcul pas après
if (QdeltaSigmaRatdt >= ConstMath::petit)
// on ne calcul beta que si la norme de deltaSigmaRatdt est significative
{if (xnp == 2.) { beta = - 2. * xmu * unsurwprimeQ0_puiss_np;}
else {beta = - 2. * xmu * unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRatdt,xnp-2.);};
};
// -- calcul du résidu
residuBH = rdelta_sigma_barre_tdt_BH
- 2. * xmu * delta_barre_alpha_epsBH - beta * phidt * rdelta_sigma_barre_BH_Ratdt;
};
// -- retour du résidu
for (int i=1;i<=9;i++)
residu(i) = residuBH(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i));
return residu;
};
// calcul de la matrice tangente de la résolution de l'équation constitutive
// l'argument test ramène
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
// fatal, qui invalide le calcul du résidu et de la dérivée
Mat_abstraite& Hysteresis3D::Mat_tangente_constitutif
(const double & alphap,const Vecteur & x, Vecteur& residu, int& test)
{ // récupération du vecteur delta_sigma
for (int i=1;i<=9;i++) // accroissement de t à tdt
rdelta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = x(i);
test = 1; // init par défaut
rdelta_sigma_barre_BH_Ratdt = sigma_i_barre_BH + rdelta_sigma_barre_tdt_BH - sigma_barre_BH_R;
delta_barre_alpha_epsBH = alphap * delta_barre_epsBH; // limitation de la charge
// -- calcul de QdeltaSigmaRatdt
QdeltaSigmaRatdt = sqrt(rdelta_sigma_barre_BH_Ratdt.II());
// -- calculs relatifs aux angles de phases
if (avecVarWprimeS)
{Cal_wprime_et_der(QdeltaSigmaRatdt,rdelta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR
,wBase,rdelta_sigma_barre_tdt_BH,wprime,d_wprime,wprime_point,d_wprime_point,trajet_neutre);
}
else
{Cal_wprime(QdeltaSigmaRatdt,rdelta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR
,wBase,rdelta_sigma_barre_tdt_BH,wprime,wprime_point,trajet_neutre);
d_wprime.Inita(0.);
d_wprime_point.Inita(0.);
};
// --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
if (avec_phaseDeltaSig)
{xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,rdelta_sigma_barre_BH_Ratdt,xmu,Qzero);
};
//---- débug
// wprime = 1.;wprime_point=0.;trajet_neutre = false;
// --- fin débug
// -- modif 7 mai 2012 ---
// dans le cas où le phidt est négatif, l'équation constitutive n'est plus valide, dans ce cas on considère
// que phi est nulle, c-a-d que l'on évolue suivant un trajet neutre ce qui n'est sans doute pas l'idéal, mais permet d'éviter l'explosion
// calcul de phidt
double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH) ; // init partie radial // -- modif 16 juin 2012 --
if (!trajet_neutre) // -- modif 16 juin 2012 --
phidt -= QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime); // -- modif 16 juin 2012 --
// if ((phidt < 0) && force_phi_zero_si_negatif) // -- modif 16 juin 2012 --
if (phidt < 0) // -- modif 9 fev 2015 --
{if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 8)) || (Permet_affichage() > 8) )
cout << "\n phidt= " << phidt << "\n Hysteresis3D::Mat_tangente_constitutif(..." << endl;
if ( Dabs(phidt) < force_phi_zero_si_negatif) // -- modif 9 fev 2015 --
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
{ cout << "\n attention on utilise une valeur de phi negative dans pour le calcul de l'equation constitutive "
<< " phidt= " << phidt;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::Mat_tangente_constitutif(..." << endl;
};
phidt = 0.; trajet_neutre=true;
test = -1; // c'est un pb mais pas fatal
};
};
// fin modif modif 7 mai 2012 ---
if (trajet_neutre)
{ // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
// -- calcul du résidu
residuBH = rdelta_sigma_barre_tdt_BH - 2. * xmu * delta_barre_alpha_epsBH;
// variation du résidu et retour du résidu
for (int ig=1;ig<=9;ig++)
{int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
residu(ig) = residuBH(i,j);
for (int jg=1;jg<=9;jg++)
{int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
derResidu(ig,jg) = IdBH3(i,k)*IdBH3(l,j);
};
};
// --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
double phitau = 0.;
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
}
else
{ // cas d'un trajet non neutre
// modif 7 mai 2012 --- on ne recalcule plus le phi qui a été calculé plus haut pour savoir s'il était négatif
// // calcul de phidt
// double phidt = (rdelta_sigma_barre_BH_Ratdt && delta_barre_alpha_epsBH)
// - QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2. * xmu * wprime);
// --- si dépendance à la déformation équivalente
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phidt,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
// -- calcul de Beta
double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
double beta=0.; // init dans le cas où ne le calcul pas après
double unsurQdeltaSigma= ConstMath::grand; // on limite la valeur par exemple pour le démarrage
double puiss = 1.; // variable intermédiaire qui resert par la suite
if (QdeltaSigmaRatdt >= ConstMath::petit)
// on ne calcul beta que si la norme de deltaSigmaRatdt est significative
// idem pour unsurQdeltaSigma
{unsurQdeltaSigma=1./QdeltaSigmaRatdt;
if (xnp == 2.) { beta = - 2. * xmu * unsurwprimeQ0_puiss_np;}
else {puiss = pow(QdeltaSigmaRatdt,xnp-2.);
beta = - 2. * xmu * unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRatdt,xnp-2.);
};
};
// calcul du résidu
residuBH = rdelta_sigma_barre_tdt_BH
- 2. * xmu * delta_barre_alpha_epsBH - beta * phidt * rdelta_sigma_barre_BH_Ratdt;
// -- calcul de la variation du résidu
// 1 # variation de QdeltaSigmaRatdt ;
// -- calcul de la variation de QdeltaSigmaRatdt par rapport à sigma_barre
Tenseur3BH d_QdeltaSigmaRatdt(unsurQdeltaSigma * rdelta_sigma_barre_BH_Ratdt);
d_QdeltaSigmaRatdt.PermuteHautBas();
// 2 # variation de beta
Tenseur3BH d_beta; // initialisé par défaut à 0
if (QdeltaSigmaRatdt >= ConstMath::petit)
// dans le cas ou la norme de deltaSigmaRatdt n'est pas significative, beta est mis à zero
// il en est de même pour le cas où xnp = 2, pour lequel beta est constant ainsi que sa variation.
// A part ces deux cas on peut calculer d_beta
{if (xnp != 2.)
{double a = 2. * xmu * unsurwprimeQ0_puiss_np * puiss;
d_beta = (a * (2.-xnp) * unsurQdeltaSigma) * d_QdeltaSigmaRatdt
+ ( a *xnp / wprime) * d_wprime;
}; // si xnp==2: cas où d_beta est nul, déjà fait dans l'initialisation de d_beta
};
// 3 # variation de phidt
Tenseur3BH d_phidt(delta_barre_alpha_epsBH);
d_phidt.PermuteHautBas();
d_phidt -= (wprime_point*QdeltaSigmaRatdt/(xmu*wprime))*d_QdeltaSigmaRatdt;
d_phidt -= (QdeltaSigmaRatdt*QdeltaSigmaRatdt/2. * xmu)*
(d_wprime_point/wprime - (wprime_point/(wprime*wprime))*d_wprime);
// 4 # variation du résidu et retour du résidu
for (int ig=1;ig<=9;ig++)
{int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
residu(ig) = residuBH(i,j);
for (int jg=1;jg<=9;jg++)
{int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
derResidu(ig,jg) = (1.-beta * phidt)*IdBH3(i,k)*IdBH3(l,j)
- (d_beta(k,l) * phidt + beta * d_phidt(k,l))
* rdelta_sigma_barre_BH_Ratdt(i,j);
};
};
};// fin du cas non neutre
calcul_dResi_dsig = true; // on indique que la matrice est calculée
return derResidu;
};
// calcul de l'opérateur tangent : dsigma/depsilon
TenseurHHHH& Hysteresis3D::Dsig_depsilon
(bool en_base_orthonormee,const Tenseur3HH & gijHH_tdt, TenseurHHHH& dsig_deps_
,TenseurHH & sigHH)
{ // -- calcul de la matrice MPC
Dsig_d_(MPC);
// -- calcul de dsigma / d_eps
// comme le tenseur des contraintes est symétrique, idem pour le tenseur de déformation
// on doit avoir un tenseur d'ordre 4 symétrique par rapport aux deux premiers indices
// et symétrique par rapport aux deux derniers indices
// donc modif à faire -> on devrait itérer de 1 à 6 et non de 1 à 9 <<<<< !!!!!
// ===> non ! car en fait il y a de la non symétrie qui traîne... voir modif plus bas
Tenseur3HHHH & dsig_deps = *((Tenseur3HHHH*) &dsig_deps_);
// TenseurQ3geneHHHH & dsig_deps = *((TenseurQ3geneHHHH*) &dsig_deps_);
// dans le cas d'une base orthonormeee, la position des indices importe peu, donc l'opérateur
// déjà calculé est ok
if (en_base_orthonormee)
{for (int ig=1;ig<=9;ig++)
{int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
for (int jg=1;jg<=9;jg++)
{int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
dsig_deps.Change(i,j,k,l,-MPC(ig,jg));
};
};
}
else // dans le cas d'une base quelconque il faut intervenir sur le premier et dernier indice
{
// for (int i=1;i<4;i++)
// for (int j=1;j<4;j++)
// for (int k=1;k<4;k++)
// for (int l=1;l<4;l++)
// {double inter=0.;
// for (int e=1;e<4;e++)
// for (int f=1;f<4;f++)
// {int ej = Tenseur3BH::OdVect(e,j);
// int kf = Tenseur3HB::OdVect(k,f); // en fait idem mixte dans un sens ou l'autre
//// inter += gijHH_tdt(i,e) * MPC(ej,kf) * gijHH_tdt(f,l);
//// essai (a virer si pas concluant ) 1 oct 2020
//// a priori c'est concluant, cela fct nettement mieux, en particulier quand on l'utilise
//// en implicit : LoiAdditiveEnSigma::Calcul_DsigmaHH_via_eps_tdt
//// Je ne comprends pas pourquoi la matrice ne conduit pas à une symétrie, mais ...
//// forcé de constater qu'il faut intégrer tous les termes... et donc on force une symétrie !!
// int je = Tenseur3BH::OdVect(j,e);
// int fk = Tenseur3HB::OdVect(f,k); // en fait idem mixte dans un sens ou l'autre
// inter += gijHH_tdt(i,e) *
// 0.25* (MPC(ej,kf)+MPC(je,kf)+MPC(ej,fk)+MPC(je,fk))
// * gijHH_tdt(f,l);
////$$$$$ il manque les variation de gijHH_tdt par rapport à epsilon !!
//
//// fin essai: on garde car a priori concluant
// }
// dsig_deps.Change(i,j,k,l,-inter);
// };
////$$$$ à voir la prise en compte de la partie der de Jauman
//// essai en suivant exactement la même procédure que dsigma/d_ddl
// nouveau calcul car il manquait des termes
for (int oj=1;oj<=9;oj++)
{int o=Tenseur3BH::idx_i(oj); int j=Tenseur3BH::idx_j(oj);
for (int ef=1;ef<=9;ef++)
{int e=Tenseur3BH::idx_i(ef); int f=Tenseur3BH::idx_j(ef);
double inter=0.;
for (int i=1;i<4;i++)
for (int k=1;k<4;k++)
for (int l=1;l<4;l++)
{ int ij = Tenseur3BH::OdVect(i,j);
int kl = Tenseur3BH::OdVect(k,l);
// ici un - car MPC(ij,kl) est = - d sig_i^j / d eps_k^l
inter -= gijHH_tdt(o,i) * MPC(ij,kl) * IdBH3(k,e) * gijHH_tdt(f,l);
};
// ici - car c'est un - dans la formule
inter -= 2.* gijHH_tdt(o,e) * sigHH(f,j);
dsig_ojef_HHHH.Change(o,j,e,f,inter);
};
};
// on symétrise
dsig_deps.TransfertDunTenseurGeneral(dsig_ojef_HHHH.Symetrise1et2_3et4());
};
// retour
return dsig_deps;
};
// calcul de l'opérateur tangent : dsigma/d_ddl
void Hysteresis3D::Dsig_d_ddl(Tableau <TenseurHH *>& d_sigHH,TenseurBB & epsBB_tdt
,TenseurHH & gijHH_tdt,Tableau <TenseurBB *>& d_gijBB_tdt
,Tableau <TenseurHH *>& d_gijHH_tdt,Tableau <TenseurBB *>& d_epsBB )
{ // -- calcul de la matrice MPC
Dsig_d_(MPC);
////----- debug
//cout << "\n debug Hysteresis3D::Dsig_d_ddl( ";
//MPC.Affiche();
//cout << endl;
////--- fin debug
// -- calcul de dsigma / d_ddl
// nombre de ddl pour la variation du tenseur des contraintes
int nbddl = d_gijBB_tdt.Taille();
Tenseur3BH depsBH_dll;Tenseur3BH dsig3BH; Tenseur_ns3HH dsig1_HH;
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // " " " "
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
for (int i = 1; i<= nbddl; i++)
{ // on fait uniquement une égalité d'adresse de manière à ne pas utiliser
// le constructeur d'ou la profusion d'* et de ()
Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i))); // passage en dim 1
const Tenseur3BB & d_gijBB = *((Tenseur3BB*)(d_gijBB_tdt(i))); // passage en dim 1
const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); // "
depsBH_dll = (depsBB * gijHH + epsBB * dgijHH );
for (int ig=1;ig<=9;ig++)
{int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
dsig3BH.Coor(i,j) = 0.;
for (int jg=1;jg<=9;jg++)
{int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
// dsig3BH(i,j) += -(dRdeps(ig,jg) * depsBH_dll(k,l));
dsig3BH.Coor(i,j) += -(MPC(ig,jg) * depsBH_dll(k,l));
};
};
// et en deux fois contra
dsig1_HH = gijHH * dsig3BH + dgijHH * sigma_t_barre_tdt ;
// et on symétrise
dsigHH = 0.5*(dsig1_HH + dsig1_HH.Transpose());
}
};
// calcul de la première phase de l'opérateur tangent, est utilisé par Dsig_d_ddl et
// Dsig_depsilon
// calcul de la matrice MPC
void Hysteresis3D::Dsig_d_(MatLapack& MPC)
{ // dans le cas où la résolution est faite avec du Runge on appelle une routine spécifique
if (type_calcul_comportement_tangent==2)
{ Dsig_d_Runge(MPC);
return;
};
// ---- sinon il s'agit de l'algo avec linéarisation
double const untiers=1./3.;
double beta=0.; // init dans le cas où ne le calcul pas après
if (!trajet_neutre)
{ // si c'est un trajet non neutre, on utilise l'ensemble de la loi constitutive
// -- calcul de Beta
// remarque: le xnp que l'on utilise ici est le dernier calculé dans l'algo
// d'avancement temporel (on ne le met pas à jour)
double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
if (QdeltaSigmaRatdt >= ConstMath::petit)
// on ne calcul beta que si la norme de deltaSigmaRatdt est significative
{if (xnp == 2.)
{ beta = - 2. * xmu * unsurwprimeQ0_puiss_np ;}
else
{ beta = - 2. * xmu * unsurwprimeQ0_puiss_np * pow(QdeltaSigmaRatdt,xnp-2.);};
};
}; // fin du calcul de beta pour le cas d'un trajet non neutre
// normalement la matrice der_at_racine (qui vient d'être calculée
// pour résoudre l'équation constitutive, contient d_R_i^j / d_sigbarre_k^l)
// mais ce serait peut-être mieux de systématiquement la recalculer pour alpha = 1
// on aurait ainsi une matrice sécante ???
if (!calcul_dResi_dsig)
{// cas où la matrice tangente n'a pas été calculée
// calcul de la matrice tangente de la résolution de l'équation constitutive
// essai d'une nouvelle méthode
/* a prioiri ne fonctionne pas
// on utilise une petite portion de tàtdt
MatLapack der_at_racine_inter(9,9);
int nb_decoup = 10;
double alphap = 1./nb_decoup; //00; // pour voir 1/100
// sigma_i_barre_BH représente le point de ref sur la même courbe d'évolution
// delta_sigma_barre_tdt_BH représente le delta sigma depuis le point de ref
Tenseur3BH sigma_courant_BH(sigma_i_barre_BH);
Tenseur3BH sigma_i_barre_sauve_BH(sigma_i_barre_BH); // sauvegarde
Tenseur3BH incr_sig_BH(delta_sigma_barre_tdt_BH);incr_sig_BH /= nb_decoup;
for (int ig=1;ig<=9;ig++)
racine(ig) = incr_sig_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
for (int i=1;i<=nb_decoup;i++)
{ // on prend une nouvelle référence ad Hoc
sigma_i_barre_BH = sigma_courant_BH;
sigma_courant_BH += incr_sig_BH;
// on calcul l'opérateur tangent pour le sigma_courant_BH
der_at_racine_inter = Mat_tangente_constitutif(alphap,racine,residu);
// on cumule l'opérateur tangent
der_at_racine += der_at_racine_inter;
};
// calcul de l'opérateur tangent finale
der_at_racine *= alphap;
*/
// ici une solution qui fonctionne bien pour alpha = 10 !! incompréhensible !!
/* double alphap = 10;
// récup de la solution
for (int ig=1;ig<=9;ig++)
// racine(ig) = alphap * delta_sigma_barre_tdt_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
racine(ig) = delta_sigma_barre_tdt_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
// on prend une nouvelle référence, proche de l'arrivée
// sigma_i_barre_BH = sigma_t_barre_tdt-alphap * delta_sigma_barre_tdt_BH;
sigma_i_barre_BH = sigma_t_barre_tdt-delta_sigma_barre_tdt_BH;
der_at_racine = Mat_tangente_constitutif(alphap,racine,residu); */
// fin essai de la nouvelle méthode
// la méthode normale !
double alphap = 1.;
// récup de la solution
for (int ig=1;ig<=9;ig++)
racine(ig) = delta_sigma_barre_tdt_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
int test; // ne sert pas pour l'instant
der_at_racine = Mat_tangente_constitutif(alphap,racine,residu,test);
// !!!!!!!! pour le débug
//cout << " norme residu= " << residu.Max_val_abs()/delta_sigma_barre_tdt_BH.MaxiComposante() ;
//if (residu.Max_val_abs() <= ConstMath::petit)
// cout << " la norme est nulle !";
// !!!!!!!! pour le débug
};
// 1-- on calcul son inverse
der_at_racine.Inverse(dRdsigMoinsun);
// ------ pour le débug
// MatLapack produit(9,9);
// produit = dRdsigMoinsun * der_at_racine;
// for (int i=1;i<=9;i++)
// { if (Dabs(produit(i,i)-1. ) > ConstMath::unpeupetit)
// cout << "\n diag("<<i<<","<<i<<") n'est pas unitaire =" << produit(i,i);
// for (int j=1;j<= 9; j++)
// { if (i!=j)
// if (Dabs(produit(i,j)) > ConstMath::unpeupetit)
// cout << "\n terme ("<<i<<","<<j<<") non nul =" << produit(i,j);
// if ((produit(i,j) != produit(j,i)) && (Dabs(produit(i,j)) > ConstMath::unpeupetit))
// cout << "\n non symetrique ("<<i<<","<<j<<") =" << produit(i,j) << " et " << produit(j,i);
//
// }
// }
//
//------ fin pour le débug
// 2-- maintenant il faut calculer la variation du résidu constitutif / aux déformations
dRdeps.Initialise(0.);
// ************ important ********************
// en fait dans la suite on calcul - dSigma/ dEps et non +
// ensuite dans l'utilisation on met un moins (ce serait peut-être plus judicieux de mettre moins ici
// et pas après !!!
if (trajet_neutre)
{ // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
for (int ig=1;ig<=9;ig++)
{int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
for (int jg=1;jg<=9;jg++)
{int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
dRdeps(ig,jg) = -2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j));
};
};
}
else
{ // cas non neutre
for (int ig=1;ig<=9;ig++)
{int i=Tenseur3BH::idx_i(ig); int j=Tenseur3BH::idx_j(ig);
for (int jg=1;jg<=9;jg++)
{int k=Tenseur3BH::idx_i(jg); int l=Tenseur3BH::idx_j(jg);
dRdeps(ig,jg) = -2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j)) // originale
// dRdeps(ig,jg) = 2. * xmu * (IdBH3(i,k)*IdBH3(l,j)-untiers*IdBH3(l,k)*IdBH3(i,j)) // essai
// // + beta * rdelta_sigma_barre_BH_Ratdt(i,j)
// // * rdelta_sigma_barre_BH_Ratdt(l,k);
// + beta * delta_sigma_barre_BH_Ratdt(i,j) // originale
// normalement suivant les formules, on devrait avoir un - ???? incompréhensible
// car avec un + ça marche mieux ??
- beta * delta_sigma_barre_BH_Ratdt(i,j) // essai
* delta_sigma_barre_BH_Ratdt(l,k);
};
};
};
// 3-- puis on calcul le produit des deux matrices
MPC = dRdsigMoinsun * dRdeps; // originale
// MPC = -dRdeps * dRdsigMoinsun; // essai
};
// calcul de l'avancement temporel sur 1 pas,
// utilisé par les 3 programmes principaux:
// Calcul_SigmaHH, Calcul_DsigmaHH_tdt, Calcul_dsigma_deps,
void Hysteresis3D::Avancement_temporel(const Tenseur3BB & delta_epsBB,const Tenseur3HH & gijHH
,const Tenseur3BB & gijBB,int cas,SaveResulHysteresis3D & save_resul
,Tenseur3HH & sigHH,const EnergieMeca & energ_t,EnergieMeca & energ)
{
// . le tenseur des contraintes initiale en mixte, qui sera considéré comme la contrainte de début
// de calcul, pour toute la suite,
// . s'il y a plusieurs passage dans la boucle while qui suit, sigma_i_barre_BH variera, à chaque passage
// . on utilise le sigmaBH sauvegardé plutôt que le sigmaHH ce qui signifie que l'on effectue un transport
// en mixte plutôt qu'en 2 fois contravariants !!
// mise en place d'un transport éventuelle sur toutes les grandeurs
Transport_grandeur(gijBB,gijHH,save_resul);
sigma_i_barre_BH = save_resul.sigma_barre_BH_t; // récup des info
// même procédé pour la def équivalente de début de calcul, initié à t qui sert pour calculer la def équi de R à tdt
def_i_equi = save_resul.def_equi_at;
// init énergie
energ = energ_t;
double W_a = save_resul.fonction_aide_t; // valeur courante de la fonction d'aide
double delta_W_a = 0.; // valeur courante du delta W
// -- récup et initialisation des variables de travail
// iatens_princ: pointeur sur la liste déjà enregistrée, iatens: pointeur courant
List_io <Tenseur3BH>::iterator iatens_princ = save_resul.sigma_barre_BH_R.begin();
List_io <Tenseur3BH>::iterator iatens = iatens_princ;
// idem pour la fonction d'aide
List_io <double>::iterator iafct_princ = save_resul.fct_aide.begin();
List_io <double>::iterator iafct = iafct_princ;
aucun_pt_inversion = true; // pour gérer le fait qu'au début il n'y a pas de pt d'inversion !
if (iatens_princ != save_resul.sigma_barre_BH_R.end())
aucun_pt_inversion = false; // cas où on a déjà des inversions enregistrées
// idem pour les déformations équivalentes, fonctionnent avec l'indicateur "aucun_pt_inversion"
List_io <double>::iterator iadef_equi_princ = save_resul.def_equi.begin();
List_io <double>::iterator iadef_equi = iadef_equi_princ;
// idem pour les centres de surfaces de coïncidence
List_io <Tenseur3BH>::iterator iaOi_princ = save_resul.oc_BH_R.begin();
List_io <Tenseur3BH>::iterator iaOi = iaOi_princ;
// on récupère la def maxi de première charge, c'est celle qui correspond au premier point d'inversion
// si on est sur la première charge on prend la dernière valeur enregistrée sinon on récupère
// la valeur du premier point d'inversion
if (aucun_pt_inversion) { defEqui_maxi = def_i_equi;}
else { List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();
iil++;defEqui_maxi = *iil; // la dernière valeur est nulle, c'est l'avant dernière qui est la bonne
};
centre_initial = true; // pour gérer le fait qu'au début le centre initial = tenseur nul
if (iaOi_princ != save_resul.oc_BH_R.end())
centre_initial = false; // cas où on a déjà de nouveaux centres enregistrés
// un pointeur qui indique si les pointeurs courants sont sur les listes principales ou
// sur les nouvelles listes
bool pt_sur_principal = true; bool oi_sur_principal = true;
wBase = save_resul.wBase_t; // initialisation du paramètre de masing
double unSur_wBaseCarre= 1./(wBase*wBase);
// -- la déformation
delta_epsBH = delta_epsBB * gijHH; // incrément de def en mixte
double trace = delta_epsBH.Trace();
delta_barre_epsBH = delta_epsBH - (trace/3.)*IdBH3;
//-- partie spécifique gestion via fct aide
Tenseur3BH delta_epsBH_totale(delta_epsBH); // sauvegarde delta epsilon totale
Tenseur3BH delta_barre_epsBH_totale(delta_barre_epsBH); // sauvegarde delta_barre epsilon totale
// a priori, on fait un pas de calcul correspondant à la totalité de la déformation
// sauf cas inversion ou coïncidence
double ratio_deformation_en_cours = 0.; // utilisée pour l'instant uniquement avec l'algo fct aide
double ratio_multiplicatif_coincidence = 1.; // ""
int phase_coincidence = 0; // pour la gestion de la coincidence
//-- fin partie spécifique gestion via fct aide
// on définit un indicateur pour le type de gestion de la mémorisation discrète
bool bascule_fct_aide=false; // init
bool premier_passage = true; // init
// indicateur pour gérer la fin du traitement
bool fin_traitement= false; bool force_fin_traitement = false;
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\n$$$ === Hysteresis3D::Avancement_temporel: debut de l'algo " << endl ;
};
#endif
//----- fin débug
// ==== 1 === calcul de la contrainte à tdt =============
// faire tant que fin_traitement n'est pas bon
int nb_coin_inver=0; // compteur pour éviter une boucle infinie
while (!fin_traitement)
{// initialisation des variables de travail qui peuvent varier à chaque passage du while
// -- 1) sigma_barre_BH_R : initialisé sur la liste principal, mais ensuite peut évoluer sur liste secondaire
if (aucun_pt_inversion) {sigma_barre_BH_R.Inita(0.);def_equi_R=0.;}
else {sigma_barre_BH_R = *iatens; // la contrainte de référence
def_equi_R = *iadef_equi; }; // def équi au pt de référence
delta_sigma_barre_BH_Rat = sigma_i_barre_BH - sigma_barre_BH_R;
// -- 2) centre et rayon de la première surface de coïncidence
if (centre_initial) {oc_BH_R.Inita(0.);}
else {oc_BH_R = *iaOi;};
delta_sigma_barre_BH_OiaR = sigma_barre_BH_R-oc_BH_R; // le tenseur rayon
// -- 3) calcul du rayon: la première fois le rayon est infini
if (aucun_pt_inversion) {Q_OiaR = ConstMath::tresgrand;}
else {Q_OiaR = sqrt(delta_sigma_barre_BH_OiaR.II());};
// -- détermination du type de gestion de la mémorisation discrète au premier passage
// restera constant pour tout le reste du traitement
// le choix est donc fait par rapport aux valeurs de départ
if (premier_passage)
{Tenseur3BH delta_sigma_oi_a_t_BH = sigma_i_barre_BH - oc_BH_R; // OA_t
double rayon_t = sqrt(delta_sigma_oi_a_t_BH.II());
if (rayon_t > niveau_bascule_fct_aide * Qzero)
bascule_fct_aide = true;
premier_passage = false;
};
// le calcul de la contrainte s'effectue par la résolution de l'équation différentielle
// du schéma constitutif en BH, et toutes les grandeurs sont relatives à BH
// bien voir que que ce n'est qu'une partie de l'équation, et que pour l'intégration de
// la dérivée de Jauman, il faut aussi la partie HB qui est différente, mais on a une symétrie
// formelle (cf théorie)
// sigma_point = 2*mu*D_barre + beta*phi*(delta_barre de t à R de Sigma)
// ---7 mai 2012 ---
// point nouveau: cette équation n'est valide "que" si phi est positif !!
// en effet, si phi est négatif, le sigma_point au lieu de saturé, a tendance a exploser (c'est normal car la tendance est inverse)
// méthode permettant le calcul de sigma à tdt par différente méthodes: linéarisation
// ou kutta
// en sortie calcul de :
// sigma_t_barre_tdt, delta_sigma_barre_tdt_BH, delta_sigma_barre_BH_Ratdt
// également les grandeurs: QdeltaSigmaRatdt, et wprime,wprime_point, calculées à tdt
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << " \nsigR_sur_principal= " << pt_sur_principal << ", oi_sur_principal= " << oi_sur_principal
<< " save_resul.modif= " << save_resul.modif
<< " \n nb_Oi= "<<save_resul.oc_BH_R.size() << " nb_sigR= "<<save_resul.sigma_barre_BH_R.size()
<< " || nb_Oi_sec= " << save_resul.oc_BH_t_a_tdt.size()
<< " nb_sigR_sec= " << save_resul.sigma_barre_BH_R_t_a_tdt.size()
<< " || wBase= " << wBase;
};
#endif
//----- fin débug
double phi_2_dt = 0.; // on initialise avec un trajet neutre par défaut
CalculContrainte_tdt(phi_2_dt,save_resul.indicateurs_resolution);
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\nCalculContrainte_tdt(phi_2_dt...";
cout << " phi_2_dt = " << setprecision(ParaGlob::NbdigdoEC()) << phi_2_dt
<< " || wBase= " << wBase
<< endl ;
};
#endif
//----- fin débug
// 2==== gestion de la mémoire discrète
Tenseur3BH delta_sigma_oi_a_tdt_BH = sigma_t_barre_tdt - oc_BH_R; // OA
double rayon_tdt = sqrt(delta_sigma_oi_a_tdt_BH.II());
double d_W_a_essai = 0.; // init
if (wprime != 0.) // cas d'un trajet non neutre
{d_W_a_essai = 2. * phi_2_dt / (wprime * wprime);};
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\n rayon_tdt = " << setprecision(ParaGlob::NbdigdoEC()) << rayon_tdt
<< " Q_OiaR= " << setprecision(ParaGlob::NbdigdoEC()) << Q_OiaR << endl ;
};
#endif
//----- fin débug
if (!bascule_fct_aide) // cas d'une gestion à l'aide des rayons
{ if ((rayon_tdt < Q_OiaR)&&(!force_fin_traitement))
// cas où le nouveau rayon est inférieur au rayon de référence (= infini si courbe de première charge)
{int pb_calcul_centre_ok = 0; // init pour une gestion locale du cas où le calcul du centre cercle pose pb,
if (phi_2_dt < 0) //===========================================================================
// ce qui correspond également au cas d_W_a_essai < O
// a priori on a détecté un point d'inversion
// le calcul effectué n'est pas valide et on doit refaire une itération en initialisant les grandeurs ad hoc
// ------------------ cas d'une inversion ----------
{// -- on vérifie la cohérence de la valeur de la fonction d'aide (qui n'a rien à voir ici
// et donc doit être inférieur au maxi)
if (!aucun_pt_inversion) // si on a déjà des pt d'inversion, on vérifie les niveaux de la fonction d'aide
if (W_a > *iafct)
{if ( (ParaGlob::NiveauImpression() > 5) || (Permet_affichage() > 5))
{cout << "\n warning dans l'algorithme de gestion de la memoire discrete, "
<< " bien que l'on ait detecte un point d'inversion, il se trouve que la fonction d'aide "
<< " est superieur au dernier maxi, ce qui signifie que le calcul de l'integrale de la puissance "
<< " non reversible (phi(t)) n'a pas ete correcte (c'est une erreur cumulee), du sans doute a des pas de temps trop grand, "
<< " le calcul aura peut-etre des pb si on gere avec la fct d'aide, essayer de le reprendre avec un pas de temps plus petit ! "
<< " \n W_a= " << W_a << " iafct= " << *iafct
<< "\n Hysteresis3D::Avancement_temporel(...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
};
};
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\ninversion " << endl ;
};
#endif
//----- fin débug
// -- calcul des centre de cercle
if (!aucun_pt_inversion) // s'il y a au moins déjà un point d'inversion valide on peut traiter
{Tenseur3BH sig0i_BH;
pb_calcul_centre_ok=CentreCercle(sig0i_BH,oc_BH_R,delta_sigma_barre_BH_OiaR
,delta_sigma_barre_BH_Rat);
if (!pb_calcul_centre_ok) // si le calcul a été ok
{save_resul.oc_BH_t_a_tdt.push_front(sig0i_BH);
iaOi = save_resul.oc_BH_t_a_tdt.begin(); // init du centre courant pour la suite
oi_sur_principal = false; // va être mis à jour au début de la prochaine boucle
centre_initial = false; // du while (!fin_traitement)
};
}
else // si auparavant il n'y avait pas de pt d'inversion, maintenant il y en a un
{aucun_pt_inversion = false; // on le signal pour la suite
// point d'inversion
wBase = 2; unSur_wBaseCarre = 0.25; //=1./(wBase*wBase);
defEqui_maxi = def_i_equi;
};
// -- s'il n'y a pas eu de pb de calcul de centre on met tout à jour ==> fonctionnement normal
if (!pb_calcul_centre_ok)
{ // gestion du paramètre "modif" de saveresult
bool inversion = true;
Gestion_para_Saveresult_Modif(pt_sur_principal,save_resul,inversion);
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\npas de probleme de calcul de centre " ;
};
#endif
//----- fin débug
// calcul du nouveau centre d'inversion et sauvegarde, sauf dans le cas du premier point
// d'inversion, pour lequel le centre reste encore le tenseur nul pour la suite
// ajout de la contrainte de référence = sigma at
save_resul.sigma_barre_BH_R_t_a_tdt.push_front(sigma_i_barre_BH);
save_resul.fct_aide_t_a_tdt.push_front(W_a);
iafct = save_resul.fct_aide_t_a_tdt.begin();
W_a=0.; delta_W_a = 0.; // on ré-initialise pour les futures coincidences éventuelles
iatens = save_resul.sigma_barre_BH_R_t_a_tdt.begin(); // pour l'init du début du while générale
pt_sur_principal = false; // on indique que maintenant on pointe sur la nouvelle liste
// on indique le point d'inversion pour les futurs traitements
sigma_barre_BH_R =save_resul.sigma_barre_BH_R_t_a_tdt.front();
// -- idem pour la déformation équivalente
save_resul.def_equi_t_a_tdt.push_front(def_i_equi);
iadef_equi = save_resul.def_equi_t_a_tdt.begin(); // pour l'init du début du while générale
// on indique la déformation équivalente d'inversion pour les futurs traitements
def_equi_R = save_resul.def_equi_t_a_tdt.front();
// on met à jour la def de début de calcul (valant t au début, puis tau= pt d'inver ensuite ) pour def_equi_atdt
def_i_equi = def_equi_R;
// le calcul d'énergie se fera ensuite dans le cas normale ou avec coïncidence
// car la valeur actuelle de sigma est erroné, et phi est négatif
fin_traitement = false; // -- il faut recalculer avec le nouveau point d'inversion
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << " fin_traitement = " << fin_traitement
<< " save_resul.modif= " << save_resul.modif
<< endl ;
};
#endif
//----- fin débug
};
// sinon on est dans le cas où il y a eu un pb dans le calcul du centre --> ce qui invalide le point d'inverion
// on indique qu'il faut traiter la suite comme si il n'y avait pas de point d'inversion
// ceci se fait au travers de l'indicateur: pb_calcul_centre_ok
};//-- fin du cas de l'inversion
// --- maintenant on traite le cas courant ou il n'y a pas de point d'inversion ---
if ((phi_2_dt >= 0) || pb_calcul_centre_ok)
//======================================================================================
// c-a-d si phi_2_dt = 0 : on se promène sur un cercle
// si phi_2_dt > 0 : progression normale
// si pb_calcul_centre_ok != 0 : on est dans un cas limite pour lequel il est préférable de ne pas
// introduire de point d'inversion
{ // --- cas d'une évolution normale sans inversion ni coincidence ---
// mais cette évolution peut-être la fin après une coincidence par exemple, donc
// on ne modifie pas save_resul.modif
fin_traitement = true;
// dans l'expression suivante, on utilise une intégration incrémentale ...
// dans l'expression suivante, on utilise une intégration incrémentale ...
// s'il y a eu un pb de calcul de centre, cela veut dire que phi_2_dt était légèrement
// négatif à la prec près, il ne faut pas le compter
if (pb_calcul_centre_ok) {delta_W_a = 0;}
else if (wprime != 0.) // cas d'un trajet non neutre
{ delta_W_a = 2. * phi_2_dt / (wprime * wprime);}
else {delta_W_a = 0.;};
W_a += delta_W_a;
save_resul.fonction_aide_tdt = W_a;
// -- uptdate des énergies
UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R,def_i_equi,def_equi_atdt);
if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
else // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
// non c'est pas bon
/* defEqui_maxi = 0; // init
// pour la liste principale, la dernière valeur est nulle, c'est l'avant dernière qui est la bonne
if (save_resul.def_equi.size())
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
// pour la liste secondaire c'est directement la dernière valeur
if (save_resul.def_equi_t_a_tdt.size())
{ List_io <double>::reverse_iterator iil = save_resul.def_equi_t_a_tdt.rbegin();
defEqui_maxi = MaX(*iil,defEqui_maxi);
};*/
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\ncas courant sans point d'inversion " ;
cout << " fin_traitement = " << fin_traitement
<< " save_resul.modif= " << save_resul.modif
<< endl ;
};
#endif
//----- fin débug
};
} // -- fin du cas où le nouveau rayon est inférieur au rayon de référence
else if (!force_fin_traitement) // sinon cas d'une coincidence (rayon_tdt >= Q_OiaR )
{ // -- on calcul exactement le point de coïncidence
// auparavant on sauvegarde le delta eps pour connaître la partie concernée pour le
Tenseur3BH delta_eps_sur_deltaTau_BH = delta_barre_epsBH; // calcul d'énergie
double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas
// directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie sur une valeur = def_equi_R)
// mais c'est pour bien marquer le fait que à la sortie de la coïncidence, le pointeur iadef_equi de def_equi à changé
fin_traitement = Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB,save_resul
,W_a,iatens,iadef_equi
,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ,iadef_equi_princ
,iaOi_princ,iafct_princ,rayon_tdt,false);
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\ncas d'une coïncidence ";
cout << " fin_traitement = " << fin_traitement
<< " save_resul.modif= " << save_resul.modif
<< " save_resul.nb_coincidence= " << save_resul.nb_coincidence
<< endl ;
};
#endif
//----- fin débug
// dans le cas d'un trajet strictement radial, on a deux coïncidences à suivre
// mais on laisse l'algorithme passer deux fois, car sinon il y a des pb dans le cas
// de grand incrément
// -- uptdate des énergies
// récup du delta eps concerné: delta_barre_epsBH est modifié dans Coïncidence, maintenant
// il contient le reste de la déformation
delta_eps_sur_deltaTau_BH -= delta_barre_epsBH;
// on recalcule la valeur de phi car ici on est après le calcul d'une coïncidence, la contrainte qu'il faut considérer est
// celle de la coïncidence, la deformation à considérer est delta_eps_sur_deltaTau_BH, il faut donc calculer le phi qui
// correspond à ces deux grandeurs.
delta_sigma_barre_BH_Ra_i = sigma_i_barre_BH - sigma_barre_BH_R;
double QdeltaSigmaRai = sqrt(delta_sigma_barre_BH_Ra_i.II());
phi_2_dt = (delta_sigma_barre_BH_Ra_i && delta_eps_sur_deltaTau_BH)
- QdeltaSigmaRai * QdeltaSigmaRai * wprime_point / (2.*xmu * wprime);
// si le phi est négatif, on le met à 0 pour limiter les ennuis dans les calculs d'énergies (ici phi ne sert
// ne sert que là par la suite), mais c'est quand même pas normal !! mais ça peut être due à une légère erreur
// au niveau de la coïncidence
if (phi_2_dt < 0.)
phi_2_dt = 0.;
UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_eps_sur_deltaTau_BH
,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
else // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
// non c'est pas bon
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
// on met à jour la def de début de calcul (valant t au début )
// il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
// au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
def_i_equi = def_equi_R_ancien; // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
// NB: 1)on pourrait également prendre def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
// 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
}
else // cas où on force la fin du traitement: force_fin_traitement = true
{ fin_traitement = true;
// comme (rayon_tdt >= Q_OiaR ) on force la coincidence
double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas
// directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie sur une valeur = def_equi_R)
// mais c'est pour bien marquer le fait que à la sortie de la coïncidence, le pointeur iadef_equi de def_equi à changé
Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB,save_resul,W_a,iatens,iadef_equi
,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ,iadef_equi_princ
,iaOi_princ,iafct_princ,rayon_tdt,true);
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\non force la fin du traitement et la coincidence ";
cout << " fin_traitement = " << fin_traitement
<< " save_resul.modif= " << save_resul.modif
<< endl ;
};
#endif
//----- fin débug
// si le phi est négatif, on le met à 0 pour limiter les ennuis dans les calculs d'énergies (ici phi ne sert
// ne sert que là par la suite)
if (phi_2_dt < 0.)
phi_2_dt = 0.;
// -- uptdate des énergies
UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
else // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
// non c'est pas bon
/* defEqui_maxi = 0; // init
// pour la liste principale, la dernière valeur est nulle, c'est l'avant dernière qui est la bonne
if (save_resul.def_equi.size())
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
// pour la liste secondaire c'est directement la dernière valeur
if (save_resul.def_equi_t_a_tdt.size())
{ List_io <double>::reverse_iterator iil = save_resul.def_equi_t_a_tdt.rbegin();
defEqui_maxi = MaX(*iil,defEqui_maxi);
};*/
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;
};
// on met à jour la def de début de calcul (valant t au début )
// il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
// au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
def_i_equi = def_equi_R_ancien; // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
// NB: 1)on pourrait également prendre def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
// 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
};
}
// --- fin de la gestion par rapport au rayon
else // sinon cas d'une gestion avec la fonction d'aide
{ int pb_calcul_centre_ok = 0; // init pour une gestion locale du cas où le calcul du centre cercle pose pb,
if (((phase_coincidence==0) && (d_W_a_essai < 0.)) && (!force_fin_traitement))
{ // on a un point d'inversion
// -- on vérifie la cohérence de la valeur de la fonction d'aide
if (!aucun_pt_inversion) // si on a déjà des pt d'inversion, on vérifie les niveaux de la fonction d'aide
if (W_a > *iafct)
{cout << "\n erreur dans l'algorithme de gestion de la memoire discrete, "
<< " bien que l'on ait detecte un point d'inversion, il se trouve que la fonction d'aide "
<< " est superieur au dernier maxi, ce qui signifie que le calcul de l'integrale de la puissance "
<< " non reversible (phi(t)) n'a pas ete correcte (c'est une erreur cumulee), du sans doute a des pas de temps trop grand, "
<< " il n'est pas possible de continuer le calcul, essayer de le reprendre avec un pas de temps plus petit ! "
<< " \n W_a= " << W_a << " iafct= " << *iafct
<< "\n Hysteresis3D::Avancement_temporel(...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
Sortie(1);
};
// -- calcul des centre de cercle
if (!aucun_pt_inversion) // s'il y a au moins déjà un point d'inversion valide on peut traiter
{Tenseur3BH sig0i_BH;
pb_calcul_centre_ok=CentreCercle(sig0i_BH,oc_BH_R,delta_sigma_barre_BH_OiaR
,delta_sigma_barre_BH_Rat);
if (!pb_calcul_centre_ok) // si le calcul a été ok
{save_resul.oc_BH_t_a_tdt.push_front(sig0i_BH);
iaOi = save_resul.oc_BH_t_a_tdt.begin(); // init du centre courant pour la suite
oi_sur_principal = false; // va être mis à jour au début de la prochaine boucle
centre_initial = false; // du while (!fin_traitement)
};
}
else // si auparavant il n'y avait pas de pt d'inversion, maintenant il y en a un
{aucun_pt_inversion = false; // on le signal pour la suite
// point d'inversion
wBase = 2; unSur_wBaseCarre = 0.25; //=1./(wBase*wBase);
defEqui_maxi = def_i_equi;
};
// -- s'il n'y a pas eu de pb de calcul de centre on met tout à jour ==> fonctionnement normal
if (!pb_calcul_centre_ok)
{ // calcul du nouveau centre d'inversion et sauvegarde, sauf dans le cas du premier point
// d'inversion, pour lequel le centre reste encore le tenseur nul pour la suite
// gestion du paramètre "modif" de saveresult
bool inversion = true;
Gestion_para_Saveresult_Modif(pt_sur_principal,save_resul,inversion);
// ajout de la contrainte de référence = sigma at
save_resul.sigma_barre_BH_R_t_a_tdt.push_front(sigma_i_barre_BH);
// on ajoute le W_a et non le W_a_essai, car ce dernier n'est pas correcte
// par contre le W_a est de dernier maxi correcte
save_resul.fct_aide_t_a_tdt.push_front(W_a);
W_a=0.; delta_W_a = 0.; // on ré-initialise pour les futures coincidences éventuelles
iatens = save_resul.sigma_barre_BH_R_t_a_tdt.begin(); // pour l'init du début du while générale
iafct = save_resul.fct_aide_t_a_tdt.begin();
pt_sur_principal = false; // on indique que maintenant on pointe sur la nouvelle liste
// on indique le point d'inversion pour les futurs traitements
sigma_barre_BH_R =save_resul.sigma_barre_BH_R_t_a_tdt.front();
// -- idem pour la déformation équivalente
save_resul.def_equi_t_a_tdt.push_front(def_i_equi);
iadef_equi = save_resul.def_equi_t_a_tdt.begin(); // pour l'init du début du while générale
// on indique la déformation équivalente d'inversion pour les futurs traitements
def_equi_R = save_resul.def_equi_t_a_tdt.front();
// on met à jour la def de début de calcul (valant t au début, puis tau= pt d'inver ensuite ) pour def_equi_atdt
def_i_equi = def_equi_R;
// le calcul d'énergie se fera ensuite dans le cas normale ou avec coïncidence
// car la valeur actuelle de sigma est erroné, et phi est négatif
fin_traitement = false; // -- il faut recalculer avec le nouveau point d'inversion
};
// sinon on est dans le cas où il y a eu un pb dans le calcul du centre --> ce qui invalide le point d'inverion
// on indique qu'il faut traiter la suite comme si il n'y avait pas de point d'inversion
// ceci se fait au travers de l'indicateur: pb_calcul_centre_ok
};//-- fin du cas de l'inversion
//-------------------------------------------------
// arrivée ici on peu éventuellement calculer une valeur de la fonction d'aide correcte
// dans l'expression suivante, on utilise une intégration incrémentale ...
// s'il y a eu un pb de calcul de centre, cela veut dire que phi_2_dt était légèrement
// négatif à la prec près, il ne faut pas le compter
if (pb_calcul_centre_ok) {delta_W_a = 0;}
else { delta_W_a = d_W_a_essai;};
if (phase_coincidence != 1)
W_a += delta_W_a;
if (((phase_coincidence==1) || (delta_W_a > 0.)) || force_fin_traitement)
{if (W_a > *iafct) // on regarde si on est en coïncidence
phase_coincidence=1;
if (((phase_coincidence==0)&&(W_a <= *iafct)) && (!force_fin_traitement))// cas normal après plusieurs inversion
//======================================================================================
// c-a-d si phi_2_dt = 0 : on se promène sur un cercle
// si phi_2_dt > 0 : progression normale
// si pb_calcul_centre_ok != 0 : on est dans un cas limite pour lequel il est préférable de ne pas
// introduire de point d'inversion
{ // --- cas d'une évolution normale sans inversion ni coincidence ---
// mais cette évolution peut-être la fin après une coincidence par exemple, donc
// on ne modifie pas save_resul.modif
fin_traitement = true;
save_resul.fonction_aide_tdt = W_a;
// -- uptdate des énergies
UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R,def_i_equi,def_equi_atdt);
if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
else // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
// non c'est pas bon
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
}
else if ((phase_coincidence==1) && (!force_fin_traitement))
//======================================================================================
// cas où on est en coïncidence: on gère toutes les coïncidence en même temps c-a-d courbe secondaire
// ou retour sur la courbe principale
{ // -- on calcul exactement le point de coïncidence
// auparavant on sauvegarde le delta eps pour connaître la partie concernée pour le
Tenseur3BH delta_eps_sur_deltaTau_BH = delta_barre_epsBH; // calcul d'énergie
double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas
// directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie
// sur une valeur = def_equi_R)
// mais c'est pour bien marquer le fait que à la sortie de la coïncidence,
// le pointeur iadef_equi de def_equi à changé
double ratio_deformation_pour_coincidence=0;
bool force_coincidence = false; // init par défaut
if (phase_coincidence ==1) // si on suit une phase de coincidence qui vient d'être repéré, à la boucle précédente
force_coincidence = true; // on force au deuxième tour la coincidence
//--- debug
//-- fin debug
bool coincidence_exacte = Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB
,save_resul,W_a
,iatens,iadef_equi,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ
,iadef_equi_princ,ratio_deformation_pour_coincidence,iaOi_princ,iafct_princ,delta_W_a,force_coincidence);
// si la coïncidence est parfaite, dans ce cas seulement elle a été validée dans l'algo Coïncidence
// et on peut updater les énergies, sinon on refait un passage avec une déformation diminuée pour
// obtenir une coïncidence exacte
if (coincidence_exacte)
{ // -- cas d'une coïncidence validée
// dans le cas d'un trajet strictement radial, on a deux coïncidences à suivre
// mais on laisse l'algorithme passer deux fois, car sinon il y a des pb dans le cas
// de grand incrément
// -- uptdate des énergies
UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
else // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
// non c'est pas bon
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;};
// on met à jour la def de début de calcul (valant t au début )
// il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
// au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
def_i_equi = def_equi_R_ancien; // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
// NB: 1)on pourrait également prendre def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
// 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
// Mise en place de la déformation restante pour le calcul qui suit
delta_epsBH_totale -= delta_epsBH; delta_epsBH = delta_epsBH_totale;
delta_barre_epsBH_totale -= delta_barre_epsBH; delta_barre_epsBH = delta_barre_epsBH_totale;
// idem pour la contrainte de départ
sigma_i_barre_BH = sigma_t_barre_tdt;
// ratio_deformation_en_cours += ratio_multiplicatif_coincidence;
// ratio_multiplicatif_coincidence = 1.;
phase_coincidence = 0; // on réinitialise pour les autres passages
}
else
{ // -- cas non valide, on prépare pour une passe suivante
phase_coincidence = 1;
// ratio_multiplicatif_coincidence *= ratio_deformation_pour_coincidence;
// if (ratio_deformation_en_cours + ratio_multiplicatif_coincidence > 1.)
// // on ne peut pas dépasser la déformation imposée on limite donc
// ratio_multiplicatif_coincidence = 1.-ratio_deformation_en_cours;
delta_barre_epsBH *= ratio_deformation_pour_coincidence;
delta_epsBH *= ratio_deformation_pour_coincidence;
}
}
else // cas où on force la fin du traitement: force_fin_traitement = true
{ // ici comme on force une coïncidence exacte, il n'y a pas de calcul de point intermédiaire
// donc on reprend l'appel du cas des rayons, mais c'est équivalent à l'autre appel de Coincidence
// avec la fonction d'aide
fin_traitement = true;
// comme (rayon_tdt >= Q_OiaR ) on force la coincidence
double def_equi_R_ancien= def_equi_R; // sauvegarde qui ne sert à rien, car Coincidence ne change pas
// directement def_equi_R, seulement le pointeur iadef_equi (qui lui ne pointe plus en sortie sur une valeur = def_equi_R)
// mais c'est pour bien marquer le fait que à la sortie de la coïncidence, le pointeur iadef_equi de def_equi à changé
Coincidence(bascule_fct_aide,unSur_wBaseCarre,gijHH,gijBB,save_resul,W_a,iatens,iadef_equi
,iafct,pt_sur_principal,iaOi,oi_sur_principal,iatens_princ,iadef_equi_princ
,iaOi_princ,iafct_princ,rayon_tdt,true);
// si le phi est négatif, on le met à 0 pour limiter les ennuis dans les calculs d'énergies (ici phi ne sert
// ne sert que là par la suite)
if (phi_2_dt < 0.)
phi_2_dt = 0.;
// -- uptdate des énergies
UpdateEnergieHyste(energ,sigma_i_barre_BH,sigma_t_barre_tdt,delta_sigma_barre_tdt_BH
,delta_sigma_barre_BH_Ratdt,phi_2_dt,delta_barre_epsBH
,QdeltaSigmaRatdt,wprime,wprime_point,trajet_neutre,def_equi_R_ancien,def_i_equi,def_equi_atdt);
if (aucun_pt_inversion) { defEqui_maxi = def_equi_atdt;} // si sur la première charge on prend la valeur courante
else // on a des pt d'inversion, on prend le maxi des listes principales ou secondaires
// non c'est pas bon
{ List_io <double>::reverse_iterator iil = save_resul.def_equi.rbegin();iil++;defEqui_maxi = *iil;
};
// on met à jour la def de début de calcul (valant t au début )
// il s'agit d'une grandeur scalaire (pas tensorielle) on peut donc reprendre la valeur qu'elle avait
// au moment du précédent pt d'inversion, ce qui évitera de cumuler des erreurs d'intégration
def_i_equi = def_equi_R_ancien; // def_equi_R n'a pas encore été modifié, il correspond à la précédente valeur d'inversion
// NB: 1)on pourrait également prendre def_equi_atdt, mais qui cumulle les erreurs d'integ entre les deux derniers pt d'integ
// 2) c'est la même formule que dans le cas d'une inversion, mais cela ne donne pas du tout le même résultat
};
};
};
// calcul de la contrainte de retour si le traitement est fini
// et mise à jour a doc
if (fin_traitement)
{save_resul.sigma_barre_BH_tdt = sigma_t_barre_tdt;
save_resul.wBase_tdt = wBase;
save_resul.def_equi_atdt = def_equi_atdt;
switch (cas)
{case 1:
{ // cas de la dérivée de jauman
Tenseur_ns3HH sig1HH(gijHH * sigma_t_barre_tdt);
sigHH = 0.5*(sig1HH+sig1HH.Transpose()); // jauman
break;
}
case 2:
{ // pas de symétrisation du tenseur
sigHH = gijHH * sigma_t_barre_tdt;
break;
}
default:
cout << "\n erreur: cas= " << cas << " non prevu pour le type de derivee de la contrainte ";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::Avancement_temporel(...)" << endl;
Sortie(1);
};
//--------débug- traçage
#ifdef MISE_AU_POINT
if (Permet_affichage() == -77)
{ cout << "\n$$$ === fin de l'algo d'avancement : calcul de la contrainte " << endl ;
};
#endif
//----- fin débug
};
// on regarde si l'on n'a pas dépassé le nombre de boucle maxi permis: ce qui permet
// d'éviter les boucles infinies
nb_coin_inver++;
if (nb_coin_inver > Abs(nb_maxInvCoinSurUnPas))
{ // dans ce cas cela signifie qu'il y a pb et on arrête en générant éventuellement
// une interruption d'erreur de convergence dans une loi de comportement
if (((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
{cout << "\n warning dans le calcul de la contrainte a t+dt, le nombre de coincidence et/ou"
<< " d'inversion (" << nb_coin_inver << ") est superieur a la limite fixee: = "
<< nb_maxInvCoinSurUnPas ;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
};
// débug---------
// Verif_coincidence(delta_epsBB,gijHH,cas,save_resul,sigHH,energ_t,energ,true);
// Sortie(1);
// fin débug -------------
if (save_resul.nb_coincidence == 0)
{ // on est dans un cas particulier où l'algo fait du sur-place: un coup il trouve une coïncidence
// puis le coup après il trouve une inversion --> il revient au même endroit --> boucle infini
// dans ce cas on décide de manière arbitraire a imposer la coïncidence ce qui garantie la cohérence de la
// succession des rayons on contraire d'imposer l'inversion, qui conduira à un rayon incohérent (ce qui conduit
// normalement tout de suite après à une coïncidence)
force_fin_traitement = true;
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
cout << " on est dans le cas d'aller-retour identiques !! on impose la coincidence " << endl;
}
else // on regarde s'il faut imposer la valeur finale obtenue ou générer une erreur
{ if (nb_maxInvCoinSurUnPas > 0)
{ throw ErrNonConvergence_loiDeComportement();break;} // envoie d'un signal a traiter en dehors
else
{ force_fin_traitement = true;
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
cout << " on impose la coincidence " << endl;
}; // cas où on impose coûte que coûte (résultat pas sûr !!)
};
};
};// -- fin du While (!fin_traitement)
};
// initialisation éventuelle des variables thermo-dépendantes
void Hysteresis3D::Init_thermo_dependance()
{ // cas d'une thermo dépendance, on calcul les grandeurs en fonction de la température
if (xnp_temperature != NULL) xnp = xnp_lue = xnp_temperature->Valeur(*temperature);
if (Qzero_temperature != NULL) Qzero = Qzero_lue = Qzero_temperature->Valeur(*temperature);
if (xmu_temperature != NULL) xmu = xmu_lue = xmu_temperature->Valeur(*temperature);
//--- pour le debug
//cout << "\n np=" << xnp << " Q0=" << Qzero << " xmu=" << xmu << " T=" << *temperature;
//--- fin debug
};
// méthode permettant le calcul de sigma à tdt par différente méthodes: linéarisation
// ou kutta
// en sortie calcul de :
// sigma_t_barre_tdt, delta_sigma_barre_tdt_BH, delta_sigma_barre_BH_Ratdt
// également les grandeurs: QdeltaSigmaRatdt, et wprime,wprime_point, calculées à tdt
// ramène en retour la valeur de phi_2_dt
void Hysteresis3D::CalculContrainte_tdt(double & phi_2_dt,Tableau<double>& indicateurs_resolution)
{ // le calcul de la contrainte s'effectue par la résolution de l'équation différentielle
// du schéma constitutif en BH, et toutes les grandeurs sont relatives à BH
// bien voir que que ce n'est qu'une partie de l'équation, et que pour l'intégration de
// la dérivée de Jauman, il faut aussi la partie HB qui est différente, mais on a une symétrie
// formelle (cf théorie)
// sigma_point = 2*mu*D_barre + beta*phi*(delta_barre de t à R de Sigma)
if ((sortie_post)&&(indicateurs_resolution.Taille()!= 5)) // dimensionnement éventuelle de la sortie d'indicateurs
indicateurs_resolution.Change_taille(5);
calcul_dResi_dsig=false; // on indique que pour l'instant la matrice tangente n'est pas calculée
// comme la matrice n'est pas forcément définie positive on utilise GAUSS avec MatLapack
der_at_racine.Change_Choix_resolution(GAUSS,RIEN_PRECONDITIONNEMENT);
bool passage_par_newton = false; // permettra de savoir si on est passé par du newton dans le cas d'un basculement de méthode
//-- choix de la méthode
switch (type_resolution_equa_constitutive)
{case 1: // cas de la linéarisation de l'équation
{// ----- pour ce faire on appelle une methode de recherche de zero
val_initiale.Zero(); // on démarre la recherche à la valeur à t
racine.Zero(); // init du résultat à 0.
der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
// on met à jour éventuellement le tableau d'indicateurs
// 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
bool conver = false; // init par défaut
try // on met le bloc sous surveillance
{
conver=alg_zero.Newton_raphson
(*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
,maxi_delta_var_sig_sur_iter_pour_Newton);
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... ) //(ErrNonConvergence_Newton erreur)
{ double absracinemax=racine.Max_val_abs();
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 0))
{cout << "\n non convergence 2 sur l'algo de la resolution du schema constitutif"
<< " abs_racine_max " << absracinemax;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::CalculContrainte_tdt(.." << endl;
};
LibereTenseur();
LibereTenseurQ();
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;
////---- debug
// cout << "\n Hysteresis3D::CalculContrainte_tdt: indicateurs_resolution(1)= "<< indicateurs_resolution(1)
// << " indicateurs_resolution(2)=" << indicateurs_resolution(2)
// << " nb_incr_total "<< nb_incr_total << " nb_iter_total "<< nb_iter_total<< endl;
////---- fin debug
};
// if (!conver)
double absracinemax=racine.Max_val_abs();
if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 0))
{cout << "\n non convergence sur l'algo de la resolution du schema constitutif"
<< " abs_racine_max " << absracinemax;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::CalculContrainte_tdt(.." << endl;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
// récup de la solution
for (int ig=1;ig<=9;ig++)
delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH;
//---- debug
if ((ParaGlob::NiveauImpression() > 5)|| (Permet_affichage() > 5))
// on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
{ double sigma_tau_II = sigma_tauBH.II();
if ( sigma_tau_II > Qzero*Qzero)
cout << "\n sigma_tau_II= "<< sigma_tau_II << " Hysteresis3D::CalculContrainte_tdt" << endl;
};
//---- fin debug
// ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et
// pour t+dt donc pas de pb
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
passage_par_newton = true; // c'est vrai compte tenu de la méthode
//modif : 11 juin 2014
// on vérifie que l'amplitude de la contrainte calculée
{double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
{if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur10 fatale dans l'algo de la resolution du schema constitutif"
<< " concernant le niveau de saturation de la contrainte "
<< "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase);
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::CalculContrainte_tdt() (..." << endl;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
};
break;
}
case 2: // cas d'une résolution par intégration explicite par du kutta
// si cela ne converge pas on bascule sur de la linéarisation
{//initialisation des variables de calcul
for (int ig=1;ig<=9;ig++)
sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
double tdeb=0.,tfi=1.;
// calcul de la dérivée initiale
int erreur=0; //init d'une erreur de calcul de Sigma_point
Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
//// --- debug ---
//cout << "\n Hysteresis3D::CalculContrainte_tdt avant kutta"
// << "\n sig_initiale= ";sig_initiale.Affiche();
//cout << "\n dersig_initiale= "; dersig_initiale.Affiche();
//cout << endl; double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
//// --- fin debug ---
if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
{ if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur0 fatale dans l'algo de la resolution du schema constitutif"
<< " au niveau du calcul de la derivee initiale "
<< "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II())
<< " au lieu de "<< (Qzero*depassement_Q0*wBase);
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
cout << "\n Hysteresis3D::CalculContrainte_tdt() (..." << endl;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
int nombreAppelF=0,nb_step=0; // " "
double erreur_maxi_global=0.; // "
// appel de la fonction kutta
int conver=alg_edp.Pilotage_kutta
(cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale
,tdeb,tfi,erreurAbsolue,erreurRelative
,sig_finale,dersig_finale,dernierTemps,dernierdeltat
,nombreAppelF,nb_step,erreur_maxi_global);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
indicateurs_resolution(5)+=erreur_maxi_global;
};
// gestion de l'erreur de retour
// if ((conver != 2) || (sig_finale.Max_val_abs() > ConstMath::tresgrand))
// --- debug piloté ---
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt "
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
<< erreur_maxi_global << endl;
// --- fin debug piloté---
double abssigmax=sig_finale.Max_val_abs();
if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
{// dans le cas où le Runge ne marche pas on essaie avec la première méthode
// linéarisation de l'équation
// ----- pour ce faire on appelle une methode de recherche de zero
val_initiale.Zero(); // on démarre la recherche à la valeur à t
racine.Zero(); // init du résultat à 0.
der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
// 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
bool converge=alg_zero.Newton_raphson
(*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
,maxi_delta_var_sig_sur_iter_pour_Newton);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;};
double absracinemax=racine.Max_val_abs();
if ((converge) && (isfinite(absracinemax)) && (!isnan(absracinemax)))
{ // on a convergé, donc on termine le calcul et on sort du switch
// récup de la solution
for (int ig=1;ig<=9;ig++)
delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH;
// ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et
// pour t+dt donc pas de pb
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
passage_par_newton = true; // on indique donc que l'on est passé par du newton
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
{ cout << "\n probleme dans la resolution de l'equation constitutive avec Runge Kutta";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
if (conver != 2)
alg_edp.Affiche_Explication_erreur_RG(conver);
if (!isfinite(abssigmax))
cout << " abssigmax n'est pas fini = " << abssigmax ;
if (isnan(abssigmax))
cout << " abssigmax = nan = " << abssigmax ;
if (Permet_affichage() > 2)
{ cout << "\n abs_racine_max " << absracinemax;}
if (Permet_affichage() > 4)
{cout << " (2) utilisation avec succes de la linearisation ";};
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt (..."<<endl;
};
// puis on met conver à 2 pour signaler que pour la suite que ça a convergé
conver=2;
break;
}
else // sinon on n'a pas convergé également avec la linéarisation
// on indique à toute fin utile que la première n'a pas convergé au cas ou cela n'aurait pas été noté
{ if (conver == 2)
conver = 0 ;
};
};
// arrivé ici cela signifie que: soit kutta a convergé -> conver=2,
// soit kutta n'a pas convergé et la méthode de linéarisation n'a pas également convergé -> conver != 2
// --- si on n'a pas convergé:
if (conver !=2) // cas d'une non convergence
{ // on appel du kutta adoc sans gestion d'erreur !!
double deltat=tfi-tdeb;
switch (cas_kutta)
{ case 3:
{ alg_edp.Runge_Kutta_step23(*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale,tdeb,deltat,sig_finale,estime_erreur);
break;
}
case 4:
{ alg_edp.Runge_Kutta_step34(*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale,tdeb,deltat,sig_finale,estime_erreur);
break;
}
case 5:
{ alg_edp.Runge_Kutta_step45(*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale,tdeb,deltat,sig_finale,estime_erreur);
break;
}
default:
cout << "\n erreur , cas Runge kutta non prevu ! "
<< "\n Hysteresis3D::CalculContrainte_tdt(..";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();
Sortie(1);
};
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
{cout << "\n probleme dans la resolution de l'equation constitutive avec Runge Kutta";
alg_edp.Affiche_Explication_erreur_RG(conver);
cout << ", appel direct de kutta-> erreur estimee= " << estime_erreur(1)
<< "\n Hysteresis3D::CalculContrainte_tdt (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
// Sortie(1);
// pour le debug
// if ((Dabs(estime_erreur(1)) > 10) || (Dabs(sig_finale(1)) > 10))
// alg_edp.Runge_Kutta_step45(*this,& Hysteresis3D::Sigma_point,sig_initiale,dersig_initiale
// ,tdeb,deltat,sig_finale,estime_erreur);
};
double abssigmax=sig_finale.Max_val_abs();
if ((!isfinite(abssigmax)) || (isnan(abssigmax)) )
// dans le cas où la valeur de sig_finale est infinie, ce n'est pas la peine d'aller plus loin
{
cout << "\n erreur , l'appel directe de Kutta donne une valeur infinie de contrainte ! "
<< "\n Hysteresis3D::CalculContrainte_tdt(..";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
// on génère une exception
throw ErrNonConvergence_loiDeComportement();
Sortie(1);
};
}
// -- arrivée ici, cela veut dire que: 1) soit le RG de départ a marché correctement,
// 2) soit la linéarisation a marché correctement, 3) soit si les deux autres n'ont pas marché,
// on a eu un appel directe à un Rg d'un ordre donné, avec un résultat non contrôlé
else if (!passage_par_newton)// sinon cela veut dire que c'est le résultat avec Runge que l'on exploite
{// récup des résultats
for (int ig=1;ig<=9;ig++)
sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
delta_sigma_barre_tdt_BH = sigma_t_barre_tdt - sigma_i_barre_BH;
// ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
// d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
};
// le dernier cas possible avec convergence est celui du passage par newton, et dans ce cas
// la sauvegarde a déjà été faite
//modif : 11 juin 2014
// on vérifie que l'amplitude de la contrainte calculée
{double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
{if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur11 fatale dans l'algo de la resolution du schema constitutif"
<< " concernant le niveau de saturation de la contrainte "
<< "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
};
break;
}
case 3: // cas d'une résolution par intégration explicite par du kutta
// mais ici si ça ne converge pas on génère une exception c'est tout , on ne cherche pas
// a utiliser une autre méthode de secours
{//initialisation des variables de calcul
for (int ig=1;ig<=9;ig++)
sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
double tdeb=0.,tfi=1.;
int erreur=0; //init d'une erreur de calcul de Sigma_point
// calcul de la dérivée initiale
Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
{ if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur1 fatale dans l'algo de la resolution du schema constitutif"
<< " au niveau du calcul de la derivee initiale "
<< "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "<< Qzero
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
int nombreAppelF=0,nb_step=0; // " "
double erreur_maxi_global=0.; // "
// appel de la fonction kutta
int conver=alg_edp.Pilotage_kutta
(cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale
,tdeb,tfi,erreurAbsolue,erreurRelative
,sig_finale,dersig_finale,dernierTemps,dernierdeltat
,nombreAppelF,nb_step,erreur_maxi_global);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
indicateurs_resolution(5)+=erreur_maxi_global;
};
// --- debug piloté ---
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt "
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
<< erreur_maxi_global << endl;
// --- fin debug piloté---
// gestion de l'erreur de retour
// if ((conver != 2) || (sig_finale.Max_val_abs() > ConstMath::tresgrand))
double abssigmax=sig_finale.Max_val_abs();
if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 0))
{cout << "\n non convergence sur l'algo de la resolution du schema constitutif"
<< "\n Hysteresis3D::CalculContrainte_tdt(.." ;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
// récup des résultats
for (int ig=1;ig<=9;ig++)
sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
delta_sigma_barre_tdt_BH = sigma_t_barre_tdt - sigma_i_barre_BH;
// ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
// d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
//modif : 11 juin 2014
// on vérifie que l'amplitude de la contrainte calculée
{if ( QdeltaSigmaRatdt > (Qzero*depassement_Q0*wBase))
{if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 0))
{cout << "\n erreur12 fatale dans l'algo de la resolution du schema constitutif"
<< " concernant le niveau de saturation de la contrainte "
<< "\n QdeltaSigma =" << QdeltaSigmaRatdt << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
};
break;
}
case 4: // cas d'une résolution par intégration explicite par du kutta
{//initialisation des variables de calcul
for (int ig=1;ig<=9;ig++)
sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
double tdeb=0.,tfi=1.;
// calcul de la dérivée initiale
int erreur=0; //init d'une erreur de calcul de Sigma_point
Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 0))
{cout << "\n erreur2 fatale dans l'algo de la resolution du schema constitutif"
<< " au niveau du calcul de la derivee initiale "
<< "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
int nombreAppelF=0,nb_step=0; // " "
double erreur_maxi_global=0.; // "
// appel de la fonction kutta
int conver=alg_edp.Pilotage_kutta
(cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale
,tdeb,tfi,erreurAbsolue,erreurRelative
,sig_finale,dersig_finale,dernierTemps,dernierdeltat
,nombreAppelF,nb_step,erreur_maxi_global);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
indicateurs_resolution(5)+=erreur_maxi_global;
};
// --- debug piloté ---
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
{cout << "\n Hysteresis3D::CalculContrainte_tdt "
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
<< erreur_maxi_global;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature(); cout << flush;
};
// --- fin debug piloté---
// gestion de l'erreur de retour
// if ((conver != 2) || (sig_finale.Max_val_abs() > ConstMath::tresgrand))
double abssigmax=sig_finale.Max_val_abs();
if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
{// dans le cas où le Runge ne marche pas on essaie avec la première méthode
// linéarisation de l'équation
// ----- pour ce faire on appelle une methode de recherche de zero
val_initiale.Zero(); // on démarre la recherche à la valeur à t
racine.Zero(); // init du résultat à 0.
der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
// 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
bool converge=alg_zero.Newton_raphson
(*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
,maxi_delta_var_sig_sur_iter_pour_Newton);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;};
double absracinemax=racine.Max_val_abs();
if ((converge) && (isfinite(absracinemax)) && (!isnan(absracinemax)))
{ // on a convergé, donc on termine le calcul et on sort du switch
// récup de la solution
for (int ig=1;ig<=9;ig++)
delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH;
// ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et
// pour t+dt donc pas de pb
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
passage_par_newton = true; // on indique donc que l'on est passé par du newton
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 2)) || (Permet_affichage() > 2))
{ if (Permet_affichage() > 4)
{cout << " -> probleme dans la resolution de l'equation constitutive avec Runge Kutta";
alg_edp.Affiche_Explication_erreur_RG(conver);
cout << ", (1) utilisation avec succes de la linearisation ";
cout << "\n abs_racine_max " << absracinemax;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt (..."<<endl;
};
// puis on met conver à 2 pour signaler que pour la suite que ça a convergé
conver=2;
break;
}
else
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 2)) || (Permet_affichage() > 2))
{cout << "\n non convergence sur l'algo de la resolution du schema constitutif"
<< " abs_racine_max " << absracinemax << " conver= " << conver << " abssigmax= " << abssigmax;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 7)) || (Permet_affichage() > 7))
cout << "\n Hysteresis3D::CalculContrainte_tdt(.." << endl;
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
}
else // sinon cela veut dire que le calcul avec Runge est correcte on sauvegarde les résutats
{// récup des résultats
for (int ig=1;ig<=9;ig++)
sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
delta_sigma_barre_tdt_BH = sigma_t_barre_tdt - sigma_i_barre_BH;
// ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
// d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
};
//modif : 11 juin 2014
// on vérifie que l'amplitude de la contrainte calculée
{double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
{if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur13 fatale dans l'algo de la resolution du schema constitutif"
<< " concernant le niveau de saturation de la contrainte "
<< "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< "\n Hysteresis3D::CalculContrainte_tdt() (..." ;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
};
break;
}
case 5: case 6:// cas d'une résolution par intégration explicite par du kutta
// l'idée est ici de diminuer la précision jusqu'à un résultat si on ne converge pas normalement
{//initialisation des variables de calcul
for (int ig=1;ig<=9;ig++)
sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
double tdeb=0.,tfi=1.;
// calcul de la dérivée initiale
int erreur=0; //init d'une erreur de calcul de Sigma_point
Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
{ if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur3 fatale dans l'algo de la resolution du schema constitutif"
<< " au niveau du calcul de la derivee initiale "
<< "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "
<< (Qzero*depassement_Q0*wBase)
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
int nombreAppelF=0,nb_step=0; // " "
double erreur_maxi_global=0.; // "
// appel de la fonction kutta
int conver=alg_edp.Pilotage_kutta
(cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale
,tdeb,tfi,erreurAbsolue,erreurRelative
,sig_finale,dersig_finale,dernierTemps,dernierdeltat
,nombreAppelF,nb_step,erreur_maxi_global);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
indicateurs_resolution(5)+=erreur_maxi_global;
};
// --- debug piloté ---
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt "
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
<< erreur_maxi_global << endl;
// --- fin debug piloté---
// gestion de l'erreur de retour
// if ((conver != 2) || (sig_finale.Max_val_abs() > ConstMath::tresgrand))
double abssigmax=sig_finale.Max_val_abs();
if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
{// dans le cas où le Runge ne marche pas on essaie avec la première méthode
// linéarisation de l'équation
// ----- pour ce faire on appelle une methode de recherche de zero
val_initiale.Zero(); // on démarre la recherche à la valeur à t
racine.Zero(); // init du résultat à 0.
der_at_racine.Initialise(0.); // init de la matrice dérivée à 0.
int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
// 1== résolution de l'équation constitutive d'avancement discrétisée en euler implicite
bool converge=alg_zero.Newton_raphson
(*this,&Hysteresis3D::Residu_constitutif,&Hysteresis3D::Mat_tangente_constitutif
,val_initiale,racine,der_at_racine,nb_incr_total,nb_iter_total
,maxi_delta_var_sig_sur_iter_pour_Newton);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(1)+=nb_incr_total;indicateurs_resolution(2)+=nb_iter_total;};
// if (converge)
double absracinemax=racine.Max_val_abs();
if ((converge) && (isfinite(absracinemax)) && (!isnan(absracinemax)))
{ // on a convergé, donc on termine le calcul et on sort du switch
// récup de la solution
for (int ig=1;ig<=9;ig++)
delta_sigma_barre_tdt_BH.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=racine(ig);
sigma_t_barre_tdt = sigma_i_barre_BH + delta_sigma_barre_tdt_BH;
// ici il n'est pas nécessaire de calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs sont calculées à l'intérieur de la boucle de Newton, et
// pour t+dt donc pas de pb
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
passage_par_newton = true; // on indique donc que l'on est passé par du newton
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 2)) || (Permet_affichage() > 2))
{ cout << "\n estime_erreur= " << estime_erreur(1);
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
if (Permet_affichage() > 4)
{cout << " -> probleme dans la resolution de l'equation constitutive avec Runge Kutta";
alg_edp.Affiche_Explication_erreur_RG(conver);
cout << ", (3) utilisation avec succes de la linearisation ";
};
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt (..."<<endl;
};
// puis on met conver à 2 pour signaler que pour la suite que ça a convergé
conver=2;
break;
}
else // sinon on n'a pas convergé également avec la linéarisation
// on indique à toute fin utile que la première n'a pas convergé au cas ou cela n'aurait pas été noté
{ if (conver == 2)
conver = 0 ;
};
};
// arrivé ici cela signifie que: soit kutta a convergé -> conver=2,
// soit kutta n'a pas convergé et la méthode de linéarisation n'a pas également convergé -> conver != 2
if (conver !=2) // cas d'une non convergence
{ // l'idée est de diminuer la précision jusqu'à une solution
int conver_interne = 0;
double erreurAsolue_inter=erreurAbsolue;
double erreurRelative_inter = erreurRelative;
int nb_boucle=1;
double abssigmax0; // init
do
{ erreurAsolue_inter *=10;
erreurRelative_inter *=10;
nb_boucle ++;
for (int ig=1;ig<=9;ig++)
sig_initiale(ig)=sigma_i_barre_BH(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig));
double tdeb=0.,tfi=1.;
// calcul de la dérivée initiale
int erreur=0; //init d'une erreur de calcul de Sigma_point
Sigma_point(tdeb,sig_initiale,dersig_initiale,erreur);
if (erreur) // cas où le calcul de la dérivée initiale n'est pas possible, on ne peut pas aller plus loin
{ if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur4 fatale dans l'algo de la resolution du schema constitutif"
<< " au niveau du calcul de la derivee initiale "
<< "\n QdeltaSigma =" << sqrt(delta_sigma_barre_R_a_tauBH.II()) << " au lieu de "<< Qzero
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
double dernierTemps=0.,dernierdeltat=0.; // valeurs de retour
int nombreAppelF=0,nb_step=0; // " "
erreur_maxi_global=0.;
// appel de la fonction kutta
conver_interne=alg_edp.Pilotage_kutta
(cas_kutta,*this,& Hysteresis3D::Sigma_point,& Hysteresis3D::Verif_integrite_Sigma
,sig_initiale,dersig_initiale
,tdeb,tfi,erreurAsolue_inter,erreurRelative_inter
,sig_finale,dersig_finale,dernierTemps,dernierdeltat
,nombreAppelF,nb_step,erreur_maxi_global);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{indicateurs_resolution(3)+=nombreAppelF;indicateurs_resolution(4)+=nb_step;
indicateurs_resolution(5)+=erreur_maxi_global;
};
// --- debug piloté ---
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n Hysteresis3D::CalculContrainte_tdt "
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
<< erreur_maxi_global << endl;
// --- fin debug piloté---
if (nb_boucle > 10) // si en * par 10^10 on n'a pas réussi on arrête
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 0))
{cout << "\n non convergence extreme !! sur l'algo de la resolution du schema constitutif"
<< "\n Hysteresis3D::CalculContrainte_tdt(.." ;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
abssigmax=sig_finale.Max_val_abs();
}
while ((conver_interne != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) );
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
{cout << "\n probleme dans la resolution de l'equation constitutive avec Runge Kutta";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
alg_edp.Affiche_Explication_erreur_RG(conver);
cout << ", decroissance de la precision -> erreur estimee= "
<< erreur_maxi_global << " nb cycle= " << nb_boucle
<< "\n Hysteresis3D::CalculContrainte_tdt (..."<< endl;
};
}
else if (!passage_par_newton)// sinon cela veut dire que le calcul avec Runge est correcte on sauvegarde les résutats
{// récup des résultats
for (int ig=1;ig<=9;ig++)
sigma_t_barre_tdt.Coor(Tenseur3BH::idx_i(ig),Tenseur3BH::idx_j(ig))=sig_finale(ig);
delta_sigma_barre_tdt_BH = sigma_t_barre_tdt - sigma_i_barre_BH;
// ici on doit calculer QdeltaSigmaRatdt, et wprime,wprime_point
// car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
// d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
delta_sigma_barre_BH_Ratdt = sigma_t_barre_tdt - sigma_barre_BH_R;
QdeltaSigmaRatdt = sqrt(delta_sigma_barre_BH_Ratdt.II());
};
// le dernier cas possible avec convergence est celui du passage par newton, et dans ce cas
// la sauvegarde a déjà été faite
//modif : 11 juin 2014
// on vérifie que l'amplitude de la contrainte calculée
{double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
{if ((ParaGlob::NiveauImpression() > 3)|| (Permet_affichage() > 0))
{cout << "\n erreur14 fatale dans l'algo de la resolution du schema constitutif"
<< " concernant le niveau de saturation de la contrainte "
<< "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< "\n Hysteresis3D::CalculContrainte_tdt() (...";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();break;
};
};
break;
}
};
//--- // fin du switch sur les différentes méthodes pour calculer la contraintes à tdt
// --- on calcul certaines grandeurs globales et locales
//------ debug------
//cout << "\n dd " << delta_sigma_barre_BH_Ratdt.MaxiComposante() << ", "
// << " fini ? " << isfinite(delta_sigma_barre_BH_Ratdt.MaxiComposante())
// << " isnan . " << isnan(delta_sigma_barre_BH_Ratdt.MaxiComposante());
//cout << "\n delta_sigma_barre_BH_Ratdt=== ";
//delta_sigma_barre_BH_Ratdt.Ecriture(cout);
//if ((!isfinite(delta_sigma_barre_BH_Ratdt.MaxiComposante()))||(isnan(delta_sigma_barre_BH_Ratdt.MaxiComposante())))
// { cout << " delta_sigma_barre_BH_Ratdt= ";
// delta_sigma_barre_BH_Ratdt.Ecriture(cout);
// };
//----- fin debug -----
if ((type_resolution_equa_constitutive != 1)&& (!passage_par_newton))
{ // sauf le cas 1, tous les autres c'est des kutta on doit calculer wprime,wprime_point
// car ces grandeurs ne sont pas calculées au temps final !! dans la méthode de Runge-Kutta
// d'ailleurs elles ne sont même pas calculées en temps que grandeurs globales !!
Cal_wprime(QdeltaSigmaRatdt,delta_sigma_barre_BH_Ratdt,Q_OiaR,delta_sigma_barre_BH_OiaR
,wBase,delta_sigma_barre_tdt_BH,wprime,wprime_point,trajet_neutre);
// --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
if (avec_phaseDeltaSig)
{xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,delta_sigma_barre_BH_Ratdt,xmu,Qzero);
};
};
// --<< calcul de phi >>---
if (trajet_neutre) // cas où l'on n'a pas un trajet neutre
{ phi_2_dt = 0.;} // pas de variation sur un trajet neutre
else
{ phi_2_dt = (delta_sigma_barre_BH_Ratdt && delta_barre_epsBH)
- QdeltaSigmaRatdt * QdeltaSigmaRatdt * wprime_point / (2.*xmu * wprime);
//------ debug------
//if ((!isfinite(phi_2_dt))||(isnan(phi_2_dt)))
// { cout << "\n";
// delta_sigma_barre_BH_Ratdt.Ecriture(cout);
// delta_barre_epsBH.Ecriture(cout);
// cout << " q= " << QdeltaSigmaRatdt;
//// Sortie(1);
// };
//// if (phi_2_dt < 0.) cout << "\n phi_2_dt = " << phi_2_dt;
// cout << " " << phi_2_dt;
//----- fin debug -----
};
// --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
ParaMateriauxAvecDefEqui(xnp,QdeltaSigmaRatdt,phi_2_dt,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
};
// calcul de l'expression permettant d'obtenir la dérivée temporelle de la contrainte
// en fait il s'agit de l'équation constitutive
// utilisée dans la résolution explicite (runge par exemple) de l'équation constitutive
// erreur : indique si le calcul est valide =0 dans ce cas
// sinon (diff de 0) indique que les valeurs calculées ne sont pas licite
// =1: la norme de sigma est supérieure à la valeur limite de saturation
Vecteur& Hysteresis3D::Sigma_point(const double & tau, const Vecteur & sigma_tau, Vecteur& sig_point,int& erreur)
{
double titi= tau; // sert à rien, c'est pour taire le compilo car tau ne sert pas directement
// récup de la contrainte à tau
for (int i=1;i<=9;i++)
sigma_tauBH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = sigma_tau(i);
//modif : 1 juin 2014
// non on fait la vérif sur le delta cf. plus bas
// // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
// double sigma_tau_II = sigma_tauBH.II();
// if ( sigma_tau_II > Qzero*Qzero)
// {erreur = 1;
// //--- pour le debug
// if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
// cout << "\n sigma_tau_II =" << sigma_tau_II << " au lieu de "<< Qzero*Qzero
// << " Hysteresis3D::Sigma_point " << endl ;
// //--- fin debug
// }
// else // sinon on continue
// { erreur = 0;
//--- pour le debug
//cout << "\n tt np=" << xnp << " Q0=" << Qzero << " xmu=" << xmu << " T=" << *temperature;
//--- fin debug
// variation de sigma de R à tau
delta_sigma_barre_R_a_tauBH = sigma_tauBH - sigma_barre_BH_R; // ne sert que dans Sigma_point
delta_sigma_barre_t_a_tauBH = sigma_tauBH - sigma_i_barre_BH; // ne sert que dans Sigma_point
// calcul de QdeltaSigma: on prend la valeur absolue car dans les petits nombres on peut avoir des presques 0 !!
// négatifs qui conduisent à une racine infinie !!
double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
//modif : 1 juin 2014
// on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
// on fait la vérification sur le delta, car si on est sur une branche non initiale, on peut très bien
// dépasser le Qzero, localement, mais ce sera invalidé ensuite avec l'algo de coincidence
if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
{erreur = 1;
//--- pour le debug
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< " Hysteresis3D::Sigma_point " << endl ;
//--- fin debug
}
else // sinon on continue
{ erreur = 0;
// -- calculs relatifs aux angles de phases : calcul de -> wprime,wprime_point,trajet_neutre
Cal_wprime(QdeltaSigma,delta_sigma_barre_R_a_tauBH,Q_OiaR,delta_sigma_barre_BH_OiaR
,wBase,delta_sigma_barre_t_a_tauBH,wprime,wprime_point,trajet_neutre);
// --- si dépendance de la phase de sigma0i_tdt (et non delta sigma) mise à jour des para matériaux
if (avec_phaseDeltaSig)
{xnp = xnp_lue; xmu = xmu_lue; Qzero = Qzero_lue;
// calcul de xnp, xmu et Qzero
ParaMateriauxAvecPhaseDeltaSig(xnp,delta_sigma_barre_BH_OiaR,delta_sigma_barre_R_a_tauBH,xmu,Qzero);
};
double phitau = (delta_sigma_barre_R_a_tauBH && delta_barre_epsBH) ; // init partie radial // -- modif 16 juin 2012 --
if (!trajet_neutre) // -- modif 16 juin 2012 --
phitau -= QdeltaSigma * QdeltaSigma * wprime_point / (2. * xmu * wprime); // -- modif 16 juin 2012 --
// -- modif 7 mai 2012 ---
// dans le cas où le phidt est négatif, l'équation constitutive n'est plus valide, dans ce cas on considère
// que phi est nulle, c-a-d que l'on évolue suivant un trajet neutre ce qui n'est sans doute pas l'idéal, mais permet d'éviter l'explosion
// -- calcul de phidt
// if ((phitau < 0) && force_phi_zero_si_negatif)
if (phitau < 0) // -- modif 9 fev 2015 --
{if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 8)) || (Permet_affichage() > 8) )
cout << "\n phidt= " << phitau << "\n Hysteresis3D::Sigma_point(..." << endl;
if ( Dabs(phitau) < force_phi_zero_si_negatif) // -- modif 9 fev 2015 --
{ if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
// if (MaX(ParaGlob::NiveauImpression(),Permet_affichage()) >= 5)
{ cout << "\n attention on utilise une valeur de phi negative dans pour le calcul de l'equation constitutive "
<< " phitau= " << phitau
<< "\n Hysteresis3D::Sigma_point(..." ;
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
phitau = 0.;trajet_neutre = true;
};
};
// fin modif modif 7 mai 2012 ---
if (trajet_neutre) // -- modif 16 juin 2012 --
{ // si c'est un trajet neutre, l'équation constitutive est beaucoup plus simple: sigma_point = 2 mu D_barre
// variables intermédiaires pour les tests
betaphideltasigHB.Inita(0.);
deuxmudeltaepsHB = (2. * xmu ) * delta_barre_epsBH;
// formule pour un trajet neutre : sigma_pointBH = 2. * xmu * delta_barre_epsBH ;
sigma_pointBH = deuxmudeltaepsHB;
// --- si dépendance à la déformation équivalente on calcul le xnp correspondant à un phi = 0
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
double phitau = 0.;
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
ParaMateriauxAvecDefEqui(xnp,QdeltaSigma,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
}
else
{ // cas d'un trajet non neutre
// modif 7 mai 2012 --- on ne recalcule plus le phi qui a été calculé plus haut pour savoir s'il était négatif
// // --- calcul de phi
// double phitau = (delta_sigma_barre_R_a_tauBH && delta_barre_epsBH)
// - QdeltaSigma * QdeltaSigma * wprime_point / (2. * xmu * wprime);
// --- si dépendance à la déformation équivalente
if (avec_defEqui)
{ if (!avec_phaseDeltaSig && !avec_phaseEps) xnp = xnp_lue; // sinon ça déjà été fait
if (Qzero_defEqui!=NULL) Qzero = Qzero_lue;
// calcul de xnp (qui est le seul paramètre de retour de la méthode)
ParaMateriauxAvecDefEqui(xnp,QdeltaSigma,phitau,wprime,def_equi_R,trajet_neutre,def_i_equi,defEqui_maxi);
};
// --- calcul de Beta
double unsurwprimeQ0_puiss_np = 1./pow(wprime*Qzero,xnp);
double beta=0.; // init dans le cas où ne le calcul pas après
if (QdeltaSigma >= ConstMath::petit)
// on ne calcul beta que si la norme de QdeltaSigma est significative
{if (xnp == 2.) { beta = - 2. * xmu * unsurwprimeQ0_puiss_np;}
else {beta = - 2. * xmu * unsurwprimeQ0_puiss_np * pow(QdeltaSigma,xnp-2.);};
};
// -- calcul de la dérivée temporelle de la contrainte
// on considère que D est constant = delta_barre_epsBH/deltat et comme deltat = 1
// ==> = delta_barre_epsBH
// variables intermédiaires pour les tests
betaphideltasigHB = beta * phitau * delta_sigma_barre_R_a_tauBH;
deuxmudeltaepsHB = (2. * xmu ) * delta_barre_epsBH;
// formule normale : sigma_pointBH = 2. * xmu * delta_barre_epsBH + beta * phitau * delta_sigma_barre_R_a_tauBH;
sigma_pointBH = deuxmudeltaepsHB + betaphideltasigHB;
}; // fin du test sur trajet_neutre
// si pb : a faire sur toutes les composantes: limitations des dérivées: spécifiques à la loi d'évolution
// on limite la variation de la dérivée de la contrainte
/* for (int i=1;i<=3;i++) for (int j=1;j<=3;j++)
if (deuxmudeltaepsHB(i,j) > 0.)
{// cas d'une contrainte positive et normalement la dérivée doit évoluer entre 2mu et 0 (saturation)
if (sigma_pointBH(i,j) < 0. ) sigma_pointBH(i,j) = 0.;
if (sigma_pointBH(i,j) > deuxmudeltaepsHB(i,j)) sigma_pointBH(i,j) = deuxmudeltaepsHB(i,j);
}
else
{// cas d'une contrainte négative et normalement la dérivée doit évoluer entre -2mu et 0 (saturation)
if (sigma_pointBH(i,j) < deuxmudeltaepsHB(i,j) ) sigma_pointBH(i,j) = deuxmudeltaepsHB(i,j);
if (sigma_pointBH(i,j) > 0.) sigma_pointBH(i,j) = 0.;
}; */
// si pb : a faire sur toutes les composantes: limitations des dérivées: spécifiques à la loi d'évolution
// on limite la variation de la dérivée de la contrainte au niveau de la dérivée à l'origine
// en fait cela introduit une erreur systématique "dramatique !!" sur le résultat qui est alors erroné !!!
// car il y a accumulation, et il n'y a pas de détection de cette erreur
// donc au lieu de regarder si l'on est juste supérieur, on met une petite marge ???
// non pas bon non plus, donc on supprime !!! à voir par la suite
// sans doute que s'il y a une erreur il faut récupérer une exception à un autre endroit
if (type_resolution_equa_constitutive==6) // cas exploratoire
{ double marge = 100.;
//// for (int i=1;i<=3;i++) for (int j=1;j<=3;j++)
//// if (Dabs(sigma_pointBH(i,j)) > marge*Dabs(deuxmudeltaepsHB(i,j))) sigma_pointBH(i,j) = deuxmudeltaepsHB(i,j);
// on tente une autre stratégie: on regarde l'intensité et on limite, parcontre on garde la même direction que celle calculée
// calcul de l'intensité de cisaillement de la vitesse de déformation
double deuxmuQdeltaeps = sqrt(deuxmudeltaepsHB.II());
// on n'intervient que si deuxmuQdeltaeps a une valeur non nulle
if (deuxmuQdeltaeps > ConstMath::petit)
if (QdeltaSigma > marge * deuxmuQdeltaeps) // normalement pb
{sigma_pointBH *= marge*deuxmuQdeltaeps/QdeltaSigma;
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 3)) || (Permet_affichage() > 3))
{cout << "\n limitation de sigma_point ";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
};
};
};
// retour de la dérivée temporelle de la contrainte
for (int i=1;i<=9;i++)
sig_point(i) = sigma_pointBH(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i));
};
// retour de la dérivée, qui n'a été calculée que si l'erreur = 0
return sig_point;
};
// vérification de l'intégrité du sigma calculé
// erreur : =0: le calcul est licite, si diff de 0, indique qu'il y a eu une erreur
// =1: la norme de sigma est supérieure à la valeur limite de saturation
void Hysteresis3D::Verif_integrite_Sigma(const double & , const Vecteur & sigma_tau,int & erreur)
{ // récup de la contrainte à tau
for (int i=1;i<=9;i++)
sigma_tauBH.Coor(Tenseur3BH::idx_i(i),Tenseur3BH::idx_j(i)) = sigma_tau(i);
//modif : 1 juin 2014
// non on fait la vérif sur le delta cf. plus bas
// // on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
// double sigma_tau_II = sigma_tauBH.II();
// if ( sigma_tau_II > Qzero*Qzero)
// {erreur = 1;
// //--- pour le debug
// if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
// cout << "\n sigma_tau_II =" << sigma_tau_II << " au lieu de "<< Qzero*Qzero
// << " Hysteresis3D::Sigma_point " << endl ;
// //--- fin debug
// }
// else // sinon on continue
// { erreur = 0;
//--- pour le debug
//cout << "\n tt np=" << xnp << " Q0=" << Qzero << " xmu=" << xmu << " T=" << *temperature;
//--- fin debug
// variation de sigma de R à tau
delta_sigma_barre_R_a_tauBH = sigma_tauBH - sigma_barre_BH_R; // ne sert que dans Sigma_point
delta_sigma_barre_t_a_tauBH = sigma_tauBH - sigma_i_barre_BH; // ne sert que dans Sigma_point
// calcul de QdeltaSigma: on prend la valeur absolue car dans les petits nombres on peut avoir des presques 0 !!
// négatifs qui conduisent à une racine infinie !!
double QdeltaSigma = sqrt(delta_sigma_barre_R_a_tauBH.II());
//modif : 1 juin 2014
// on vérifie que l'amplitude de la contrainte transmise, n'est pas supérieure à la saturation
// on fait la vérification sur le delta, car si on est sur une branche non initiale, on peut très bien
// dépasser le Qzero, localement, mais ce sera invalidé ensuite avec l'algo de coincidence
if ( QdeltaSigma > (Qzero*depassement_Q0*wBase))
{erreur = 1;
//--- pour le debug
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
cout << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
<< " Hysteresis3D::Sigma_point " << endl ;
//--- fin debug
}
else // sinon on continue
{ erreur = 0;
};
};
// update des différentes énergies sur le pas de temps ou la partie de pas de temps effectué
// et calcul de defEqui finalisée
void Hysteresis3D::UpdateEnergieHyste(EnergieMeca & energ
,const Tenseur3BH & sigma_deb_BH,const Tenseur3BH & sigma_fin_BH
,const Tenseur3BH & delta_sigma_sur_tau_BH
,const Tenseur3BH & delta_sigma_R_tau_BH
,const double& phi_2_dt
,const Tenseur3BH & delta_eps_sur_deltaTau_BH,double& QdeltaSigmaRatdt
,const double& wprime,const double& w_prime_point,const bool& trajet_neutre
,const double& def_equi_R,const double& def_i_equi, double& def_equi)
{ // delta tau = 1 ici, en fait c'est un paramètre formel
// double delta_ener_interne = - 0.5*((sigma_deb_BH + sigma_fin_BH) && delta_eps_sur_deltaTau_BH);
double delta_ener_interne = - (sigma_fin_BH && delta_eps_sur_deltaTau_BH);
// -------- débug au cas où -----------
//if (phi_2_dt < 0)
// cout << "\n attention, phi_2_dt = " << phi_2_dt << " est negatif ??, normalement, a l'issue d'une convergence"
// << " phi est positif ou nul , a regarder en debug "
// << "\n Hysteresis3D::UpdateEnergieHyste(.... ";
// ------- fin débug --------
if (trajet_neutre)
{// -- calcul du delta énergie interne
// a priori ici tout est élastique: phi = 0, tau_C=0,
// -- calcul du taux fondamental C (formule thèse Mayssa Chache p. 71: formule 2.51, avec les formules
// 2.52, )
// calcul de la déformation équivalente
def_equi = def_i_equi ;
// -- delta énergies : on considère un tau C constant sur l'incrément de temps ainsi que phi
double delta_energie_elastique = - delta_ener_interne ;
double delta_energie_frotteurs = 0.;
EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
energ += delta_ener;
}
else
{// -- calcul du delta énergie interne
// -- calcul du taux fondamental C (formule thèse Mayssa Chache p. 71: formule 2.51, avec les formules
// 2.52, )
double unsurwprime = 1./wprime;
// si QdeltaSigmaRatdt est très petit on considère qu'il n'y a qu'un calcul élastique, pas de dissipation
// on plus simple on considère que la def équivalente ne bouge pas
if (QdeltaSigmaRatdt >= ConstMath::petit )
{ double unsurQdeltaSigma=1./QdeltaSigmaRatdt;
def_equi = def_i_equi + unsurwprime * phi_2_dt * unsurQdeltaSigma;
double def_equi_Ratdt = def_equi-def_equi_R;
// calcul de QpointDeltaSigma
double QpointDeltaSigma = unsurQdeltaSigma * (delta_sigma_R_tau_BH && delta_sigma_sur_tau_BH);
// calcul du tau_C qui représente en fait : delta_R^tdt eps : sigma_point
double tau_C = (wprime * QpointDeltaSigma - w_prime_point * QdeltaSigmaRatdt) * def_equi_Ratdt;
// -- delta énergies : on considère un tau C constant sur l'incrément de temps ainsi que phi
double delta_energie_elastique = - delta_ener_interne - unsurwprime * phi_2_dt + unsurwprime * tau_C;
double delta_energie_frotteurs = unsurwprime * ( phi_2_dt - tau_C);
// si l'énergie des frotteurs est négative, on considère que c'est due aux différentes imprécisions,
// comme je ne vois pas de solution, donc
// on la met à 0, et on modifie le tau_C en conséquence, et donc
// on recalcule l'énergie élastique en conséquence, ce qui conduit, tout calcul fait à un trajet
// purement élastique ( finalement idem que neutre ??)
if (delta_energie_frotteurs < 0.)
{ delta_energie_frotteurs = 0.;
delta_energie_elastique = - delta_ener_interne;
}
//-----pour le debug
// if (delta_energie_frotteurs < 0.)
// cout << "\n (1) delta_energie_frotteurs= " << delta_energie_frotteurs << endl;
// ----- fin débug
// si l'énergie frotteur est négative, on considère que c'est due à une erreur d'arrondie, on met 0 à la place
// if (delta_energie_frotteurs < 0.)
// delta_energie_frotteurs = 0.;
EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
energ += delta_ener;
}
else
{ // cas où l'accroissement de sigma est trop faible donc considéré nul
// on a sigma_point = delta_sigma/delta_t, avec ici delta_t=1 donc delta_sigma = 0 => sigma_point=0
// comme en fait : tau_C= delta_R^(t+deltat) eps : sigma_point, on a tau_C = 0
// on considère donc tau_C = 0
double delta_energie_elastique = - delta_ener_interne - unsurwprime * phi_2_dt ;
double delta_energie_frotteurs = unsurwprime * ( phi_2_dt );
//-----pour le debug
if ( ((Permet_affichage()==0) && (ParaGlob::NiveauImpression() > 5)) || (Permet_affichage() > 5))
if (delta_energie_frotteurs < 0)
cout << "\n (2) delta_energie_frotteurs= " << delta_energie_frotteurs << endl;
// ----- fin débug
EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
energ += delta_ener;
};
//-------- ancien calcul en commentaire -------
// double unsurQdeltaSigma= ConstMath::grand; // on limite la valeur par exemple pour le démarrage
// if (QdeltaSigmaRatdt >= ConstMath::petit) unsurQdeltaSigma=1./QdeltaSigmaRatdt;
// // calcul de la déformation équivalente
// def_equi = def_equi_R + unsurwprime * phi_2_dt * unsurQdeltaSigma;
// // calcul de QpointDeltaSigma
// double QpointDeltaSigma = unsurQdeltaSigma * (delta_sigma_R_tau_BH && delta_sigma_sur_tau_BH);
// calcul du tau_C qui représente en fait : delta_R^tdt eps : sigma_point
// double tau_C = (wprime * QpointDeltaSigma - w_prime_point * QdeltaSigmaRatdt) * def_equi;
// -- delta énergies : on considère un tau C constant sur l'incrément de temps ainsi que phi
// double delta_energie_elastique = - delta_ener_interne - unsurwprime * phi_2_dt + unsurwprime * tau_C;
// double delta_energie_frotteurs = unsurwprime * ( phi_2_dt - tau_C);
// EnergieMeca delta_ener(delta_energie_elastique,delta_energie_frotteurs,0.);
// energ += delta_ener;
//-------- fin ancien calcul en commentaire -------
};
};
// vérif de la coïncidence: prog de mise au point
// premiere_fois: comme il y a une récursion, pour arrêter la boucle infinie
void Hysteresis3D::Verif_coincidence(const Tenseur3BB & delta_epsBB,const Tenseur3HH & gijHH
,const Tenseur3BB & gijBB,int cas,SaveResulHysteresis3D & save_resul
,Tenseur3HH & sigHH,const EnergieMeca & energ_t,EnergieMeca & energ
,bool premiere_fois)
{
// l'objectif est ici de vérifier que l'on a bien la position de la nouvelle contrainte cohérente avec
// le rayon de référence
// --- 1ere étape, récupérer les infos de position de référence qui vont-être enregistrées
Tenseur3BH sigBH_de_ref; Tenseur3BH centreBH_de_ref; // init à 0 par défaut
int aucun_pt_inv=false;
switch (save_resul.modif)
{ case 0: // rien ne change sur le pas en terme de coincidence et d'inversion
{ if(save_resul.sigma_barre_BH_R.begin() != save_resul.sigma_barre_BH_R.end())
{sigBH_de_ref = *(save_resul.sigma_barre_BH_R.begin());}
else aucun_pt_inv=true;
if(save_resul.oc_BH_R.begin() != save_resul.oc_BH_R.end())
centreBH_de_ref = *(save_resul.oc_BH_R.begin());
break;
}
case 1: // cas où l'on a une ou plusieurs coincidences
{ // a chaque coincidence il faudra supprimer 1 pt de ref et un centre de cercle
// sauf si l'on arrive au centre nul ou au point de référence null
List_io <Tenseur3BH>::iterator itens;List_io <Tenseur3BH>::iterator ioc_i;
itens = save_resul.sigma_barre_BH_R.begin(); ioc_i = save_resul.oc_BH_R.begin();
for (int i=1;i<=save_resul.nb_coincidence;i++)
{ if (itens != save_resul.sigma_barre_BH_R.end()) // si itens == sigma_barre_BH_R.end(), cela signifie
itens++; // que la référence est le tenseur nulle, on arrête de dépiler
if (ioc_i != save_resul.oc_BH_R.end()) // si ioc_i == oc_BH_R.end(), cela signifie que le centre de référence
ioc_i++; // est le tenseur nulle
};
if(itens != save_resul.sigma_barre_BH_R.end())
{sigBH_de_ref = *(itens);}
else aucun_pt_inv=true;
if (ioc_i != save_resul.oc_BH_R.end())
centreBH_de_ref = *(ioc_i);
break;
}
case 2: // cas où l'on a une ou plusieurs inversions
{ // on traite les centres différement des points d'inversions
// --- cas des centres ---
if (save_resul.oc_BH_t_a_tdt.size() != 0)
{centreBH_de_ref = *(save_resul.oc_BH_t_a_tdt.begin());} // on récupère le premier
else if (save_resul.oc_BH_R.size() != 0)
{centreBH_de_ref = *(save_resul.oc_BH_R.begin());}; // on récupère le premier
// sinon tous les centres sont nuls on utilise le centre de départ initialisé à 0
// --- cas des ref ---
if (save_resul.sigma_barre_BH_R_t_a_tdt.size() != 0)
{sigBH_de_ref = *(save_resul.sigma_barre_BH_R_t_a_tdt.begin());}
else if (save_resul.sigma_barre_BH_R.size() != 0)
{sigBH_de_ref = *(save_resul.sigma_barre_BH_R.begin());}
else aucun_pt_inv=true;
break;
}
case 3: // cas où l'on a une ou plusieurs coincidence(s) et inversion (s)
{ // on traite les centres différement des points d'inversions
// --- cas des centres ---
if (save_resul.oc_BH_t_a_tdt.size() != 0) // cas on ajoutera des inversions
{centreBH_de_ref = *(save_resul.oc_BH_t_a_tdt.begin());} // on récupère le premier
else // cas où on a déjà une coincidence et pas d'inversion
{ List_io <Tenseur3BH>::iterator ioc_i = save_resul.oc_BH_R.begin();
for (int i=1;i<=save_resul.nb_coincidence;i++)
if (ioc_i != save_resul.oc_BH_R.end()) // si ioc_i == oc_BH_R.end(), cela signifie que le centre de référence
ioc_i++; // est le tenseur nulle
if (ioc_i != save_resul.oc_BH_R.end())
centreBH_de_ref = *(ioc_i);
};
// --- cas des références ----
if (save_resul.sigma_barre_BH_R_t_a_tdt.size() != 0)
{sigBH_de_ref = *(save_resul.sigma_barre_BH_R_t_a_tdt.begin());}
else // cas où on a déjà une coincidence et pas d'inversion
{ List_io <Tenseur3BH>::iterator itens = save_resul.sigma_barre_BH_R.begin();
for (int i=1;i<=save_resul.nb_coincidence;i++)
if (itens != save_resul.sigma_barre_BH_R.end()) // si itens == sigma_barre_BH_R.end(), cela signifie
itens++; // que la référence est le tenseur nulle, on arrête de dépiler
if(itens != save_resul.sigma_barre_BH_R.end())
{sigBH_de_ref = *(itens);}
else aucun_pt_inv=true;
};
break;
}
default:
cout << "\n erreur modif= " << save_resul.modif
<< "\n Hysteresis3D::Verif_coincidence(... ";
// on génère une exception
throw ErrNonConvergence_loiDeComportement();
Sortie(1);
}
// --- calcul du rayon de Oi à R sigBH_de_ref
Tenseur3BH ddelta_sigma_barre_BH_OiaR = sigBH_de_ref - centreBH_de_ref; // le tenseur rayon
// calcul du rayon: la première fois le rayon est infini
double QQ_OiaR=0.;
if (aucun_pt_inv) {QQ_OiaR = ConstMath::tresgrand;}
else {QQ_OiaR = sqrt(ddelta_sigma_barre_BH_OiaR.II());};
// --- calcul du nouveau rayon ----
Tenseur3BH delta_sigma_oi_a_tdt_BH = save_resul.sigma_barre_BH_tdt - centreBH_de_ref; // OA
double rayon_tdt = sqrt(delta_sigma_oi_a_tdt_BH.II());
// ---- vérif ---
if (rayon_tdt > QQ_OiaR) // cas d'un inchoérence
{ cout << "\n *********** inchoerence de coincidence ";
if (Ptintmeca_en_cours() != NULL) Ptintmeca_en_cours()->Signature();cout << flush;
// re-appel de l'algo pour passer en débug dedans
// initialisation des variables de travail
save_resul.Init_debut_calcul();
// initialisation éventuelle des variables thermo-dépendantes
Init_thermo_dependance();
// vérif de l'emboitement
save_resul.Verif_centre_reference(Permet_affichage());
// re-calcul pour vérif pas à pas
Avancement_temporel(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ);
if (premiere_fois)
{premiere_fois = false;
Verif_coincidence(delta_epsBB,gijHH,gijBB,cas,save_resul,sigHH,energ_t,energ,premiere_fois);
};
// on génère une exception
throw ErrNonConvergence_loiDeComportement();
Sortie(1);
};
};
// vérif de la cohérence des centres et références: les rayons associés doivent être de taille
// décroissante
void Hysteresis3D::SaveResulHysteresis3D::Verif_centre_reference(const int &permet_affichage)
{ // on va balayer les différents centres et points de références, pour construire les différents rayons
// qui devront être de tailles croissantes à partir du premier
double rayon_precedent = 0.; List_io <double> les_rayons;
List_io <Tenseur3BH>::iterator itens,itens_fin = sigma_barre_BH_R.end();
List_io <Tenseur3BH>::iterator ioc_i,ioc_i_fin = oc_BH_R.end();
for (itens=sigma_barre_BH_R.begin(), ioc_i = oc_BH_R.begin();
itens != itens_fin; itens++,ioc_i++)
{double rayon = 0;
if (ioc_i != ioc_i_fin)
{ Tenseur3BH dsigmaBH = *itens - *ioc_i;
rayon = sqrt(dsigmaBH.II());
}
else
{ Tenseur3BH dsigmaBH = *itens;
rayon = sqrt(dsigmaBH.II());
}
les_rayons.push_front(rayon);
if (rayon >= rayon_precedent)
{ rayon_precedent = rayon;}
else
{ cout << "\n ---- ******* attention la croissance des rayons n'est pas verifiee !! " ;
if ( ((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5)) || (permet_affichage > 5))
{// les listes
cout << "\n save_resul.sigma_barre_BH_R_t_a_tdt= " << sigma_barre_BH_R_t_a_tdt;
cout << "\n save_resul.sigma_barre_BH_R= " << sigma_barre_BH_R;
cout << "\n save_resul.oc_BH_t_a_tdt= " << oc_BH_t_a_tdt;
cout << "\n save_resul.oc_BH_R= " << oc_BH_R;
cout << "\n les rayons= " << les_rayons;
// on vide le buffer
cout << endl;
};
// Sortie(1);
};
};
};
Reference in a new issue View git blame Copy permalink