917 lines
41 KiB
C++
917 lines
41 KiB
C++
// FICHIER : Hyper_W_gene_3D.cc
|
|
// CLASSE : Hyper_W_gene_3D
|
|
|
|
|
|
// This file is part of the Herezh++ application.
|
|
//
|
|
// The finite element software Herezh++ is dedicated to the field
|
|
// of mechanics for large transformations of solid structures.
|
|
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
|
|
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
|
|
//
|
|
// Herezh++ is distributed under GPL 3 license ou ultérieure.
|
|
//
|
|
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
|
|
// AUTHOR : Gérard Rio
|
|
// E-MAIL : gerardrio56@free.fr
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License,
|
|
// or (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty
|
|
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
// See the GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
//
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
|
|
|
|
|
#include "Hyper_W_gene_3D.h"
|
|
|
|
//#include "Debug.h"
|
|
|
|
# include <iostream>
|
|
using namespace std; //introduces namespace std
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include "Sortie.h"
|
|
#include "TypeConsTens.h"
|
|
#include "ParaGlob.h"
|
|
#include "Enum_TypeQuelconque.h"
|
|
#include "TypeQuelconqueParticulier.h"
|
|
#include "CharUtil.h"
|
|
|
|
// CONSTRUCTEURS :
|
|
// ---------- classe de stockage des grandeurs spécifiques pour la loi ---------
|
|
|
|
// constructeur par défaut
|
|
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::SaveResulHyper_W_gene_3D() :
|
|
invP(NULL),invP_t(NULL),map_type_quelconque()
|
|
,para_loi(NULL),para_loi_t(NULL)
|
|
{ // mise à jour de la liste des grandeurs quelconques internes
|
|
Mise_a_jour_map_type_quelconque();
|
|
};
|
|
// avec_para : le nombre indique si on sauvegarde des paramètres matériaux qui varient
|
|
// et également des invariants.
|
|
// si == 0 : pas de stockage (joue en fait le role de sortie_post)
|
|
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::SaveResulHyper_W_gene_3D(int avec_para) :
|
|
invP(NULL),invP_t(NULL)
|
|
,para_loi(NULL),para_loi_t(NULL)
|
|
,map_type_quelconque()
|
|
{ if (avec_para)
|
|
{invP = new Invariantpost3D();
|
|
invP_t = new Invariantpost3D();
|
|
};
|
|
if (avec_para)
|
|
{para_loi = new Vecteur(avec_para,-ConstMath::trespetit);
|
|
para_loi_t = new Vecteur(avec_para,-ConstMath::trespetit);
|
|
};
|
|
// mise à jour de la liste des grandeurs quelconques internes
|
|
Mise_a_jour_map_type_quelconque();
|
|
};
|
|
|
|
// constructeur de copie
|
|
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::SaveResulHyper_W_gene_3D(const SaveResulHyper_W_gene_3D& sav):
|
|
invP(NULL),invP_t(NULL),para_loi(NULL),para_loi_t(NULL)
|
|
,map_type_quelconque(sav.map_type_quelconque)
|
|
{ if (sav.invP != NULL)
|
|
{invP = new Invariantpost3D(*(sav.invP));
|
|
invP_t = new Invariantpost3D(*(sav.invP_t));
|
|
};
|
|
if (sav.para_loi != NULL)
|
|
{para_loi = new Vecteur(*sav.para_loi);
|
|
para_loi_t = new Vecteur(*sav.para_loi_t);
|
|
};
|
|
// mise à jour de la liste des grandeurs quelconques internes
|
|
Mise_a_jour_map_type_quelconque();
|
|
|
|
};
|
|
|
|
// destructeur
|
|
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::~SaveResulHyper_W_gene_3D()
|
|
{ if (invP != NULL)
|
|
{ delete invP;
|
|
delete invP_t;
|
|
};
|
|
if (para_loi != NULL) delete para_loi;
|
|
if (para_loi_t != NULL) delete para_loi_t;
|
|
};
|
|
|
|
// affectation
|
|
Loi_comp_abstraite::SaveResul & Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::operator = ( const Loi_comp_abstraite::SaveResul & a)
|
|
{ Hyper_W_gene_3D::SaveResulHyper_W_gene_3D& sav = *((Hyper_W_gene_3D::SaveResulHyper_W_gene_3D*) &a);
|
|
if (sav.invP != NULL)
|
|
{ if (invP == NULL)
|
|
{invP = new Invariantpost3D(*(sav.invP));
|
|
invP_t = new Invariantpost3D(*(sav.invP_t));
|
|
}
|
|
else
|
|
{ *invP = (*(sav.invP)); *invP_t = (*(sav.invP_t));
|
|
};
|
|
}
|
|
else
|
|
// sinon cela veut dire que invP et invP_t ne doivent pas être affecté
|
|
{ if (invP != NULL)
|
|
{delete invP;invP=NULL;
|
|
delete invP_t;invP_t=NULL;
|
|
};
|
|
};
|
|
map_type_quelconque = sav.map_type_quelconque;
|
|
if (sav.para_loi != NULL)
|
|
{if (para_loi == NULL)
|
|
{para_loi = new Vecteur(sav.para_loi->Taille(),-ConstMath::trespetit);
|
|
para_loi_t = new Vecteur(sav.para_loi->Taille(),-ConstMath::trespetit);
|
|
};
|
|
*para_loi = *(sav.para_loi);
|
|
*para_loi_t = *(sav.para_loi_t);
|
|
};
|
|
|
|
return *this;
|
|
};
|
|
|
|
//------- lecture écriture dans base info -------
|
|
// cas donne le niveau de la récupération
|
|
// = 1 : on récupère tout
|
|
// = 2 : on récupère uniquement les données variables (supposées comme telles)
|
|
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Lecture_base_info (ifstream& ent,const int )
|
|
{ string toto; int presence=0;
|
|
ent >> toto >> presence;
|
|
#ifdef MISE_AU_POINT
|
|
if (toto != "dat_sp_hyperelas:")
|
|
{ cout << "\n erreur dans la lecture des donnees specifiques hyperelastiques, on ne trouve pas le mot cle: dat_sp_hyperelas: "
|
|
<< "\n Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Lecture_base_info (... ";
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
if (presence == 1)
|
|
{if (invP == NULL)
|
|
{invP = new Invariantpost3D();
|
|
invP_t = new Invariantpost3D();
|
|
};
|
|
// lecture
|
|
ent >> (*invP); (*invP_t)= (*invP); //toto >> invP->V >> toto >> invP->Qeps >> toto >> invP->cos3phi;
|
|
}
|
|
else
|
|
{ if (invP != NULL)
|
|
{ delete invP;invP=NULL;
|
|
delete invP_t;invP_t=NULL;
|
|
};
|
|
};
|
|
// cas des paramètres de lois éventuelles
|
|
ent >> toto;
|
|
if (toto == "sans_para")
|
|
{if (para_loi != NULL)
|
|
{delete para_loi;para_loi=NULL;
|
|
delete para_loi_t;para_loi_t=NULL;
|
|
};
|
|
}
|
|
else if (toto == "avec_para_taille=")
|
|
{int taille;
|
|
ent >> taille;
|
|
if (para_loi == NULL)
|
|
{para_loi = new Vecteur(taille);
|
|
para_loi_t = new Vecteur(taille);
|
|
};
|
|
ent >> *para_loi;
|
|
*para_loi_t=*para_loi;
|
|
};
|
|
// mise à jour de la liste des grandeurs quelconques internes
|
|
Mise_a_jour_map_type_quelconque();
|
|
|
|
};
|
|
|
|
// cas donne le niveau de sauvegarde
|
|
// = 1 : on sauvegarde tout
|
|
// = 2 : on sauvegarde uniquement les données variables
|
|
//(supposées comme telles)
|
|
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Ecriture_base_info(ofstream& sort,const int )
|
|
{ sort << "\n dat_sp_hyperelas: ";
|
|
if (invP != NULL)
|
|
{ sort << " 1 " << (*invP); }
|
|
else { sort << " 0 ";};
|
|
// cas des paramètres de lois éventuelles
|
|
if (para_loi == NULL)
|
|
{sort << " sans_para ";}
|
|
else
|
|
{sort << " avec_para_taille= "<<para_loi->Taille() << " ";
|
|
sort << *para_loi;
|
|
};
|
|
|
|
};
|
|
|
|
// mise à jour des informations transitoires
|
|
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::TdtversT()
|
|
{if (invP != NULL)
|
|
(*invP_t) = (*invP);
|
|
if (para_loi != NULL)
|
|
{*para_loi_t = *para_loi;};
|
|
// mise à jour de la liste des grandeurs quelconques internes
|
|
Mise_a_jour_map_type_quelconque();
|
|
};
|
|
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::TversTdt()
|
|
{if (invP != NULL)
|
|
(*invP) = (*invP_t);
|
|
if (para_loi != NULL)
|
|
{*para_loi = *para_loi_t;};
|
|
// mise à jour de la liste des grandeurs quelconques internes
|
|
Mise_a_jour_map_type_quelconque();
|
|
};
|
|
|
|
|
|
// affichage à l'écran des infos
|
|
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Affiche() const
|
|
{ cout << "\n data specifiques hyperelastiques : " ;
|
|
if (invP != NULL)
|
|
{ cout << (*invP); };
|
|
if (para_loi == NULL)
|
|
{cout << " sans_para ";}
|
|
else
|
|
{cout << " avec_para ";
|
|
cout << *para_loi;
|
|
};
|
|
};
|
|
|
|
|
|
// mise à jour de la liste des grandeurs quelconques internes
|
|
void Hyper_W_gene_3D::SaveResulHyper_W_gene_3D::Mise_a_jour_map_type_quelconque()
|
|
{ // on parcours la map
|
|
map < EnumTypeQuelconque , TypeQuelconque, std::less < EnumTypeQuelconque> >::iterator
|
|
il,ilfin=map_type_quelconque.end();
|
|
for (il=map_type_quelconque.begin();il != ilfin;il++)
|
|
|
|
{ switch ((*il).first)
|
|
{ case POTENTIEL:
|
|
{ Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) (*il).second.Grandeur_pointee()); // pour simplifier
|
|
(*tyTQ.ConteneurDouble()) = invP->potentiel;
|
|
break;
|
|
}
|
|
// -----cas des paramètres de la loi
|
|
case PARA_LOI_TYPE_MOONEY:
|
|
{ Grandeur_Vecteur& gr= *((Grandeur_Vecteur*) map_type_quelconque[PARA_LOI_TYPE_MOONEY].Grandeur_pointee()); // pour simplifier
|
|
if (para_loi != NULL)
|
|
{*(gr.ConteneurVecteur()) = *(para_loi);}
|
|
else {gr.ConteneurVecteur()->Zero();};
|
|
break;
|
|
}
|
|
|
|
default:
|
|
cout << "\n *** erreur on demande l'acces a : "
|
|
<< NomTypeQuelconque((*il).first)
|
|
<< " or celui-ci n'est pas dispo pour la loi ";
|
|
this->Affiche();
|
|
cout << " revoir la mise en donnees ! " << endl;
|
|
Sortie(1);
|
|
};
|
|
|
|
};
|
|
|
|
};
|
|
|
|
// ---------- classe Hyper_W_gene_3D ---------
|
|
|
|
// Constructeur par defaut
|
|
Hyper_W_gene_3D::Hyper_W_gene_3D ():
|
|
Loi_comp_abstraite()
|
|
,I_B(0.),I_BB(0.),II_B(0.),III_B(0.)
|
|
,sortie_post(0),nb_para_loi(0)
|
|
// I_B(0.),I_BB(0.),I_BBB(0.),II_B(0.),III_B(0.)
|
|
,J_r(3)
|
|
,d_I_B_epsBB_HH(),d_II_B_epsBB_HH(),d_III_B_epsBB_HH()
|
|
,d_J_r_epsBB_HH(3)
|
|
// ,d_I_B_eps2BB_HHHH(),d_II_B_eps2BB_HHHH(),d_III_B_eps2BB_HHHH()
|
|
,d_J_1_eps2BB_HHHH(),d_J_2_eps2BB_HHHH(),d_J_3_eps2BB_HHHH()
|
|
,V(0.),BB_HH(),BB_HB()
|
|
,IxI_HHHH(),IxbarreI_HHHH(),IxB_HHHH(),BxI_HHHH()
|
|
{
|
|
|
|
};
|
|
|
|
|
|
// Constructeur utile si l'identificateur du nom de la loi
|
|
// de comportement et la dimension sont connus
|
|
// vit_def indique si oui ou non la loi utilise la vitesse de déformation
|
|
Hyper_W_gene_3D::Hyper_W_gene_3D (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp
|
|
,int dimension,bool vit_def):
|
|
Loi_comp_abstraite(id_compor,categorie_comp,dimension,vit_def)
|
|
,I_B(0.),I_BB(0.),II_B(0.),III_B(0.)
|
|
,sortie_post(0),nb_para_loi(0)
|
|
// I_B(0.),I_BB(0.),I_BBB(0.),II_B(0.),III_B(0.)
|
|
,J_r(3)
|
|
,d_I_B_epsBB_HH(),d_II_B_epsBB_HH(),d_III_B_epsBB_HH()
|
|
,d_J_r_epsBB_HH(3)
|
|
// ,d_I_B_eps2BB_HHHH(),d_II_B_eps2BB_HHHH(),d_III_B_eps2BB_HHHH()
|
|
,d_J_1_eps2BB_HHHH(),d_J_2_eps2BB_HHHH(),d_J_3_eps2BB_HHHH()
|
|
,V(0.),BB_HH(),BB_HB()
|
|
,IxI_HHHH(),IxbarreI_HHHH(),IxB_HHHH(),BxI_HHHH()
|
|
{
|
|
|
|
};
|
|
|
|
// Constructeur de copie
|
|
Hyper_W_gene_3D::Hyper_W_gene_3D (const Hyper_W_gene_3D& loi):
|
|
Loi_comp_abstraite(loi)
|
|
,nb_para_loi(loi.nb_para_loi)
|
|
,I_B(loi.I_B),I_BB(loi.I_BB),II_B(II_B),III_B(III_B)
|
|
// I_B(loi.I_B),I_BB(loi.I_BB),I_BBB(loi.I_BBB),II_B(II_B),III_B(III_B)
|
|
,J_r(loi.J_r)
|
|
,d_I_B_epsBB_HH(loi.d_I_B_epsBB_HH),d_II_B_epsBB_HH(loi.d_II_B_epsBB_HH)
|
|
// ,d_III_B_epsBB_HH(loi.d_III_B_epsBB_HH)
|
|
,d_J_r_epsBB_HH(loi.d_J_r_epsBB_HH)
|
|
// ,d_I_B_eps2BB_HHHH(loi.d_I_B_eps2BB_HHHH)
|
|
,d_J_1_eps2BB_HHHH(loi.d_J_1_eps2BB_HHHH)
|
|
,d_J_2_eps2BB_HHHH(loi.d_J_2_eps2BB_HHHH)
|
|
,d_J_3_eps2BB_HHHH(loi.d_J_3_eps2BB_HHHH)
|
|
// ,d_II_B_eps2BB_HHHH(loi.d_II_B_eps2BB_HHHH)
|
|
// ,d_III_B_eps2BB_HHHH(loi.d_III_B_eps2BB_HHHH)
|
|
,V(loi.V),BB_HH(loi.BB_HH),BB_HB(loi.BB_HB)
|
|
,IxI_HHHH(),IxbarreI_HHHH(),IxB_HHHH(),BxI_HHHH()
|
|
,sortie_post(loi.sortie_post)
|
|
{
|
|
|
|
};
|
|
|
|
// DESTRUCTEUR :
|
|
|
|
Hyper_W_gene_3D::~Hyper_W_gene_3D ()
|
|
{
|
|
|
|
};
|
|
|
|
// récupération des grandeurs particulière (hors ddl )
|
|
// correspondant à liTQ
|
|
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
|
|
void Hyper_W_gene_3D::Grandeur_particuliere
|
|
(bool ,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list<int>& decal) const
|
|
{ // ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas
|
|
// on passe en revue la liste
|
|
List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end();
|
|
list<int>::iterator idecal=decal.begin();
|
|
for (itq=liTQ.begin();itq!=itqfin;itq++,idecal++)
|
|
{TypeQuelconque& tipParticu = (*itq); // pour simplifier
|
|
if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
|
|
switch (tipParticu.EnuTypeQuelconque().EnumTQ())
|
|
{ case POTENTIEL:
|
|
{ SaveResulHyper_W_gene_3D & save_resul = *((SaveResulHyper_W_gene_3D*) saveDon);
|
|
Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
|
|
if (sortie_post)
|
|
tyTQ(1+(*idecal)) = (save_resul.invP)->potentiel;
|
|
else
|
|
tyTQ(1+(*idecal)) = 0.;
|
|
(*idecal)++; break;
|
|
}
|
|
case PARA_LOI_TYPE_MOONEY:
|
|
// ----- cas des coefficients de loi
|
|
{ SaveResulHyper_W_gene_3D & save_resul = *((SaveResulHyper_W_gene_3D*) saveDon);
|
|
Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
|
|
if (save_resul.para_loi != NULL) {tyTQ(1+(*idecal))= *(save_resul.para_loi);}
|
|
else {tyTQ(1+(*idecal)).Zero();};
|
|
(*idecal)++; break;
|
|
}
|
|
|
|
default: ;// on ne fait rien
|
|
};
|
|
|
|
};
|
|
};
|
|
|
|
// récupération et création de la liste de tous les grandeurs particulières
|
|
// ces grandeurs sont ajoutées à la liste passées en paramètres
|
|
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
|
|
void Hyper_W_gene_3D::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
|
|
{// ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas
|
|
Tableau <double> tab_1(1);
|
|
Tab_Grandeur_scalaire_double grand_courant(tab_1);
|
|
// on ne propose des grandeurs que si elles ont été stockée
|
|
// non car le choix est un choix global, par contre les valeurs ne seront différentes de 0 que si sortie_post
|
|
// est activé !!
|
|
// if (sortie_post)
|
|
{ // def d'un type quelconque représentatif à chaque grandeur
|
|
// a priori ces grandeurs sont défini aux points d'intégration identique à la contrainte par exemple
|
|
// enu_ddl_type_pt est définit dans la loi Abtraite générale
|
|
// on n'ajoute que si sortie_post est vraie, sinon aucune grandeur n'est sauvegardé, donc on ne peut
|
|
// plus y accèder
|
|
//on regarde si ce type d'info existe déjà: si oui on augmente la taille du tableau, si non on crée
|
|
// cas du potentiel
|
|
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
|
|
for (itq=liTQ.begin();itq!=itqfin;itq++)
|
|
if ((*itq).EnuTypeQuelconque() == POTENTIEL)
|
|
{ Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
|
|
int taille = tyTQ.Taille()+1;
|
|
tyTQ.Change_taille(taille); nexistePas = false;
|
|
};
|
|
if (nexistePas)
|
|
{TypeQuelconque typQ1(POTENTIEL,enu_ddl_type_pt,grand_courant);
|
|
liTQ.push_back(typQ1);
|
|
};
|
|
};
|
|
};
|
|
// $$$ cas des paramètres de la loi de comportement
|
|
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
|
|
for (itq=liTQ.begin();itq!=itqfin;itq++)
|
|
if ((*itq).EnuTypeQuelconque() == PARA_LOI_TYPE_MOONEY)
|
|
{Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
|
|
int taille = tyTQ.Taille()+1;
|
|
tyTQ.Change_taille(taille); nexistePas = false;
|
|
tyTQ(taille).Change_taille(nb_para_loi);
|
|
};
|
|
if (nexistePas)
|
|
{Vecteur tens(nb_para_loi); //
|
|
Tab_Grandeur_Vecteur gr_vec_para(tens,1);
|
|
// def d'un type quelconque représentatif
|
|
TypeQuelconque typQ(PARA_LOI_TYPE_MOONEY,EPS11,gr_vec_para);
|
|
liTQ.push_back(typQ);
|
|
};
|
|
};
|
|
|
|
};
|
|
|
|
// activation du stockage de grandeurs quelconques qui pourront ensuite être récupéré
|
|
// via le conteneur SaveResul, si la grandeur n'existe pas ici, aucune action
|
|
void Hyper_W_gene_3D::Activation_stockage_grandeurs_quelconques(list <EnumTypeQuelconque >& listEnuQuelc)
|
|
{ // récup de la liste de stockage
|
|
list <EnumTypeQuelconque >& listlocale = ListQuelc_mis_en_acces_localement();
|
|
|
|
// on parcours la liste des grandeurs à activer
|
|
// et on remplit la liste locale
|
|
list <EnumTypeQuelconque >::iterator il, ilfin = listEnuQuelc.end();
|
|
for (il = listEnuQuelc.begin();il != ilfin; il++)
|
|
// for (EnumTypeQuelconque enu : listEnuQuelc)
|
|
// on ne remplit que s'il s'agit d'une grandeur qui peut-être accessible
|
|
{switch (*il)
|
|
{case POTENTIEL: case PARA_LOI_TYPE_MOONEY:
|
|
listlocale.push_back(*(il));
|
|
break;
|
|
default: ; // pour les autres cas on ne fait rien
|
|
};
|
|
};
|
|
|
|
// on supprime les doublons localement
|
|
listlocale.sort(); // on ordonne la liste
|
|
listlocale.unique(); // suppression des doublons
|
|
|
|
};
|
|
|
|
// insertion des conteneurs ad hoc concernant le stockage de grandeurs quelconques
|
|
// passée en paramètre, dans le save result: ces conteneurs doivent être valides
|
|
// c-a-d faire partie de listdeTouslesQuelc_dispo_localement
|
|
void Hyper_W_gene_3D::Insertion_conteneur_dans_save_result(SaveResul * a)
|
|
{
|
|
// -- autre stockage éventuel en fonction des grandeurs quelconques demandées par d'autres lois
|
|
// on va regarder s'il faut activer ou non invP et invP_t, pour le cas
|
|
// de la récupération de grandeur quelconque au moment de l'utilisation de la loi
|
|
// récup de la liste de stockage
|
|
list <EnumTypeQuelconque >& listlocale = ListQuelc_mis_en_acces_localement();
|
|
// def
|
|
Hyper_W_gene_3D::SaveResulHyper_W_gene_3D& sav = *((Hyper_W_gene_3D::SaveResulHyper_W_gene_3D*) a);
|
|
|
|
List_io <EnumTypeQuelconque>::iterator jk,jkfin = listlocale.end();
|
|
for (jk=listlocale.begin();jk != jkfin;jk++)
|
|
// for (EnumTypeQuelconque enu : listQuelc_dispo_localement)
|
|
// { switch (enu)
|
|
{EnumTypeQuelconque enu = *jk;
|
|
switch (enu)
|
|
{ case POTENTIEL: // cas on des grandeurs sont demandées
|
|
{ if (sav.invP == NULL) // def si ils n'existent pas
|
|
{sav.invP = new Invariantpost3D();
|
|
sav.invP_t = new Invariantpost3D();
|
|
if ((Permet_affichage() > 0 ) && (sortie_post == 0))
|
|
{cout << "\n *** warning: on demande pour une loi hyper elastique, l'utilisation de la grandeur "
|
|
<< NomTypeQuelconque(enu) << " or on n'a pas demande sortie_post_ different de 0 "
|
|
<< " on met arbitrairement sortie_post_ = 1 ";
|
|
};
|
|
if (sortie_post == 0)
|
|
sortie_post = 1;
|
|
};
|
|
// on crée le conteneur ad hoc pour le passage d'info
|
|
// def d'un conteneur de grandeurs quelconques, initialisée à 0
|
|
Grandeur_scalaire_double grand_courant(0.);
|
|
TypeQuelconque typQ1(enu,EPS11,grand_courant);
|
|
sav.map_type_quelconque[enu]=(typQ1);
|
|
break;
|
|
}
|
|
case PARA_LOI_TYPE_MOONEY: // cas on des grandeurs sont demandées
|
|
{ if (sav.para_loi == NULL) // def si ils n'existent pas
|
|
{sav.para_loi = new Vecteur; // la dimension sera ensuite précisée
|
|
sav.para_loi_t = new Vecteur; // dans la loi dérivée
|
|
};
|
|
// on crée le conteneur ad hoc pour le passage d'info
|
|
// def d'un conteneur de grandeurs quelconques, initialisée à 0
|
|
Grandeur_Vecteur grand_courant(3);grand_courant.ConteneurVecteur()->Zero();
|
|
TypeQuelconque typQ1(enu,EPS11,grand_courant);
|
|
sav.map_type_quelconque[enu]=(typQ1);
|
|
break;
|
|
}
|
|
|
|
default:
|
|
cout << "\n *** erreur on demande l'acces a : "
|
|
<< NomTypeQuelconque(enu)
|
|
<< " or celui-ci n'est pas dispo pour la loi ";
|
|
this->Affiche();
|
|
cout << " revoir la mise en donnees ! " << endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
};
|
|
|
|
|
|
// --------------- méthodes internes -------------
|
|
// affichage et definition interactive des commandes particulières à chaques lois
|
|
void Hyper_W_gene_3D::Info_commande_LoisDeComp_hyper3D(UtilLecture& entreePrinc)
|
|
{ ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
|
|
sort << "\n# il est egalement possible de stocker et donc de recuperer en post-traitement, differentes grandeurs "
|
|
<< "\n# specifiquement utilises : ici la valeur du potentiel ,"
|
|
<< "\n# pour cela sur la derniere ligne des donnees de la loi il "
|
|
<< "\n# faut la presence du mot cle: sortie_post_ suivi de 1, par defaut sa valeur est 0 "
|
|
<< "\n# ex: sortie_post_ 1 ";
|
|
};
|
|
|
|
|
|
// calcul des invariants et de leurs variations premières
|
|
void Hyper_W_gene_3D::InvariantsEtVar1(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
|
|
,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
|
|
,const double& jacobien_0,const double& jacobien)
|
|
{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
|
|
const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); // "
|
|
const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); // "
|
|
const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // "
|
|
// tout d'abord on définit deux tenseurs qui nous seront utiles par la suite
|
|
const Tenseur3HH & B_HH = gijHH_0; // le tenseur B
|
|
B_HB = B_HH * gijBB_tdt; // le tenseur B en mixte
|
|
BB_HB = B_HB * B_HB;
|
|
BB_HH = BB_HB * gijHH_tdt; // le tenseur B . B en deux fois contra
|
|
// BB_HH = B_HB * B_HH ;
|
|
// calcul des trois invariants traces primaires de B
|
|
I_B= B_HB.Trace(); //gijBB_tdt && B_HH ; //gijHH_0;
|
|
I_BB= BB_HB.Trace(); //gijBB_tdt && BB_HH; //((gijBB_tdt * gijHH_0) * gijBB_tdt) && gijHH_0;
|
|
// I_BBB= gijBB_tdt && ((BB_HH * gijBB_tdt) * B_HH); // ((((gijBB_tdt * gijHH_0) * gijBB_tdt) * gijHH_0) * gijBB_tdt) && gijHH_0;
|
|
// calcul des trois invariants en I de B: I_B, II_B, III_B
|
|
II_B=0.5*(I_B*I_B-I_BB);
|
|
const double untiers=1./3.;
|
|
// III_B=untiers*I_BBB - 0.5*I_B*I_BB + I_B*I_B*I_B/6;
|
|
V= jacobien/jacobien_0; // variation relative de volume
|
|
III_B=V*V; // calcul plus rapide que le précédent
|
|
// calcul des trois invariants en J
|
|
if (III_B <= ConstMath::petit)
|
|
{ if (Permet_affichage() >= 1)
|
|
cout << "\n **** erreur, le volume final sur volume initial est quasi-nulle ! III_B=" << (III_B)
|
|
<< " on limite (V**(2)) a : " << ConstMath::petit;
|
|
if (Permet_affichage() >= 5)
|
|
cout << "\n Hyper_W_gene_3D::Invariants_et_var1( ";
|
|
III_B = ConstMath::petit;
|
|
};
|
|
// variables intermédiaires
|
|
J3_puiss_untiers=pow(III_B,untiers);
|
|
unSurJ3_puissuntiers=1./J3_puiss_untiers;
|
|
unSurJ3_puissuntiers2=unSurJ3_puissuntiers*unSurJ3_puissuntiers;
|
|
unSurJ3_puissuntiers4=unSurJ3_puissuntiers2*unSurJ3_puissuntiers2;
|
|
// les J
|
|
J_r(1)= I_B*unSurJ3_puissuntiers;
|
|
J_r(2)=II_B*unSurJ3_puissuntiers2;
|
|
J_r(3)=III_B;
|
|
// V = sqrt(III_B); // autre calcul de la variation relative de volume
|
|
// calcul des variations premières des I_r par rapports aux composantes epsBB d'Almansi
|
|
// ( en fait on ne les calcul pas pour optimiser le temps)
|
|
// d_I_B_epsBB_HH=2.* B_HH ; //gijHH_0;
|
|
// d_II_B_epsBB_HH= 2.* (I_B * B_HH - BB_HH);
|
|
// // - 2. * ((gijHH_0 * gijBB_tdt) * gijHH_0);
|
|
// d_III_B_epsBB_HH=(2.*III_B)*gijHH_tdt;
|
|
// calcul des variation premières des J_r par rapports aux composantes epsBB d'Almansi
|
|
|
|
d_J_r_epsBB_HH(1)=(2.*unSurJ3_puissuntiers)*(B_HH - (untiers*I_B)*gijHH_tdt);
|
|
d_J_r_epsBB_HH(2)=(2.*unSurJ3_puissuntiers2)*(I_B * B_HH - BB_HH - (2.*untiers*II_B)*gijHH_tdt);
|
|
d_J_r_epsBB_HH(3)=(2.*III_B)*gijHH_tdt;
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 4)
|
|
{ cout << "\nHyper_W_gene_3D::InvariantsEtVar1(..."
|
|
<< "\n I_B= " << I_B << ", II_B= "<< II_B <<", III_B= "<< III_B
|
|
<< "\n J_r(1)= " << J_r(1) <<", J_r(2)= " << J_r(2) << ", J_r(3)= " << J_r(3)
|
|
<< "\n d_J_r_epsBB_HH(1)= " << d_J_r_epsBB_HH(1)
|
|
<< ", d_J_r_epsBB_HH(2)= " << d_J_r_epsBB_HH(2)
|
|
<< ", d_J_r_epsBB_HH(3)= " << d_J_r_epsBB_HH(3);
|
|
};
|
|
#endif
|
|
};
|
|
|
|
// calcul des invariants et de leurs variations premières et secondes
|
|
void Hyper_W_gene_3D::InvariantsEtVar2(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
|
|
,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
|
|
,const double& jacobien_0,const double& jacobien)
|
|
{ // calcul des invariants et de leurs variations premières
|
|
Invariants_et_var1(gijBB_0_,gijHH_0_,gijBB_tdt_,gijHH_tdt_,jacobien_0,jacobien);
|
|
|
|
const double untiers=1./3.;
|
|
// calcul des variations secondes
|
|
const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
|
|
const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); // "
|
|
const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); // "
|
|
const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // "
|
|
const Tenseur3HH & B_HH = gijHH_0; // récup du tenseur B
|
|
// le tenseur BB est déjà calculé
|
|
// calcul de variables intermédiaires
|
|
IxI_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,gijHH_tdt);
|
|
IxbarreI_HHHH=Tenseur3HHHH::Prod_tensoriel_barre(gijHH_tdt,gijHH_tdt);
|
|
IxB_HHHH=Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,B_HH);
|
|
BxI_HHHH=Tenseur3HHHH::Prod_tensoriel(B_HH,gijHH_tdt);
|
|
|
|
// variations secondes des I_r
|
|
//(finalement on ne les calcule pas pour optimiser)
|
|
//// d_I_B_eps2BB_HHHH.Inita(0.);
|
|
//// d_II_B_eps2BB_HHHH = 4.*(Tenseur3HHHH::Prod_tensoriel(B_HH,B_HH)
|
|
//// - Tenseur3HHHH::Prod_tensoriel_barre(B_HH,B_HH));
|
|
//// d_III_B_eps2BB_HHHH = (4. * III_B) *
|
|
//// ( Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,gijHH_tdt)
|
|
//// - Tenseur3HHHH::Prod_tensoriel_barre(gijHH_tdt,gijHH_tdt));
|
|
|
|
// variations secondes des J_r
|
|
// d_J_1_eps2BB_HHHH = (-2.*untiers*unSurJ3_puissuntiers4)*
|
|
// Tenseur3HHHH::Prod_tensoriel(d_I_B_epsBB_HH,d_III_B_epsBB_HH)
|
|
// - (untiers * I_B * unSurJ3_puissuntiers4) * d_III_B_eps2BB_HHHH;
|
|
|
|
d_J_1_eps2BB_HHHH = (-4.*untiers*unSurJ3_puissuntiers)*
|
|
( BxI_HHHH + IxB_HHHH - I_B*(untiers*IxI_HHHH + IxbarreI_HHHH) );
|
|
|
|
d_J_2_eps2BB_HHHH = (-8. * untiers * unSurJ3_puissuntiers2 )*
|
|
(I_B* BxI_HHHH
|
|
- Tenseur3HHHH::Prod_tensoriel(BB_HH,gijHH_tdt)
|
|
- (2.*untiers)*II_B * IxI_HHHH
|
|
)
|
|
+(4. * unSurJ3_puissuntiers2 )*
|
|
(Tenseur3HHHH::Prod_tensoriel(B_HH,B_HH)
|
|
- Tenseur3HHHH::Prod_tensoriel_barre(B_HH,B_HH)
|
|
-(2.*untiers )*(I_B * IxB_HHHH
|
|
- Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,BB_HH))
|
|
+(2.*untiers * II_B)*IxbarreI_HHHH
|
|
);
|
|
|
|
d_J_3_eps2BB_HHHH = (4. * III_B) * ( IxI_HHHH - IxbarreI_HHHH);
|
|
|
|
};
|
|
|
|
//-------------- vérification de la dérivée avec une dérivée numérique -------------
|
|
|
|
|
|
// calcul des dérivées numériques premières
|
|
void Hyper_W_gene_3D::Calcul_derivee_numerique(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
|
|
,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
|
|
,const double& jacobien_0,const double& jacobien)
|
|
{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
|
|
const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); // "
|
|
const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); // "
|
|
const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // "
|
|
|
|
Tenseur3BB gijBBtdt_N; // tenseur modifié
|
|
Tenseur3HH gijHHtdt_N; // idem_0
|
|
const Tenseur3HH & B_HH = gijHH_0;
|
|
// on choisit un incrément
|
|
double delta = ConstMath::unpeupetit;///10.;
|
|
double unSurDelta = 1./delta;
|
|
// sauvegarde des invariants actuelles et de leur variations
|
|
Vecteur J_r_sauve(J_r) ;
|
|
double I_B_sauve=I_B; double II_B_sauve=II_B; double III_B_sauve=III_B;
|
|
Tableau <Tenseur3HH> d_J_r_epsBB_HH_sauve(d_J_r_epsBB_HH);
|
|
Tenseur3HH d_I_B_epsBB_HH_sauve=d_I_B_epsBB_HH;
|
|
Tenseur3HH d_II_B_epsBB_HH_sauve = d_II_B_epsBB_HH;
|
|
Tenseur3HH d_III_B_epsBB_HH_sauve = d_III_B_epsBB_HH;
|
|
// dimensionnement
|
|
Tableau <Tenseur3HH> d_J_r_epsBB_HH_num(d_J_r_epsBB_HH);
|
|
// on va boucler sur les composantes de gijBB
|
|
for (int i=1;i<=3;i++)
|
|
for (int j=1;j<=3;j++)
|
|
{ gijBBtdt_N = gijBB_tdt;
|
|
gijBBtdt_N.Coor(i,j) += delta;
|
|
// en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i)
|
|
// car le tenseur est symétrique
|
|
// on a donc en variation numérique la somme des deux dérivées
|
|
// on définit un coeff multiplicatif qui vaut 1 ou 0.5
|
|
double coef=1.; if (i != j) coef = 0.5;
|
|
gijHHtdt_N = gijBBtdt_N.Inverse();
|
|
double jacobien_N=sqrt(gijBBtdt_N.Det());
|
|
// calcul des invariants
|
|
InvariantsEtVar1(gijBB_0_,gijHH_0_,(const TenseurBB &)gijBBtdt_N,(const TenseurHH &)gijHHtdt_N
|
|
,jacobien_0,jacobien_N);
|
|
for (int r=1;r<=3;r++)
|
|
{double& der_num = d_J_r_epsBB_HH_num(r).Coor(i,j);
|
|
// comme on a la somme des dérivées (i,j) et (j,i) qui doivent être égale
|
|
// on divise par 2 lorsque i != j
|
|
der_num=coef*2.*(J_r(r)-J_r_sauve(r))*unSurDelta; // dérivée numérique
|
|
double& der = d_J_r_epsBB_HH(r).Coor(i,j); // dérivée analytique
|
|
// comparaison avec les valeurs de dérivées analytiques
|
|
bool erreur = false;
|
|
// if (r!=2)
|
|
if (diffpourcent(der_num,der,MaX(Dabs(der_num),Dabs(der)),0.1))
|
|
if (MaX(Dabs(der_num),Dabs(der)) > 1.e-5) //1000.*delta)
|
|
{if (MiN(Dabs(der_num),Dabs(der)) == 0.)
|
|
{if ( MaX(Dabs(der_num),Dabs(der)) > 5.e-6 ) erreur = true;} //50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
if (erreur)
|
|
{
|
|
// calcul des dérivées des I pour voir
|
|
double derI_B_num = coef * 2. * (I_B-I_B_sauve)*unSurDelta;
|
|
double derII_B_num = coef * 2. * (II_B-II_B_sauve)*unSurDelta;
|
|
double derIII_B_num = coef * 2. * (III_B-III_B_sauve)*unSurDelta;
|
|
double& derI_B = d_I_B_epsBB_HH_sauve.Coor(i,j);
|
|
double& derII_B = d_II_B_epsBB_HH_sauve.Coor(i,j);
|
|
double& derIII_B = d_III_B_epsBB_HH_sauve.Coor(i,j);
|
|
|
|
// calcul de dérivées partielles pour voir
|
|
// double dernumpartielle = pow(III_B,-1./3.)*2.*(I_B-I_B_sauve)*unSurDelta;
|
|
// double deranapartielle = 2.*J_r(3)*B_HH(i,j);
|
|
// double dernumpartielle = I_B * 2. *(pow(III_B,-1./3.)-pow(III_B_sauve,-1./3.))*unSurDelta;
|
|
// double deranapartielle = -2./3.*pow(III_B_sauve,-1./3.)*gijHH_tdt(i,j);
|
|
double dernumpartielle = 2. *(III_B-III_B_sauve)*unSurDelta;
|
|
double deranapartielle = 2.*III_B_sauve*gijHH_tdt(i,j);
|
|
cout << "\n dernumpartielle= " << dernumpartielle << ", deranapartielle= " << deranapartielle;
|
|
//
|
|
cout << "\n erreur dans le calcul analytique des derivees de l'invariants: " << r
|
|
<< "\n der_num= " << der_num << " der= " << der <<" der_N=" << d_J_r_epsBB_HH(r)(i,j)
|
|
<< " ietj= " << i << " " << j
|
|
<< "\n derI_B_num= " << derI_B_num << " derI_B= " << derI_B
|
|
<< " derII_B_num= " << derII_B_num << " derII_B= " << derII_B
|
|
<< " derIII_B_num= " << derIII_B_num << " derIII_B= " << derIII_B ;
|
|
cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
|
|
cout << "\n un caractere ";
|
|
string toto;
|
|
toto=lect_chaine();
|
|
// Sortie(1);
|
|
};
|
|
};
|
|
};
|
|
// retour aux valeurs sauvegardées
|
|
J_r = J_r_sauve;
|
|
I_B=I_B_sauve; II_B = II_B_sauve; III_B = III_B_sauve;
|
|
for (int e=1;e<=3;e++)
|
|
d_J_r_epsBB_HH(e)=d_J_r_epsBB_HH_sauve(e);
|
|
};
|
|
|
|
// calcul des dérivées numériques premières et seconde
|
|
void Hyper_W_gene_3D::Calcul_derivee_numerique2(const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
|
|
,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
|
|
,const double& jacobien_0,const double& jacobien)
|
|
{ const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) &gijBB_0_); // passage en dim 3 explicit
|
|
const Tenseur3BB & gijBB_tdt = *((Tenseur3BB*) &gijBB_tdt_); // "
|
|
const Tenseur3HH & gijHH_0 = *((Tenseur3HH*) &gijHH_0_); // "
|
|
const Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // "
|
|
|
|
Tenseur3BB gijBBtdt_N; // tenseur modifié
|
|
Tenseur3HH gijHHtdt_N; // idem_0
|
|
const Tenseur3HH & B_HH = gijHH_0;
|
|
double delta = ConstMath::unpeupetit*10.;
|
|
double unSurDelta = 1./delta;
|
|
// sauvegarde des invariants actuelles et de leur variations
|
|
Vecteur J_r_sauve(J_r) ;
|
|
double I_B_sauve=I_B; double II_B_sauve=II_B; double III_B_sauve=III_B;
|
|
Tableau <Tenseur3HH> d_J_r_epsBB_HH_sauve(d_J_r_epsBB_HH);
|
|
// idem pour la partie variation seconde
|
|
Tenseur3HHHH d_J_1_eps2BB_HHHH_sauve = d_J_1_eps2BB_HHHH;
|
|
Tenseur3HHHH d_J_2_eps2BB_HHHH_sauve = d_J_2_eps2BB_HHHH;
|
|
Tenseur3HHHH d_J_3_eps2BB_HHHH_sauve = d_J_3_eps2BB_HHHH;
|
|
|
|
// dimensionnement
|
|
Tableau <Tenseur3HH> d_J_r_epsBB_HH_num(d_J_r_epsBB_HH);
|
|
// on va boucler sur les composantes de gijBB
|
|
for (int i=1;i<=3;i++)
|
|
for (int j=1;j<=3;j++)
|
|
{ gijBBtdt_N = gijBB_tdt;
|
|
gijBBtdt_N.Coor(i,j) += delta;
|
|
// en fait dans l'opération précédente on a modifier les termes (i,j) et (j,i)
|
|
// car le tenseur est symétrique
|
|
// on a donc en variation numérique la somme des deux dérivées
|
|
// on définit un coeff multiplicatif qui vaut 1 ou 0.5
|
|
double coef=1.; if (i != j) coef = 0.5;
|
|
gijHHtdt_N = gijBBtdt_N.Inverse();
|
|
double jacobien_N=sqrt(gijBBtdt_N.Det());
|
|
// calcul des invariants
|
|
InvariantsEtVar1(gijBB_0_,gijHH_0_,(const TenseurBB &)gijBBtdt_N,(const TenseurHH &)gijHHtdt_N
|
|
,jacobien_0,jacobien_N);
|
|
for (int r=1;r<=3;r++)
|
|
{double& der_num = d_J_r_epsBB_HH_num(r).Coor(i,j);
|
|
der_num=coef*2.*(J_r(r)-J_r_sauve(r))*unSurDelta; // dérivée numérique
|
|
double& der = d_J_r_epsBB_HH(r).Coor(i,j); // dérivée analytique
|
|
// comparaison avec les valeurs de dérivées analytiques
|
|
bool erreur = false;
|
|
// if (r!=2)
|
|
if (diffpourcent(der_num,der,MaX(Dabs(der_num),Dabs(der)),0.01))
|
|
if (MaX(Dabs(der_num),Dabs(der)) > 1000.*delta)
|
|
{if (MiN(Dabs(der_num),Dabs(der)) == 0.)
|
|
{if ( MaX(Dabs(der_num),Dabs(der)) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
if (erreur)
|
|
{
|
|
// calcul des dérivées des I pour voir
|
|
cout << "\n erreur dans le calcul analytique des derivees de l'invariants: " << r
|
|
<< "\n der_num= " << der_num << " der= " << der << " ietj= " << i << " " << j
|
|
<< " J_r= " << J_r(r);
|
|
cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
|
|
cout << "\n un caractere ";
|
|
string toto;
|
|
toto=lect_chaine();
|
|
};
|
|
};
|
|
// cas des variations secondes
|
|
// calcul des dérivées secondes analytiques
|
|
InvariantsEtVar2(gijBB_0_,gijHH_0_,(const TenseurBB &)gijBBtdt_N,(const TenseurHH &)gijHHtdt_N
|
|
,jacobien_0,jacobien_N);
|
|
// calcul des dérivées numériques et comparaisons
|
|
for (int k=1;k<=3;k++)
|
|
for (int l=1;l<=3;l++)
|
|
{ // cas de d_J_1_eps2BB_HHHH
|
|
double der2num = coef * 2.*(d_J_r_epsBB_HH(1)(k,l) - d_J_r_epsBB_HH_sauve(1)(k,l) )*unSurDelta;
|
|
double der2ana = 0.5*(d_J_1_eps2BB_HHHH(k,l,i,j) + d_J_1_eps2BB_HHHH(k,l,j,i));
|
|
bool erreur = false;
|
|
if ((Dabs(d_J_r_epsBB_HH(1)(k,l))> 0.8) && (Dabs(der2ana)<30))
|
|
if (diffpourcent(der2num,der2ana,MaX(Dabs(der2num),Dabs(der2ana)),0.01))
|
|
if (MaX(Dabs(der2num),Dabs(der2ana)) > 1000.*delta)
|
|
{if (MiN(Dabs(der2num),Dabs(der2ana)) == 0.)
|
|
{if ( MaX(Dabs(der2num),Dabs(der2ana)) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
if (erreur)
|
|
{
|
|
// calcul des dérivées d'éléments intermédiaires pour voir
|
|
double der2anapartielle=-4./3.*pow(III_B,-1./3.)*(B_HH(k,l)-1./3.*I_B*gijHH_tdt(k,l))*gijHH_tdt(i,j)
|
|
+2.*pow(III_B,-1./3.)*(-2./3.*gijHH_tdt(k,l)*B_HH(i,j)+2./3.*I_B*gijHH_tdt(k,i)*gijHH_tdt(l,j));
|
|
double der2numpartielle=coef*2.*(2.*pow(III_B,-1./3.)*(B_HH(k,l)-1./3.*I_B*gijHHtdt_N(k,l))
|
|
-(2.*pow(III_B_sauve,-1./3.)*(B_HH(k,l)-1./3.*I_B_sauve*gijHH_tdt(k,l)) ))*unSurDelta;
|
|
cout << "\n der2numpartielle= " << der2numpartielle << " der2anapartielle= " << der2anapartielle ;
|
|
|
|
|
|
|
|
//
|
|
cout << "\n erreur dans le calcul analytique des derivees secondes de l'invariants: " << 1
|
|
<< "\n der2num= " << der2num << " der2ana= " << der2ana
|
|
<< " klij= " << k << " " << l << " " << i << " " << j
|
|
<< " d_J_r_epsBB_HH(1)(k,l)= " << d_J_r_epsBB_HH(1)(k,l);
|
|
cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
|
|
cout << "\n un caractere ";
|
|
string toto;
|
|
toto=lect_chaine();
|
|
};
|
|
// d_J_1_eps2BB_HHHH.Change(k,l,i,j,der2num); // a virer
|
|
// cas de d_J_2_eps2BB_HHHH
|
|
der2num = coef * 2.*(d_J_r_epsBB_HH(2)(k,l) - d_J_r_epsBB_HH_sauve(2)(k,l) )*unSurDelta;
|
|
der2ana = 0.5*(d_J_2_eps2BB_HHHH(k,l,i,j) + d_J_2_eps2BB_HHHH(k,l,j,i));
|
|
erreur = false;
|
|
if (diffpourcent(der2num,der2ana,MaX(Dabs(der2num),Dabs(der2ana)),0.01))
|
|
if (MaX(Dabs(der2num),Dabs(der2ana)) > 1000.*delta)
|
|
{if (MiN(Dabs(der2num),Dabs(der2ana)) == 0.)
|
|
{if ( MaX(Dabs(der2num),Dabs(der2ana)) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
if (erreur)
|
|
{
|
|
// calcul des dérivées des I pour voir
|
|
cout << "\n erreur dans le calcul analytique des derivees secondes de l'invariants: " << 2
|
|
<< "\n der2num= " << der2num << " der2ana= " << der2ana
|
|
<< " klij= " << k << " " << l << " " << i << " " << j
|
|
<< " d_J_r_epsBB_HH(2)(k,l)= " << d_J_r_epsBB_HH(2)(k,l);
|
|
cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
|
|
cout << "\n un caractere ";
|
|
string toto;
|
|
// cin >> toto;
|
|
};
|
|
// d_J_2_eps2BB_HHHH.Change(k,l,i,j,der2num); // a virer
|
|
// cas de d_J_3_eps2BB_HHHH
|
|
der2num = coef * 2.*(d_J_r_epsBB_HH(3)(k,l) - d_J_r_epsBB_HH_sauve(3)(k,l) )*unSurDelta;
|
|
der2ana = 0.5*(d_J_3_eps2BB_HHHH(k,l,i,j) + d_J_3_eps2BB_HHHH(k,l,j,i));
|
|
erreur = false;
|
|
if (diffpourcent(der2num,der2ana,MaX(Dabs(der2num),Dabs(der2ana)),0.1))
|
|
if (MaX(Dabs(der2num),Dabs(der2ana)) > 1000.*delta)
|
|
{if (MiN(Dabs(der2num),Dabs(der2ana)) == 0.)
|
|
{if ( MaX(Dabs(der2num),Dabs(der2ana)) > 50.*delta) erreur = true;}
|
|
else erreur = true;
|
|
};
|
|
if (erreur)
|
|
{
|
|
// calcul des dérivées des I pour voir
|
|
cout << "\n erreur dans le calcul analytique des derivees secondes de l'invariants: " << 3
|
|
<< "\n der2num= " << der2num << " der2ana= " << der2ana
|
|
<< " klij= " << k << " " << l << " " << i << " " << j
|
|
<< " d_J_r_epsBB_HH(3)(k,l)= " << d_J_r_epsBB_HH(3)(k,l);
|
|
cout << "\n Hyper_W_gene_3D::Calcul_derivee_numerique(..";
|
|
cout << "\n un caractere ";
|
|
string toto;
|
|
// cin >> toto;
|
|
};
|
|
// d_J_3_eps2BB_HHHH.Change(k,l,i,j,der2num); // a virer
|
|
|
|
};
|
|
};
|
|
// retour aux valeurs sauvegardées
|
|
J_r = J_r_sauve;
|
|
I_B=I_B_sauve; II_B = II_B_sauve; III_B = III_B_sauve;
|
|
for (int e=1;e<=3;e++)
|
|
d_J_r_epsBB_HH(e)=d_J_r_epsBB_HH_sauve(e);
|
|
// pour les dérivées secondes
|
|
d_J_1_eps2BB_HHHH = d_J_1_eps2BB_HHHH_sauve;
|
|
d_J_2_eps2BB_HHHH = d_J_2_eps2BB_HHHH_sauve;
|
|
d_J_3_eps2BB_HHHH = d_J_3_eps2BB_HHHH_sauve;
|
|
};
|