// 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) .
//
// 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 .
//
// For more information, please consult: .
#include "HyperD.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"
// ---------- classe de stockage des grandeurs spécifiques pour la loi ---------
// constructeur par défaut
HyperD::SaveResulHyperD::SaveResulHyperD() :
invP(NULL),invP_t(NULL),map_type_quelconque()
{ // mise à jour de la liste des grandeurs quelconques internes
Mise_a_jour_map_type_quelconque();
// // def du conteneur de grandeurs quelconques, initialisée à 0
// Grandeur_scalaire_double grand_courant(0.);
// {TypeQuelconque typQ1(V_vol,EPS11,grand_courant);
// map_type_quelconque[V_vol]=(typQ1);
// };
// {TypeQuelconque typQ1(COS3PHI_EPS,EPS11,grand_courant);
// map_type_quelconque[COS3PHI_EPS]=(typQ1);
// };
// {TypeQuelconque typQ1(SPHERIQUE_EPS,EPS11,grand_courant);
// map_type_quelconque[SPHERIQUE_EPS]=(typQ1);
// };
// {TypeQuelconque typQ1(Q_EPS,EPS11,grand_courant);
// map_type_quelconque[Q_EPS]=(typQ1);
// };
};
// avec init ou pas
HyperD::SaveResulHyperD::SaveResulHyperD(int sortie_post) :
invP(new Invariantpost3D()),invP_t(new Invariantpost3D())
,map_type_quelconque()
{ // mise à jour de la liste des grandeurs quelconques internes
Mise_a_jour_map_type_quelconque();
// // def du conteneur de grandeurs quelconques, initialisée à 0
// Grandeur_scalaire_double grand_courant(0.);
// {TypeQuelconque typQ1(V_vol,EPS11,grand_courant);
// map_type_quelconque[V_vol]=(typQ1);
// };
// {TypeQuelconque typQ1(COS3PHI_EPS,EPS11,grand_courant);
// map_type_quelconque[COS3PHI_EPS]=(typQ1);
// };
// {TypeQuelconque typQ1(SPHERIQUE_EPS,EPS11,grand_courant);
// map_type_quelconque[SPHERIQUE_EPS]=(typQ1);
// };
// {TypeQuelconque typQ1(Q_EPS,EPS11,grand_courant);
// map_type_quelconque[Q_EPS]=(typQ1);
// };
};
// constructeur de copie
HyperD::SaveResulHyperD::SaveResulHyperD(const SaveResulHyperD& sav):
invP(NULL),invP_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));
};
};
// destructeur
HyperD::SaveResulHyperD::~SaveResulHyperD()
{ if (invP != NULL)
{ delete invP;
delete invP_t;
};
};
// affectation
Loi_comp_abstraite::SaveResul & HyperD::SaveResulHyperD::operator = ( const Loi_comp_abstraite::SaveResul & a)
{ HyperD::SaveResulHyperD& sav = *((HyperD::SaveResulHyperD*) &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;
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 HyperD::SaveResulHyperD::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 HyperD::SaveResulHyperD::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;
};
};
// mise à jour éventuelle 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 HyperD::SaveResulHyperD::Ecriture_base_info(ofstream& sort,const int )
{ sort << "\n dat_sp_hyperelas: ";
if (invP != NULL)
{ sort << " 1 " << (*invP); }// V= " << invP->V << " Qeps= "<< invP->Qeps << " cos3phi= "<< invP->cos3phi << " "; }// V,Qeps,cos3phi;
else { sort << " 0 ";};
};
// affichage à l'écran des infos
void HyperD::SaveResulHyperD::Affiche() const
{ cout << "\n data specifiques hyperelastiques : " ;
if (invP != NULL)
{ cout << (*invP); } //" V= " << invP->V << " Qeps= "<< invP->Qeps << " cos3phi= "<< invP->cos3phi << " ";}
else
{ cout << " aucun data sauvegarde "; };
};
// mise à jour des informations transitoires
void HyperD::SaveResulHyperD::TdtversT()
{ if (invP != NULL) { (*invP_t) = (*invP);};
// mise à jour de la liste des grandeurs quelconques internes
Mise_a_jour_map_type_quelconque();
};
void HyperD::SaveResulHyperD::TversTdt()
{ if (invP != NULL) { (*invP) = (*invP_t);};
// mise à jour de la liste des grandeurs quelconques internes
Mise_a_jour_map_type_quelconque();
};
// mise à jour de la liste des grandeurs quelconques internes
void HyperD::SaveResulHyperD::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++)
// for (auto ent1 : map_type_quelconque) // ent1.first est la première clée
// { switch (ent1.first)
{ switch ((*il).first)
{ // ----- cos_3phi
case COS3PHI_EPS:
{ Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) (*il).second.Grandeur_pointee()); // pour simplifier
(*tyTQ.ConteneurDouble()) = invP->cos3phi;
break;
}
// ----- Q_EPS
case Q_EPS:
{ Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) (*il).second.Grandeur_pointee()); // pour simplifier
(*tyTQ.ConteneurDouble()) = invP->Qeps;
break;
}
// ----- V_vol
case V_vol:
{ Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) (*il).second.Grandeur_pointee()); // pour simplifier
(*tyTQ.ConteneurDouble()) = invP->V;
break;
}
case POTENTIEL:
{ Grandeur_scalaire_double& tyTQ= *((Grandeur_scalaire_double*) (*il).second.Grandeur_pointee()); // pour simplifier
(*tyTQ.ConteneurDouble()) = invP->potentiel;
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 HyperD ---------
// Constructeur par defaut
//++template
//++HyperD::HyperD () :
HyperD::HyperD () :
avec_phase(false) // par defaut sans phase
,avec_regularisation(false),fact_regularisation(ConstMath::pasmalpetit)
,I_xbarre_I_HHHH()
,ex_impli_hyper(NULL),ex_expli_tdt_hyper(NULL),ex_umat_hyper(NULL)
,sortie_post(0)
// ,invariant_t()
{// --- on remplit avec les grandeurs succeptible d'être utilisées
// acces à listdeTouslesQuelc_dispo_localement
list & list_tousQuelc = ListdeTouslesQuelc_dispo_localement();
list_tousQuelc.push_back(COS3PHI_EPS);
list_tousQuelc.push_back(Q_EPS);
list_tousQuelc.push_back(V_vol);
list_tousQuelc.push_back(POTENTIEL);
// on supprime les doublons localement
list_tousQuelc.sort(); // on ordonne la liste
list_tousQuelc.unique(); // suppression des doublons
}
// Constructeur utile si l'identificateur du nom de la loi
// de comportement, la dimension et le paramètre phase sont connus
//++template
HyperD::HyperD (Enum_comp id_compor,Enum_categorie_loi_comp categorie_comp,int dimension,bool avec_ph) :
Loi_comp_abstraite(id_compor,categorie_comp,dimension),avec_phase(avec_ph),avec_regularisation(false)
,fact_regularisation(ConstMath::pasmalpetit),I_xbarre_I_HHHH()
,ex_impli_hyper(NULL),ex_expli_tdt_hyper(NULL),ex_umat_hyper(NULL)
,sortie_post(0)
// ,invariant_t()
{// --- on remplit avec les grandeurs succeptible d'être utilisées
// acces à listdeTouslesQuelc_dispo_localement
list & list_tousQuelc = ListdeTouslesQuelc_dispo_localement();
list_tousQuelc.push_back(COS3PHI_EPS);
list_tousQuelc.push_back(Q_EPS);
list_tousQuelc.push_back(V_vol);
list_tousQuelc.push_back(POTENTIEL);
// on supprime les doublons localement
list_tousQuelc.sort(); // on ordonne la liste
list_tousQuelc.unique(); // suppression des doublons
}
// Constructeur utile si l'identificateur du nom de la loi
// de comportement, la dimension et le paramètre phase sont connus
//++template
//++HyperD::HyperD (char* nom,int dimension,bool avec_ph) :
HyperD::HyperD (char* nom,Enum_categorie_loi_comp categorie_comp,int dimension,bool avec_ph) :
Loi_comp_abstraite(nom,categorie_comp,dimension),avec_phase(avec_ph),avec_regularisation(false)
,fact_regularisation(ConstMath::pasmalpetit),I_xbarre_I_HHHH()
,ex_impli_hyper(NULL),ex_expli_tdt_hyper(NULL),ex_umat_hyper(NULL)
,sortie_post(0)
// ,invariant_t()
{// --- on remplit avec les grandeurs succeptible d'être utilisées
// acces à listdeTouslesQuelc_dispo_localement
list & list_tousQuelc = ListdeTouslesQuelc_dispo_localement();
list_tousQuelc.push_back(COS3PHI_EPS);
list_tousQuelc.push_back(Q_EPS);
list_tousQuelc.push_back(V_vol);
list_tousQuelc.push_back(POTENTIEL);
// on supprime les doublons localement
list_tousQuelc.sort(); // on ordonne la liste
list_tousQuelc.unique(); // suppression des doublons
}
// DESTRUCTEUR :
//++template
//++HyperD::~HyperD ()
HyperD::~HyperD ()
{}
// constructeur de copie
//++template
//++HyperD::HyperD (const HyperD & a) :
HyperD::HyperD (const HyperD & a) :
Loi_comp_abstraite (a),avec_phase(a.avec_phase),avec_regularisation(a.avec_regularisation)
,fact_regularisation(a.fact_regularisation),I_xbarre_I_HHHH()
,ex_impli_hyper(NULL),ex_expli_tdt_hyper(NULL),ex_umat_hyper(NULL)
,sortie_post(a.sortie_post)
// ,invariant_t(a.invariant_t)
{}
HyperD::SaveResul * HyperD::New_et_Initialise()
{ // création des stockages internes éventuels en fonction de sortie_post
SaveResulHyperD* pt = new SaveResulHyperD(this->sortie_post);
// insertion éventuelle de conteneurs de grandeurs quelconque
Insertion_conteneur_dans_save_result(pt);
// retour
return pt;
};
// 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 HyperD::Grandeur_particuliere
(bool ,List_io& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list& 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::iterator itq,itqfin=liTQ.end();
list::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())
{ // ----- cos_3phi
case COS3PHI_EPS:
{ SaveResulHyperD & save_resul = *((SaveResulHyperD*) saveDon);
Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
if (sortie_post)
tyTQ(1+(*idecal)) = (save_resul.invP)->cos3phi;
else
tyTQ(1+(*idecal)) = 0.;
(*idecal)++; break;
////------- debug
//if (Dabs((save_resul.invP)->cos3phi) > 1.)
// cout << "\n debug HyperD::Grandeur_particuliere *** (save_resul.invP)->cos3phi= "<< (save_resul.invP)->cos3phi << endl;
////------ fin debug
}
// ----- Q_EPS
case Q_EPS:
{ SaveResulHyperD & save_resul = *((SaveResulHyperD*) saveDon);
Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
if (sortie_post)
tyTQ(1+(*idecal)) = (save_resul.invP)->Qeps;
else
tyTQ(1+(*idecal)) = 0.;
(*idecal)++; break;
}
// ----- V_vol
case V_vol:
{ SaveResulHyperD & save_resul = *((SaveResulHyperD*) saveDon);
Tab_Grandeur_scalaire_double& tyTQ= *((Tab_Grandeur_scalaire_double*) (*itq).Grandeur_pointee()); // pour simplifier
if (sortie_post)
tyTQ(1+(*idecal)) = (save_resul.invP)->V;
else
tyTQ(1+(*idecal)) = 0.;
(*idecal)++; break;
}
case POTENTIEL:
{ SaveResulHyperD & save_resul = *((SaveResulHyperD*) 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;
}
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 HyperD::ListeGrandeurs_particulieres(bool absolue,List_io& liTQ) const
{// ici on est en 3D et les grandeurs sont par principe en absolue, donc la variable absolue ne sert pas
Tableau 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 de l'angle de phase
{List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == COS3PHI_EPS)
{ 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(COS3PHI_EPS,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ1);
};
};
// $$$ cas de Qeps
{List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == Q_EPS)
{ 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(Q_EPS,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ1);
};
};
// cas de V
{List_io::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == V_vol)
{ 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(V_vol,enu_ddl_type_pt,grand_courant);
liTQ.push_back(typQ1);
};
};
// cas du potentiel
{List_io::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);
};
};
};
};
// 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 HyperD::Activation_stockage_grandeurs_quelconques(list & listEnuQuelc)
{ // récup de la liste de stockage
list & listlocale = ListQuelc_mis_en_acces_localement();
// on parcours la liste des grandeurs à activer
// et on remplit la liste locale
list ::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 COS3PHI_EPS : case Q_EPS: case V_vol: case POTENTIEL:
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 HyperD::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 & listlocale = ListQuelc_mis_en_acces_localement();
// def
HyperD::SaveResulHyperD& sav = *((HyperD::SaveResulHyperD*) a);
List_io ::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: case V_vol: case COS3PHI_EPS: case SPHERIQUE_EPS: case Q_EPS: // 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 ((ParaGlob::NiveauImpression() > 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;
}
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);
};
};
};
// ========== codage des METHODES VIRTUELLES protegees:================
// calcul des contraintes a t+dt
//++template
//++void HyperD::Calcul_SigmaHH
void HyperD::Calcul_SigmaHH
( TenseurHH & ,TenseurBB& ,DdlElement & tab_ddl,
TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& , TenseurBB & epsBB_,
TenseurBB & ,
TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau & d_gijBB_,
double& jacobien_0,double& jacobien_,TenseurHH & sigHH
,EnergieMeca & energ,const EnergieMeca & ,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Expli_t_tdt& ex )
{
#ifdef MISE_AU_POINT
if (tab_ddl.NbDdl() != d_gijBB_.Taille())
{ cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_ !\n";
cout << " HyperD::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
// attribution des pointeurs de métriques pour l'utilisation par les potentiels
ex_impli_hyper=NULL;ex_expli_tdt_hyper=&ex;ex_umat_hyper=NULL;
int dim_tens = gijBB_.Dimension(); // récup de la dimension des tenseurs
//++ TensBH epsBH; // la déformation en mixte
//++ TensBH * IdGBH; // l'identitée
TenseurBH* epBH; // la déformation en mixte
epBH = NevezTenseurBH(dim_tens);TenseurBH & epsBH = *(epBH);
TenseurBH * IdGBH; // l'identitée
A_i a_i; // définition des a_i
Invariant invariant; // def des invariants
// Calcul des invariants, de epsBH,
// retour de IdGBH qui pointe sur le bon tenseur
IdGBH = Invariants (epsBB_,gijBB_,gijHH_,jacobien_0,jacobien_,invariant,epsBH);
energ.Inita(0.); // init des énergies
if (!avec_phase)
{ // cas sans phase
// calcul du potentiel
PotenSansPhaseSansVar potenSansPhaseSansVar = Potentiel(invariant,jacobien_0);
// calcul des coefficients a_i
AAA_i(potenSansPhaseSansVar,invariant,jacobien_,a_i);
// en fait Potentiel_et_var calcul E*racine(g), mais pour l'énergie il faut E
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenSansPhaseSansVar.E/(invariant.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenSansPhaseSansVar.Ks;
}
else
{ // cas avec la phase
// calcul du potentiel
PotenAvecPhaseSansVar potenAvecPhaseSansVar=PotentielPhase(invariant,jacobien_0);
// calcul des coefficients a_i
AAA_iPhase(potenAvecPhaseSansVar,invariant,jacobien_,a_i);
// en fait Potentiel_et_var calcul E*racine(g), mais pour l'énergie il faut E
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenAvecPhaseSansVar.E/(invariant.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenAvecPhaseSansVar.Ks;
}
// calcul des contraintes finales
//++ TensBH sigBH = a_i.a_0 * (*IdGBH) + a_i.a_1 * epsBH
TenseurBH * siBH = NevezTenseurBH(dim_tens);TenseurBH & sigBH = *(siBH);
sigBH = a_i.a_0 * (*IdGBH) + a_i.a_1 * epsBH
+ a_i.a_2 * (epsBH * epsBH);
sigHH = gijHH_ * sigBH; // en deux fois contravariants
//++ }
// on détruit les tenseurs intermédiaires
delete epBH; delete siBH;
LibereTenseur();
};
// calcul des contraintes et de ses variations a t+dt
//++template
//++ void HyperD::Calcul_DsigmaHH_tdt
void HyperD::Calcul_DsigmaHH_tdt
(TenseurHH & ,TenseurBB& ,DdlElement & tab_ddl
,BaseB& ,TenseurBB & gijBB_t,TenseurHH & gijHH_t,BaseB& ,Tableau & ,BaseH& ,Tableau &
,TenseurBB & epsBB_tdt,Tableau &
,TenseurBB & delta_epsBB,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt
,Tableau & d_gijBB_tdt
,Tableau & d_gijHH_tdt,double& jacobien_0,double& jacobien_tdt
,Vecteur& d_jacobien_tdt,TenseurHH& sigHH,Tableau & d_sigHH
,EnergieMeca & energ,const EnergieMeca & ,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Impli& ex)
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
cout << "\n calcul HyperD::Calcul_DsigmaHH_tdt(... "< depsBH(nbddl);
Tableau depsBH(nbddl);
for (int i=1;i<=nbddl;i++) depsBH(i) = NevezTenseurBH(dim_tens);
//++ TensBH * IdGBH; // l'identitée
TenseurBH * IdGBH; // l'identitée
// retour de IdGBH qui pointe sur le bon tenseur
InvariantVarDdl invariantVarDdl(nbddl); // def des invariants
// Calcul des invariants, de epsBH et de depsBH
// retour de IdGBH qui pointe sur le bon tenseur
IdGBH = Invariants_et_var (epsBB_tdt,gijBB_tdt,d_gijBB_tdt,
gijHH_tdt,d_gijHH_tdt,jacobien_0,jacobien_tdt,d_jacobien_tdt,
invariantVarDdl,epsBH,depsBH);
//------ pour le debug
/*cout << "\n HyperD::Calcul_DsigmaHH_tdt - \n";
cout << "\n epsBB_tdt: "; epsBB_tdt.Ecriture(cout);
cout << "\n gijBB_tdt: ";gijBB_tdt.Ecriture(cout);
cout << "\n gijHH_tdt: ";gijHH_tdt.Ecriture(cout);
for (int i=1;i<=12;i++)
{cout << "\n d_gijBB_tdt("<Ecriture(cout);};
//cout << "\n d_gijBB_tdt(2): ";d_gijBB_tdt(2)->Ecriture(cout);
//cout << "\n d_gijBB_tdt(3): ";d_gijBB_tdt(2)->Ecriture(cout);
cout << "\n d_gijHH_tdt(2): ";d_gijHH_tdt(2)->Ecriture(cout);
cout << "\n jacobien_0: " << jacobien_0 << " jacobien_tdt: " << jacobien_tdt << " d_jacobien_tdt: " << d_jacobien_tdt << endl(cout) ;
*/
// ------- fin pour le debug
energ.Inita(0.); // init des énergies
if (!avec_phase)
{ // cas sans phase
// calcul du potentiel
PotenSansPhaseAvecVar potenSansPhaseAvecVar = Potentiel_et_var(invariantVarDdl,jacobien_0);
// calcul des coefficients a_i
AAA_i_var(potenSansPhaseAvecVar,jacobien_0,invariantVarDdl,jacobien_tdt,d_jacobien_tdt,a_i);
// en fait Potentiel_et_var calcul E*racine(g), mais pour l'énergie il faut E
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenSansPhaseAvecVar.E/(invariantVarDdl.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenSansPhaseAvecVar.Ks;
}
else
{ // cas avec phase
double norm_i_epsBH = epsBH.MaxiComposante(); // norme infinie de la déformation
if ((norm_i_epsBH < ConstMath::petit)&&(!avec_regularisation))
// on est dans le cas on la déformation est très petite, on ne considère pas dans ce cas la phase
{ // calcul du potentiel
if (ParaGlob::NiveauImpression() >= 5)
cout << "\n calcul sans phase ";
PotenSansPhaseAvecVar potenSansPhaseAvecVar = Potentiel_et_var(invariantVarDdl,jacobien_0);
// calcul des coefficients a_i
AAA_i_var(potenSansPhaseAvecVar,jacobien_0,invariantVarDdl,jacobien_tdt,d_jacobien_tdt,a_i);
// en fait Potentiel_et_var calcul E*racine(g), mais pour l'énergie il faut E
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenSansPhaseAvecVar.E/(invariantVarDdl.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenSansPhaseAvecVar.Ks;
}
else // pour la suite il s'agit d'un tenseur de déformation pas "petit"
{ // on regarde s'il faut faire un scale du Qeps pour le calcul de la phase
// double scale_eps = 0.;
// if (invariantVarDdl.bIIb < ConstMath::petit)
// { // dans ce cas, on fait un scale de la norme infinie d'epsilon, qui normalement n'est pas nulle
// ce scale peut ne pas agir sur Qeps, cas d'un eps spérique, dans ce cas, le cos3phi est indéterminé
// on le prendra nulle, mais dans les autres cas normalement le scale agit: et sur Qeps et sur bIIIb
// d'où normalement un calcul correcte de l'angle
// scale_eps = 1./norm_i_epsBH;
// }
//----------------- partie exploratoire normalement n'ajoute pas de pb --------------------
// calcul des invariants pour le temps t: sera utilisé par les méthodes dérivées : ici potentiel de Laurent
// Calcul des invariants à t
// on récup la def initiale, ce n'est pas parfait car on est dans un nouveau repère, mais j'espère que c'est une bonne approximation
// TenseurBB* epsBB_t = NevezTenseurBB(epsBB_tdt); (*epsBB_t) -= delta_epsBB;
// Invariants ((*epsBB_t),gijBB_t,gijHH_t,jacobien_0,jacobien_t,invariant_t,epsBH_t);
//----------------- fin partie exploratoire ----------------------------------------------
// calcul du potentiel
PotenAvecPhaseAvecVar potenAvecPhaseAvecVar = PotentielPhase_et_var(invariantVarDdl,jacobien_0);
// calcul des coefficients a_i
AAA_iPhase_var(potenAvecPhaseAvecVar,jacobien_0,invariantVarDdl,jacobien_tdt,d_jacobien_tdt,a_i);
// en fait Potentiel_et_var calcul E*racine(g), mais pour l'énergie il faut E
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenAvecPhaseAvecVar.E/(invariantVarDdl.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenAvecPhaseAvecVar.Ks;
};
}
// calcul des contraintes finales
//* TensBH epsepsBH = epsBH && epsBH; // epsBH:epsBH
//* TensBH sigBH = alpha_0 * (*IdGBH) + alpha_1 * epsBH
TenseurBH* epsepBH = NevezTenseurBH(dim_tens);
TenseurBH& epsepsBH = *(epsepBH); epsepsBH = epsBH * epsBH; // epsBH * epsBH
TenseurBH * siBH = NevezTenseurBH(dim_tens);TenseurBH & sigBH = *(siBH);
sigBH = a_i.a_0 * (*IdGBH) + a_i.a_1 * epsBH
+ a_i.a_2 * epsepsBH;
sigHH = gijHH_tdt * sigBH; // en deux fois contravariants
// cas le la variation du tenseur des contraintes
// déclaration des variables tenseurs
TenseurBB *depBB = NevezTenseurBB(dim_tens);TenseurBB& depsBB =*(depBB);
TenseurBH *dsiBH = NevezTenseurBH(dim_tens);TenseurBH& dsigBH =*(dsiBH);
for (int i = 1; i<= nbddl; i++)
//++ { const TensHH & dgijHH = (*(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
//++ TensHH& dsigHH = *((TensHH*) (d_sigHH(i))); //
//++ TensBB depsBB = 0.5 * (*(d_gijBB_tdt(i)));
//++ TensBH depsBH = epsBB_tdt * dgijHH + depsBB * gijHH_tdt ;
//++ TensBH dsigBH = dalpha_0(i) * (*IdGBH) + dalpha_1(i) * epsBH
{ const TenseurHH & dgijHH = (*(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
TenseurHH& dsigHH = *((TenseurHH*) (d_sigHH(i))); //
depsBB = 0.5 * (*(d_gijBB_tdt(i)));
dsigBH = a_i.da_0(i) * (*IdGBH) + a_i.da_1(i) * epsBH
+ a_i.da_2(i) * epsepsBH + a_i.a_1 * (*depsBH(i))
+ a_i.a_2 * ((*depsBH(i)) *epsBH + epsBH*(*depsBH(i)));
dsigHH = dgijHH * sigBH + gijHH_tdt * dsigBH;
}
// on détruit les tenseurs intermédiaires
delete epBH; delete siBH; delete epsepBH;
delete depBB; delete dsiBH;
for (int i=1;i<=nbddl;i++)
delete depsBH(i);
LibereTenseur();
//++ }
};
// 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, sinon il s'agit de la base orthonormee fixe
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void HyperD::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & ,TenseurBB&
,TenseurBB & epsBB_tdt,TenseurBB &, double& jacobien_0,double& jacobien_tdt
,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 (Permet_affichage() > 4)
cout << "\n calcul HyperD::Calcul_dsigma_deps(... "< 1)&& (jacobien_tdt<0))
{cout << "\n *** attention le jacobien "<< jacobien_tdt <<" est negatif "
<< " cela va entrainer des erreurs dans le calcul du potentiel "
<< " \nHyperD::Calcul_dsigma_deps(... \n"
<< endl ;
};
#endif
// Calcul des invariants, et de epsBH et de depsBH
Invariants_et_varEps (epsBB_tdt,*ex.gijBB_tdt,*ex.gijHH_tdt,jacobien_0,jacobien_tdt,
invariantVarEps,*epBH);
Tenseur3BH & epsBH = *((Tenseur3BH*) epBH); // passage en dim 3
energ.Inita(0.); // init des énergies
if (!avec_phase)
{ // cas sans phase
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{cout << "\n calcul sans phase ";
if (avec_regularisation) {cout << "avec_regularisation= "<< fact_regularisation << flush;}
else {cout << "sans_regularisation= "< 7)
{cout << "\n calcul avec phase ";
if (avec_regularisation) {cout << "avec_regularisation= "<< fact_regularisation < passage en calcul sans phase ";
PotenSansPhaseAvecVar potenSansPhaseAvecVar = Potentiel_et_var(invariantVarEps,jacobien_0);
// calcul des coefficients a_i
AAA_i_varEps(potenSansPhaseAvecVar,jacobien_0,invariantVarEps,jacobien_tdt,a_i);
// en fait Potentiel_et_var calcul E*racine(g), mais pour l'énergie il faut E
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenSansPhaseAvecVar.E/(invariantVarEps.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenSansPhaseAvecVar.Ks;
}
else // pour la suite il s'agit d'un tenseur de déformation pas "petit"
{ // on regarde s'il faut faire un scale du Qeps pour le calcul de la phase
// double scale_eps = 0.;
// if (invariantVarDdl.bIIb < ConstMath::petit)
// { // dans ce cas, on fait un scale de la norme infinie d'epsilon, qui normalement n'est pas nulle
// ce scale peut ne pas agir sur Qeps, cas d'un eps spérique, dans ce cas, le cos3phi est indéterminé
// on le prendra nulle, mais dans les autres cas normalement le scale agit: et sur Qeps et sur bIIIb
// d'où normalement un calcul correcte de l'angle
// scale_eps = 1./norm_i_epsBH;
// }
// calcul du potentiel
PotenAvecPhaseAvecVar potenAvecPhaseAvecVar = PotentielPhase_et_var(invariantVarEps,jacobien_0);
// calcul des coefficients a_i
AAA_iPhase_varEps(potenAvecPhaseAvecVar,jacobien_0,invariantVarEps,jacobien_tdt,a_i);
// a priori le jacobien n'est pas nul sinon pb de distorsion d'élément ou d'intégration
energ.ChangeEnergieElastique(potenAvecPhaseAvecVar.E/(invariantVarEps.V*jacobien_0)); // traitement des énergies
// module de compressibilité
module_compressibilite = potenAvecPhaseAvecVar.Ks;
};
};
// calcul des contraintes finales
Tenseur3BH epsepsBH = epsBH * epsBH; // epsBH * epsBH
Tenseur3BH sigBH = a_i.a_0 * IdBH3 + a_i.a_1 * epsBH + a_i.a_2 * epsepsBH;
Tenseur3HH sig_localeHH = gijHH_tdt * sigBH; // en deux fois contravariants
// passage dans la base I_a si nécessaire
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{if (en_base_orthonormee)
cout << ", en base orthonormee "< Ia de l'opérateur calculer précédemment
dSigdepsHHHH.ChangeBase(d_sigma_depsFinHHHH,*(ex.giB_tdt));
}
else // cas de l'utilisation des bases curvilignes
{ I_xbarre_I_HHHH=Tenseur3HHHH::Prod_tensoriel_barre(gijHH_tdt,gijHH_tdt);
d_sigma_depsFinHHHH=
(-2.*a_i.a_0) * I_xbarre_I_HHHH // partie: d_gijHH_tdt / d_eps
+ Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,a_i.da_0_deps_HH) // partie: gijHH_tdt * d_sigHB / d_eps
+ a_i.a_1 * ( I_xbarre_I_HHHH - 2.* d1_HHHH) + Tenseur3HHHH::Prod_tensoriel(epsHH,a_i.da_1_deps_HH)
+ a_i.a_2 * (d1_HHHH -2. * d2_HHHH) + Tenseur3HHHH::Prod_tensoriel(epsepsHH,a_i.da_2_deps_HH);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n a_i.a_0= "<< a_i.a_0 << ", a_i.a_1= "<< a_i.a_1 << ", a_i.a_2= "<< a_i.a_2;
cout << "\n a_i.da_0_deps_HH= "<< a_i.da_0_deps_HH
<< ", a_i.da_1_deps_HH= " << a_i.da_1_deps_HH
<< ", a_i.da_2_deps_HH= " << a_i.da_2_deps_HH;
cout << "\n I_xbarre_I_HHHH= "<< I_xbarre_I_HHHH
<< "\n d1_HHHH= " << d1_HHHH
<< "\n d2_HHHH= " << d2_HHHH;
Tenseur3HHHH toto = Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,a_i.da_0_deps_HH);
cout << "\n Tenseur3HHHH::Prod_tensoriel(gijHH_tdt,a_i.da_0_deps_HH)= " ; toto.Affiche_bidim(cout);
toto = ( I_xbarre_I_HHHH - 2.* d1_HHHH);
cout << "\n ( I_xbarre_I_HHHH - 2.* d1_HHHH)= "; toto.Affiche_bidim(cout);
toto = Tenseur3HHHH::Prod_tensoriel(epsHH,a_i.da_1_deps_HH);
cout << "\n Tenseur3HHHH::Prod_tensoriel(epsHH,a_i.da_1_deps_HH)= "; toto.Affiche_bidim(cout);
toto = (d1_HHHH -2. * d2_HHHH);
cout << "\n (d1_HHHH -2. * d2_HHHH)= "; toto.Affiche_bidim(cout);
toto = Tenseur3HHHH::Prod_tensoriel(epsepsHH,a_i.da_2_deps_HH);
cout << "\n Tenseur3HHHH::Prod_tensoriel(epsepsHH,a_i.da_2_deps_HH)= "; toto.Affiche_bidim(cout);
};
#endif
};
// on détruit les tenseurs intermédiaires
delete epBH;
LibereTenseur();
LibereTenseurQ();
//++ }
};
//-------------------------------------------------------------------------------------------------------
// METHODES PROTEGEES :
// affichage et definition interactive des commandes particulières à chaques lois
void HyperD::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, les invariants"
<< "\n# specifiquement utilises : V, Qeps et cos3phi, 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 coefficients alpha dans le cas sans la phase c'est à dire ici sans le troisième invariant
//++template
//++ inline void HyperD::Alpha
void HyperD::AAA_i (const PotenSansPhaseSansVar& pot,const Invariant& inv,const double& jaco,
A_i& a_i)
{
// dans le cas où on ne dépend pas de la phase on n'a pas de dépendance avec Ieps
// dans les formules les terme b1,c1,d1 n'ont pas a être utilisé
double unsurjaco = 1./jaco;
// calcul des grandeurs ci di
double b2 = inv.V; double b3 = -inv.Ieps / 3.;
double c2 = 0.; double c3 = 2. * inv.Ieps / 3. + 1.;
double d2 = 0.; double d3 = -2.;
// dans le cas où on dépend de la phase calcul de a_i
a_i.a_0 = unsurjaco * (b2 * pot.EV + b3 * pot.EbIIb) ;
a_i.a_1 = unsurjaco * (c2 * pot.EV + c3 * pot.EbIIb) ;
a_i.a_2 = unsurjaco * (d2 * pot.EV + d3 * pot.EbIIb) ;
};
// calcul des coefficients alpha dans le cas avec la phase
//++template
//++ inline void AAA_iPhase(PotenAvecPhaseSansVar& potenAvecPhaseSansVar,Invariant& invariant
void HyperD::AAA_iPhase(const PotenAvecPhaseSansVar& pot,const Invariant& inv
,const double& jaco,A_i& a_i)
{ double unsurjaco = 1./jaco;
// calcul des grandeurs bi ci di
double b1 = 1.; double b2 = inv.V; double b3 = -inv.Ieps / 3.;
double c1 = -2.; double c2 = 0.; double c3 = 2. * inv.Ieps / 3. + 1.;
double d1 = 0.; double d2 = 0.; double d3 = -2.;
// dans le cas où on dépend de la phase calcul de a_i
a_i.a_0 = unsurjaco * (b1 * pot.EIeps + b2 * pot.EV + b3 * pot.EbIIb) ;
a_i.a_1 = unsurjaco * (c1 * pot.EIeps + c2 * pot.EV + c3 * pot.EbIIb) ;
a_i.a_2 = unsurjaco * (d1 * pot.EIeps + d2 * pot.EV + d3 * pot.EbIIb) ;
};
// calcul des coefficients alpha dans le cas sans la phase et de leurs variations / ddl
void HyperD::AAA_i_var(const PotenSansPhaseAvecVar& pot,const double& jaco0
,const InvariantVarDdl & inv
,const double& jaco,const Vecteur& ,A_iAvecVarDdl& a_i)
{ // dans le cas où on ne dépend pas de la phase on n'a pas de dépendance avec Ieps
// dans les formules les terme b1,c1,d1 n'ont pas a être utilisé
double unsurjaco = 1./jaco;
// calcul des grandeurs ci di
double b2 = inv.V; double b3 = -inv.Ieps / 3.;
double c2 = 0.; double c3 = 2. * inv.Ieps / 3. + 1.;
double d2 = 0.; double d3 = -2.;
// dans le cas où on dépend de la phase calcul de a_i
a_i.a_0 = unsurjaco * (b2 * pot.EV + b3 * pot.EbIIb) ;
a_i.a_1 = unsurjaco * (c2 * pot.EV + c3 * pot.EbIIb) ;
a_i.a_2 = unsurjaco * (d2 * pot.EV + d3 * pot.EbIIb) ;
// calcul des variations
int nbddl = a_i.da_0.Taille();
for (int i=1;i<= nbddl;i++)
{ // variation des termes bi, ci, di
double db2 = inv.dV(i); double db3 = -inv.dIeps(i) /3.;
double dc3 = 2. * inv.dIeps(i) / 3.;
// variation de 1/racine de g
double dunsurjaco = - inv.dV(i) / (inv.V * inv.V * jaco0) ;
// variation mixte du potentiel
double d2EdbIIbdddl = pot.EbIIb2 * inv.dbIIb(i) + pot.EbIIbV * inv.dV(i);
double d2EdVdddl = pot.EbIIbV * inv.dbIIb(i) + pot.EVV * inv.dV(i);
// variation des a_i
a_i.da_0(i) = a_i.a_0 * jaco * dunsurjaco
+ unsurjaco * (db2 * pot.EV + db3 * pot.EbIIb
+ b2 * d2EdVdddl + b3 * d2EdbIIbdddl);
a_i.da_1(i) = a_i.a_1 * jaco * dunsurjaco
+ unsurjaco * (dc3 * pot.EbIIb
+ c2 * d2EdVdddl + c3 * d2EdbIIbdddl);
a_i.da_2(i) = a_i.a_2 * jaco * dunsurjaco
+ unsurjaco * ( d2 * d2EdVdddl + d3 * d2EdbIIbdddl);
};
};
// calcul des coefficients alpha dans le cas avec la phase et de leurs variations / ddl
void HyperD::AAA_iPhase_var (const PotenAvecPhaseAvecVar& pot,const double& jaco0
,const InvariantVarDdl & inv,const double& jaco
,const Vecteur& ,A_iAvecVarDdl& a_i)
{ double unsurjaco = 1./jaco;
// calcul des grandeurs bi ci di
double b1 = 1.; double b2 = inv.V; double b3 = -inv.Ieps / 3.;
double c1 = -2.; double c2 = 0.; double c3 = 2. * inv.Ieps / 3. + 1.;
double d1 = 0.; double d2 = 0.; double d3 = -2.;
// dans le cas où on dépend de la phase calcul de a_i
a_i.a_0 = unsurjaco * (b1 * pot.EIeps + b2 * pot.EV + b3 * pot.EbIIb) ;
a_i.a_1 = unsurjaco * (c1 * pot.EIeps + c2 * pot.EV + c3 * pot.EbIIb) ;
a_i.a_2 = unsurjaco * (d1 * pot.EIeps + d2 * pot.EV + d3 * pot.EbIIb) ;
// calcul des variations
int nbddl = a_i.da_0.Taille();
for (int i=1;i<= nbddl;i++)
{ // variation des termes bi, ci, di
double db2 = inv.dV(i); double db3 = -inv.dIeps(i) /3.;
double dc3 = 2. * inv.dIeps(i) / 3.;
// variation de 1/racine de g
double dunsurjaco = - inv.dV(i) / (inv.V * inv.V * jaco0) ;
// variation mixte du potentiel
double d2EdIepsdddl = pot.EIeps2 * inv.dIeps(i) + pot.EbIIbIeps * inv.dbIIb(i)
+ pot.EIepsV * inv.dV(i);
double d2EdbIIbdddl = pot.EbIIbIeps * inv.dIeps(i) + pot.EbIIb2 * inv.dbIIb(i)
+ pot.EbIIbV * inv.dV(i);
double d2EdVdddl = pot.EIepsV * inv.dIeps(i) + pot.EbIIbV * inv.dbIIb(i)
+ pot.EVV * inv.dV(i);
// variation des a_i
a_i.da_0(i) = a_i.a_0 * jaco * dunsurjaco
+ unsurjaco * (db2 * pot.EV + db3 * pot.EbIIb +
b1 * d2EdIepsdddl + b2 * d2EdVdddl + b3 * d2EdbIIbdddl);
a_i.da_1(i) = a_i.a_1 * jaco * dunsurjaco
+ unsurjaco * (dc3 * pot.EbIIb +
c1 * d2EdIepsdddl + c2 * d2EdVdddl + c3 * d2EdbIIbdddl);
a_i.da_2(i) = a_i.a_2 * jaco * dunsurjaco
+ unsurjaco * ( d1 * d2EdIepsdddl + d2 * d2EdVdddl + d3 * d2EdbIIbdddl);
};
};
// calcul des coefficients alpha dans le cas sans la phase et de leurs variations / eps
void HyperD::AAA_i_varEps(const PotenSansPhaseAvecVar& pot,const double& jaco0
,const InvariantVarEps & inv,const double& jaco,A_iAvecVarEps& a_i)
{ // dans le cas où on ne dépend pas de la phase on n'a pas de dépendance avec Ieps
// dans les formules les terme b1,c1,d1 n'ont pas a être utilisé
double unsurjaco = 1./jaco;
// calcul des grandeurs ci di
double b2 = inv.V; double b3 = -inv.Ieps / 3.;
double c2 = 0.; double c3 = 2. * inv.Ieps / 3. + 1.;
double d2 = 0.; double d3 = -2.;
// dans le cas où on dépend de la phase calcul de a_i
a_i.a_0 = unsurjaco * (b2 * pot.EV + b3 * pot.EbIIb) ;
a_i.a_1 = unsurjaco * (c2 * pot.EV + c3 * pot.EbIIb) ;
a_i.a_2 = unsurjaco * (d2 * pot.EV + d3 * pot.EbIIb) ;
// calcul des variations
// variation des termes bi, ci, di
const Tenseur3HH& db2 = inv.dV_deps_HH;
Tenseur3HH db3 = (-1./3.)*inv.dIeps_deps_HH ;
Tenseur3HH dc3 = (2./3.) * inv.dIeps_deps_HH ;
// variation de 1/racine de g
Tenseur3HH dunsurjaco = - inv.dV_deps_HH / (inv.V * inv.V * jaco0) ;
// variation mixte du potentiel
Tenseur3HH d2EdbIIbdEps = pot.EbIIb2 * inv.dbIIb_deps_HH + pot.EbIIbV * inv.dV_deps_HH;
Tenseur3HH d2EdVdEps = pot.EbIIbV * inv.dbIIb_deps_HH + pot.EVV * inv.dV_deps_HH;
// variation des a_i
a_i.da_0_deps_HH = a_i.a_0 * jaco * dunsurjaco
+ unsurjaco * (db2 * pot.EV + db3 * pot.EbIIb
+ b2 * d2EdVdEps + b3 * d2EdbIIbdEps);
a_i.da_1_deps_HH = a_i.a_1 * jaco * dunsurjaco
+ unsurjaco * (dc3 * pot.EbIIb
+ c2 * d2EdVdEps + c3 * d2EdbIIbdEps);
a_i.da_2_deps_HH = a_i.a_2 * jaco * dunsurjaco
+ unsurjaco * ( d2 * d2EdVdEps + d3 * d2EdbIIbdEps);
};
// calcul des coefficients alpha dans le cas avec la phase et de leurs variations / eps
void HyperD::AAA_iPhase_varEps (const PotenAvecPhaseAvecVar& pot,const double& jaco0
,const InvariantVarEps & inv
,const double& jaco,A_iAvecVarEps& a_i)
{ double unsurjaco = 1./jaco;
// calcul des grandeurs bi ci di
double b1 = 1.; double b2 = inv.V; double b3 = -inv.Ieps / 3.;
double c1 = -2.; double c2 = 0.; double c3 = 2. * inv.Ieps / 3. + 1.;
double d1 = 0.; double d2 = 0.; double d3 = -2.;
// dans le cas où on dépend de la phase calcul de a_i
a_i.a_0 = unsurjaco * (b1 * pot.EIeps + b2 * pot.EV + b3 * pot.EbIIb) ;
a_i.a_1 = unsurjaco * (c1 * pot.EIeps + c2 * pot.EV + c3 * pot.EbIIb) ;
a_i.a_2 = unsurjaco * (d1 * pot.EIeps + d2 * pot.EV + d3 * pot.EbIIb) ;
// calcul des variations
// variation des termes bi, ci, di
const Tenseur3HH& db2 = inv.dV_deps_HH;
Tenseur3HH db3 = (-1./3.)*inv.dIeps_deps_HH ;
Tenseur3HH dc3 = (2./3.) * inv.dIeps_deps_HH ;
// variation de 1/racine de g
Tenseur3HH dunsurjaco = - inv.dV_deps_HH / (inv.V * inv.V * jaco0) ;
// variation mixte du potentiel
Tenseur3HH d2EdIepsdEps = pot.EIeps2 * inv.dIeps_deps_HH + pot.EbIIbIeps * inv.dbIIb_deps_HH
+ pot.EIepsV * inv.dV_deps_HH;
Tenseur3HH d2EdbIIbdEps = pot.EbIIbIeps * inv.dIeps_deps_HH + pot.EbIIb2 * inv.dbIIb_deps_HH
+ pot.EbIIbV * inv.dV_deps_HH;
Tenseur3HH d2EdVdEps = pot.EIepsV * inv.dIeps_deps_HH + pot.EbIIbV * inv.dbIIb_deps_HH
+ pot.EVV * inv.dV_deps_HH;
// variation des a_i
a_i.da_0_deps_HH = a_i.a_0 * jaco * dunsurjaco
+ unsurjaco * (db2 * pot.EV + db3 * pot.EbIIb +
b1 * d2EdIepsdEps + b2 * d2EdVdEps + b3 * d2EdbIIbdEps);
a_i.da_1_deps_HH = a_i.a_1 * jaco * dunsurjaco
+ unsurjaco * (dc3 * pot.EbIIb +
c1 * d2EdIepsdEps + c2 * d2EdVdEps + c3 * d2EdbIIbdEps);
a_i.da_2_deps_HH = a_i.a_2 * jaco * dunsurjaco
+ unsurjaco * ( d1 * d2EdIepsdEps + d2 * d2EdVdEps + d3 * d2EdbIIbdEps);
};