rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/comportement/Hyper_elastique/Hart_Smith3D.cc
Gérard Rio 0460624e94 mise à jour 7.010
2023-05-03 17:23:49 +02:00

2044 lines
95 KiB
C++
Raw Permalink Blame History

// FICHIER : Hart_Smith3D.cc
// CLASSE : Hart_Smith3D
// 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 "Debug.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "CharUtil.h"
#include "Hart_Smith3D.h"
#include "MathUtil.h"
Hart_Smith3D::Hart_Smith3D () : // Constructeur par defaut
Hyper_W_gene_3D(HART_SMITH3D,CAT_MECANIQUE,3)
,C1(-ConstMath::trespetit),C2(-ConstMath::trespetit),C3(-ConstMath::trespetit),K(-ConstMath::trespetit)
,C1_use(0.),C2_use(0.),C3_use(0.),K_use(0.) // recalculés avant chaque usage
,C1_temperature(NULL),C2_temperature(NULL),C3_temperature(NULL),K_temperature(NULL),type_pot_vol(1)
,K_nD(NULL),C1_nD(NULL),C2_nD(NULL),C3_nD(NULL)
,W_d(0.),W_v(0.)
,W_d_J1(0.),W_d_J2(0.),W_v_J3(0.),W_v_J3J3(0.)
,W_d_J1_2(0.),W_d_J1_J2(0.),W_d_J2_2(0.)
,int_J1(),W_r(),W_rs()
,avec_courbure(false),a_courbure(0.),r_courbure(0),a_temperature(NULL),r_temperature(NULL)
,W_c(0.),W_c_J1(0.),W_c_J3(0.),W_c_J1_2(0.),W_c_J3_2(0.),W_c_J1_J3(0.)
{ // définition de la courbe représentant l'évolution de l'énergie en fonction de J1
Calcul_courbe_evolW_J1();
nb_para_loi = 4; // nombre de paramètre de la loi
};
// Constructeur de copie
Hart_Smith3D::Hart_Smith3D (const Hart_Smith3D& loi) :
Hyper_W_gene_3D(loi)
,C1(loi.C1),C2(loi.C2),C3(loi.C3),K(loi.K)
,C1_use(0.),C2_use(0.),C3_use(0.),K_use(0.) // recalculés avant chaque usage
,C1_temperature(loi.C1_temperature),C2_temperature(loi.C2_temperature),C3_temperature(loi.C3_temperature)
,K_temperature(loi.K_temperature),type_pot_vol(loi.type_pot_vol)
,K_nD(loi.K_nD),C1_nD(loi.C1_nD),C2_nD(loi.C2_nD),C3_nD(loi.C3_nD)
,W_d(loi.W_d),W_v(loi.W_v)
,W_d_J1(loi.W_d_J1),W_d_J2(loi.W_d_J2),W_v_J3(loi.W_v_J3),W_v_J3J3(loi.W_v_J3J3)
,W_d_J1_2(loi.W_d_J1_2),W_d_J1_J2(loi.W_d_J1_J2),W_d_J2_2(loi.W_d_J2_2)
,int_J1(loi.int_J1),W_r(loi.W_r),W_rs(loi.W_rs)
,avec_courbure(loi.avec_courbure),a_courbure(loi.a_courbure),r_courbure(loi.r_courbure)
,a_temperature(loi.a_temperature),r_temperature(loi.r_temperature)
,W_c(loi.W_c),W_c_J1(loi.W_c_J1),W_c_J3(loi.W_c_J3),W_c_J1_2(loi.W_c_J1_2)
,W_c_J3_2(loi.W_c_J3_2),W_c_J1_J3(loi.W_c_J1_J3)
{// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale
if (C1_temperature != NULL)
if (C1_temperature->NomCourbe() == "_")
C1_temperature = Courbe1D::New_Courbe1D(*(loi.C1_temperature));
if (C2_temperature != NULL)
if (C2_temperature->NomCourbe() == "_")
C2_temperature = Courbe1D::New_Courbe1D(*(loi.C2_temperature));
if (C3_temperature != NULL)
if (C3_temperature->NomCourbe() == "_")
C3_temperature = Courbe1D::New_Courbe1D(*(loi.C3_temperature));
if (K_temperature != NULL)
if (K_temperature->NomCourbe() == "_")
K_temperature = Courbe1D::New_Courbe1D(*(loi.K_temperature));;
if(avec_courbure)
{ if (a_temperature != NULL)
if (a_temperature->NomCourbe() == "_")
a_temperature = Courbe1D::New_Courbe1D(*(loi.a_temperature));
if (r_temperature != NULL)
if (r_temperature->NomCourbe() == "_")
r_temperature = Courbe1D::New_Courbe1D(*(loi.r_temperature));
};
// idem pour les fonctions nD
if (K_nD != NULL)
if (K_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
K_nD = Fonction_nD::New_Fonction_nD(*loi.K_nD);
};
if (C1_nD != NULL)
if (C1_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
C1_nD = Fonction_nD::New_Fonction_nD(*loi.C1_nD);
};
if (C2_nD != NULL)
if (C2_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
C2_nD = Fonction_nD::New_Fonction_nD(*loi.C2_nD);
};
if (C3_nD != NULL)
if (C3_nD->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_courbe("_");
C3_nD = Fonction_nD::New_Fonction_nD(*loi.C3_nD);
};
nb_para_loi = 4; // nombre de paramètre de la loi
};
Hart_Smith3D::~Hart_Smith3D ()
// Destructeur
{ if (C1_temperature != NULL)
if (C1_temperature->NomCourbe() == "_") delete C1_temperature;
if (C2_temperature != NULL)
if (C2_temperature->NomCourbe() == "_") delete C2_temperature;
if (C3_temperature != NULL)
if (C3_temperature->NomCourbe() == "_") delete C3_temperature;
if (K_temperature != NULL)
if (K_temperature->NomCourbe() == "_") delete K_temperature;
if (a_temperature != NULL)
if (a_temperature->NomCourbe() == "_") delete a_temperature;
if (r_temperature != NULL)
if (r_temperature->NomCourbe() == "_") delete r_temperature;
// idem pour les fonctions nD
if (K_nD != NULL) if (K_nD->NomFonction() == "_") delete K_nD;
if (C1_nD != NULL) if (C1_nD->NomFonction() == "_") delete C1_nD;
if (C2_nD != NULL) if (C2_nD->NomFonction() == "_") delete C2_nD;
if (C3_nD != NULL) if (C3_nD->NomFonction() == "_") delete C3_nD;
};
// Lecture des donnees de la classe sur fichier
void Hart_Smith3D::LectureDonneesParticulieres
(UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D,LesFonctions_nD& lesFonctionsnD)
{ string toto,nom;
string nom_class_methode = "Hart_Smith3D::LectureDonneesParticulieres";
// lecture du coefficient C1
*(entreePrinc->entree) >> nom;
if(nom != "C1=")
{ cout << "\n erreur en lecture du parametre C1, on attendait la chaine C1= et on a lu: " << nom;
cout << "\n Hart_Smith3D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
entreePrinc->MessageBuffer("erreur 1 ");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde si le coefficient est thermo dépendante
if(strstr(entreePrinc->tablcar,"C1_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "C1_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de C1, on aurait du lire le mot cle C1_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
entreePrinc->MessageBuffer("**erreur2 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ C1_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
C1_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
C1_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else
{ // lecture de C1
*(entreePrinc->entree) >> C1 ;
};
// ---- lecture du coefficient C2 -----
*(entreePrinc->entree) >> nom;
if(nom != "C2=")
{ cout << "\n erreur en lecture du parametre C2, on attendait la chaine C2= et on a lu: " << nom;
cout << "\n Hart_Smith3D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
entreePrinc->MessageBuffer("erreur 2 ");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde si le coefficient est thermo dépendante
if(strstr(entreePrinc->tablcar,"C2_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "C2_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de C2, on aurait du lire le mot cle C2_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
entreePrinc->MessageBuffer("**erreur4 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ C2_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
C2_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
C2_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else
{ // lecture de C2
*(entreePrinc->entree) >> C2 ;
};
// ---- lecture du coefficient C3 -----
*(entreePrinc->entree) >> nom;
if(nom != "C3=")
{ cout << "\n erreur en lecture du parametre C3, on attendait la chaine C3= et on a lu: " << nom;
cout << "\n Hart_Smith3D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
entreePrinc->MessageBuffer("erreur 2 ");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde si le coefficient est thermo dépendante
if(strstr(entreePrinc->tablcar,"C3_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "C3_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de C3, on aurait du lire le mot cle C3_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
entreePrinc->MessageBuffer("**erreur4 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ C3_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
C3_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
C3_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else
{ // lecture de C3
*(entreePrinc->entree) >> C3 ;
};
// ---- lecture du coefficient K (dernier coefficient obligatoire) -----
*(entreePrinc->entree) >> nom;
if(nom != "K=")
{ cout << "\n erreur en lecture du parametre K, on attendait la chaine K= et on a lu: " << nom;
cout << "\n Hart_Smith3D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
entreePrinc->MessageBuffer("erreur 5 ");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// on regarde si le coefficient est thermo dépendante
if(strstr(entreePrinc->tablcar,"K_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "K_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de K, on aurait du lire le mot cle K_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
entreePrinc->MessageBuffer("**erreur6 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ K_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
K_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
K_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else
{ // lecture de K
*(entreePrinc->entree) >> K ;
// s'il n'y a plus rien n'a lire, il faut passer un enregistrement
if((strstr(entreePrinc->tablcar,"type_potvol_")==0)
&& (strstr(entreePrinc->tablcar,"avec_courbure_")==0)
&& (strstr(entreePrinc->tablcar,"avec_nD_")==0)
&& (strstr(entreePrinc->tablcar,"sortie_post_")==0)
)
entreePrinc->NouvelleDonnee();
};
// lecture éventuelle du type de potentiel
if(strstr(entreePrinc->tablcar,"type_potvol_")!=0)
{ *(entreePrinc->entree) >> nom;
if (nom != "type_potvol_")
{ cout << "\n erreur en lecture du type de variation volumique, on aurait du lire le mot cle type_potvol_"
<< " suivi d'un nombre entier et on a lue: "<< nom ;
entreePrinc->MessageBuffer("**erreur7 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture du type de variation volumique
*(entreePrinc->entree) >> type_pot_vol;
// on regarde si ce type est possible
switch (type_pot_vol)
{case 1: case 2: case 3: case 4 : break; // ok
default:
{ cout << "\n erreur en lecture du type de variation volumique: valeur lue: " << type_pot_vol
<< " , actuellement seule les types 1, 2 ,3 et 4 sont implantes ";
entreePrinc->MessageBuffer("**erreur8 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// s'il n'y a plus rien n'a lire, il faut passer un enregistrement
if(strstr(entreePrinc->tablcar,"avec_courbure_")==0) entreePrinc->NouvelleDonnee();
};
// ---- lecture éventuelle d'un terme de courbure
if(strstr(entreePrinc->tablcar,"avec_courbure_")!=0)
{ *(entreePrinc->entree) >> nom;
if (nom != "avec_courbure_")
{ cout << "\n erreur en lecture du mot cle, on aurait du lire le mot cle avec_courbure_"
<< " et on a lue: " << nom ;
entreePrinc->MessageBuffer("**erreur9 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
*(entreePrinc->entree) >> nom;
if (nom != "a_courbure=")
{ cout << "\n erreur en lecture du parametre a, on aurait du lire le mot cle a_courbure="
<< " suivi d'un nombre reel et on a lue: " << nom ;
entreePrinc->MessageBuffer("**erreur10 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
avec_courbure=true;
// on regarde si le coefficient est thermo dépendante
if(strstr(entreePrinc->tablcar,"a_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "a_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de a, on aurait du lire le mot cle a_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme, et on a lue: " << nom ;
entreePrinc->MessageBuffer("**erreur10 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ a_temperature = lesCourbes1D.Trouve(nom);}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
a_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
a_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
};
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else
{ // lecture de a
*(entreePrinc->entree) >> a_courbure ;
};
//-- on lit maintenant le paramètre r
*(entreePrinc->entree) >> nom;
if (nom != "r_courbure=")
{ cout << "\n erreur en lecture du parametre r, on aurait du lire le mot cle r_courbure="
<< " suivi d'un nombre reel et on a lue: " << nom ;
entreePrinc->MessageBuffer("**erreur11 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde si le coefficient est thermo dépendante
if(strstr(entreePrinc->tablcar,"r_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "r_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de r, on aurait du lire le mot cle r_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme, et on a lue: " << nom ;
entreePrinc->MessageBuffer("**erreur12 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ r_temperature = lesCourbes1D.Trouve(nom);}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
r_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
r_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else
{ // lecture de a
*(entreePrinc->entree) >> r_courbure ;
};
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}; //-- fin de la lecture éventuelle d'un terme de courbure
// -------- on regarde s'il y a une dépendance à une fonction nD
if (strstr(entreePrinc->tablcar,"avec_nD_")!=NULL)
{ string nom1,nom2;
// lecture des fonctions nD
entreePrinc->NouvelleDonnee(); // passage d'une ligne
// on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_avec_nD_"
// ou un nouveau mot clé global auquel cas il y a pb !!
MotCle motCle; // ref aux mots cle
while (strstr(entreePrinc->tablcar,"fin_parametres_avec_nD_")==0)
{
// si on a un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur
if ( motCle.SimotCle(entreePrinc->tablcar))
{ cout << "\n erreur de lecture des parametre de dependance a une fonction nD: on n'a pas trouve le mot cle "
<< " fin_parametres_avec_nD_ et par contre la ligne courante contient un mot cle global ";
entreePrinc->MessageBuffer
("** erreur51 des parametres de reglage de la loi de comportement Hart_Smith3D**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// lecture d'un mot clé
*(entreePrinc->entree) >> nom1;
if ((entreePrinc->entree)->rdstate() == 0)
{} // lecture normale
#ifdef ENLINUX
else if ((entreePrinc->entree)->fail())
// on a atteind la fin de la ligne et on appelle un nouvel enregistrement
{ entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement
*(entreePrinc->entree) >>nom1;
}
#else
else if ((entreePrinc->entree)->eof())
// la lecture est bonne mais on a atteind la fin de la ligne
{ if(nom1 != "fin_parametres_avec_nD_")
{entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom1;};
}
#endif
else // cas d'une erreur de lecture
{ cout << "\n erreur de lecture inconnue ";
entreePrinc->MessageBuffer
("** erreur61 des parametres de reglage de la loi de comportement Hart_Smith3D**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
if(nom1 == "K_nD=")
{string nom_fonct;
// lecture du nom de la fonction
*(entreePrinc->entree) >> nom_fonct;
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{K_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
K_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la fonction
K_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (K_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << K_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur031** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// on regarde si la fonction nD intègre la température
const Tableau <Ddl_enum_etendu>& tab_enu = K_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else if(nom1 == "C1_nD=")
{string nom_fonct;
// lecture du nom de la fonction
*(entreePrinc->entree) >> nom_fonct;
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{C1_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
C1_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la fonction
C1_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (C1_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << C1_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur031** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// on regarde si la fonction nD intègre la température
const Tableau <Ddl_enum_etendu>& tab_enu = C1_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else if(nom1 == "C2_nD=")
{string nom_fonct;
// lecture du nom de la fonction
*(entreePrinc->entree) >> nom_fonct;
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{C2_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
C2_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la fonction
C2_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (C2_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << C2_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur031** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// on regarde si la fonction nD intègre la température
const Tableau <Ddl_enum_etendu>& tab_enu = C2_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else if(nom1 == "C3_nD=")
{string nom_fonct;
// lecture du nom de la fonction
*(entreePrinc->entree) >> nom_fonct;
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{C3_nD = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
C3_nD = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la fonction
C3_nD->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (C3_nD->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << C3_nD->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur031** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// on regarde si la fonction nD intègre la température
const Tableau <Ddl_enum_etendu>& tab_enu = C3_nD->Tab_enu_etendu();
if (tab_enu.Contient(TEMP))
thermo_dependant=true;
// prepa du flot de lecture
entreePrinc->NouvelleDonnee();
}
else // cas d'une erreur
{ cout << "\n erreur en lecture de la dependance a une fonction nD des coefficients, on aurait du lire "
<< " le mot cle K_nD= ou C1_nD= ou C2_nD= ou C3_nD= et on a lu: "
<< nom1;
entreePrinc->MessageBuffer("**erreur7 Hart_Smith3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
}; //-- fin du while
entreePrinc->NouvelleDonnee(); // passage du mot clé "fin_parametres_avec_nD_"
};
// lecture de l'indication éventuelle du post traitement
string le_mot_cle = "sortie_post_";
if (entreePrinc->Lecture_un_parametre_int(0,nom_class_methode,0,1,le_mot_cle,sortie_post))
entreePrinc->NouvelleDonnee(); // passage éventuel du mot clé sortie_post_;
// --- appel au niveau de la classe mère
// ici il n'y a pas de type de déformation associé
// mais on prend la def standart d'almansi, pour les fonctions associées éventuelles
bool avec_passage_nouvelle_donnee=false;
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD,avec_passage_nouvelle_donnee);
};
// affichage de la loi
void Hart_Smith3D::Affiche() const
{ cout << " \n loi de comportement hyper elastique 3D de type Hart Smith ";
if ( C1_temperature != NULL) { cout << " C1 thermo dependant "
<< " courbe C1=f(T): " << C1_temperature->NomCourbe() <<" ";}
else { cout << " C1= " << C1 ;}
if ( C2_temperature != NULL) { cout << " C2 thermo dependant "
<< " courbe C2=f(T): " << C2_temperature->NomCourbe() <<" ";}
else { cout << " C2= " << C2 << " ";}
if ( C3_temperature != NULL) { cout << " C3 thermo dependant "
<< " courbe C3=f(T): " << C3_temperature->NomCourbe() <<" ";}
else { cout << " C3= " << C3 << " ";}
if ( K_temperature != NULL) { cout << " K thermo dependant "
<< " courbe K=f(T): " << K_temperature->NomCourbe() <<" ";}
else { cout << " K= " << K << " ";};
cout << " type de potentiel= " << type_pot_vol << " ";
if (avec_courbure)
{ if ( a_temperature != NULL) { cout << " a thermo dependant "
<< " courbe a=f(T): " << a_temperature->NomCourbe() <<" ";}
else { cout << " a= " << a_courbure ;}
if ( r_temperature != NULL) { cout << " r thermo dependant "
<< " courbe r=f(T): " << r_temperature->NomCourbe() <<" ";}
else { cout << " r= " << r_courbure ;}
};
// les dépendances à une fonction nD
if (K_nD != NULL)
{cout << " fonction nD K_nD: ";
if (K_nD->NomFonction() != "_")
cout << K_nD->NomFonction();
else
K_nD->Affiche();
};
if (C1_nD!= NULL)
{cout << " fonction nD C1_nD: ";
if (C1_nD->NomFonction() != "_")
cout << C1_nD->NomFonction();
else
C1_nD->Affiche();
};
if (C2_nD!= NULL)
{cout << " fonction nD C2_nD: ";
if (C2_nD->NomFonction() != "_")
cout << C2_nD->NomFonction();
else
C2_nD->Affiche();
};
if (C3_nD!= NULL)
{cout << " fonction nD C3_nD: ";
if (C3_nD->NomFonction() != "_")
cout << C3_nD->NomFonction();
else
C3_nD->Affiche();
};
// appel de la classe mère
Loi_comp_abstraite::Affiche_don_classe_abstraite();
};
// affichage et definition interactive des commandes particulières à chaques lois
void Hart_Smith3D::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
{ ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? ";
string rep = "_";
// procédure de lecture avec prise en charge d'un retour chariot
rep = lect_return_defaut(true,"o");
sort << "\n# ----------------------------------------------------------------------------------"
<< "\n# |...... loi de comportement hyper elastique 3D de type Hart-Smith ....... |"
<< "\n# | .. quatre coefficients C1, C2, C3 et K .. |"
<< "\n# ----------------------------------------------------------------------------------"
<< "\n\n# exemple de definition de loi"
<< "\n C1= 1. C2= 0.145 C3= 1.e-4 K= 100 "
<< "\n# .. fin de la definition de la loi Hart Smith \n" ;
if ((rep != "o") && (rep != "O" ) && (rep != "0") )
{ sort << "\n#------------------------------------------------------------------------------------"
<< "\n# il est possible de definir des parametres thermo-dependants (1 a 4 parametres)"
<< "\n# par exemple pour les quatres parametres on ecrit: "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# C1= C1_thermo_dependant_ courbe1 "
<< "\n# C2= C2_thermo_dependant_ courbe2 "
<< "\n# C3= C3_thermo_dependant_ courbe3 "
<< "\n# K= K_thermo_dependant_ courbe4 "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# noter qu'apres la definition de chaque courbe, on change de ligne, a l'inverse "
<< "\n# si la valeur du parametre est fixe, on poursuit sur la meme ligne. "
<< "\n# par exemple supposons C1 et C2 fixes et K thermo-dependant, on ecrit: "
<< "\n#------------------------------------------------------------------------------------"
<< "\n# C1= 1. C2= 0.145 C3= C3_thermo_dependant_ courbe4 "
<< "\n# K= K_thermo_dependant_ courbe4 "
<< "\n#---------------------------------- type_potvol_ ----------------------------------------"
<< "\n# un parametre facultatif permet d'indiquer le type de variation volumique "
<< "\n# que l'on desire: par defaut il s'agit de : K(1-(1+log(V))/V) qui correspond au type 1"
<< "\n# mais on peut choisir: K/2(V-1) qui correspond au type 2 "
<< "\n# ou: K/2(log(V))^2 qui correspond au type 3 "
<< "\n# ou: K/2(V-1)^2 qui correspond au type 4 "
<< "\n# en indiquant (en fin de ligne) le mot cle: type_potvol_ suivi du type"
<< "\n# par defaut type_potvol_ a la valeur 1 "
<< "\n# "
<< "\n#---------------------------------- avec_courbure_ -------------------------------------"
<< "\n# apres le type de variation volumique on peut indiquer facultativement l'ajout au potentiel "
<< "\n# d'un terme permettant de raidir le comportement a partir d'un certain niveau de chargement "
<< "\n# pour cela on indique le mot cle: avec_courbure_ puis on change de ligne "
<< "\n# le terme additionnelle depend de deux parametres: a et r, a positionne la valeur de J1 "
<< "\n# a partir de laquelle il y a durcissement, r controle la courbure du changement de regime "
<< "\n# exemple de declaration : "
<< "\n# C1= 1. C2= 0.145 C3= 1.e-4 K= 100 type_potvol_ 2 avec_courbure_ "
<< "\n# a_courbure= 94 r_courbure= 1.24 "
<< "\n# ces deux parametres peuvent etre l'un et/ou l'autre dependant de la temperature "
<< "\n# dans ce cas la declaration de dependance suit les regles habituelles "
<< "\n# exemple de declaration dans le cas d'une dependance a la temperature : "
<< "\n# C1= 1. C2= 0.145 C3= 1.e-4 K= 100 type_potvol_ 2 avec_courbure_ "
<< "\n# a_courbure= a_thermo_dependant_ "
<< "\n# r_courbure= r_thermo_dependant_ "
<< "\n# "
<< "\n#--------------------- dependance explicite a une fonction nD -----------------------------------"
<< "\n# il est egalement possible de definir une dependance des parametres a une fonction nD "
<< "\n# dans ce cas, a la suite de la dependance explicite a la temperature via une courbr 1D "
<< "\n# ou a la suite du dernier parametre fixe, on met le mot cle: avec_nD_ "
<< "\n# puis sur les lignes qui suivent on met des fonctions nD multiplicatives des parametres "
<< "\n# initiaux, puis un mot cle signalant la fin des fonctions: fin_parametres_avec_nD_ "
<< "\n# suit un exemple de declaration (il est possible d'ommettre certaine fonction ! "
<< "\n# la presence de fct nD est facultative), de meme il est possible soit d'utiliser une fonction "
<< "\n# deja defini soit de donner directement la fonction nD (comme pour les autres lois) "
<< "\n# Enfin, apres chaque definition de fonction il faut passer a une nouvelle ligne "
<< "\n# "
<< "\n# exemple de definition de loi"
<< "\n# C1= 1. C2= 0.145 C3= 1.e-4 K= 100 avec_nD_ "
<< "\n# K_nD= fct_nD_1 "
<< "\n# C1_nD= fct_nD_2 "
<< "\n# C3_nD= fct_nD_3 "
<< "\n# fin_parametres_avec_nD_ "
<< "\n# "
<< "\n# Remarques: "
<< "\n# - les fonction nD sont des fonctions multiplicatives "
<< "\n# des parametres initiaux (ou de ceux obtenus apres dependance explicite a la temperature "
<< "\n# ( bien noter que la loi obtenue peut-etre quelconque, en particulier plus du tout "
<< "\n# hyperelastique, tout depend des choix des fcts nD !) "
<< "\n# "
<< "\n#----------------------------------- sortie_post_ -------------------------------------"
<< "\n# il est possible de recuperer differentes grandeurs de travail par exemple "
<< "\n# l'intensite du potentiel, comme ces grandeurs sont calculees au moment de la resolution "
<< "\n# et ne sont pas stockees. On peut egalement recuperer la valeur des coefficients de la loi "
<< "\n# (utile s'ils varient) sous forme d'un tableau de quadruplets : K C1 C2 C3 "
<< "\n# Il faut indiquer le mot cle sortie_post_ suivi de 1 (par defaut = 0) "
<< "\n# ensuite au moment de la constitution du .CVisu on aura acces aux grandeurs de travail "
<< "\n# ex: "
<< "\n# sortie_post_ 1 "
<< "\n# ce mot cle est le dernier des parametres specifiques de la loi il doit se situe "
<< "\n# a la fin de la derniere ligne de donnees "
<< "\n#"
<< "\n#------------------------------------------------------------------------------------"
<< "\n#------------------------------------------------------------------------------------"
<< endl;
};
sort << endl;
// appel de la classe Hyper_W_gene_3D
Hyper_W_gene_3D::Info_commande_LoisDeComp_hyper3D(entreePrinc);
// appel de la classe mère
Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);
};
// test si la loi est complete
int Hart_Smith3D::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
if ((C1_temperature == NULL) && (C1 == (-ConstMath::trespetit)))
{ cout << " \n le coefficient C1 de la loi de Hart Smith n'est pas defini pour la loi "
<< Nom_comp(id_comp)
<< '\n' << endl;
ret = 0;
};
if ((C2_temperature == NULL) && (C2 == (-ConstMath::trespetit)))
{ cout << " \n le coefficient C2 de la loi de Hart Smith n'est pas defini pour la loi "
<< Nom_comp(id_comp)
<< '\n' << endl;
ret = 0;
};
if ((C3_temperature == NULL) && (C3 == (-ConstMath::trespetit)))
{ cout << " \n le coefficient C3 de la loi de Hart Smith n'est pas defini pour la loi "
<< Nom_comp(id_comp)
<< '\n' << endl;
ret = 0;
};
if ((K_temperature == NULL) && (K == (-ConstMath::trespetit)))
{ cout << " \n le coefficient K de la loi de Hart Smith n'est pas defini pour la loi "
<< Nom_comp(id_comp)
<< '\n' << endl;
ret = 0;
};
if (avec_courbure)
{ if ((a_temperature == NULL) && (a_courbure == (-ConstMath::trespetit)))
{ cout << " \n le coefficient de courbure a (loi de Hart Smith + raidissement) n'est pas defini pour la loi "
<< Nom_comp(id_comp)
<< '\n' << endl;
ret = 0;
};
if ((r_temperature == NULL) && (r_courbure == (-ConstMath::trespetit)))
{ cout << " \n le coefficient de courbure r (loi de Hart Smith + raidissement) n'est pas defini pour la loi "
<< Nom_comp(id_comp)
<< '\n' << endl;
ret = 0;
};
};
//
if (ret == 0)
{this-> Affiche();
ret = 0;
};
return ret;
};
//----- lecture écriture de restart -----
// 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 Hart_Smith3D::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string nom; bool test;
if (cas == 1)
{ // les coefficients
// C1
ent >> nom >> test;
// vérification
#ifdef MISE_AU_POINT
if (nom != "C1=")
{ cout << "\n erreur en lecture du parametres C1 de la loi de Hart Smith 3D"
<< " on devait lire C1= avant le premier parametre "
<< " et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
}
#endif
if (!test)
{ ent >> C1;
if (C1_temperature != NULL) {if (C1_temperature->NomCourbe() == "_")
delete C1_temperature; C1_temperature = NULL;};
}
else
{ ent >> nom; C1_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,C1_temperature); };
// on regarde si il y a une fonction nD multiplicative
ent >> nom;
if (nom == "C1_nD:")
C1_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,C1_nD);
// C2
ent >> nom >> test;
// vérification
#ifdef MISE_AU_POINT
if (nom != "C2=")
{ cout << "\n erreur en lecture du parametres C2 de la loi de Hart Smith 3D"
<< " on devait lire C2= avant le second parametre "
<< " et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
};
#endif
if (!test)
{ ent >> C2;
if (C2_temperature != NULL) {if (C2_temperature->NomCourbe() == "_")
delete C2_temperature; C2_temperature = NULL;};
}
else
{ ent >> nom; C2_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,C2_temperature); };
// on regarde si il y a une fonction nD multiplicative
ent >> nom;
if (nom == "C2_nD:")
C2_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,C1_nD);
// C3
ent >> nom >> test;
// vérification
#ifdef MISE_AU_POINT
if (nom != "C3=")
{ cout << "\n erreur en lecture du parametres C3 de la loi de Hart Smith 3D"
<< " on devait lire C3= avant le second parametre "
<< " et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
}
#endif
if (!test)
{ ent >> C3;
if (C3_temperature != NULL) {if (C3_temperature->NomCourbe() == "_")
delete C3_temperature; C3_temperature = NULL;};
}
else
{ ent >> nom; C3_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,C3_temperature); };
// on regarde si il y a une fonction nD multiplicative
ent >> nom;
if (nom == "C3_nD:")
C3_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,C3_nD);
// K
ent >> nom >> test;
// vérification
#ifdef MISE_AU_POINT
if (nom != "K")
{ cout << "\n erreur en lecture du parametres K de la loi de Hart Smith 3D"
<< " on devait lire K= avant le second parametre "
<< " et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
};
#endif
if (!test)
{ ent >> K;
if (K_temperature != NULL) {if (K_temperature->NomCourbe() == "_")
delete K_temperature; K_temperature = NULL;};
}
else
{ ent >> nom; K_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,K_temperature); };
// on regarde si il y a une fonction nD multiplicative
ent >> nom;
if (nom == "K_nD:")
K_nD = lesFonctionsnD.Lecture_pour_base_info(ent,cas,K_nD);
// le type de potentiel
ent >> nom >> type_pot_vol ;
// --- potentiel de gestion de courbure, éventuellement
ent >> nom ;
#ifdef MISE_AU_POINT
if (nom != "avec_courbure=")
{ cout << "\n erreur en lecture de l'indicateur de courbure de la loi de Hart Smith 3D"
<< " on devait lire avec_courbure= et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
}
#endif
ent >> avec_courbure;
if (avec_courbure)
{// a_courbure
ent >> nom >> test;
// vérification
#ifdef MISE_AU_POINT
if (nom != "a=")
{ cout << "\n erreur en lecture du parametres a de la gestion eventuelle de courbure de la loi de Hart Smith 3D"
<< " on devait lire a= avant le premier parametre "
<< " et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
};
#endif
if (!test)
{ ent >> a_courbure;
if (a_temperature != NULL) {if (a_temperature->NomCourbe() == "_")
delete a_temperature; a_temperature = NULL;};
}
else
{ ent >> nom; a_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,a_temperature); };
// r_courbure
ent >> nom >> test;
// vérification
#ifdef MISE_AU_POINT
if (nom != "r=")
{ cout << "\n erreur en lecture du parametres r de la gestion eventuelle de courbure de la loi de Hart Smith 3D"
<< " on devait lire r= avant le premier parametre "
<< " et on a lue: " << nom << " "
<< "\n Hart_Smith3D::Lecture_base_info_loi(..."
<< endl;
Sortie(1);
};
#endif
if (!test)
{ ent >> r_courbure;
if (r_temperature != NULL) {if (r_temperature->NomCourbe() == "_")
delete r_temperature; r_temperature = NULL;};
}
else
{ ent >> nom; r_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,r_temperature); };
};
}
// appel de la classe mère
Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Hart_Smith3D::Ecriture_base_info_loi(ofstream& sort,const int cas)
{ if (cas == 1)
{ // les coefficients
// C1
sort << "\n C1= ";
if (C1_temperature == NULL)
{ sort << false << " " << C1 << " ";}
else
{ sort << true << " fonction_C1_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,C1_temperature);
};
if (C1_nD != NULL)
{sort << " 2 C1_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,C1_nD);
}
else {sort << " non_C1_nD "; };
// C2
sort << " C2= ";
if (C2_temperature == NULL)
{ sort << false << " " << C2 << " ";}
else
{ sort << true << " fonction_C2_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,C2_temperature);
};
if (C2_nD != NULL)
{sort << " 2 C2_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,C2_nD);
}
else {sort << " non_C2_nD "; };
// C3
sort << " C3= ";
if (C3_temperature == NULL)
{ sort << false << " " << C3 << " ";}
else
{ sort << true << " fonction_C3_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,C3_temperature);
};
if (C3_nD != NULL)
{sort << " 2 C3_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,C3_nD);
}
else {sort << " non_C3_nD "; };
// K
sort << " K= ";
if (K_temperature == NULL)
{ sort << false << " " << K << " ";}
else
{ sort << true << " fonction_K_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,K_temperature);
};
if (K_nD != NULL)
{sort << " K_nD: " << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,K_nD);
}
else {sort << " non_K_nD "; };
sort << " typ_pot= " << type_pot_vol << " "; // le type de potentiel
// --- potentiel de gestion de courbure, éventuellement
sort << " avec_courbure= " << avec_courbure;
if (avec_courbure)
{ sort << " a= ";
if (a_temperature == NULL)
{ sort << false << " " << a_courbure << " ";}
else
{ sort << true << " fonction_a_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,a_temperature);
};
sort << "\n r= ";
if (r_temperature == NULL)
{ sort << false << " " << r_courbure << " ";}
else
{ sort << true << " fonction_r_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,r_temperature);
};
};
}
// appel de la classe mère
Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);
};
// ========== codage des METHODES VIRTUELLES protegees:================
// calcul des contraintes a t+dt
void Hart_Smith3D::Calcul_SigmaHH (TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl,
TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB & epsBB_,
TenseurBB & ,
TenseurBB & gijBB_,TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,
double& jacobien_0,double& jacobien,TenseurHH & sigHH_
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Expli_t_tdt& ex)
{
#ifdef MISE_AU_POINT
if (epsBB_.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hart_Smith3D::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 << " Hart_Smith3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage explicite en tenseur dim 3
// init coefficient de la loi
C1_use = C1; C2_use = C2; C3_use = C3; K_use = K;
// cas de la thermo dépendance, on calcul les grandeurs
if (C1_temperature != NULL) C1_use = C1_temperature->Valeur(*temperature);
if (C2_temperature != NULL) C2_use = C2_temperature->Valeur(*temperature);
if (C3_temperature != NULL) C3_use = C3_temperature->Valeur(*temperature);
if (K_temperature != NULL) K_use = K_temperature->Valeur(*temperature);
if (avec_courbure)
{ if (a_temperature != NULL) a_courbure = a_temperature->Valeur(*temperature);
if (r_temperature != NULL) r_courbure = r_temperature->Valeur(*temperature);
};
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
// dans le cas de l'utilisation de fct nD
if ( (K_nD != NULL) || (C1_nD != NULL)
|| (C2_nD != NULL) || (C3_nD != NULL)
)
{ // opération de transmission de la métrique
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = &ex;
const Met_abstraite::Umat_cont* ex_expli = NULL;
if (K_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
K_use *= tab_val(1);
};
if (C1_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C1_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C1_use *= tab_val(1);
};
if (C2_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C2_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C2_use *= tab_val(1);
};
if (C3_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C3_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C3_use *= tab_val(1);
};
};
if (sortie_post)
{// sauvegarde des paramètres matériau
(*save_resulHyper_W.para_loi)(1) = K_use;
(*save_resulHyper_W.para_loi)(2) = C1_use;
(*save_resulHyper_W.para_loi)(3) = C2_use;
(*save_resulHyper_W.para_loi)(4) = C3_use;
};
// calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D)
Invariants_et_var1(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_,gijHH_,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
Potentiel_et_var(module_compressibilite);
// stockage éventuel pour du post-traitement
if (sortie_post)
{ // récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
save_resulHyper_W.invP->potentiel= W_d + W_v + W_c;
};
// calcul du tenseur des contraintes
sigHH = ((W_d_J1+W_c_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3+W_c_J3)/V) * d_J_r_epsBB_HH(3);
// // calcul du module de compressibilité
// module_compressibilite = 2. * V * (W_v_J3+W_c_J3) / (MaX(ConstMath::petit, log(V) ));
// pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
module_cisaillement = 0.;
// traitement des énergies
energ.Inita(0.);
energ.ChangeEnergieElastique((W_d+W_v+W_c)/V);
LibereTenseur();
};
// calcul des contraintes a t+dt et de ses variations
void Hart_Smith3D::Calcul_DsigmaHH_tdt
( TenseurHH& ,TenseurBB& ,DdlElement & tab_ddl
,BaseB& ,TenseurBB & ,TenseurHH & ,
BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> & ,
TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,TenseurBB & ,
TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt_,
Tableau <TenseurBB *>& d_gijBB_tdt,
Tableau <TenseurHH *>& ,double& jacobien_0,double& jacobien,
Vecteur& ,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Impli& ex
)
{
#ifdef MISE_AU_POINT
if (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Hart_Smith3D::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 << " Hart_Smith3D::Calcul_SDsigmaHH_tdt\n";
Sortie(1);
};
#endif
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // passage en dim 3 explicite
Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) &gijHH_tdt_); // passage en dim 3 explicite
// init coefficient de la loi
C1_use = C1; C2_use = C2; C3_use = C3; K_use = K;
// cas de la thermo dépendance, on calcul les grandeurs
if (C1_temperature != NULL) C1_use = C1_temperature->Valeur(*temperature);
if (C2_temperature != NULL) C2_use = C2_temperature->Valeur(*temperature);
if (C3_temperature != NULL) C3_use = C3_temperature->Valeur(*temperature);
if (K_temperature != NULL) K_use = K_temperature->Valeur(*temperature);
if (avec_courbure)
{ if (a_temperature != NULL) a_courbure = a_temperature->Valeur(*temperature);
if (r_temperature != NULL) r_courbure = r_temperature->Valeur(*temperature);
};
//--- pour le debug
//cout << "\n C1_use=" << C1_use << " C2_use=" << C2_use << " C3_use=" << C3_use << " K_use=" << K_use << " T=" << *temperature;
//--- fin debug
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
// dans le cas de l'utilisation de fct nD
if ( (K_nD != NULL) || (C1_nD != NULL)
|| (C2_nD != NULL) || (C3_nD != NULL)
)
{ // opération de transmission de la métrique
const Met_abstraite::Impli* ex_impli = &ex;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_expli = NULL;
if (K_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
K_use *= tab_val(1);
};
if (C1_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C1_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C1_use *= tab_val(1);
};
if (C2_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C2_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C2_use *= tab_val(1);
};
if (C3_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C3_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C3_use *= tab_val(1);
};
};
if (sortie_post)
{// sauvegarde des paramètres matériau
(*save_resulHyper_W.para_loi)(1) = K_use;
(*save_resulHyper_W.para_loi)(2) = C1_use;
(*save_resulHyper_W.para_loi)(3) = C2_use;
(*save_resulHyper_W.para_loi)(4) = C3_use;
};
// calcul des invariants et de leurs variations premières et secondes
Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
Potentiel_et_var2(module_compressibilite);
// stockage éventuel pour du post-traitement
if (sortie_post)
{ // récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
save_resulHyper_W.invP->potentiel= W_d + W_v + W_c;
};
// calcul du tenseur des contraintes
double unSurV=1./V;
sigHH = ((W_d_J1+W_c_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3+W_c_J3)/V) * d_J_r_epsBB_HH(3);
// cout << "\n sigHH " << sigHH
// << "\n d_J_r_epsBB_HH1 " << d_J_r_epsBB_HH(1)
// << "\n d_J_r_epsBB_HH2 " << d_J_r_epsBB_HH(2)
// << "\n d_J_r_epsBB_HH3 " << d_J_r_epsBB_HH(3)
// << "\n (W_d_J1*unSurV) " << (W_d_J1*unSurV)
// << " (W_d_J2*unSurV) " << (W_d_J2*unSurV)
// << " (W_d_J3*unSurV) " << (W_v_J3*unSurV);
// calcul de la variation seconde du potentiel par rapport à epsij epskl
Tenseur3HHHH d2W_d2epsHHHH
= // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
(W_d_J1_2 + W_c_J1_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
+ W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
+ W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
// puis les dérivées premières du potentiel déviatoire + courbure éventuellement
+ (W_d_J1 + W_c_J1) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
// enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
+ (W_v_J3 + W_c_J3) * d_J_3_eps2BB_HHHH
+ (W_v_J3J3 + W_c_J3_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
if (avec_courbure)
d2W_d2epsHHHH += W_c_J1_J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(3));
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV);
// Tenseur3HHHH dSigdepsHHHH(1,interHH,d_J_r_epsBB_HH(3));
Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt);
// dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH;
Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH); // cas des tenseurs généraux
dSigdepsHHHH += interHHHH; // cas des tenseurs généraux
//---------------------------------------------------------------------------------
// vérif numérique de l'opérateur tangent
// Cal_dsigma_deps_num (*(ex.gijBB_0),*(ex.gijHH_0),gijBB_tdt,gijHH_tdt,jacobien_0,jacobien,dSigdepsHHHH);
//---------------------------------------------------------------------------------
// calcul des variations / aux ddl
int nbddl = d_gijBB_tdt.Taille();
for (int i = 1; i<= nbddl; i++)
{ // on fait uniquement une égalité d'adresse et de ne pas utiliser
// le constructeur d'ou la profusion d'* et de ()
Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i))); // passage en dim 3
const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); // "
dsigHH = dSigdepsHHHH && depsBB;
};
// // calcul du module de compressibilité
// module_compressibilite = 2. * V * (W_v_J3+W_c_J3) / (MaX(ConstMath::petit, log(V) ));
// pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
module_cisaillement = 0.;
// traitement des énergies
energ.Inita(0.);
energ.ChangeEnergieElastique((W_d+W_v+W_c)/V);
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 orthonormes
// 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 Hart_Smith3D::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & ,TenseurBB&
,TenseurBB & epsBB_tdt,TenseurBB &, double& jacobien_0,double& jacobien
,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 << " Hart_Smith3D::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // // passage en dim 3 explicite
Tenseur3HH & gijHH_tdt = *((Tenseur3HH*) ex.gijHH_tdt); // passage en dim 3 explicite
// init coefficient de la loi
C1_use = C1; C2_use = C2; C3_use = C3; K_use = K;
// cas de la thermo dépendance, on calcul les grandeurs
if (C1_temperature != NULL) C1_use = C1_temperature->Valeur(*temperature);
if (C2_temperature != NULL) C2_use = C2_temperature->Valeur(*temperature);
if (C3_temperature != NULL) C3_use = C3_temperature->Valeur(*temperature);
if (K_temperature != NULL) K_use = K_temperature->Valeur(*temperature);
if (avec_courbure)
{ if (a_temperature != NULL) a_courbure = a_temperature->Valeur(*temperature);
if (r_temperature != NULL) r_courbure = r_temperature->Valeur(*temperature);
};
// récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
// dans le cas de l'utilisation de fct nD
if ( (K_nD != NULL) || (C1_nD != NULL)
|| (C2_nD != NULL) || (C3_nD != NULL)
)
{ // opération de transmission de la métrique
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_expli = &ex;
if (K_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(K_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
K_use *= tab_val(1);
};
if (C1_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C1_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C1_use *= tab_val(1);
};
if (C2_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C2_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C2_use *= tab_val(1);
};
if (C3_nD != NULL)
// là il faut calculer la fonction nD
{ // on utilise la méthode générique de loi abstraite
list <SaveResul*> list_save; // inter pour l'appel de la fonction
list_save.push_back(saveResul);
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(C3_nD,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,&list_save
);
// on récupère le premier élément du tableau uniquement
C3_use *= tab_val(1);
};
};
if (sortie_post)
{// sauvegarde des paramètres matériau
(*save_resulHyper_W.para_loi)(1) = K_use;
(*save_resulHyper_W.para_loi)(2) = C1_use;
(*save_resulHyper_W.para_loi)(3) = C2_use;
(*save_resulHyper_W.para_loi)(4) = C3_use;
};
// calcul des invariants et de leurs variations premières et secondes
Invariants_et_var2(*(ex.gijBB_0),*(ex.gijHH_0),*(ex.gijBB_tdt),gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
Potentiel_et_var2(module_compressibilite);
// stockage éventuel pour du post-traitement
if (sortie_post)
{ // récup du conteneur spécifique du point, pour sauvegarde éventuelle
SaveResulHyper_W_gene_3D & save_resulHyper_W = *((SaveResulHyper_W_gene_3D*) saveResul);
save_resulHyper_W.invP->potentiel= W_d + W_v + W_c;
};
// calcul du tenseur des contraintes, on travaille ici dans le repère matériel finale correspondant
// aux coordonnées initiales X_(0)^a, on obtient donc un tenseur dans la base naturelle finale
double unSurV=1./V;
Tenseur3HH sig_localeHH = ((W_d_J1+W_c_J1)*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2)
+ ((W_v_J3+W_c_J3)*unSurV) * d_J_r_epsBB_HH(3);
// passage éventuelle dans la base I_a
if (en_base_orthonormee)
{sig_localeHH.BaseAbsolue(sigHH,*(ex.giB_tdt));}
else {sigHH = sig_localeHH;}; // sinon la base locale est la bonne
// calcul de la variation seconde du potentiel par rapport à epsij epskl
// calcul de la variation seconde du potentiel par rapport à epsij epskl
Tenseur3HHHH d2W_d2epsHHHH
= // tout d'abord les dérivées secondes du potentiel déviatoire + courbure éventuellement
(W_d_J1_2 + W_c_J1_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(1))
+ W_d_J1_J2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(1),d_J_r_epsBB_HH(2))
+ W_d_J2_2 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(2),d_J_r_epsBB_HH(2))
// puis les dérivées premières du potentiel déviatoire + courbure éventuellement
+ (W_d_J1 + W_c_J1) * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
// enfin les dérivées seconde et première du potentiel sphérique + courbure éventuellement
+ (W_v_J3 + W_c_J3) * d_J_3_eps2BB_HHHH
+ (W_v_J3J3 + W_c_J3_2) * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HH interHH = -sig_localeHH ; //* (-0.5*unSurV*unSurV);
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HHHH dSigdepsHHHH(1,interHH,gijHH_tdt);
// dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH;
Tenseur3HHHH interHHHH((unSurV) * d2W_d2epsHHHH); // cas des tenseurs généraux
dSigdepsHHHH += interHHHH; // cas des tenseurs généraux
// transfert des informations: on pas d'un tenseur de 81 composantes à 36
// avec des symétries par rapport aux deux premiers indices et par rapport aux deux derniers
/// Tenseur3HHHH d_sigma_depsHHHH; d_sigma_depsHHHH.TransfertDunTenseurGeneral(dSigdepsHHHH.Symetrise1et2_3et4());
// calcul de la première partie de l'opérateur tangent (correspond au changement de repère
// gi_tdt -> Ia de l'opérateur calculer précédemment
Tenseur3HHHH & d_sigma_depsFinHHHH = *((Tenseur3HHHH*) &d_sigma_deps_); // pour accés directe
// passage éventuelle dans la base I_a
if (en_base_orthonormee)
{dSigdepsHHHH.ChangeBase(d_sigma_depsFinHHHH,*(ex.giB_tdt));}
else
{d_sigma_depsFinHHHH = dSigdepsHHHH;};
/* Tenseur3HHHH dSigdepsHHHH(1,interHH,d_J_r_epsBB_HH(3));
dSigdepsHHHH += (unSurV) * d2W_d2epsHHHH;
// transfert des informations: on pas d'un tenseur de 81 composantes à 36
// avec des symétries par rapport aux deux premiers indices et par rapport aux deux derniers
Tenseur3HHHH & d_sigma_depsHHHH = *((Tenseur3HHHH*) &d_sigma_deps_); // pour accés directe
d_sigma_depsHHHH.TransfertDunTenseurGeneral(dSigdepsHHHH.Symetrise1et2_3et4());
// d_sigma_depsHHHH.TransfertDunTenseurGeneral(dSigdepsHHHH);
// for (int i=1;i<=3;i++) for (int j=1;j<=3;j++) for (int k=1;k<=3;k++) for (int l=1;l<=3;l++)
// d_sigma_depsHHHH.Change(i,j,k,l,0.25*(dSigdepsHHHH(i,j,k,l)+dSigdepsHHHH(j,i,k,l)+dSigdepsHHHH(i,j,l,k)+dSigdepsHHHH(j,i,l,k)));
*/
// // calcul du module de compressibilité
// module_compressibilite = 2. * V * (W_v_J3+W_c_J3) / (MaX(ConstMath::petit, log(V) ));
// pour le module de cisaillement, pour l'instant je ne fais rien !! à voir ***
module_cisaillement = 0.;
energ.Inita(0.);
energ.ChangeEnergieElastique((W_d+W_v+W_c)/V);
LibereTenseur();
LibereTenseurQ();
};
//---------------------- méthodes privées -------------------------------
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
void Hart_Smith3D::Potentiel_et_var(double & module_compressibilite)
{ // calcul de grandeurs intermédiaires
double unSurV=1./V; // le cas V=0 a normalement été traité dans Hyper_W_gene_3D
double unSurV2=unSurV*unSurV;
// calcul du potentiel
double logV = log(V);
//---- ici la valeur du potentiel est à calculer par intégration
W_d = C1_use*int_J1.Valeur(J_r(1)) + C2_use * log(J_r(2)/3.); // partie déviatorique
W_d_J1 = C1_use*exp(C3_use*Sqr(J_r(1)-3.)); // variation / J1
W_d_J2 = C2_use/J_r(2); // variation / J2
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{ cout << "\n Hart_Smith3D::Potentiel_et_var(..."
<< "\n W_d= " << W_d << ", W_d_J1= "<< W_d_J1 <<", W_d_J2= "<< W_d_J2 <<" ";
};
#endif
// puis partie volumique
double VmoinsUn = V-1.;
switch (type_pot_vol)
{case 1: {W_v = K_use * (1- (1+logV)*unSurV); // potentiel 1
W_v_J3 = 0.5 * K_use * logV * unSurV * unSurV2; // variation / J3
module_compressibilite = K_use * unSurV2;
break;
}
case 2: {W_v = 0.5 * K_use * VmoinsUn; // potentiel 2
W_v_J3 = 0.25 * K_use * unSurV; // variation / J3
if (Dabs(VmoinsUn) > ConstMath::petit )
{module_compressibilite = 0.5 * K_use / logV;}
else // on fait un développement limité du log pour V proche de 1.
{module_compressibilite = 0.5 * K_use / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
break;
}
case 3: {W_v = K_use * 0.5 * logV * logV; // potentiel 3
W_v_J3 = K_use * 0.5 * unSurV2 * logV; // variation / J3
module_compressibilite = K_use / V ;
break;
}
case 4: {W_v = 0.5 * K_use * VmoinsUn * VmoinsUn; // potentiel
W_v_J3 = 0.5 * K_use * unSurV * VmoinsUn; // variation / J3
if (Dabs(VmoinsUn) > ConstMath::petit )
{module_compressibilite = VmoinsUn * K_use /logV ;}
else // on fait un développement limité du log pour V proche de 1.
{module_compressibilite = VmoinsUn * K_use / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
break;
}
};
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{switch (type_pot_vol)
{case 1: case 2: case 3: case 4:
cout << "\n W_v= " << W_v << ", W_v_J3= "<< W_v_J3 << ", K_s= " << module_compressibilite << " ";
break;
};
};
#endif
//--- partie courbure éventuelle
if (avec_courbure)
{double unSurV3=unSurV2*unSurV;
double aP2r = pow(a_courbure,2.*r_courbure);
double J_rMoins3P2r = pow((J_r(1)-3.),(2.*r_courbure));
W_c = unSurV * ( J_rMoins3P2r * (J_r(1)-3.) / aP2r / (2.*r_courbure+1) );
W_c_J1 = unSurV * ( J_rMoins3P2r/aP2r );
W_c_J3 = -0.5*unSurV3 * ( J_rMoins3P2r * (J_r(1)-3.) /aP2r/(2.*r_courbure+1) );
module_compressibilite += 2. * V * W_c_J3 / (MaX(ConstMath::petit, log(V) ));
}
else
{W_c = W_c_J1 = W_c_J3 = 0.;};
};
// calcul du potentiel et de ses dérivées premières et secondes / aux invariants J_r
void Hart_Smith3D::Potentiel_et_var2(double & module_compressibilite)
{ // calcul de grandeurs intermédiaires
double unSurV=1./V; // le cas V=0 a normalement été traité dans Hyper_W_gene_3D
// calcul du potentiel
double logV = log(V);
//---- ici la valeur du potentiel est à calculer par intégration
W_d = C1_use*int_J1.Valeur(J_r(1)) + C2_use * log(J_r(2)/3.); // partie déviatorique
// calcul des variations premières non nulles du potentiel
// les variations secondes sont nulles
W_d_J1 = C1_use*exp(C3_use*Sqr(J_r(1)-3.)); // variation / J1
W_d_J2 = C2_use/J_r(2); // variation / J2
// dérivées secondes
W_d_J1_2 = 2.*C1_use*C3_use*(J_r(1)-3.)*exp(C3_use*Sqr(J_r(1)-3.));
W_d_J1_J2 = 0.;
W_d_J2_2 = -C2_use/(Sqr(J_r(2)));
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{ cout << "\n Hart_Smith3D::Potentiel_et_var2(..."
<< "\n W_d= " << W_d << ", W_d_J1= "<< W_d_J1 <<", W_d_J2= "<< W_d_J2
<< ", W_d_J1_2= "<< W_d_J1_2 <<", W_d_J1_J2= "<< W_d_J1_J2
<<", W_d_J2_2= "<< W_d_J2_2<< " ";
};
#endif
// puis partie volumique
double unSurV2=unSurV*unSurV;
double VmoinsUn = V-1.;
switch (type_pot_vol)
{case 1: {W_v = K_use * (1- (1+logV)*unSurV); // potentiel 1
W_v_J3 = 0.5 * K_use * logV * unSurV * unSurV2; // variation / J3
// calcul des variations secondes / (J3)^2
W_v_J3J3 = 0.25 * K_use * unSurV2 * unSurV2 * unSurV *(1.-3.*logV);
module_compressibilite = K_use * unSurV2;
break;
}
case 2: {W_v = 0.5 * K_use * VmoinsUn; // potentiel 2
W_v_J3 = 0.25 * K_use * unSurV; // variation / J3
// calcul des variations secondes / (J3)^2
W_v_J3J3 = -K_use * 0.125 * unSurV2 * unSurV;
if (Dabs(VmoinsUn) > ConstMath::petit )
{module_compressibilite = 0.5 * K_use / logV;}
else // on fait un développement limité du log pour V proche de 1.
{module_compressibilite = 0.5 * K_use / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
break;
}
case 3: {W_v = K_use * 0.5 * logV * logV; // potentiel 3
W_v_J3 = K_use * 0.5 * unSurV2 * logV; // variation / J3
// calcul des variations secondes / (J3)^2
W_v_J3J3 = K_use * 0.25 * unSurV2 * unSurV2 * (1.-2.*logV);
module_compressibilite = K_use / V ;
break;
}
case 4: {W_v = 0.5 * K_use * VmoinsUn * VmoinsUn; // potentiel
W_v_J3 = 0.5 * K_use * unSurV * VmoinsUn; // variation / J3
// calcul des variations secondes / (J3)^2
W_v_J3J3 = K_use * 0.25 * unSurV2 * unSurV;
if (Dabs(VmoinsUn) > ConstMath::petit )
{module_compressibilite = VmoinsUn * K_use /logV ;}
else // on fait un développement limité du log pour V proche de 1.
{module_compressibilite = VmoinsUn * K_use / (VmoinsUn-VmoinsUn*VmoinsUn/2.+VmoinsUn*VmoinsUn*VmoinsUn/3.);}
break;
}
};
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{switch (type_pot_vol)
{case 1: case 2: case 3: case 4:
cout << "\n W_v= " << W_v << ", W_v_J3= "<< W_v_J3 << ", W_v_J3J3= "<< W_v_J3J3
<< ", K_s= " << module_compressibilite << " ";
break;
};
};
#endif
// -- partie courbure éventuelle ---
if (avec_courbure)
{double unSurV3=unSurV2*unSurV;
double aP2r = pow(a_courbure,2.*r_courbure);
double J_rMoins3P2rMoins1 = pow((J_r(1)-3.),(2.*r_courbure-1));
double J_rMoins3P2r = J_rMoins3P2rMoins1 * (J_r(1)-3.);
W_c = unSurV * ( J_rMoins3P2r * (J_r(1)-3.) / aP2r / (2.*r_courbure+1) );
W_c_J1 = unSurV * ( J_rMoins3P2r/aP2r );
W_c_J3 = -0.5*unSurV3 * ( J_rMoins3P2r * (J_r(1)-3.) /aP2r/(2.*r_courbure+1) );
W_c_J1_2 = unSurV * r_courbure * J_rMoins3P2rMoins1 / aP2r;
W_c_J3_2 = 3./4. * unSurV3 * unSurV2 * ( J_rMoins3P2r * (J_r(1)-3.) /aP2r/(2.*r_courbure+1) );
W_c_J1_J3 = -0.5*unSurV3 * ( J_rMoins3P2r/aP2r );
module_compressibilite += 2. * V * W_c_J3 / (MaX(ConstMath::petit, log(V) ));
}
else
{W_c = W_c_J1 = W_c_J3 = W_c_J1_2 = W_c_J3_2 = W_c_J1_J3 =0.;};
// cout << "\n Jr= " << J_r(1) << " " << J_r(2) << " " << J_r(3) << " pot " << W_d << " " << W_v ;
};
// définition de la courbe représentant l'évolution de l'énergie en fonction de J1
void Hart_Smith3D::Calcul_courbe_evolW_J1()
{ // le potentiel pour la partie ne dépendant pas de la variation de volume:
// W(J1,J2)=C1*int_3^{J1(finale)) [exp(C3*(J1-3)**2)] dJ1 + C2 * log(J2/3)
// donc la fontion qui nous intéresse c'est l'intégrale, qui ne dépend que de J1
// integ = int_3^{J1(finale)) [exp(C3*(J1-3)**2)] dJ1
//
// dans une première étape on constitue un tableau de n points de 3 à 10
// -> gregory: et pourquoi pas inf à 3
// ensuite on interpolera sur les points pour une valeurs quelconque
// de 3 à 6 => 30 points
int nb1=30;
double x_1=3.,x_2=6.;
// de 6 à 10 => 20 points
int nb2=20;
double x_3=10;
// l'intégration est faite par trapèze basique avec 5 points intermédiaires
int nbtra=5;
Tableau <Coordonnee2> points(nb1+nb2+1);
double integral = 0.; // la première valeur est nulle
double x=3.; // coordonnée courante
double delta_x = (x_2-x_1)/(nb1*nbtra); // le pas de la première vague de points
double foncti=0.; // valeur courante de la fonction à intégrer, initialement nulle
// calcul des nb1 premiers points
for (int i=1; i<=nb1; i++)
{ for (int j=1;j<=nbtra;j++)
{ integral += 0.5 * foncti;
x += delta_x;
foncti = exp(C3_use*Sqr(x-3.));
integral += 0.5 * foncti;
};
points(i+1)(1)=x; points(i+1)(2)=integral;
};
// calcul des nb2 seconds points
delta_x = (x_3-x_2)/(nb2*nbtra); // le pas de la seconde vague de point
for (int i=1; i<=nb2; i++)
{ for (int j=1;j<=nbtra;j++)
{ integral += 0.5 * foncti;
x += delta_x;
foncti = exp(C3_use*Sqr(x-3.));
integral += 0.5 * foncti;
};
points(i+nb1+1)(1)=x; points(i+nb1+1)(2)=integral;
};
// on crée la fonction linéaire équivalente
int_J1.Change_tabPoints(points);
};
// calcul de la dérivée numérique de la contrainte
void Hart_Smith3D::Cal_dsigma_deps_num (const TenseurBB & gijBB_0_,const TenseurHH & gijHH_0_
,const TenseurBB & gijBB_tdt_,const TenseurHH & gijHH_tdt_
,const double& jacobien_0,const double& jacobien
,Tenseur3HHHH& dSigdepsHHHH)
{ 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
double delta = ConstMath::unpeupetit*10.;
double unSurDelta = 1./delta;
// init coefficient de la loi
C1_use = C1; C2_use = C2; C3_use = C3; K_use = K;
// cas de la thermo dépendance, on calcul les grandeurs
if (C1_temperature != NULL) C1_use = C1_temperature->Valeur(*temperature);
if (C2_temperature != NULL) C2_use = C2_temperature->Valeur(*temperature);
if (C3_temperature != NULL) C3_use = C3_temperature->Valeur(*temperature);
if (K_temperature != NULL) K_use = K_temperature->Valeur(*temperature);
if (avec_courbure)
{ if (a_temperature != NULL) a_courbure = a_temperature->Valeur(*temperature);
if (r_temperature != NULL) r_courbure = r_temperature->Valeur(*temperature);
};
// cas des contraintes et de ses variations analytiques
// Tenseur3HHHH dSigdepsHHHH; // le tenseur contenant les dérivées analytiques
Tenseur3HH SigmaHH_deb;
Cal_sigmaEtDer_pour_num(gijBB_0_,gijHH_0_,gijBB_tdt_,gijHH_tdt_
,jacobien_0,jacobien,SigmaHH_deb,dSigdepsHHHH);
// dimensionnement pour la matrice numérique
Tenseur3HHHH dSigdepsHHHH_num;
// 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());
// cas des contraintes
Tenseur3HH SigmaHH_N;
Cal_sigma_pour_num(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
,jacobien_0,jacobien_N,SigmaHH_N);
// calcul des dérivées numériques et comparaisons
for (int k=1;k<=3;k++)
for (int l=1;l<=3;l++)
{ //
double derSigNum = coef * 2.*(SigmaHH_N(k,l) - SigmaHH_deb(k,l) )*unSurDelta;
dSigdepsHHHH_num.Change(k,l,i,j,derSigNum);
double derSigAna = dSigdepsHHHH(k,l,i,j);//0.5*(dSigdepsHHHH(k,l,i,j) + dSigdepsHHHH(k,l,j,i));
bool erreur = false;
if (diffpourcent(derSigNum,derSigAna,MaX(Dabs(derSigNum),Dabs(derSigAna)),0.1))
if (MaX(Dabs(derSigNum),Dabs(derSigAna)) > 200.)
{if (MiN(Dabs(derSigNum),Dabs(derSigAna)) == 0.)
{if ( MaX(Dabs(derSigNum),Dabs(derSigAna)) > 50.*delta) erreur = true;}
else erreur = true;
};
// erreur = false; // a virer
if (erreur)
{
// calcul des dérivées d'éléments intermédiaires pour voir
//
cout << "\n erreur dans le calcul analytique de l'operateur tangent "
<< "\n derSigNum= " << derSigNum << " derSigAna= " << derSigAna
<< " klij= " << k << " " << l << " " << i << " " << j
<< " SigmaHH_N(k,l)= " << SigmaHH_N(k,l);
cout << "\n Hart_Smith3D::Calcul_derivee_numerique(..";
cout << "\n un caractere ";
// --- pour le débug ----
// calcul des invariants et de leurs variations premières en numérique
Invariants_et_var1_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
,jacobien_0,jacobien_N);
// calcul des invariants et de leurs variations premières et secondes
Invariants_et_var2_deb(gijBB_0_,gijHH_0_,(const TenseurBB &) gijBBtdt_N,(const TenseurHH &) gijHHtdt_N
,jacobien_0,jacobien_N);
string toto;
toto=lect_chaine();
};
};
};
// passage des dérivées numériques aux dérivées finales
dSigdepsHHHH= dSigdepsHHHH_num;
};
// calcul de la contrainte avec le minimum de variable de passage, utilisé pour le numérique
void Hart_Smith3D::Cal_sigma_pour_num(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0
,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt
,const double& jacobien_0,const double& jacobien,TenseurHH & sigHH_)
{
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage en dim 3 explicite
// calcul des invariants et de leurs variations premières (méthode de Hyper_W_gene_3D)
// Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// pour vérif on appelle var2, mais c'est à virer
Invariants_et_var1(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
double module_compressibilite; // ne sert pas ici
Potentiel_et_var(module_compressibilite);
// calcul du tenseur des contraintes
sigHH = ((W_d_J1+W_c_J1)/V) * d_J_r_epsBB_HH(1) + (W_d_J2/V) * d_J_r_epsBB_HH(2)
+ ((W_v_J3+W_c_J3)/V) * d_J_r_epsBB_HH(3);
};
// idem avec la variation
void Hart_Smith3D::Cal_sigmaEtDer_pour_num(const TenseurBB & gijBB_0,const TenseurHH & gijHH_0
,const TenseurBB & gijBB_tdt,const TenseurHH & gijHH_tdt
,const double& jacobien_0,const double& jacobien
,TenseurHH & sigHH_,Tenseur3HHHH& dSigdepsHHHH)
{
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // passage en dim 3 explicite
// calcul des invariants et de leurs variations premières et seconde (méthode de Hyper_W_gene_3D)
Invariants_et_var2(gijBB_0,gijHH_0,gijBB_tdt,gijHH_tdt,jacobien_0,jacobien);
// calcul du potentiel et de ses dérivées premières / aux invariants J_r
double module_compressibilite; // ne sert pas ici
Potentiel_et_var2(module_compressibilite);
// calcul du tenseur des contraintes
double unSurV=1./V;
sigHH = ((W_d_J1+W_c_J1)*unSurV) * d_J_r_epsBB_HH(1) + (W_d_J2*unSurV) * d_J_r_epsBB_HH(2)
+ ((W_v_J3+W_c_J3)*unSurV) * d_J_r_epsBB_HH(3);
// calcul de la variation seconde du potentiel par rapport à epsij epskl
// !!!! je pense que la formule qui suit est fausse donc si on s'en sert il faut y re-regarder !!
// revoir par rapport aux formules générales ??
cout << "\n revoir la formule de vérif !! Hart_Smith3D::Cal_sigmaEtDer_pour_num " << endl ;
Sortie(1);
Tenseur3HHHH d2W_d2epsHHHH
= W_d_J1 * d_J_1_eps2BB_HHHH + W_d_J2 * d_J_2_eps2BB_HHHH
+ W_v_J3 * d_J_3_eps2BB_HHHH
+ W_v_J3J3 * Tenseur3HHHH::Prod_tensoriel(d_J_r_epsBB_HH(3),d_J_r_epsBB_HH(3));
// cout << "\n d2W_d2epsHHHH(1,1,1,1)= " << d2W_d2epsHHHH(1,1,1,1);
// calcul de la variation du tenseur des contraintes par rapports aux déformations
// on tient compte du fait que V*sigHH = d W/ d epsij
Tenseur3HH interHH = -sigHH ; //* (-0.5*unSurV*unSurV);
Tenseur3HHHH d_igdepsHHHH(1,interHH,gijHH_tdt);
d_igdepsHHHH += (unSurV) * d2W_d2epsHHHH;
dSigdepsHHHH = d_igdepsHHHH;
};
Reference in a new issue View git blame Copy permalink