2021-09-23 11:21:15 +02:00
|
|
|
|
|
|
|
|
// 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.
|
|
|
|
|
//
|
2023-05-03 17:23:49 +02:00
|
|
|
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
|
2021-09-23 11:21:15 +02:00
|
|
|
// AUTHOR : Gérard Rio
|
|
|
|
|
// E-MAIL : gerardrio56@free.fr
|
|
|
|
|
//
|
|
|
|
|
// This program is free software: you can redistribute it and/or modify
|
|
|
|
|
// it under the terms of the GNU General Public License as published by
|
|
|
|
|
// the Free Software Foundation, either version 3 of the License,
|
|
|
|
|
// or (at your option) any later version.
|
|
|
|
|
//
|
|
|
|
|
// This program is distributed in the hope that it will be useful,
|
|
|
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty
|
|
|
|
|
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
|
|
|
// See the GNU General Public License for more details.
|
|
|
|
|
//
|
|
|
|
|
// You should have received a copy of the GNU General Public License
|
|
|
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
|
|
|
//
|
|
|
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
//#include "Debug.h"
|
|
|
|
|
|
|
|
|
|
# include <iostream>
|
|
|
|
|
using namespace std; //introduces namespace std
|
|
|
|
|
#include <math.h>
|
|
|
|
|
#include <stdlib.h>
|
|
|
|
|
#include "Sortie.h"
|
|
|
|
|
#include "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 --
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 8)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n phidt= " << phidt << "\n Hysteresis3D::Residu_constitutif(..." << endl;
|
|
|
|
|
if ( Dabs(phidt) < force_phi_zero_si_negatif) // -- modif 9 fev 2015 --
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ 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 --
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 8)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n phidt= " << phidt << "\n Hysteresis3D::Mat_tangente_constitutif(..." << endl;
|
|
|
|
|
if ( Dabs(phidt) < force_phi_zero_si_negatif) // -- modif 9 fev 2015 --
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ 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)
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
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;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
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();
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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)) )
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
// 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;
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
// 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))
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n Hysteresis3D::CalculContrainte_tdt "
|
|
|
|
|
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
|
|
|
|
|
<< erreur_maxi_global << endl;
|
|
|
|
|
|
|
|
|
|
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
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ 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 ;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ cout << "\n abs_racine_max " << absracinemax;}
|
|
|
|
|
if (Permet_affichage() > 4)
|
|
|
|
|
{cout << " (2) utilisation avec succes de la linearisation ";};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
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);
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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))
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n Hysteresis3D::CalculContrainte_tdt "
|
|
|
|
|
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
|
|
|
|
|
<< erreur_maxi_global << endl;
|
|
|
|
|
// gestion de l'erreur de retour
|
|
|
|
|
double abssigmax=sig_finale.Max_val_abs();
|
|
|
|
|
if ((conver != 2) || (!isfinite(abssigmax)) || (isnan(abssigmax)) )
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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))
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
|
|
|
|
};
|
|
|
|
|
// gestion de l'erreur de retour
|
|
|
|
|
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
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ 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;
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
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
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
|
|
|
|
};
|
|
|
|
|
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 7)
|
2021-09-23 11:21:15 +02:00
|
|
|
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))
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n Hysteresis3D::CalculContrainte_tdt "
|
|
|
|
|
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
|
|
|
|
|
<< erreur_maxi_global << endl;
|
|
|
|
|
// gestion de l'erreur de retour
|
|
|
|
|
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
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ 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 ";
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
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
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
|
|
|
|
};
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n Hysteresis3D::CalculContrainte_tdt "
|
|
|
|
|
<< "nombreAppelF= " << nombreAppelF << " nb_step= " << nb_step << " erreur_maxi_global= "
|
|
|
|
|
<< erreur_maxi_global << endl;
|
|
|
|
|
if (nb_boucle > 10) // si en * par 10^10 on n'a pas réussi on arrête
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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)) );
|
|
|
|
|
|
|
|
|
|
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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))
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 0)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
|
|
|
|
|
<< " Hysteresis3D::Sigma_point " << endl ;
|
|
|
|
|
}
|
|
|
|
|
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 --
|
2023-05-03 17:23:49 +02:00
|
|
|
{if (Permet_affichage() > 8)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n phidt= " << phitau << "\n Hysteresis3D::Sigma_point(..." << endl;
|
|
|
|
|
if ( Dabs(phitau) < force_phi_zero_si_negatif) // -- modif 9 fev 2015 --
|
2023-05-03 17:23:49 +02:00
|
|
|
{ if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
{ 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;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 3)
|
2021-09-23 11:21:15 +02:00
|
|
|
{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;
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
cout << "\n QdeltaSigma =" << QdeltaSigma << " au lieu de "<< (Qzero*depassement_Q0*wBase)
|
2023-05-03 17:23:49 +02:00
|
|
|
<< " Hysteresis3D::Verif_integrite_Sigma " << endl ;
|
2021-09-23 11:21:15 +02:00
|
|
|
}
|
|
|
|
|
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 );
|
2023-05-03 17:23:49 +02:00
|
|
|
if (Permet_affichage() > 5)
|
2021-09-23 11:21:15 +02:00
|
|
|
if (delta_energie_frotteurs < 0)
|
|
|
|
|
cout << "\n (2) delta_energie_frotteurs= " << delta_energie_frotteurs << endl;
|
|
|
|
|
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);
|
|
|
|
|
};
|
|
|
|
|
};
|
|
|
|
|
};
|
|
|
|
|
|
