rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/comportement/anisotropie/Projection_anisotrope_3D.cc

3537 lines
153 KiB
C++
Raw Permalink Normal View History

ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// FICHIER : Projection_anisotrope_3D.cp
// CLASSE : Projection_anisotrope_3D
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
//#include "Debug.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "ConstMath.h"
#include "LesLoisDeComp.h"
#include "Projection_anisotrope_3D.h"
#include "NevezTenseur.h"
#include "MathUtil.h"
#include "Util.h"
#include "Enum_TypeQuelconque.h"
#include "TypeQuelconqueParticulier.h"
#include "TenseurQ3gene.h"
#include "NevezTenseurQ.h"
#include "MathUtil.h"
#include "MathUtil2.h"
#include "CharUtil.h"
// ========== fonctions pour la classe de sauvegarde des résultats =========
// constructeur par défaut à ne pas utiliser
Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::SaveResulProjection_anisotrope_3D():
O_B(NULL), O_H(NULL),Op_H(3,3),Op_H_t(3,3),eps_loc_HH(NULL),sig_loc_HH(NULL)
,para_loi(NULL)
{
// non on ne fait pas de définition par défaut
// if (type_transport == 0)
// {O_H = new BaseH(); sig_loc_HH = NevezTenseurHH(3); }
// else
// {O_B = new BaseB(); eps_loc_HH = NevezTenseurHH(3); };
cout << "\n erreur, le constructeur par defaut ne doit pas etre utilise !"
<< "\n Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::SaveResulProjection_anisotrope_3D()";
cout << "\n Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::SaveResulProjection_anisotrope_3D()";
Sortie(1);
};
// le constructeur courant
Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::SaveResulProjection_anisotrope_3D
(SaveResul* l_des_SaveResul,TenseurHH* l_siHH
,TenseurHH* l_siHH_t
,EnergieMeca l_energ_,EnergieMeca l_energ_t_
,bool avec_ponderation):
O_B(NULL), O_H(NULL),Op_H(3,3),Op_H_t(3,3),eps_loc_HH(NULL),sig_loc_HH(NULL)
,para_loi(NULL)
,SaveResul_interne(NULL) ,l_sigoHH(),l_sigoHH_t(),l_energ(l_energ_),l_energ_t(l_energ_t_)
,f_ponder(1.),f_ponder_t(1.)
{ SaveResul * nevez_save_result=NULL;
if (l_des_SaveResul != NULL) nevez_save_result = l_des_SaveResul->Nevez_SaveResul();
SaveResul_interne = nevez_save_result;
TenseurHH * interHH=NULL;
// dans le cas où le tenseur passé en paramètre est non null on s'en sert
if (l_siHH_t != NULL) {interHH=NevezTenseurHH(*l_siHH_t);};
l_sigoHH= interHH;
TenseurHH * interHH_t=NULL;
// idem interHH, dans le cas où le tenseur passé en paramètre est non null on s'en sert
if (l_siHH_t != NULL) {interHH_t =NevezTenseurHH(*l_siHH_t);};
l_sigoHH_t=interHH_t;
};
// constructeur de copie
Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D
::SaveResulProjection_anisotrope_3D(const SaveResulProjection_anisotrope_3D& sav):
O_B(NULL), O_H(NULL),Op_H(3,3),Op_H_t(3,3),eps_loc_HH(NULL),sig_loc_HH(NULL)
,para_loi(NULL)
,SaveResul_interne(NULL) ,l_sigoHH(NULL),l_sigoHH_t(NULL)
,l_energ(sav.l_energ),l_energ_t(sav.l_energ_t)
,f_ponder(sav.f_ponder),f_ponder_t(sav.f_ponder_t)
{
// partie repère d'anisotropie
if (sav.O_B != NULL)
O_B = new BaseB(*sav.O_B);
if (sav.O_H != NULL)
O_H = new BaseH(*sav.O_H);
if (sav.eps_loc_HH != NULL)
eps_loc_HH = NevezTenseurHH(*sav.eps_loc_HH);
if (sav.sig_loc_HH != NULL)
sig_loc_HH = NevezTenseurHH(*sav.sig_loc_HH);
if (sav.para_loi != NULL)
para_loi = new Vecteur (*sav.para_loi);
// loi interne
if (sav.SaveResul_interne != NULL) SaveResul_interne = sav.SaveResul_interne->Nevez_SaveResul();
// dans le cas où le tenseur passé en paramètre est non null on s'en sert
if (sav.l_sigoHH != NULL) {l_sigoHH=NevezTenseurHH(*(sav.l_sigoHH));};
// idem interHH, dans le cas où le tenseur passé en paramètre est non null on s'en sert
if (sav.l_sigoHH_t != NULL) {l_sigoHH_t =NevezTenseurHH(*(sav.l_sigoHH_t));};
// pour les pondérations
};
// destructeur
Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::~SaveResulProjection_anisotrope_3D()
{
// partie repère d'anisotropie
if (O_B != NULL)
delete O_B;
if (O_H != NULL)
delete O_H;
if (eps_loc_HH != NULL)
delete eps_loc_HH;
if (sig_loc_HH != NULL)
delete sig_loc_HH;
if (para_loi != NULL)
delete para_loi;
// loi interne
if (SaveResul_interne != NULL) delete SaveResul_interne;
if(l_sigoHH != NULL) delete l_sigoHH; if(l_sigoHH_t != NULL) delete l_sigoHH_t;
};
// affectation
Loi_comp_abstraite::SaveResul & Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::operator = ( const Loi_comp_abstraite::SaveResul & a)
{ Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D& sav = *((Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D*) &a);
// partie repère d'anisotropie
if (sav.O_B != NULL)
{if (O_B == NULL) {O_B = new BaseB(*sav.O_B);}
else {*O_B = *sav.O_B;};
}
else // sinon cas NULL
{if (O_B != NULL) {delete O_B;O_B=NULL;}} ;
if (sav.O_H != NULL)
{if (O_H == NULL) {O_H = new BaseH(*sav.O_H);}
else {*O_H = *sav.O_H;};
}
else // sinon cas NULL
{if (O_H != NULL) {delete O_H;O_H=NULL;}} ;
if (sav.eps_loc_HH != NULL)
{if (eps_loc_HH == NULL) {eps_loc_HH = NevezTenseurHH(*sav.eps_loc_HH);}
else {*eps_loc_HH = *sav.eps_loc_HH;};
}
else // sinon cas NULL
{if (eps_loc_HH != NULL) {delete eps_loc_HH;eps_loc_HH=NULL;}} ;
if (sav.sig_loc_HH != NULL)
{if (sig_loc_HH == NULL) {sig_loc_HH = NevezTenseurHH(*sav.sig_loc_HH);}
else {*sig_loc_HH = *sav.sig_loc_HH;};
}
else // sinon cas NULL
{if (sig_loc_HH != NULL) {delete sig_loc_HH; sig_loc_HH= NULL;};
} ;
if (sav.para_loi != NULL)
{if (para_loi == NULL) para_loi = new Vecteur (*sav.para_loi);
else {*para_loi = *sav.para_loi; };
}
else // sinon cas NULL
{if (para_loi != NULL) {delete para_loi;para_loi=NULL;}} ;
// et la base
Op_H = sav.Op_H;
Op_H_t = sav.Op_H_t;
// loi interne
if (SaveResul_interne != NULL)
{if (sav.SaveResul_interne == NULL) {delete SaveResul_interne;SaveResul_interne=NULL;}
else // on considère que les deux sont de même type
{*SaveResul_interne = *(sav.SaveResul_interne);};
}
else
{if (sav.SaveResul_interne != NULL)
{ SaveResul_interne = sav.SaveResul_interne->Nevez_SaveResul();
};
// sinon il n'y a rien à faire car les deux sont nulles
};
if (l_sigoHH != NULL)
{if (sav.l_sigoHH == NULL) {delete l_sigoHH;l_sigoHH=NULL;}
else // on considère que les deux sont de même type
{*l_sigoHH = *(sav.l_sigoHH);};
}
else
{if (sav.l_sigoHH != NULL)
{ *l_sigoHH = *(sav.l_sigoHH);
};
// sinon il n'y a rien à faire car les deux sont nulles
};
if (l_sigoHH_t != NULL)
{if (sav.l_sigoHH_t == NULL) {delete l_sigoHH_t;l_sigoHH_t=NULL;}
else // on considère que les deux sont de même type
{*l_sigoHH_t = *(sav.l_sigoHH_t);};
}
else
{if (sav.l_sigoHH_t != NULL)
{ *l_sigoHH_t = *(sav.l_sigoHH_t);
};
// sinon il n'y a rien à faire car les deux sont nulles
};
// // pour les autres conteneurs internes, affectation directe, car a priori pas de pb
l_energ = sav.l_energ;
l_energ_t = sav.l_energ_t;
f_ponder = sav.f_ponder;
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 Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Lecture_base_info(ifstream& ent,const int cas)
{ string nom;
ent >> nom ;
#ifdef MISE_AU_POINT
if (nom != "S_Proj_aniso_3D")
{ cout << "\nErreur : on attendait le mot cle: S_Proj_aniso_3D "
<< " et on a lue "<<nom;
cout << " Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Lecture_base_info\n";
Sortie(1);
};
#endif
// tout d'abord ce qui varie tout le temps
// lecture pondération
string toto;
ent >> toto >> f_ponder; f_ponder_t = f_ponder;
// données de la loi
bool avec_save_result = false;
ent >> toto >> avec_save_result;
if (avec_save_result)
{ if(SaveResul_interne != NULL)
{SaveResul_interne->Lecture_base_info(ent,cas);}
else
{ cout << "\n *** erreur en lecture des infos specifiques a la loi de comportement interne "
<< " il n'y a pas de conteneur definit ! on ne peut pas continuer "
<< "\n Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Lecture_base_info(..." << flush;
Sortie(1);
};
};
// lecture de la contrainte sauvegardée
ent >> toto; // le mot clé
(l_sigoHH_t)->Lecture(ent); // lecture du tenseur
*l_sigoHH = *l_sigoHH_t;
// lecture des énergies
ent >> l_energ_t;l_energ =l_energ_t;
// puis on lit que ce qui est à 0 et à t, ce qui est à tdt est considéré une zone de travail
switch (cas)
{ case 1 :
{
// la base initiale
if (O_B != NULL)
{ ent >> nom >> *O_B;}
else if (O_H != NULL)
{ ent >> nom >> *O_H;}
// la base convectée
ent >> nom >> Op_H_t ;
Op_H = Op_H_t;
// les stockages conditionnels
int test;
// -- // cas d'eps_loc
ent >> nom >> test;
if (test == 0)
{ if (eps_loc_HH != NULL)
{delete eps_loc_HH;eps_loc_HH=NULL;}
}
else // sinon cas où on a des données additionnelles
{ if (eps_loc_HH == NULL) eps_loc_HH = NevezTenseurHH(3);
eps_loc_HH->Lecture(ent);
};
// -- // cas de sig_loc
ent >> nom >> test;
if (test == 0)
{ if (sig_loc_HH != NULL)
{delete sig_loc_HH;sig_loc_HH=NULL;}
}
else // sinon cas où on a des données additionnelles
{ if (sig_loc_HH == NULL) sig_loc_HH = NevezTenseurHH(3);
sig_loc_HH->Lecture(ent);
};
// -- // cas des para loi
ent >> nom >> test;
if (test == 0)
{ if (para_loi != NULL)
{delete para_loi;para_loi=NULL;}
}
else // sinon cas où on a des données additionnelles
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ if (para_loi == NULL) para_loi = new Vecteur (6);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
ent >> *para_loi;
};
break;
}
case 2 :
{
// la base convectée
ent >> nom >> Op_H_t ;
// les stockages conditionnels
int test;
ent >> nom >> test;
if (test == 0)
{ if (eps_loc_HH != NULL)
{delete eps_loc_HH;eps_loc_HH=NULL;}
}
else // sinon cas où on a des données additionnelles
{ if (eps_loc_HH == NULL) eps_loc_HH = NevezTenseurHH(3);
eps_loc_HH->Lecture(ent);
};
ent >> nom >> test;
if (test == 0)
{ if (sig_loc_HH != NULL)
{delete sig_loc_HH;sig_loc_HH=NULL;}
}
else // sinon cas où on a des données additionnelles
{ if (sig_loc_HH == NULL) sig_loc_HH = NevezTenseurHH(3);
sig_loc_HH->Lecture(ent);
};
ent >> nom >> test;
if (test == 0)
{ if (para_loi != NULL)
{delete para_loi;para_loi=NULL;}
}
else // sinon cas où on a des données additionnelles
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ if (para_loi == NULL) para_loi = new Vecteur(6);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
ent >> *para_loi;
};
break;
}
default:
cout << "\n cas non considere !!: cas= " << cas
<< "\n Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Ecriture_base_info(...";
Sortie(1);
};
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables
//(supposées comme telles)
void Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Ecriture_base_info(ofstream& sort,const int cas )
{ sort << "\n S_Proj_aniso_3D ";
// tout d'abord ce qui varie tout le temps
// pondération ?
sort << " Pond= "<< f_ponder_t << " ";
// on regarde s'il y a un save_result local
sort << " avec_save_result "<< (SaveResul_interne != NULL)<<" ";
if (SaveResul_interne != NULL)
SaveResul_interne->Ecriture_base_info(sort, cas);
// écriture de la contrainte sauvegardée
sort << " sig_intHH "; (l_sigoHH)->Ecriture(sort);
// écriture des énergies
sort << l_energ;
// puis on ne sauvegarde que ce qui est à 0 et à t, ce qui est à tdt est considéré une zone de travail
switch (cas)
{ case 1 :
{
// la base initiale
if (O_B != NULL)
{ sort << " O_B: "<< *O_B;}
else if (O_H != NULL)
{ sort << " O_H: "<< *O_H;}
// la base convectée
sort << " Op_H_t: "<< Op_H_t << " ";
// les stockages conditionnels
if (eps_loc_HH != NULL)
{ sort << "\n addi_eps_loc: 1 ";
eps_loc_HH->Ecriture(sort);
}
else
{ sort << " addi_eps_loc: 0 ";};
if (sig_loc_HH != NULL)
{ sort << "\n addi_sig_loc: 1 " ;
sig_loc_HH->Ecriture(sort);
}
else
{ sort << " addi_sig_loc: 0 ";};
if (para_loi != NULL)
{ sort << " addi_para_loi: 1 ";
sort << *para_loi << " ";
}
else
{ sort << " addi_para_loi: 0 ";};
break;
}
case 2 :
{
// la base convectée
sort << " Op_H_t: "<< Op_H_t << " ";
// les stockages conditionnels
if (eps_loc_HH != NULL)
{ sort << "\n addi_eps_loc: 1 ";
eps_loc_HH->Ecriture(sort);
}
else
{ sort << " addi_eps_loc: 0 ";};
if (sig_loc_HH != NULL)
{ sort << "\n addi_sig_loc: 1 " ;
sig_loc_HH->Ecriture(sort);
}
else
{ sort << " addi_sig_loc: 0 ";};
if (para_loi != NULL)
{ sort << " addi_para_loi: 1 ";
sort << *para_loi << " ";
}
else
{ sort << " addi_para_loi: 0 ";};
break;
}
default:
cout << "\n cas non considere !!: cas= " << cas
<< "\n Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Ecriture_base_info(...";
Sortie(1);
};
};
// mise à jour des informations transitoires
void Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::TdtversT()
{ Op_H_t = Op_H;
f_ponder_t = f_ponder;
(*l_sigoHH_t) = (*l_sigoHH);
l_energ_t = l_energ;
if (SaveResul_interne != NULL)
SaveResul_interne->TdtversT();
};
void Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::TversTdt()
{ Op_H = Op_H_t;
f_ponder = f_ponder_t;
(*l_sigoHH) = (*l_sigoHH_t);
l_energ = l_energ_t;
if (SaveResul_interne != NULL)
SaveResul_interne->TversTdt();
};
// affichage à l'écran des infos
void Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Affiche() const
{ cout << "\n SaveResulProjection_anisotrope_3D: " ;
cout << "\n SaveResul_loi_interne: ";
if (SaveResul_interne != NULL)
SaveResul_interne->Affiche();
cout << "\n l_sigoHH: ";
l_sigoHH->Ecriture(cout);
cout << "\n l_sigoHH_t: ";
l_sigoHH_t->Ecriture(cout);
cout << "\n l_energ: " << l_energ;
cout << "\n l_energ_t: " << l_energ_t;
cout << "\n f_ponder: " << f_ponder;
cout << "\n f_ponder_t: " << f_ponder_t;
// la partie base d'anisotropie
// la base initiale
if (O_B != NULL)
{ cout << " O_B: "<< *O_B;}
else if (O_H != NULL)
{ cout << " O_H: "<< *O_H;}
// la base convectée
cout << " Op_H_t: "<< Op_H_t << " ";
// les stockages conditionnels
if (eps_loc_HH != NULL)
{ cout << " addi: 1 " << "\n eps_loc_HH: ";
eps_loc_HH->Ecriture(cout);
cout << "\n sig_loc_HH: ";
sig_loc_HH->Ecriture(cout);
cout << "\n para_loi: " << *para_loi << " ";
}
else
{ cout << " addi: 0 ";};
// la base convectée
cout << " Op_H_t: "<< Op_H_t << " ";
// les stockages conditionnels
if (eps_loc_HH != NULL)
{ cout << "\n addi: 1 "<< "\n eps_loc_HH: " ;
eps_loc_HH->Ecriture(cout);
cout << "\n sig_loc_HH: " ;
sig_loc_HH->Ecriture(cout);
cout << "\n para_loi: " << *para_loi << " ";
}
else
{ cout << "\n addi: 0 ";};
cout << "\n .. fin SaveResulProjection_anisotrope_3D ... ";
};
//changement de base de toutes les grandeurs internes tensorielles stockées
// b(i,j) represente les coordonnees de la nouvelle base naturelle gpB dans l'ancienne gB
// gpB(i) = beta(i,j) * gB(j), i indice de ligne, j indice de colonne
// gpH(i) = gamma(i,j) * gH(j)
void Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::ChBase_des_grandeurs(const Mat_pleine& beta,const Mat_pleine& gamma)
{ // on ne s'intéresse qu'aux grandeurs tensorielles
// il faut que l'on exprime le repère d'anisotropie dans le nouveau repère
if (O_B != NULL)
{for (int i=1;i<4;i++)
{ CoordonneeB Ap_B(3); // inter
MathUtil2::ChBase(O_B->CoordoB(i),beta,Ap_B);
O_B->CoordoB(i)=Ap_B;
}
};
if (O_H != NULL)
{for (int i=1;i<4;i++)
{ CoordonneeH Ap_H(3); // inter
MathUtil2::ChBase(O_H->CoordoH(i),gamma,Ap_H);
O_H->CoordoH(i)=Ap_H;
}
};
// les autres grandeurs sont issues du calcul
// normalement elles sont calculées en même temps que les contraintes
// mais on les change de repère néanmoins, au cas on on ferait une sortie
// d'info sans calcul (c'est pénalisant pour le cas avec calcul, mais cela
// évite peut-être des pb potentiels ??)
// Op_H est toujours calculé
{for (int i=1;i<4;i++)
{ CoordonneeH Ap_H(3); // inter
MathUtil2::ChBase(Op_H_t.CoordoH(i),gamma,Ap_H);
Op_H_t.CoordoH(i)=Ap_H;
}
};
// cas des grandeurs locales si elles sont sauvegardées
if (eps_loc_HH != NULL)
{ eps_loc_HH->ChBase(gamma);};
if (sig_loc_HH != NULL)
{ sig_loc_HH->ChBase(gamma);};
// partie relative à la loi interne
if (SaveResul_interne != NULL)
SaveResul_interne->ChBase_des_grandeurs(beta,gamma);
// les tenseurs
l_sigoHH->ChBase(gamma);
l_sigoHH_t->ChBase(gamma);
};
// procedure permettant de completer éventuellement les données particulières
// de la loi stockées
// au niveau du point d'intégration par exemple: exemple: un repère d'anisotropie
// completer est appelé apres sa creation avec les donnees du bloc transmis
// peut etre appeler plusieurs fois
Loi_comp_abstraite::SaveResul* Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D
::Complete_SaveResul(const BlocGen & bloc, const Tableau <Coordonnee>& tab_coor
,const Loi_comp_abstraite* loi)
{// on regarde s'il s'agit d'un repère d'anisotropie
if (bloc.Nom(1) == "repere_anisotropie_")
{// ensuite on vérifie le nom de l'identificateur
Projection_anisotrope_3D* loiOrtho = (Projection_anisotrope_3D*) loi;
// dimensionnement des bases
if (loiOrtho->Type_transport() == 0)
{O_H = new BaseH(); sig_loc_HH = NevezTenseurHH(3); }
else
{O_B = new BaseB(); eps_loc_HH = NevezTenseurHH(3); };
// // on crée sig et eps loc, les deux pour éviter des pb de sortie par la suite
// sig_loc_HH = NevezTenseurHH(3);
// eps_loc_HH = NevezTenseurHH(3);
// récupération du repère
if (bloc.Nom(2) == loiOrtho->NomRepere())
{// c'est le bon, récupération du repère
if (O_B != NULL)
{for (int i=1;i<4;i++)
O_B->CoordoB(i).Change_val(tab_coor(i));
};
if (O_H != NULL)
{for (int i=1;i<4;i++)
O_H->CoordoH(i).Change_val(tab_coor(i));
};
};
};
// puis idem pour la loi interne
// on transmet au conteneur 3D interne
const Projection_anisotrope_3D * loi_CPtt = (const Projection_anisotrope_3D*) loi;
// loi_CPtt : est un nom quelconque
if (SaveResul_interne != NULL)
SaveResul_interne->Complete_SaveResul(bloc,tab_coor, loi_CPtt->loi_interne);
return this;
};
// initialise les informations de travail concernant le pas de temps en cours
void Projection_anisotrope_3D::SaveResulProjection_anisotrope_3D::Init_debut_calcul()
{ };
// ========== fin des fonctions pour la classe de sauvegarde des résultats =========
Projection_anisotrope_3D::Projection_anisotrope_3D () : // Constructeur par defaut
Loi_comp_abstraite(PROJECTION_ANISOTROPE_3D,CAT_MECANIQUE,3)
,hij(3,3),fct_para(6),cas_calcul(0),type_transport(0)
,A1(ConstMath::tresgrand),A2(ConstMath::tresgrand),A3(ConstMath::tresgrand)
,B12(ConstMath::tresgrand),B13(ConstMath::tresgrand),B23(ConstMath::tresgrand)
,ratio_inf_module_compressibilite(0.01)
,sortie_post(0)
,nom_repere("")
,inv_loi(3,3),Op_B(3,3),d_Op_B(3,3),d_Op_H(3,3),pO_B(3,3),pO_H(3,3)
,beta_inv(3,3),beta(3,3),gamma(3,3),beta_transpose(3,3),gamma_transpose(3,3)
,alpha_H(3,3)
,I_x_I_HHHH(),I_xbarre_I_HHHH(),I_x_eps_HHHH(),Ixbarre_eps_HHHH()
,dsig_ijkl_HHHH()
// loi interne
,loi_interne(NULL)
,completude_calcul(CONTRAINTE_ET_TANGENT),d_sigma_deps_inter(NULL),d_sig_deps_3D_HHHH()
,ponderation_nD_quelconque(NULL)
,d_sigRef_HH()
// -- grandeurs de travail
,H_HHBB()
// un conteneur d'un point d'intégration courant
,ptintmeca(3)
{// def de hij
B12=ConstMath::tresgrand;B13=ConstMath::tresgrand;B23=ConstMath::tresgrand;
hij(1,1) = &A1;hij(2,2) = &A2;hij(3,3) = &A3;
hij(1,2) = &B12;hij(1,3) = &B13;hij(2,3) = &B23;
hij(2,1) = &B12;hij(3,1) = &B13;hij(3,2) = &B23;
// init fnD
for (int i=1;i<=6;i++)
fct_para(i) = NULL;
Bll_fct_para=1; // pour l'instant pas de fonction
// on ajoute les invariants au pt integ courant
ptintmeca.Change_statut_Invariants_contrainte (true);
};
// Constructeur de copie
Projection_anisotrope_3D::Projection_anisotrope_3D (const Projection_anisotrope_3D& loi) :
Loi_comp_abstraite(loi)
,hij(3,3),fct_para(loi.fct_para),cas_calcul(loi.cas_calcul)
,A1(loi.A1),A2(loi.A2),A3(loi.A3)
,ratio_inf_module_compressibilite(loi.ratio_inf_module_compressibilite)
,inv_loi(loi.inv_loi),Op_B(3,3),d_Op_B(3,3),d_Op_H(3,3),pO_B(3,3),pO_H(3,3)
,beta_inv(3,3),beta(3,3),gamma(3,3),beta_transpose(3,3),gamma_transpose(3,3)
,alpha_H(3,3)
,type_transport(loi.type_transport)
,sortie_post(loi.sortie_post)
,nom_repere(loi.nom_repere)
,I_x_I_HHHH(),I_xbarre_I_HHHH(),I_x_eps_HHHH(),Ixbarre_eps_HHHH()
,dsig_ijkl_HHHH()
// loi interne
,completude_calcul(loi.completude_calcul),d_sigma_deps_inter(NULL)
,d_sig_deps_3D_HHHH(loi.d_sig_deps_3D_HHHH)
,ponderation_nD_quelconque(NULL)
,d_sigRef_HH()
// grandeurs de travail
,H_HHBB(loi.H_HHBB)
// un conteneur d'un point d'intégration courant
,ptintmeca(loi.ptintmeca)
{// def de hij
B12=loi.B12;B13=loi.B13;B23=loi.B23;
hij(1,1) = &A1;hij(2,2) = &A2;hij(3,3) = &A3;
hij(1,2) = &B12;hij(1,3) = &B13;hij(2,3) = &B23;
hij(2,1) = &B12;hij(3,1) = &B13;hij(3,2) = &B23;
// init fnD
// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale
for (int i=1;i<=6;i++)
{if (fct_para(i) != NULL)
{Bll_fct_para=0; // dans tous les cas on indique qu'il y a des fonctions
if (fct_para(i)->NomFonction() == "_")
{// comme il s'agit d'une fonction locale on la redéfinie (sinon pb lors du destructeur de loi)
string non_fonction("_");
fct_para(i) = Fonction_nD::New_Fonction_nD(*fct_para(i));
};
};
};
// loi interne
if (loi.loi_interne != NULL)
loi_interne = loi.loi_interne->Nouvelle_loi_identique();
// pondération éventuelle
if (loi.ponderation_nD_quelconque != NULL)
ponderation_nD_quelconque = new Ponderation_TypeQuelconque(*(loi.ponderation_nD_quelconque));
if (loi.d_sigma_deps_inter != NULL)
d_sigma_deps_inter = NevezTenseurHHHH(*loi.d_sigma_deps_inter);
};
Projection_anisotrope_3D::~Projection_anisotrope_3D ()
// Destructeur
{ for (int i=1;i<=6;i++)
{if (fct_para(i) != NULL)
if (fct_para(i)->NomFonction() == "_") delete fct_para(i);
};
// loi interne
if (loi_interne != NULL) delete loi_interne;
if (d_sigma_deps_inter != NULL)
delete d_sigma_deps_inter;
if (d_sigRef_HH.Taille() != 0)
{ int taille = d_sigRef_HH.Taille();
for (int i=1; i<= taille; i++)
if (d_sigRef_HH(i) != NULL) delete d_sigRef_HH(i);
};
// pondération nD quelconques
if (ponderation_nD_quelconque != NULL)
delete ponderation_nD_quelconque;
};
// def d'une instance de données spécifiques, et initialisation
// valable une fois que la loi interne est définie
Projection_anisotrope_3D::SaveResul * Projection_anisotrope_3D::New_et_Initialise()
{ // petite vérif
if (loi_interne == NULL)
{cout << "\n *** erreur la loi interne n'est pas definie, on ne peut pas "
<< " creer le conteneur au pti associe !"<< flush;
Sortie(1);
};
// on crée les éléments nécessaires
SaveResul * nevez_save_result = loi_interne->New_et_Initialise();
// pour l'instant on définit des pointeurs de tenseurs qui ont la dimension de la loi,
// ensuite au moment de l'utilisation on les initialisera correctement si besoin
TenseurHH * l_sigoHH = NevezTenseurHH(dim,0.);
TenseurHH * l_sigoHH_t = NevezTenseurHH(dim,0.);
EnergieMeca l_energ_inter,l_energ_t_inter;
bool avec_ponderation = (ponderation_nD_quelconque!=NULL);
// on ramène la bonne instance
Projection_anisotrope_3D::SaveResul * retour = new SaveResulProjection_anisotrope_3D
(nevez_save_result,l_sigoHH,l_sigoHH_t
,l_energ_inter,l_energ_t_inter
,avec_ponderation);
// on supprime les grandeurs locales qui ont été crées par des new
delete nevez_save_result;delete l_sigoHH; delete l_sigoHH_t;
// retour
return retour;
};
// Lecture des donnees de la classe sur fichier
void Projection_anisotrope_3D::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{
// on lit les coefficients dans l'ordre
string nom_class_methode("Projection_anisotrope_3D::LectureDonneesParticulieres");
double val_defaut=0.;
double min = 0.; double max = -1; // max < min => la condition n'est pas prise en compte
string nom_pour_erreur("**erreur en lecture** Projection_anisotrope_3D:");
// -- lecture du type de projection
// on lit le nom du type de projection
{string nom_proj;
string mot_cle("TYPE_PROJ_ANISO_");
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_proj);
if (!lec )
{ entreePrinc->MessageBuffer(nom_pour_erreur+" mot cle "+mot_cle);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
if ( Type_Enum_proj_aniso_existe(nom_proj))
{ type_projection = Id_Nom_proj_aniso(nom_proj);
}
else
{ entreePrinc->MessageBuffer(nom_pour_erreur+" type de projection inexistant: "+mot_cle);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
// --- suivant le type de projection on lit les paramètres ad hoc ---
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
switch (type_projection)
{case PROJ_ORTHO:
{// pour faire une boucle de lecture on constitue un tableau de mots clés
Tableau <string > tab_mot_cle(6);
tab_mot_cle(1) = "A1";tab_mot_cle(2) = "A2";tab_mot_cle(3) = "A3";
tab_mot_cle(4) = "B12";tab_mot_cle(5) = "B13";tab_mot_cle(6) = "B23";
// puis un tableau pour les valeurs
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// on boucle sur les 6 coefficients
for (int i=1;i < 7;i++)
{string mot_cle1=tab_mot_cle(i)+"=";
string mot_cle2=tab_mot_cle(i)+"_fonction_nD:";
if(strstr(entreePrinc->tablcar,mot_cle2.c_str())==0)
{// lecture du paramètre
min = 0.; max = -1; // max < min => la condition n'est pas prise en compte
if (!entreePrinc->Lecture_un_parametre_double(val_defaut,nom_class_methode
,min,max,mot_cle1, *coef(i) ))
{ entreePrinc->MessageBuffer(nom_pour_erreur+" mot cle : "+mot_cle1);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
////----- debug ---
//cout << "\n debug Projection_anisotrope_3D::LectureDonneesParticulieres "
// << " mot cle : "+mot_cle1 << " = "<< *coef(i)
// << ", i= "<< i
// << flush;
////--- fin debug ---
// //----- debug ---
// cout << "\n debug Projection_anisotrope_3D::LectureDonneesParticulieres ";
// this->Affiche();
// //--- fin debug ---
}
else // on lit une fonction
{// on passe le mot clé générique
bool lec = entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle1);
// on lit le nom de la fonction
string nom_fonct;
lec = lec && entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle2,nom_fonct);
if (!lec )
{ entreePrinc->MessageBuffer(nom_pour_erreur+mot_cle2);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
Bll_fct_para=0; // on indique qu'il y a des fonctions
// maintenant on définit la fonction
if (lesFonctionsnD.Existe(nom_fonct))
{fct_para(i) = lesFonctionsnD.Trouve(nom_fonct);
}
else
{// sinon il faut la lire maintenant
string non("_");
fct_para(i) = Fonction_nD::New_Fonction_nD(non, Id_Nom_Fonction_nD(nom_fonct));
// lecture de la courbe
fct_para(i)->LectDonnParticulieres_Fonction_nD (non,entreePrinc);
// maintenant on vérifie que la fonction est utilisable
if (fct_para(i)->NbComposante() != 1 )
{ cout << "\n erreur en lecture, la fonction " << nom_fonct
<< " est une fonction vectorielle a " << fct_para(i)->NbComposante()
<< " composante alors qu'elle devrait etre scalaire ! "
<< " elle n'est donc pas utilisable !! ";
string message("\n**erreur08** \n"+nom_class_methode+"(...");
entreePrinc->MessageBuffer(message);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
};
if (i != 6)
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
};
}; // fin de la boucle for (int i=1;i < 7;i++)
}
break;
default:
cout << "\n *** le type de projection: " << Nom_proj_aniso(type_projection)
<< " n'est pas encore implemente, se plaindre !! " << flush;
entreePrinc->MessageBuffer("**erreur en lecture** Projection_anisotrope_3D ");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
break;
};
// --- lecture du repère d'anisotropie entraîné, associé ---
string mot_cle("nom_repere_associe_");
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
bool lec = entreePrinc->Lecture_mot_cle_et_string(nom_class_methode,mot_cle,nom_repere);
if (!lec )
{ entreePrinc->MessageBuffer(nom_pour_erreur+mot_cle);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// --- lecture de la loi associée ---
{// lecture du nom de la loi
string st2;
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
*(entreePrinc->entree) >> st2;
// definition de la loi
LoiAbstraiteGeneral * pt = LesLoisDeComp::Def_loi(st2);
// enregistrement de la loi
loi_interne = ((Loi_comp_abstraite*)pt);
// lecture des informations particulières propres à la loi
entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
pt->LectureDonneesParticulieres(entreePrinc,lesCourbes1D,lesFonctionsnD);
// si la loi est thermo dépendante on indique que la loi projection l'est aussi
if (((Loi_comp_abstraite*)pt)->ThermoDependante()) this->thermo_dependant = true;
// on lit le mot cle "fin_loi_interne"
string mot_cle("fin_loi_interne");
if (entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle))
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
};
// --- lecture éventuelle des paramètres de réglage ----
cas_calcul = 0; // par défaut
if(strstr(entreePrinc->tablcar,"avec_parametres_de_reglage_")!=0)
{string nom;
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
// on lit tant que l'on ne rencontre pas la ligne contenant "fin_parametres_reglage_"
// ou un nouveau mot clé global auquel cas il y a pb !!
MotCle motCle; // ref aux mots cle
while (strstr(entreePrinc->tablcar,"fin_parametres_reglage_")==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 reglage : on n'a pas trouve le mot cle "
<< " fin_parametres_reglage_ et par contre la ligne courante contient un mot cle global ";
entreePrinc->MessageBuffer(nom_pour_erreur);
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// lecture d'un mot clé
*(entreePrinc->entree) >> nom;
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) >>nom;
}
#else
else if ((entreePrinc->entree)->eof())
// la lecture est bonne mais on a atteind la fin de la ligne
{ if(nom != "fin_parametres_reglage_")
{entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom;};
}
#endif
else // cas d'une erreur de lecture
{ cout << "\n erreur de lecture inconnue ";
entreePrinc->MessageBuffer("** erreur4 des parametres de reglage de la loi de comportement de Projection_anisotrope_3D **");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// type de transport
if (nom == "type_transport_")
{ // lecture du type
*(entreePrinc->entree) >> type_transport;
if ((type_transport!=0)&&(type_transport!=1))
{ cout << "\n le type de transport lue pour la loi de Projection_anisotrope_3D: "<< type_transport
<< " n'est pas acceptable (uniquement 0 ou 1), on utilise le type par defaut (0)"
<< " qui correspond a un transport de type contravariant ";
type_transport = 0;
};
}
// forcer un affichage particulier pour les méthodes
else if (nom == "permet_affichage_")
{Lecture_permet_affichage(entreePrinc,lesFonctionsnD);
}
// lecture du ratio mini du module de compressibilité
else if (nom == "ratio_inf_module_compressibilite_")
{*(entreePrinc->entree) >> ratio_inf_module_compressibilite;
}
// on regarde si le calcul est éventuellement uniquement déviatorique
else if (nom == "seule_deviatorique")
{if (cas_calcul == 2) {cas_calcul=0;} else {cas_calcul = 1;};}
// idem pour la partie sphérique
else if (nom == "seule_spherique")
{if (cas_calcul == 1) {cas_calcul=0;} else {cas_calcul = 2;};}
// forcer un stockage pour des sorties
else if (nom == "sortie_post_")
{*(entreePrinc->entree) >> sortie_post;
}
// sinon ce n'est pas un mot clé connu, on le signale
else if (nom != "fin_parametres_reglage_")
{ cout << "\n erreur en lecture d'un parametre, le mot cle est inconnu "
<< " on a lu : " << nom << endl;
if (ParaGlob::NiveauImpression()>3)
cout << "\n Projection_anisotrope_3D::LectureDonneesParticulieres(UtilLecture * entreePrinc) " << endl ;
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
}; //-- fin du while
}; //-- fin de la lecture des paramètres de réglage
// --- appel au niveau de la classe mère ---
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD);
{ string mot_cle("fin_loi_projection_anisotrope");
if (entreePrinc->Lecture_et_verif_mot_cle(nom_class_methode,mot_cle))
entreePrinc->NouvelleDonnee(); // on se positionne sur un nouvel enreg
};
// //----- debug ---
// cout << "\n debug Projection_anisotrope_3D::LectureDonneesParticulieres ";
// this->Affiche();
// //--- fin debug ---
};
// affichage de la loi
void Projection_anisotrope_3D::Affiche() const
{ cout << "\n ....... loi de comportement Projection anisotrope 3D ........";
// pour faire une boucle on constitue un tableau de mots clés
Tableau <string > tab_mot_cle(6);
tab_mot_cle(1) = "A1";tab_mot_cle(2) = "A2";tab_mot_cle(3) = "A3";
tab_mot_cle(4) = "B12";tab_mot_cle(5) = "B13";tab_mot_cle(6) = "B23";
// puis un tableau pour les valeurs
Tableau < const double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// on boucle sur les 6 coefficients
for (int i=1;i<=6;i++)
{string mot_cle1=tab_mot_cle(i)+"= ";
string mot_cle2=tab_mot_cle(i)+"_fonction_nD:";
cout << mot_cle1 ;
if (fct_para(i) == NULL)
cout << *coef(i) << " ";
else
{cout << mot_cle2 << " ";
if (fct_para(i)->NomFonction() != "_")
cout << fct_para(i)->NomFonction();
else
fct_para(i)->Affiche();
cout << "\n";
};
};
// le nom du repère associé
cout << "\n nom_repere_associe " << nom_repere;
cout << "\n ratio_inf_module_compressibilite_ "<<ratio_inf_module_compressibilite <<" ";
// indicateur de cas de calcul
if (cas_calcul != 0)
{ if (cas_calcul == 1)
{cout << "\n calcul uniquement deviatorique ";}
else if (cas_calcul == 2)
{cout << " calcul uniquement spherique ";}
else
{cout << " cas de calcul mal defini !! ";};
};
// affichage du type de transport
cout << " type_transport: " << type_transport;
// niveau d'affichage
cout << " niveau_affichage_local: ";
Affiche_niveau_affichage();
cout << " sortie_post: "<< sortie_post;
cout << endl;
// appel de la classe mère
Loi_comp_abstraite::Affiche_don_classe_abstraite();
};
// affichage et definition interactive des commandes particulières à chaques lois
void Projection_anisotrope_3D::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");
// pour faire une boucle on constitue un tableau de mots clés
Tableau <string > tab_mot_cle(6);
tab_mot_cle(1) = "A1";tab_mot_cle(2) = "A2";tab_mot_cle(3) = "A3";
tab_mot_cle(4) = "B12";tab_mot_cle(5) = "B13";tab_mot_cle(6) = "B23";
// puis un tableau pour les valeurs
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
A1= 1.; A2= 2.; A3= 3.;
B12= 1.; B13 = 1.; B23 = 1.;
sort << "\n# ....... loi de comportement Projection anisotrope 3D ........"
<< "\n# on commence par indiquer le type de projection, par exemple: "
<< "\n# TYPE_PROJ_ANISO_ PROJ_ORTHO "
<< "\n# On definit ensuite les informations necessaires pour calculer"
<< "\n# la projection. Par exemple pour PROJ_ORTHO "
<< "\n# on definit les 6 coefficients diagonaux du tenseur H exprimee "
<< "\n# dans le repere d'anisotropie "
<< "\n# A1= 1. A2= 2. A3= 3. B12= 1. B13= 1. B23= 1. "
<< "\n# ensuite on definit le repere d'anisotropie ex: "
<< "\n# nom_repere_associe_ repere1 "
<< "\n# puis le nom de la loi 3D interne associee"
<< "\n# ISOELAS "
<< "\n# 200000 0.3 "
<< "\n# fin_loi_interne # --- fin def loi interne"
<< "\n# fin_loi_projection_anisotrope # mot cle de fin de declaration "
<< "\n# "
;
// cas avec plus d'information
if ((rep != "o") && (rep != "O" ) && (rep != "0") )
{ // cas d'une loi thermo dépendante
sort << "\n# .... infos complementaires ...."
<<"\n# 1) Pour chaque parametre materiau, individuellement, il est possible "
<<"\n# d'utiliser une fonction nD a la place d'une valeur numerique. "
<<"\n# La fonction nD peut alors dependre de toutes les grandeurs disponibles "
<<"\n# localement, en particulier la position, les donnees imposees par l'utilisateur "
<<"\n# par exe: une temperature, ou tout ddl obtenu directement ou par interpolation "
<<"\n# suivant sa disponibilite, sachant que l'on regarde d'abord si la grandeur est "
<<"\n# directement disponible au point d'integration, sinon on regarde si elle est "
<<"\n# disponible par interpolation. "
<<"\n# Par exemple pour PROJ_ORTHO: "
<<"\n# supposons que l'on veuille que A3 soit une fonction nD on ecrira: "
<<"\n# soit a: "
<<"\n# A3= A3_fonction_nD: un_nom_de_fonction_existant "
<<"\n# soit b: "
<<"\n# A3= A3_fonction_nD: un_nom_de_type_de_fonction_existant "
<<"\n# suivit sur la ligne suivante de la definition specifique de la fonction "
<<"\n# "
<<"\n# exemple d'un cas a: "
<<"\n# A3= A3_fonction_nD: f1_temperature "
<<"\n# fi_temperature doit alors avoir ete definie dans les fonctions nD"
<<"\n# "
<<"\n# exemple d'un cas b: "
<<"\n# A3= A3_fonction_nD: FONCTION_EXPRESSION_LITTERALE_nD "
<<"\n# deb_list_var_ TEMP fin_list_var_ "
<<"\n# fct= (100726-101325)/50*TEMP + 101325 "
<<"\n# fin_parametres_fonction_expression_litterale_ "
<<"\n# "
<<"\n# Remarques: "
<<"\n# a) apres chaque definition d'une fonction nD on change de ligne "
<<"\n# b) pour les autres coefficients, on remplace A3 par A1, ou A2, ou"
<<"\n# B12 ou B13 ou B23 "
<<"\n# "
<<"\n# ====================================================================="
<<"\n# Apres avoir defini les parametres de la loi, il est possible de preciser et/ou modifier"
<<"\n# le fonctionnement de la loi via des parametres de reglage qui sont encapsules entre les 2 mots"
<<"\n# clef: avec_parametres_de_reglage_ et fin_parametres_reglage_"
<<"\n# le premier se situe sur la fin de la ligne definissant les para de la loi, le second doit"
<<"\n# se situer sur une ligne seule"
<<"\n# L'ordre d'apparition des parametres n'a pas d'importance, ni le fait qu'ils soient sur une"
<<"\n# ou plusieurs lignes. la liste des parametres de reglage disponible et leurs significations sont"
<<"\n# les suivantes: "
<< "\n# 2)-------------- type de transport de la base d'orthotropie initiale ---------- "
<<"\n# Par defaut, la base initiale d'orthotropie est transportee dans l'etat courant "
<<"\n# via une methode de type transport contravariant (cf. doc) "
<<"\n# Il est possible d'indiquer un transport de type covariant. Pour ce faire "
<<"\n# on utilise le mot cle: type_transport_ suivi de 0 ou 1 "
<<"\n# 0 : pour un transport de type contravariant (valeur par defaut) "
<<"\n# 1 : pour un transport de type covariant "
<<"\n# d'une maniere generale, les parametres de reglage relatifs"
<<"\n# au repere d'anisotropie et a son evolution, sont indiques entre"
<<"\n# le mot cle: avec_parametres_de_reglage_ qui doit se trouver "
<<"\n# sur une ligne seule, apres la definition de la loi interne, "
<<"\n# et le mot cle: fin_parametres_reglage_ qui doit etre sur une seule ligne"
<<"\n# "
<<"\n# exemple de declaration: "
<< "\n# fin_loi_interne # --- fin def loi interne"
<< "\n# avec_parametres_de_reglage_ "
<< "\n# sortie_post_ 1 "
<< "\n# type_transport_ 0 "
<< "\n# fin_parametres_reglage_ "
<<"\n# "
<<"\n# 3)"
<< "\n# -------------- affichage des erreurs et des warning ---------- "
<< "\n# - l'affichage normale est fonction du parametre global d'affichage gerer"
<< "\n# par le niveau d'affichage. Cependant pour des raisons par exemple de mise au point,"
<< "\n# il est possible de permettre l'affichage local a un niveau particulier "
<< "\n# (mot cle : permet_affichage_ suivi d'un nombre entier) en plus de l'affichage normal. "
<< "\n# l'affichage s'effectuera donc en fonction de l'affichage normale et de l'affichage particulier."
<< "\n# Le fonctionnement de l'affichage particulier suit les mêmes règles que l'affichage globale"
<< "\n# soit permet_affichage_ est nulle (cas par defaut), dans ce cas l'affichage est fonction du niveau global"
<< "\n# soit permet_affichage_ vaut n , (!= 0), dans ce cas l'affichage est fonction uniquement de n "
<< "\n# "
<< "\n# ex: permet_affichage_ 5 "
<< "\n# ce mot cle ainsi que le type de deformation (cf. doc) doit etre le dernier "
<< "\n# enregistrement de la loi avant le mot cle fin_loi_projection_anisotrope "
<< "\n# Remarque: la loi interne peut egalement avoir un permet_affichage_ independant "
<< "\n# il doit etre indique juste avant le mot cle: fin_loi_interne "
<< "\n# exemple de declaration avec deux niveaux differents: "
<< "\n# "
<< "\n# ISOELAS "
<< "\n# 200000 0.3 "
<< "\n# permet_affichage_ 5 # = le niveau utilise pour la loi elastique"
<< "\n# fin_loi_interne # --- fin def loi interne"
<< "\n# permet_affichage_ 3 # = le niveau utilise pour la loi de projection"
<< "\n# fin_loi_projection_anisotrope # mot cle de fin de declaration "
<< "\n# "
<< "\n# Remarque: de maniere equivalente, il est possible de mettre "
<< "\n# le mot permet_affichage_ dans la liste des parametres de reglage de la loi "
<<"\n# cf. exemple plus bas "
<<"\n# "
<<"\n# "
<< "\n# 4) --------------- calcul spherique, deviatorique ou total -------"
<<"\n# Il est possible de calculer que la partie spherique du tenseur de contrainte "
<<"\n# ou que la partie deviatorique "
<<"\n# Pour ce faire on indique a la suite des autres parametres:"
<<"\n# soit le mot clef : seule_deviatorique "
<<"\n# soit le mot clef : seule_spherique "
<<"\n# "
<<"\n# "
<<"\n# "
<<"\n# "
<< "\n# 5) --------------- calcul compressibilite -------"
<<"\n# La compressibilite est calculee a partir de la variation de volume constatee"
<<"\n# et de la trace du tenseur des contraintes. Lorsque la variation est trop faible "
<<"\n# c'est la compressibilite initiale qui est utilisee. Lorsque la compressibilite"
<<"\n# est inferieure a 1/100 de la compressibilite initiale, on considere que le calcul"
<<"\n# de cette compressibilite n'est pas correcte, et on utilise la compressibilite initiale."
<<"\n# Le coefficient 1/100 est arbitraire, on peut le changer avec le mot cle: "
<<"\n# ratio_inf_module_compressibilite_ <un reel> "
<<"\n# ce qui permet de moduler la borne inferieure de la compressibilite"
<<"\n# "
<< "\n# 6) --------------- acces en sortie a des grandeurs intermediaires de calcul -------"
<< "\n# A chaque resolution, il est possible de stocker: "
<< "\n# - le tenseur des contraintes exprime dans le repere d'anisotropie "
<< "\n# - les parametres de la projection: interessant s'ils varient "
<< "\n# par exemple pour: TYPE_PROJ_ANISO_ PROJ_ORTHO "
<< "\n# on aura acces aux 6 parametres A1 A2 A3 B12 B13 B23 dans cet ordre"
<< "\n# "
<< "\n# le mot cle est sortie_post_ , par defaut il vaut 0, dans ce cas aucun indicateur n'est stoke"
<< "\n# s'il est different de 0, on peut acceder aux grandeurs "
<< "\n# seules les grandeurs en cours sont disponibles, il n'y a pas de stockage sur plusieurs increment "
<< "\n# "
<<"\n# exemple de declaration: "
<< "\n# fin_loi_interne # --- fin def loi interne"
<< "\n# avec_parametres_de_reglage_ "
<< "\n# sortie_post_ 1 "
<< "\n# type_transport_ 0 "
<< "\n# permet_affichage_ 3 "
<< "\n# fin_parametres_reglage_ "
<< "\n# "
<<"\n# "
;
};
// appel de la classe mère
Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);
};
// test si la loi est complete
int Projection_anisotrope_3D::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
// pour faire une boucle on constitue un tableau de mots clés
Tableau <string > tab_mot_cle(6);
tab_mot_cle(1) = "A1";tab_mot_cle(2) = "A2";tab_mot_cle(3) = "A3";
tab_mot_cle(4) = "B12";tab_mot_cle(5) = "B13";tab_mot_cle(6) = "B23";
// puis un tableau pour les valeurs
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
for (int i=1;i<=6;i++)
if ((*coef(i) == -ConstMath::trespetit) && (fct_para(i) == NULL))
{string mot_cle1=tab_mot_cle(i)+"= ";
cout << mot_cle1 << " n'est pas defini \n " ;
ret = 0;
};
// test du cas de calcul
if ((cas_calcul < 0) || (cas_calcul > 2))
{ cout << "\n l'indicateur de calcul cas_calcul= " << cas_calcul << " n'est pas correcte "
<< "\n ceci pour la Projection_anisotrope_3D";
ret = 0;
};
// le nom du repère associé
if (nom_repere.length() == 0)
{ cout << "\n le nom du repere assoce n'est pas defini "
<< "\n ceci pour la Projection_anisotrope_3D";
ret = 0;
};
// info globale disponible actuellement
if (!ret)
{ cout << "\n loi Projection_anisotrope_3D incomplete : ";
Affiche();
};
// retour
return ret;
};
// activation des données des noeuds et/ou elements nécessaires au fonctionnement de la loi
// exemple: mise en service des ddl de température aux noeuds
// ici la grandeur qui sert de proportion entre la première loi et la seconde
void Projection_anisotrope_3D::Activation_donnees(Tableau<Noeud *>& tabnoeud,bool dilatation,LesPtIntegMecaInterne& lesPtMecaInt)
{ // activation pour la loi projetée
loi_interne->Activ_donnees(tabnoeud,dilatation,lesPtMecaInt);
// appel de la méthode de la classe mère
Loi_comp_abstraite::Activ_donnees(tabnoeud,dilatation,lesPtMecaInt);
};
// 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 Projection_anisotrope_3D::Grandeur_particuliere
(bool ,List_io<TypeQuelconque>& liTQ,Loi_comp_abstraite::SaveResul * saveDon,list<int>& decal) const
{ // ici on est en 3D et les grandeurs sont par principe en locale, donc la variable absolue ne sert pas
// on passe en revue la liste
List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end();
list<int>::iterator idecal=decal.begin();
for (itq=liTQ.begin();itq!=itqfin;itq++,idecal++)
{TypeQuelconque& tipParticu = (*itq); // pour simplifier
if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
switch (tipParticu.EnuTypeQuelconque().EnumTQ())
{ case REPERE_D_ANISOTROPIE:
// a) ----- cas du repère d'anisotropie
{ SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveDon);
Tab_Grandeur_BaseH& tyTQ= *((Tab_Grandeur_BaseH*) (*itq).Grandeur_pointee()); // pour simplifier
tyTQ(1+(*idecal))=save_resul.Op_H;
(*idecal)++; break;
}
case EPS_TRANSPORTEE_ANISO:
// ----- cas de la déformation transportée dans le repère d'anisotropie
{ SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveDon);
Tab_Grandeur_TenseurHH& tyTQ= *((Tab_Grandeur_TenseurHH*) (*itq).Grandeur_pointee()); // pour simplifier
tyTQ(1+(*idecal))= *(save_resul.eps_loc_HH);
(*idecal)++; break;
}
case SIGMA_DANS_ANISO:
// ----- cas de la contrainte calculée dans le repère d'anisotropie
{ SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveDon);
Tab_Grandeur_TenseurHH& tyTQ= *((Tab_Grandeur_TenseurHH*) (*itq).Grandeur_pointee()); // pour simplifier
tyTQ(1+(*idecal))= *(save_resul.sig_loc_HH);
(*idecal)++; break;
}
case PARA_ORTHO:
// ----- cas des paramètres de la projection
{ SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveDon);
Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
tyTQ(1+(*idecal))= *(save_resul.para_loi);
(*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 Projection_anisotrope_3D::ListeGrandeurs_particulieres(bool absolue,List_io<TypeQuelconque>& liTQ) const
{
// $$$ cas du repère local d'anisotropie
int dim_espace = 3;
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == REPERE_D_ANISOTROPIE)
{ Tab_Grandeur_BaseH& tyTQ= *((Tab_Grandeur_BaseH*) (*itq).Grandeur_pointee()); // pour simplifier
int taille = tyTQ.Taille()+1;
tyTQ.Change_taille(taille); nexistePas = false;
};
if (nexistePas)
{Grandeur_BaseH v_rep(dim_espace,3);
Tab_Grandeur_BaseH grand5(v_rep,1); // def d'une grandeur courante
TypeQuelconque typQ6(REPERE_D_ANISOTROPIE,EPS11,grand5);
liTQ.push_back(typQ6);
};
};
// ---- la suite dépend de l'indicateur : sortie_post
if (sortie_post)
{ // $$$ cas de la déformation transportée dans le repère d'anisotropie
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == EPS_TRANSPORTEE_ANISO)
{Tab_Grandeur_TenseurHH& tyTQ= *((Tab_Grandeur_TenseurHH*) (*itq).Grandeur_pointee()); // pour simplifier
int taille = tyTQ.Taille()+1;
tyTQ.Change_taille(taille); nexistePas = false;
};
if (nexistePas)
{TenseurHH* tens = NevezTenseurHH(3); // un tenseur typique
Tab_Grandeur_TenseurHH eps_loc_HH(*tens,1);
// def d'un type quelconque représentatif
TypeQuelconque typQ(EPS_TRANSPORTEE_ANISO,EPS11,eps_loc_HH);
liTQ.push_back(typQ);
delete tens; // car on n'en a plus besoin
};
};
// $$$ cas de la contrainte calculée dans le repère d'anisotropie
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == SIGMA_DANS_ANISO)
{Tab_Grandeur_TenseurHH& tyTQ= *((Tab_Grandeur_TenseurHH*) (*itq).Grandeur_pointee()); // pour simplifier
int taille = tyTQ.Taille()+1;
tyTQ.Change_taille(taille); nexistePas = false;
};
if (nexistePas)
{TenseurHH* tens = NevezTenseurHH(3); // un tenseur typique
Tab_Grandeur_TenseurHH eps_loc_HH(*tens,1);
// def d'un type quelconque représentatif
TypeQuelconque typQ(SIGMA_DANS_ANISO,SIG11,eps_loc_HH);
liTQ.push_back(typQ);
delete tens; // car on n'en a plus besoin
};
};
// $$$ cas des paramètres de la loi de comportement
{List_io<TypeQuelconque>::iterator itq,itqfin=liTQ.end(); bool nexistePas = true;
for (itq=liTQ.begin();itq!=itqfin;itq++)
if ((*itq).EnuTypeQuelconque() == PARA_ORTHO)
{Tab_Grandeur_Vecteur& tyTQ= *((Tab_Grandeur_Vecteur*) (*itq).Grandeur_pointee()); // pour simplifier
int taille = tyTQ.Taille()+1;
tyTQ.Change_taille(taille); nexistePas = false;
};
if (nexistePas)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{Vecteur tens(6); // les 6 paramètres
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
Tab_Grandeur_Vecteur gr_vec_para(tens,1);
// def d'un type quelconque représentatif
TypeQuelconque typQ(PARA_ORTHO,EPS11,gr_vec_para);
liTQ.push_back(typQ);
};
};
}; // fin du cas ou sortie_post est actif, c-a-d que l'on veut des infos sur les indicateurs
// de résolution
};
// indique le type Enum_comp_3D_CP_DP_1D correspondant à une loi de comportement
// la fonction est simple dans le cas d'une loi basique, par contre dans le cas
// d'une loi combinée, la réponse dépend des lois internes donc c'est redéfini
// dans les classes dérivées
Enum_comp_3D_CP_DP_1D Projection_anisotrope_3D::Comportement_3D_CP_DP_1D()
{ // si la loi interne n'est pas nulle, on ramène le cas de la loi
if (loi_interne != NULL)
{return (loi_interne->Comportement_3D_CP_DP_1D());}
else // sinon on ramène le cas non défini
{return RIEN_COMP_3D_CP_DP_1D;};
};
//----- 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 Projection_anisotrope_3D::Lecture_base_info_loi
(ifstream& ent,const int cas,LesReferences& lesRef
,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string toto,nom;
if (cas == 1)
{ ent >> nom ;
if (nom != "PROJECTION_ANISOTROPE_3D")
{cout <<"\n *** erreur en lecture du type de la loi de comportement, on attendait"
<< " la chaine de caracteres: PROJECTION_ANISOTROPE_3D et on a lu "
<< nom <<" !!! on ne peut pas continuer "
<< "\n Projection_anisotrope_3D::Lecture_base_info_loi(.. "
<< flush;
Sortie(1);
};
// lecture du type de projection
ent >> toto >> nom;
type_projection = Id_Nom_proj_aniso(nom);
bool existe_fctnD=false; // pour la gestion des grandeurs locales éventuelles
// pour faire une boucle on constitue un tableau de mots clés
Tableau <string > tab_mot_cle(6);
tab_mot_cle(1) = "A1";tab_mot_cle(2) = "A2";tab_mot_cle(3) = "A3";
tab_mot_cle(4) = "B12";tab_mot_cle(5) = "B13";tab_mot_cle(6) = "B23";
// puis un tableau pour les valeurs
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// cela va permettre de choisir entre une valeur fixe et une fonction nD
for (int i=1;i<7;i++)
{ ent >> toto >> nom;
string mot_cle2=tab_mot_cle(i)+"_fonction_nD:";
if (nom == mot_cle2)
// cas d'une fonction nD
{fct_para(i) = lesFonctionsnD.Lecture_pour_base_info(ent,cas,fct_para(i));
existe_fctnD = true;
}
else // cas d'une valeur fixe
{(*coef(i)) = ChangeReel(nom);};
};
// lecture du repère associé
ent >> nom >> nom_repere;
// indicateur pour les calculs partielles
string type_de_calcul;
ent >> nom >> type_de_calcul ;
if (type_de_calcul == "seul_deviatorique")
{cas_calcul=1;}
else if (type_de_calcul == "seul_spherique")
{cas_calcul=2;}
else
{cas_calcul=0;}
// ratio sur la compressibilité
ent >> nom >> ratio_inf_module_compressibilite;
// le type de transport
ent >> nom >> type_transport;
// le niveau d'affichage
Lecture_permet_affichage(ent,cas,lesFonctionsnD);
// sortie_post
ent >> nom >> sortie_post;
// pondération éventuelle
int avec_nD=0;
ent >> toto >> nom >> avec_nD;
if (nom != "fct_nD")
{cout <<"\n *** erreur en lecture du type de ponderation, on attendait"
<< " la chaine de caracteres: fct_nD et on a lu "
<< nom <<" !!! on ne peut pas continuer "
<< "\n Projection_anisotrope_3D::Lecture_base_info_loi(.. "
<< flush;
Sortie(1);
};
if (avec_nD == 1 )
{existe_fctnD=true;
if (ponderation_nD_quelconque != NULL)
{ponderation_nD_quelconque->Lecture_base_info(ent, cas, lesFonctionsnD);
}
else
{ ponderation_nD_quelconque = new Ponderation_TypeQuelconque();
ponderation_nD_quelconque->Lecture_base_info(ent, cas, lesFonctionsnD);
};
}
else
{if (ponderation_nD_quelconque != NULL)
{delete ponderation_nD_quelconque;
ponderation_nD_quelconque=NULL;
};
};
// la loi de comportement
string nom_completude_calcul;
ent >> toto >> nom >> nom_completude_calcul ; // lecture du nom de la loi et du type d'action
if (nom_completude_calcul == "CONTRAINTE_ET_TANGENT")
{ completude_calcul = CONTRAINTE_ET_TANGENT;}
else if (nom_completude_calcul == "CONTRAINTE_UNIQUEMENT")
{ completude_calcul=CONTRAINTE_UNIQUEMENT;}
else if (nom_completude_calcul == "TANGENT_UNIQUEMENT")
{ completude_calcul = TANGENT_UNIQUEMENT;}
else
{ cout << "\n *** erreur en lecture du type d'action a faire avec la loi, on a lue " << nom_completude_calcul
<< " au lieu d'un des noms suivants CONTRAINTE_ET_TANGENT CONTRAINTE_UNIQUEMENT TANGENT_UNIQUEMENT "
<< "\n Projection_anisotrope_3D::Lecture_base_info_loi(.. "
<< flush;
Sortie(1);
};
// definition de la loi
if (loi_interne != NULL)
delete loi_interne;
LoiAbstraiteGeneral * pt = LesLoisDeComp::Def_loi(nom);
// lecture des informations propres à la loi
pt->Lecture_base_info_loi(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
// récup de la loi
loi_interne = (Loi_comp_abstraite*) pt;
// s'il y a pondération avec des grandeurs locales, on va vérifier que ces grandeurs
// sont dispo et on va préparer l'accès à ces grandeurs
if (existe_fctnD)
Verif_et_preparation_acces_grandeurs_locale();
}
else
{// lecture des informations propres à la loi
loi_interne->Lecture_base_info_loi(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
};
// 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 Projection_anisotrope_3D::Ecriture_base_info_loi(ofstream& sort,const int cas)
{ if (cas == 1)
{ sort << " PROJECTION_ANISOTROPE_3D " ;
// on indique le type de projection
sort << "\n TYPE_PROJ_ANISO_ "<< Nom_proj_aniso(type_projection) << "\n" ;
// pour faire une boucle on constitue un tableau de mots clés
Tableau <string > tab_mot_cle(6);
tab_mot_cle(1) = "A1";tab_mot_cle(2) = "A2";tab_mot_cle(3) = "A3";
tab_mot_cle(4) = "B12";tab_mot_cle(5) = "B13";tab_mot_cle(6) = "B23";
// puis un tableau pour les valeurs
Tableau < const double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// on boucle sur les 6 coefficients
for (int i=1;i<7;i++)
{string mot_cle1=tab_mot_cle(i)+"= ";
string mot_cle2=tab_mot_cle(i)+"_fonction_nD:";
sort << mot_cle1 ;
if (fct_para(i) == NULL)
sort << *coef(i) << " ";
else
{sort << mot_cle2 << " ";
LesFonctions_nD::Ecriture_pour_base_info(sort, cas,fct_para(i));
};
};
// le nom du repère associé
sort << "\n nom_repere_associe_ "<< nom_repere;
// indicateur de cas de calcul
if (cas_calcul != 0)
{ if (cas_calcul == 1)
{sort << "\n seul_deviatorique ";}
else if (cas_calcul == 2)
{sort << " seul_spherique ";}
else
{sort << " cas_de_calcul_mal_defini ";};
};
// ratio sur la compressibilité
sort << "\n ratio_inf_module_compressibilite_ "<< ratio_inf_module_compressibilite;
// affichage du type de transport
sort << " type_transport: " << type_transport;
// niveau d'affichage
Affiche_niveau_affichage(sort,cas);
sort << " sortie_post: "<<sortie_post;
// pondération éventuelle
if (ponderation_nD_quelconque!=NULL)
{bool sans_courbe = false;
sort << " ponderation: fct_nD 1 ";
ponderation_nD_quelconque->Ecriture_base_info(sort,cas,sans_courbe);
}
else {sort << " ponderation: fct_nD 0 ";};
// la loi de comportement
sort << "\n loi_interne: " << loi_interne->Nom_comport() << " ";
switch (completude_calcul)
{case CONTRAINTE_ET_TANGENT: sort << " CONTRAINTE_ET_TANGENT "; break;
case CONTRAINTE_UNIQUEMENT: sort << " CONTRAINTE_UNIQUEMENT "; break;
case TANGENT_UNIQUEMENT: sort << " TANGENT_UNIQUEMENT "; break;
};
loi_interne->Ecriture_base_info_loi(sort,cas);
sort << endl;
}
else
{// on demande à la loi
loi_interne->Ecriture_base_info_loi(sort,cas);
};
// appel de la classe mère
Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);
};
// calcul d'un module d'young équivalent à la loi pour un chargement nul
double Projection_anisotrope_3D::Module_young_equivalent(Enum_dure temps,const Deformation & def,SaveResul * )
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ cout << "\n *** methode non encore implante "
<< " Projection_anisotrope_3D::Module_young_equivalent(Enum_dure temps,const Deformation & def,SaveResul * ) "
<< endl;
Sortie(1);
return 0.; // taire le compilo
/*if (!thermo_dependant)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ return E;}
else
{ temperature_tdt = def.DonneeInterpoleeScalaire(TEMP,temps);
return E_temperature->Valeur(temperature_tdt);
};
*/ };
// récupération d'un module de compressibilité équivalent à la loi pour un chargement nul
double Projection_anisotrope_3D::Module_compressibilite_equivalent(Enum_dure temps,const Deformation & def,SaveResul * av)
{ if (!Bll_fct_para)
{cout << "\n **** petit probleme: dans le cas de fonction nd decrivant les parametres de projection "
<< " pour l'instant le calcul du module de compressibilite equivalent n'est pas possible "
<< "\n Projection_anisotrope_3D::Module_compressibilite_equivalent(... " << flush;
Sortie(1);
mise à jour 7.010
2023-05-03 17:23:49 +02:00
return 0.;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
}
else
{ // récup du conteneur spécifique
SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) av);
double compress_inter = loi_interne->Module_compressibilite_equivalent(temps, def, save_resul.SaveResul_interne);
double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
double module_compressibilite = 1./3. * coeff * compress_inter;
return module_compressibilite;
};
};
// ========== codage des METHODES VIRTUELLES protegees:================
void Projection_anisotrope_3D::Calcul_SigmaHH (TenseurHH& sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl
,TenseurBB & gijBB_t,TenseurHH & gijHH_t,BaseB& giB,BaseH& gi_H, TenseurBB & epsBB_
,TenseurBB& delta_epsBB ,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 << " Projection_anisotrope_3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
mise à jour 7.010
2023-05-03 17:23:49 +02:00
bool affichage = (Permet_affichage() > 3);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n --- loi Projection_anisotrope_3D Calcul_SigmaHH --- ";
Signature_pti_encours(cout);
};
#endif
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_); // " " " "
const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_); // " " " "
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // " " " "
Tenseur3BH epsBH = epsBB * gijHH; // deformation en mixte
// récup du conteneur spécifique
SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveResul);
// --- on commence par calculer le repère d'anisotropie transporté
// Op_H_i sont les coordonnées contravariantes de la nouvelle base
// donc par rapport à g_i,
BaseH& Op_H = save_resul.Op_H;
Tableau <double> tab_norme(3);
if (type_transport == 0)
// transport de type contravariant
{ // on calcule les coordonnées de la base O' dans la base naturelle
#ifdef MISE_AU_POINT
if (save_resul.O_H == NULL)
{cout << "\n *** erreur, le repere d'anisotropie n'est pas defini "
<<"on ne peut pas continuer " << flush ;
cout << " Projection_anisotrope_3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
BaseH& O_H = (*save_resul.O_H); // coordonnées en absolu de la base d'anisotropie
// calcul des coordonnées alpha_a^{.i} = O_H(a) * g^i à t= 0
for (int a = 1;a < 4; a++)
{CoordonneeH& inter = alpha_H.CoordoH(a);
for (int i=1;i < 4;i++)
inter(i)= O_H(a).ScalHH((*ex.giH_0)(i));
};
// calcule des beta_a^{.i} tel que O'_a = beta_a^{.i} \hat{g}_i donc au temps t+dt
for (int a=1; a<4;a++)
{ // tout d'abord la base non normalisée
// Op_H(a) non normé a les mêmes coordonnées alpha_a^{.j}
// mais exprimés dans le repère actuel \hat \vec g_j
CoordonneeH& Op_H_a = Op_H.CoordoH(a);
Op_H_a = alpha_H(a);
// calcul de la norme du vecteur Op_H_a
double norme = sqrt(Op_H_a * gijBB * Op_H_a);
tab_norme(a) = norme;
// coordonnées finales
Op_H_a /= norme;
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n Op_H("<<a<<"): ";
Op_H_a.Affiche();
cout << "\n alpha_H("<<a<<"): ";
alpha_H(a).Affiche();
cout << "\n giB_tdt"; (*ex.giB_tdt)(a).Affiche();
};
#endif
};
}
else
// transport de type covariant
{ // on calcule les coordonnées de la base O' dans la base duale
#ifdef MISE_AU_POINT
if (save_resul.O_B == NULL)
{cout << "\n *** erreur, le repere d'anisotropie n'est pas defini "
<<"on ne peut pas continuer " << flush ;
cout << " Projection_anisotrope_3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// BaseB& O_B = (*save_resul.O_B); // pour simplifier
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
for (int a=1; a<4;a++)
{ // tout d'abord la base non normalisée
CoordonneeB& Op_B_a = Op_B.CoordoB(a);
// calcul de la norme du vecteur Op_B_a
double norme = sqrt(Op_B_a * gijHH * Op_B_a);
// coordonnées finales
Op_B_a /= norme;
};
// maintenant on calcule la base Op_H correspondante
for (int i=1; i<4;i++)
Op_H.CoordoH(i) = Op_B.CoordoB(i) * gijHH;
};
// on calcul la matrice de passage de la base g_i vers la base O'_i
// O'_i = beta_i^{.j} * g_j
// et on a également \hat{\vec g}^i = {\beta}_{a}^{.i}~\hat{\vec O'}^a
// comme O'_i est déjà exprimé dans g_j, ses coordonnées sont directement béta
for (int i=1; i<4;i++)
for (int j=1; j<4;j++)
{ beta(i,j) = Op_H(i)(j);};
// puis on calcul les coordonnées de la base duale
beta_transpose = beta.Transpose();
gamma = beta_transpose.Inverse();
gamma_transpose = gamma.Transpose();
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n beta: ";beta.Affiche();
cout << "\n beta_transpose: ";beta_transpose.Affiche();
cout << "\n gamma: ";gamma.Affiche();
cout << "\n gamma_transpose: ";gamma_transpose.Affiche();
};
#endif
// changement de base (cf. théorie) : la matrice beta est telle que:
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// calcul des coordonnées de la déformation dans le repère O'_i
#ifdef MISE_AU_POINT
if (affichage)
{ Tenseur3BB eps_p_BB(epsBB);
cout <<"\n eps_p_BB: "; eps_p_BB.Ecriture(cout);
};
#endif
double untiers = 1./3.;
// dans le cas de fonctions nD récupération des valeurs
if (!Bll_fct_para)
// Tableau <Fonction_nD* > fct_para; // fonction nD éventuelle d'évolution des paramètres
{// un tableau de travail pour les valeurs sous forme indicée
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// 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;
// on passe en revue les fonctions nD
for (int i=1;i<7;i++)
{if (fct_para(i) != NULL)
{ Fonction_nD* pt_fonct = fct_para(i); // pour simplifier
// on utilise la méthode générique de loi abstraite
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(pt_fonct,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,NULL
);
/*
// ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
);
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 1)
{ cout << "\nErreur : la fonction nD relative au parametre materiau " << i
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " Projection_anisotrope_3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
*/
// on récupère le premier élément du tableau uniquement
(*coef(i)) = tab_val(1);
};
};
};
// --- on calcule le tenseur de contrainte de référence
Tenseur3HH sig_ref_HH; // def et init
double compress_inter=0.; double cisaill_inter=0; // init
{ // passage des informations spécifique à la loi liste_des_SaveResul
loi_interne->IndiqueSaveResult(save_resul.SaveResul_interne);
loi_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
loi_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
// calcul de la contrainte résultante,
if (completude_calcul != TANGENT_UNIQUEMENT)
{ loi_interne->Calcul_SigmaHH(*(save_resul.l_sigoHH_t),DepsBB,tab_ddl,gijBB_t,gijHH_t,giB,gi_H
,epsBB_,delta_epsBB,gijBB_,gijHH_,d_gijBB_,jacobien_0,jacobien
,*(save_resul.l_sigoHH),save_resul.l_energ,save_resul.l_energ_t
,compress_inter,cisaill_inter,ex);
if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{
double coef = save_resul.f_ponder;
compress_inter *= coef;
cisaill_inter *= coef;
energ = save_resul.l_energ * coef; // update des énergies totales
sig_ref_HH = coef*(*(save_resul.l_sigoHH)); // on coefficiente la contrainte
}
else
{energ = save_resul.l_energ; // update des énergies totales
sig_ref_HH = (*(save_resul.l_sigoHH)); // récup de la contrainte
};
};
};
// --- on calcule la projection du tenseur de contrainte dans le repère de travail
{ sigHH.Inita(0.);
for (int i=1;i<4;i++) for (int j=i;j<4;j++)
{for (int k=1;k<4;k++) for (int l=1;l<4;l++) for (int a=1;a<4;a++)for (int b=1;b<4;b++)
sigHH.Coor(i,j) += (*hij(a,b)) * beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)
* sig_ref_HH(k,l);
};
};
// dans le cas où on veut une sortie des grandeurs on sauvegarde
if (sortie_post)
{ if (save_resul.eps_loc_HH == NULL) {save_resul.eps_loc_HH = NevezTenseurHH(3);};
if (save_resul.sig_loc_HH == NULL) {save_resul.sig_loc_HH = NevezTenseurHH(3);};
if (save_resul.para_loi == NULL ) {save_resul.para_loi = new Vecteur (6);};
// calcul des coordonnées de la déformation dans le repère O'^i
// c-a-d les coordonnées dans le dual de O'_i
Tenseur3BB eps_p_BB(epsBB);
//cout <<"\n eps_p_BB: "; eps_p_BB.Ecriture(cout);
// il faut passer en 2 fois contravariant
(*save_resul.sig_loc_HH) = gijHH * eps_p_BB * gijHH;
save_resul.sig_loc_HH->ChBase(gamma); // res = (gamma * res) * gamma.Transpose();
// idem pour les contraintes
(*save_resul.sig_loc_HH) = sigHH;
save_resul.sig_loc_HH->ChBase(gamma);
// paramètres de la loi
(*save_resul.para_loi)(1) = A1; (*save_resul.para_loi)(2) = A2;
(*save_resul.para_loi)(3) = A3;
(*save_resul.para_loi)(4) = B12; (*save_resul.para_loi)(5) = B13;
(*save_resul.para_loi)(6) = B23;
};
// passage dans la bonne variance
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Tenseur3BH sigBH = gijBB * sigHH ;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
double trace_sig = sigBH.Trace();
switch (cas_calcul)
{ case 0: // calcul normal (tous les termes)
{ // on ne fait rien de spécial
break;
}
case 1: // calcul de la partie déviatorique seule
{ sigBH -= (untiers * trace_sig) * IdBH3;
sigHH -= (untiers * trace_sig) * gijHH;
break;
}
case 2: // calcul de la partie sphérique seule
{ sigBH = (untiers * trace_sig) * IdBH3;
sigHH = (untiers * trace_sig) * gijHH;
break;
}
default:
{ cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! "
<< "\n Projection_anisotrope_3D::Calcul_SigmaHH (.... ";
Sortie(1);
}
};
//----------------- debug
// const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_tdt); // passage en dim 3
//cout << "\n epsBB="<<epsBB << " gijBB_tdt= " << gijBB << " gijHH_tdt= " << gijHH<<" ";
//cout << "\n sighh="<<sigHH<<" "<< endl ;
//----------- fin debug
// traitement des énergies
{ // on commence par calculer l'incrément de l'énergie de référence
EnergieMeca delta_ener = save_resul.l_energ - save_resul.l_energ_t;
// on calcule le rapport des incréments linéarisés des énergies globales
double puiss_deltat = sigHH && delta_epsBB;
double puiss_ref_deltat = (*(save_resul.l_sigoHH)) && delta_epsBB;
double r = puiss_deltat / (puiss_ref_deltat + ConstMath::trespetit);
delta_ener *= r ; // l'incrément dilaté
energ = energ_t + delta_ener; // l'énergie finale
};
// -- calcul des modules
{
// on n'utilise plus la forme linéaire, mais à la place la variation relative de volume
// constaté, au sens d'une mesure logarithmique, sauf dans le cas où cette variation est trop petite
// calcul de la valeur de la variation relative de volume en log
double log_var_vol = log((*(ex.jacobien_tdt))/(*(ex.jacobien_0)));
// pour le module de compressibilité, choix entre les différents cas
if ((cas_calcul == 0) || (cas_calcul == 2))
{// calcul du module initiale : cas particulier
// celui d'une variation purement volumique (cf. théorie)
//K_{s} = \frac{1}{3} {g"}_{ab}~\left (\lambda(a,b) ~{g"}^{ab}\right ) K_{(ref)s}
// avec {g"}_{ab} et {g"}^{ab} les composantes de la métriques dans O'
// dans le cas ou la def est null
// le repère initial est orthonormé d'où un calcul très simple
double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
double module_compressibilite_initial = untiers * coeff * compress_inter;
double module_courant=module_compressibilite_initial; // init
if (log_var_vol > ConstMath::petit)
{module_courant = untiers * sigBH.Trace() / (log_var_vol);};
// maintenant on va essayer de faire un choix
if (module_courant >= module_compressibilite_initial * ratio_inf_module_compressibilite)
{ module_compressibilite = module_courant;}
else
// cas où le module n'est pas recevable (a priori)
{module_compressibilite = module_compressibilite_initial; };
// if (log_var_vol > ConstMath::petit)
// {module_compressibilite = untiers * sigBH.Trace() / (log_var_vol);}
// else // si la variation de volume est trop faible on passe par un cas particulier
// // celui d'une variation purement volumique (cf. théorie)
// {//K_{s} = \frac{1}{3} {g"}_{ab}~\left (\lambda(a,b) ~{g"}^{ab}\right ) K_{(ref)s}
// // avec {g"}_{ab} et {g"}^{ab} les composantes de la métriques dans O'
// // dans le cas ou la def est null
// // le repère initial est orthonormé d'où un calcul très simple
// double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
// module_compressibilite = untiers * coeff * compress_inter;
// };
}
else
// dans le cas d'une loi purement déviatorique, le module de compressibilité = 0 (approximativement)
{module_compressibilite = 0.;};
// pour la partie cisaillement
if ((cas_calcul == 0) || (cas_calcul == 1))
{// on commence par calculer l'intensité du déviateur des déformations
Tenseur3BH eps_barre_BH(epsBH); // init
double I_eps = epsBH.Trace();
eps_barre_BH -= (untiers * I_eps) * IdBH3;
double II_eps = eps_barre_BH.II();
if (II_eps > ConstMath::petit)
{Tenseur3BH sig_barre_BH(sigBH);
sig_barre_BH -= (untiers * trace_sig) * IdBH3;
module_cisaillement = sig_barre_BH.II()/II_eps;
}
else // sinon c'est le cas à déformation nulle
{module_cisaillement = 0.5*untiers * (B12+B13+B23)
*cisaill_inter;
};
}
else
// en purement sphérique, le module est supposé nul
{module_cisaillement = 0.; };
}
LibereTenseur();
};
// calcul des contraintes a t+dt et de ses variations
void Projection_anisotrope_3D::Calcul_DsigmaHH_tdt
(TenseurHH & sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl
,BaseB& giB_t,TenseurBB & gijBB_t,TenseurHH & gijHH_t
,BaseB& giB_tdt,Tableau <BaseB> & d_giB_tdt,BaseH& giH_tdt,Tableau <BaseH> & d_giH_tdt
,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB
,TenseurBB & delta_epsBB,TenseurBB & gijBB_tdt,TenseurHH & gijHH_tdt
,Tableau <TenseurBB *>& d_gijBB_tdt
,Tableau <TenseurHH *>& d_gijHH_tdt,double& jacobien_0,double& jacobien
,Vecteur& d_jacobien_tdt,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 << " Projection_anisotrope_3D::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 << " Projection_anisotrope_3D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
#endif
mise à jour 7.010
2023-05-03 17:23:49 +02:00
bool affichage = (Permet_affichage() > 3);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n --- loi Projection_anisotrope_3D Calcul_DsigmaHH_tdt --- ";
Signature_pti_encours(cout);
};
#endif
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // " " " "
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " "
const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_tdt); // " " " "
Tenseur3BH epsBH = epsBB * gijHH; // deformation en mixte
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Tenseur3HH epsHH = gijHH * epsBH; // en deuxfois contra
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
int dim_tens = sigHH_tdt.Dimension(); // ici 3, mais on laisse le cas général
TenseurHH* zeroHH = (NevezTenseurHH(dim_tens,0.));
// récup du conteneur spécifique
SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveResul);
// on commence par calculer le repère d'anisotropie transporté
// Op_H_i sont les coordonnées contravariantes de la nouvelle base
// donc par rapport à g_i,
BaseH& Op_H = save_resul.Op_H;
Tableau <double> tab_norme(3);
if (type_transport == 0)
// transport de type contravariant
{ // on calcule les coordonnées de la base O' dans la base naturelle
#ifdef MISE_AU_POINT
if (save_resul.O_H == NULL)
{cout << "\n *** erreur, le repere d'anisotropie n'est pas defini "
<<"on ne peut pas continuer " << flush ;
cout << " Projection_anisotrope_3D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
#endif
BaseH& O_H = (*save_resul.O_H); // coordonnées en absolu de la base d'anisotropie
// calcul des coordonnées alpha_a^{.i} = O_H(a) * g^i à t= 0
for (int a = 1;a < 4; a++)
{CoordonneeH& inter = alpha_H.CoordoH(a);
for (int i=1;i < 4;i++)
inter(i)= O_H(a).ScalHH((*ex.giH_0)(i));
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n O_H("<<a<<"): ";O_H(a).Affiche();
cout << "\n giH_0("<<a<<"): "; (*ex.giH_0)(a).Affiche();
};
#endif
};
// calcule des beta_a^{.i} tel que O'_a = beta_a^{.i} \hat{g}_i donc au temps t+dt
for (int a=1; a<4;a++)
{ // tout d'abord la base non normalisée
// Op_H(a) non normé a les mêmes coordonnées alpha_a^{.j}
// mais exprimés dans le repère actuel \hat \vec g_j
CoordonneeH& Op_H_a = Op_H.CoordoH(a); // pour simplifier
Op_H_a = alpha_H(a);
// calcul de la norme du vecteur Op_H_a
double norme = sqrt(Op_H_a * gijBB * Op_H_a);
tab_norme(a) = norme;
// coordonnées finales
Op_H_a /= norme;
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n Op_H("<<a<<"): ";
Op_H_a.Affiche();
cout << "\n alpha_H("<<a<<"): ";
alpha_H(a).Affiche();
cout << "\n giB_tdt"; (*ex.giB_tdt)(a).Affiche();
cout << "\n giH_tdt("<<a<<"): "; (*ex.giH_tdt)(a).Affiche();
// affichage du repère tournée
CoordonneeB V(3);
for (int j=1;j<4;j++) V += (*ex.giB_tdt)(j) * Op_H_a(j);
cout << "\n op_H("<<a<<"): en absolu " << V ;
};
#endif
};
}
else
// transport de type covariant
{ // on calcule les coordonnées de la base O' dans la base duale
#ifdef MISE_AU_POINT
if (save_resul.O_B == NULL)
{cout << "\n *** erreur, le repere d'anisotropie n'est pas defini "
<<"on ne peut pas continuer " << flush ;
cout << " Projection_anisotrope_3D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
#endif
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// BaseB& O_B = (*save_resul.O_B); // pour simplifier
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
for (int a=1; a<4;a++)
{ // tout d'abord la base non normalisée
CoordonneeB& Op_B_a = Op_B.CoordoB(a);
// calcul de la norme du vecteur Op_B_a
double norme = sqrt(Op_B_a * gijHH * Op_B_a);
tab_norme(a) = norme;
// coordonnées finales
Op_B_a /= norme;
};
// maintenant on calcule la base Op_H correspondante
for (int i=1; i<4;i++)
Op_H.CoordoH(i) = Op_B.CoordoB(i) * gijHH;
};
// on calcul la matrice de passage de la base g_i vers la base O'_i
// O'_i = beta_i^{.j} * g_j
// et on a également \hat{\vec g}^i = {\beta}_{a}^{.i}~\hat{\vec O'}^a
// comme O'_i est déjà exprimé dans g_j, ses coordonnées sont directement béta
for (int i=1; i<4;i++)
for (int j=1; j<4;j++)
{ beta(i,j) = Op_H(i)(j);};
// puis on calcul les coordonnées de la base duale
beta_transpose = beta.Transpose();
gamma = beta_transpose.Inverse();
gamma_transpose = gamma.Transpose();
#ifdef MISE_AU_POINT
if ( affichage)
{cout << "\n beta: ";beta.Affiche();
cout << "\n beta_transpose: ";beta_transpose.Affiche();
cout << "\n gamma: ";gamma.Affiche();
cout << "\n gamma_transpose: ";gamma_transpose.Affiche();
};
#endif
// changement de base (cf. théorie) : la matrice beta est telle que:
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// calcul des coordonnées de la déformation dans le repère O'_i
#ifdef MISE_AU_POINT
if (affichage)
{ Tenseur3BB eps_p_BB(epsBB);
cout <<"\n eps_p_BB: "; eps_p_BB.Ecriture(cout);
Tenseur3BB tiutiu;
epsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
cout << "\n eps en absolu :";tiutiu.Ecriture(cout);
};
#endif
// il faut passer en 2 fois contravariant
Tenseur3HH eps_p_HH(epsHH); // init
#ifdef MISE_AU_POINT
if (affichage)
{cout <<"\n eps_p_HH en g_i: "; eps_p_HH.Ecriture(cout);
};
#endif
eps_p_HH.ChBase(gamma); // res = (gamma * res) * gamma.Transpose();
#ifdef MISE_AU_POINT
if (affichage)
{ cout <<"\n eps_p_HH en O_p: "; eps_p_HH.Ecriture(cout);
cout <<"\n gijHH: "; gijHH.Ecriture(cout);
Tenseur3HH toto(eps_p_HH);
toto.ChBase(beta_transpose);
cout << "\n retour dans la base g^i de eps^{ij} :";toto.Ecriture(cout);
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Tenseur3BB titi = gijBB * toto * gijBB;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n def eps_ij :";titi.Ecriture(cout);
};
#endif
double untiers = 1./3.;
// dans le cas de fonctions nD récupération des valeurs
if (!Bll_fct_para)
// Tableau <Fonction_nD* > fct_para; // fonction nD éventuelle d'évolution des paramètres
{// un tableau de travail pour les valeurs sous forme indicée
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// 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;
// on passe en revue les fonctions nD
for (int i=1;i<7;i++)
{if (fct_para(i) != NULL)
{ Fonction_nD* pt_fonct = fct_para(i); // pour simplifier
// on utilise la méthode générique de loi abstraite
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(pt_fonct,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,NULL
);
/* // ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
);
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 1)
{ cout << "\nErreur : la fonction nD relative au parametre materiau " << i
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " Projection_anisotrope_3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
*/
// on récupère le premier élément du tableau uniquement
(*coef(i)) = tab_val(1);
};
};
};
// --- on calcule le tenseur de contrainte de référence et l'opérateur tangent ----
Tenseur3HH sig_ref_HH; // def et init
double compress_inter=0.; double cisaill_inter=0; // init
int taille = d_sigHH.Taille();
{
// on vérifie que le tenseur de travail intermédiaire est correctement
// dimensionné sinon on le modifie
{ if (d_sigRef_HH.Taille() == 0)
{ d_sigRef_HH.Change_taille(taille,NULL);}
else if ( taille != d_sigRef_HH.Taille())
{ int ancienne_taille = d_sigRef_HH.Taille();
d_sigRef_HH.Change_taille(taille);
if (ancienne_taille < taille)
for (int i=ancienne_taille+1;i<=taille;i++) d_sigRef_HH(i)=NULL;
};
for (int i=1;i<=taille;i++)
{ // mise à zéro des tenseurs du tableau avec création si nécessaire
if (d_sigRef_HH(i) == NULL) {d_sigRef_HH(i)=NevezTenseurHH(dim_tens,0.);}
else if ( d_sigRef_HH(i)->Dimension() != dim_tens)
{ delete d_sigRef_HH(i); d_sigRef_HH(i)=NevezTenseurHH(dim_tens,0.);}
else // mise à zéro simple
*(d_sigRef_HH(i))= (*zeroHH);
};
};
// passage des informations spécifique à la loi liste_des_SaveResul
loi_interne->IndiqueSaveResult(save_resul.SaveResul_interne);
loi_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
loi_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
// calcul de la contrainte résultante,
loi_interne->Calcul_DsigmaHH_tdt(*(save_resul.l_sigoHH_t),DepsBB,tab_ddl
,giB_t,gijBB_t,gijHH_t,giB_tdt,d_giB_tdt,giH_tdt,d_giH_tdt
,epsBB_tdt,d_epsBB,delta_epsBB,gijBB_tdt,gijHH_tdt
,d_gijBB_tdt,d_gijHH_tdt,jacobien_0,jacobien
,d_jacobien_tdt,*(save_resul.l_sigoHH),d_sigRef_HH
,save_resul.l_energ,save_resul.l_energ_t
,compress_inter,cisaill_inter
,ex);
switch (completude_calcul)
{ case CONTRAINTE_ET_TANGENT: // cas normal
{if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{
double coef = save_resul.f_ponder;
compress_inter *= coef;
cisaill_inter *= coef;
energ = save_resul.l_energ * coef; // update des énergies totales
sig_ref_HH = coef*(*(save_resul.l_sigoHH)); // on coefficiente la contrainte
// update de l'opérateur tangent
for (int j=1;j<=taille;j++)
(*d_sigRef_HH(j)) *= coef;
}
else
{energ = save_resul.l_energ; // update des énergies totales
sig_ref_HH = (*(save_resul.l_sigoHH)); // récup de la contrainte
};
break;
};
case CONTRAINTE_UNIQUEMENT:
{if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{
double coef = save_resul.f_ponder;
compress_inter *= coef;
cisaill_inter *= coef;
energ = save_resul.l_energ * coef; // update des énergies totales
sig_ref_HH = coef*(*(save_resul.l_sigoHH)); // on coefficiente la contrainte
}
else
{energ = save_resul.l_energ; // update des énergies totales
sig_ref_HH = (*(save_resul.l_sigoHH)); // récup de la contrainte
};
break;
};
case TANGENT_UNIQUEMENT:
{ // on remet à 0 l'énergie de ref car pas de prise en compte de la contrainte
save_resul.l_energ.Inita(0.);energ = save_resul.l_energ;
sig_ref_HH.Inita(0.);// on remet à 0 la contrainte de référence
(*(save_resul.l_sigoHH)) = sig_ref_HH;
// on ne tiens pas en compte du module de compressibilité ni de cisaillement
module_compressibilite = module_cisaillement = 0.;
if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{ double coef = save_resul.f_ponder;
// update de l'opérateur tangent
for (int j=1;j<=taille;j++)
(*d_sigRef_HH(j)) *= coef;
};
// sinon rien, on utilisera l'opérateur tangent déjà calculé
break;
};
};
};
//----------------- debug
#ifdef MISE_AU_POINT
if (affichage)
{ const Tenseur3BB & gijBB = *((Tenseur3BB*) &gijBB_tdt); // passage en dim 3
cout << "\n epsBB="<<epsBB << " gijBB_tdt= " << gijBB << " gijHH_tdt= " << gijHH<<" ";
Tenseur3BB titi;
epsBB.BaseAbsolue(titi, giH_tdt);
cout << "\n epsBB= en absolue "; titi.Ecriture(cout);
Tenseur3HH toto;
sig_ref_HH.BaseAbsolue(toto,giB_tdt);
cout << "\n sighh_ref= en absolue: ";toto.Ecriture(cout);
};
#endif
//----------- fin debug
// --- on calcule la projection du tenseur de contrainte dans le repère de travail
{ sigHH.Inita(0.);
if (completude_calcul != TANGENT_UNIQUEMENT)
{ for (int i=1;i<4;i++) for (int j=i;j<4;j++)
{double inter = 0.;
for (int k=1;k<4;k++) for (int l=1;l<4;l++) for (int a=1;a<4;a++)for (int b=1;b<4;b++)
inter += (*hij(a,b)) * beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)
* sig_ref_HH(k,l);
sigHH.Coor(i,j) = inter;
};
};
};
// dans le cas où on veut une sortie des grandeurs on sauvegarde
if (sortie_post)
{ if (save_resul.eps_loc_HH == NULL) {save_resul.eps_loc_HH = NevezTenseurHH(3);};
if (save_resul.sig_loc_HH == NULL) {save_resul.sig_loc_HH = NevezTenseurHH(3);};
if (save_resul.para_loi == NULL ) {save_resul.para_loi = new Vecteur (6);};
(*save_resul.eps_loc_HH) = eps_p_HH;
// idem pour les contraintes
(*save_resul.sig_loc_HH) = sigHH;
save_resul.sig_loc_HH->ChBase(gamma); // res = (gamma * res) * gamma.Transpose();
// paramètres de la loi
(*save_resul.para_loi)(1) = A1; (*save_resul.para_loi)(2) = A2;
(*save_resul.para_loi)(3) = A3;
(*save_resul.para_loi)(4) = B12; (*save_resul.para_loi)(5) = B13;
(*save_resul.para_loi)(6) = B23;
};
// passage dans la bonne variance
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Tenseur3BH sigBH = gijBB * sigHH ;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
double trace_sig = sigBH.Trace();
switch (cas_calcul)
{ case 0: // calcul normal (tous les termes)
{ // on ne fait rien de spécial
break;
}
case 1: // calcul de la partie déviatorique seule
{ sigBH -= (untiers * trace_sig) * IdBH3;
sigHH -= (untiers * trace_sig) * gijHH;
break;
}
case 2: // calcul de la partie sphérique seule
{ sigBH = (untiers * trace_sig) * IdBH3;
sigHH = (untiers * trace_sig) * gijHH;
break;
}
default:
{ cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! "
<< "\n Projection_anisotrope_3D::Calcul_DsigmaHH (.... ";
Sortie(1);
}
};
#ifdef MISE_AU_POINT
if (affichage)
{ Tenseur3HH toto;
sigHH.BaseAbsolue(toto,giB_tdt);
cout << "\n sighh= en absolue ";toto.Ecriture(cout);
};
#endif
// traitement des énergies
{ // on commence par calculer l'incrément de l'énergie de référence
EnergieMeca delta_ener = save_resul.l_energ - save_resul.l_energ_t;
// on calcule le rapport des incréments linéarisés des énergies globales
double puiss_deltat = sigHH && delta_epsBB;
double puiss_ref_deltat = (*(save_resul.l_sigoHH)) && delta_epsBB;
double r = puiss_deltat / (puiss_ref_deltat + ConstMath::trespetit);
delta_ener *= r ; // l'incrément dilaté
energ = energ_t + delta_ener; // l'énergie finale
};
// -- calcul des modules
{
// on n'utilise plus la forme linéaire, mais à la place la variation relative de volume
// constaté, au sens d'une mesure logarithmique, sauf dans le cas où cette variation est trop petite
// calcul de la valeur de la variation relative de volume en log
double log_var_vol = log((*(ex.jacobien_tdt))/(*(ex.jacobien_0)));
// pour le module de compressibilité, choix entre les différents cas
if ((cas_calcul == 0) || (cas_calcul == 2))
{// calcul du module initiale : cas particulier
// celui d'une variation purement volumique (cf. théorie)
//K_{s} = \frac{1}{3} {g"}_{ab}~\left (\lambda(a,b) ~{g"}^{ab}\right ) K_{(ref)s}
// avec {g"}_{ab} et {g"}^{ab} les composantes de la métriques dans O'
// dans le cas ou la def est null
// le repère initial est orthonormé d'où un calcul très simple
double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
double module_compressibilite_initial = untiers * coeff * compress_inter;
double module_courant=module_compressibilite_initial; // init
if (log_var_vol > ConstMath::petit)
{module_courant = untiers * sigBH.Trace() / (log_var_vol);};
// maintenant on va essayer de faire un choix
if (module_courant >= module_compressibilite_initial * ratio_inf_module_compressibilite)
{ module_compressibilite = module_courant;}
else
// cas où le module n'est pas recevable (a priori)
{module_compressibilite = module_compressibilite_initial; };
// if (log_var_vol > ConstMath::petit)
// {module_compressibilite = untiers * sigBH.Trace() / (log_var_vol);}
// else // si la variation de volume est trop faible on passe par un cas particulier
// // celui d'une variation purement volumique (cf. théorie)
// {//K_{s} = \frac{1}{3} {g"}_{ab}~\left (\lambda(a,b) ~{g"}^{ab}\right ) K_{(ref)s}
// // avec {g"}_{ab} et {g"}^{ab} les composantes de la métriques dans O'
// // dans le cas ou la def est null
// // le repère initial est orthonormé d'où un calcul très simple
// double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
// module_compressibilite = untiers * coeff * compress_inter;
// };
}
else
// dans le cas d'une loi purement déviatorique, le module de compressibilité = 0 (approximativement)
{module_compressibilite = 0.;};
// pour la partie cisaillement
if ((cas_calcul == 0) || (cas_calcul == 1))
{// on commence par calculer l'intensité du déviateur des déformations
Tenseur3BH eps_barre_BH(epsBH); // init
double I_eps = epsBH.Trace();
eps_barre_BH -= (untiers * I_eps) * IdBH3;
double II_eps = eps_barre_BH.II();
if (II_eps > ConstMath::petit)
{Tenseur3BH sig_barre_BH(sigBH);
sig_barre_BH -= (untiers * trace_sig) * IdBH3;
double II_sig = sig_barre_BH.II();
module_cisaillement = II_sig /II_eps;
}
else // sinon c'est le cas à déformation nulle
{module_cisaillement = 0.5*untiers * (B12+B13+B23)
*cisaill_inter;
};
}
else
// en purement sphérique, le module est supposé nul
{module_cisaillement = 0.; };
};
// cas le la variation du tenseur des contraintes
int nbddl = d_gijBB_tdt.Taille();
// des tenseurs de travail
Tenseur3BH d_sigBH; Tenseur3HH d_sig_HH;
Tenseur3HH d_epsHH;Mat_pleine mat_d_epsHH(3,3);
Vecteur d_eps_ii(3);Vecteur d_sig_ii(3);
Mat_pleine d_beta(3,3);Mat_pleine d_beta_inv(3,3);
Mat_pleine d_beta_transpose(3,3);
Mat_pleine d_gamma(3,3);
Mat_pleine mat_d_eps_p_HH(3,3);
// on récupère la matrice des composantes de déformation dans O_p_a
Mat_pleine mat_epsHH(3,3);epsHH.Matrice_composante(mat_epsHH);
//cout << "\n mat_epsHH: ";mat_epsHH.Affiche();
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 Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
const Tenseur3BB & dgijBB = *((Tenseur3BB*)(d_gijBB_tdt(i))) ; // pour simplifier l'ecriture
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); // "
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
Tenseur3HH & dsig_ref_HH = *((Tenseur3HH*) (d_sigRef_HH(i))); // passage en dim 3
//Tableau <BaseB> * d_giB_tdt
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// BaseB & d_giB = (*ex.d_giB_tdt)(i);
// BaseH & d_giH = (*ex.d_giH_tdt)(i);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// on calcul la variation de matrice de passage de la base g_i vers la base O'_a
// O'_a = beta_a^{.i} * g_i
if (type_transport == 0)
// transport de type contravariant
{for (int a=1; a<4;a++)
{const CoordonneeH& alpha_H_a = alpha_H.CoordoH(a);
// on calcul d'abord le produit alpha_a^l * alpha_a^m * d_g_ij(l,m)
double inter = 0.;
for (int l =1;l<4;l++)
for (int m =1;m<4;m++)
inter += alpha_H_a(l) * alpha_H_a(m) * dgijBB(l,m);
for (int i=1; i<4;i++)
{double & d_beta_ai = d_beta(a,i); // pour simplifier
// ajout des termes or boucle
d_beta_ai *= - 0.5 * alpha_H_a(i) * inter / PUISSN(tab_norme(a), 3);
};
};
// on calcule ensuite la variation de gamma
d_beta_transpose = d_beta.Transpose();
d_gamma = - gamma * d_beta_transpose * gamma;
// on calcule maintenant la variation de gamma
for (int b=1;b<4;b++)
for (int j=1;j<4;j++)
{for (int a=1;a<4;a++)
for (int i=1;i<4;i++)
d_gamma(b,j) -= (gamma(a,j) * gamma(b,i)) * d_beta(a,i);
};
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n d_beta: ";d_beta.Affiche();
cout << "\n d_beta_transpose: ";d_beta_transpose.Affiche();
cout << "\n d_gamma: ";d_gamma.Affiche();
d_gamma = -d_gamma;
cout << "\n d_gamma: ";d_gamma.Affiche();
};
#endif
}
else
{cout << "\n **** cas en attente "
<< "\n Projection_anisotrope_3D::Calcul_DsigmaHH_tdt (.... "<< endl;
Sortie(1);
};
// -- 3) calcul de la variation des contraintes dans le repère g_i
// ici on utilise directement la relation de calcul de contrainte
{ dsigHH.Inita(0.);
//cout << "\n *** debug Projection_anisotrope_3D::Calcul_DsigmaHH_tdt ";
for (int i=1;i<4;i++) for (int j=i;j<4;j++)
// on boucle à partir de j==i jusqu'à 3 inclus -> on ne calcul que la moitié du tenseur
{double inter = 0.;
for (int k=1;k<4;k++) for (int l=1;l<4;l++) for (int a=1;a<4;a++)for (int b=1;b<4;b++)
inter += (*hij(a,b)) *(
(d_beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)* sig_ref_HH(k,l))
+(beta(a,i) * d_beta(b,j) * gamma(a,k) * gamma(b,l)* sig_ref_HH(k,l))
+(beta(a,i) * beta(b,j) * d_gamma(a,k) * gamma(b,l)* sig_ref_HH(k,l))
+(beta(a,i) * beta(b,j) * gamma(a,k) * d_gamma(b,l)* sig_ref_HH(k,l))
+(beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)* dsig_ref_HH(k,l))
) ;
dsigHH.Coor(i,j) = inter;
//if (i != j)
// cout << "\n dsigHH.Coor("<<i<<","<<j<<")= "<<dsigHH.Coor(i,j);
};
};
switch (cas_calcul)
{ case 0: // calcul normal (tous les termes)
{ // on ne fait rien de spécial
break;
}
case 1: // calcul de la partie déviatorique seule
{// passage en mixte pour le calcul de la trace
d_sigBH = dgijBB * sigHH + gijBB * dsigHH;
double d_trace_sig = d_sigBH.Trace();
d_sigBH -= (untiers * d_trace_sig) * IdBH3;
// passage en deux fois contravariant
dsigHH = dgijHH * sigBH + gijHH * d_sigBH;
break;
}
case 2: // calcul de la partie sphérique seule
{// passage en mixte pour le calcul de la trace
d_sigBH = dgijBB * sigHH + gijBB * dsigHH;
double d_trace_sig = d_sigBH.Trace();
d_sigBH = (untiers * d_trace_sig) * IdBH3;
// passage en deux fois contravariant
dsigHH = dgijHH * sigBH + gijHH * d_sigBH;
break;
}
default: break;//l'erreur a déjà été traitée dans le calcul de la contrainte
};
};
delete zeroHH;
LibereTenseur();
};
// void Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB& DepsBB
// ,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
// ,TenseurHH& sigHH,TenseurHHHH& d_sigma_deps
// ,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
// ,const Met_abstraite::Umat_cont& ex) ; //= 0;
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// en_base_orthonormee:
// si true: 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 orthonormeee fixe
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void Projection_anisotrope_3D::Calcul_dsigma_deps
(bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB& DepsBB
,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,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 << " Projection_anisotrope_3D::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
mise à jour 7.010
2023-05-03 17:23:49 +02:00
bool affichage = (Permet_affichage() > 3);
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n --- loi de comportement Projection anisotrope 3D Calcul_dsigma_deps --- ";
Signature_pti_encours(cout);
};
#endif
const Tenseur3HH & gijHH = *((Tenseur3HH*) ex.gijHH_tdt); // " " " "
const Tenseur3BB & gijBB = *((Tenseur3BB*) ex.gijBB_tdt); // " " " "
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " "
Tenseur3HHHH& d_sigma_deps = *((Tenseur3HHHH*) &d_sigma_deps_);
Tenseur3BH epsBH = epsBB * gijHH; // deformation en mixte
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Tenseur3HH epsHH = gijHH * epsBH; // en deuxfois contra
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// int dim_tens = sigHH_tdt.Dimension(); // ici 3, mais on laisse le cas général
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// TenseurHH* zeroHH = (NevezTenseurHH(dim_tens,0.));
// récup du conteneur spécifique
SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveResul);
// on commence par calculer le repère d'anisotropie transporté
// Op_H_i sont les coordonnées contravariantes de la nouvelle base
// donc par rapport à g_i,
BaseH& Op_H = save_resul.Op_H;
Tableau <double> tab_norme(3);
if (type_transport == 0)
// transport de type contravariant
{ // on calcule les coordonnées de la base O' dans la base naturelle
#ifdef MISE_AU_POINT
if (save_resul.O_H == NULL)
{cout << "\n *** erreur, le repere d'anisotropie n'est pas defini "
<<"on ne peut pas continuer " << flush ;
cout << " Projection_anisotrope_3D::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
BaseH& O_H = (*save_resul.O_H); // coordonnées en absolu de la base d'anisotropie
// calcul des coordonnées alpha_a^{.i} = O_H(a) * g^i à t= 0
for (int a = 1;a < 4; a++)
{CoordonneeH& inter = alpha_H.CoordoH(a);
for (int i=1;i < 4;i++)
inter(i)= O_H(a).ScalHH((*ex.giH_0)(i));
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n O_H("<<a<<"): ";O_H(a).Affiche();
cout << "\n giH_0("<<a<<"): "; (*ex.giH_0)(a).Affiche();
};
#endif
};
// calcule des beta_a^{.i} tel que O'_a = beta_a^{.i} \hat{g}_i donc au temps t+dt
for (int a=1; a<4;a++)
{ // tout d'abord la base non normalisée
// récup du conteneur et affectation des coordonnées initiales
// Op_H(a) non normé a les mêmes coordonnées alpha_a^{.j}
// que O_H(a), mais exprimés dans le repère actuel \hat \vec g_j
CoordonneeH& Op_H_a = Op_H.CoordoH(a); // pour simplifier
Op_H_a = alpha_H(a);
// calcul de la norme du vecteur Op_H_a
double norme = sqrt(Op_H_a * gijBB * Op_H_a);
tab_norme(a) = norme;
// coordonnées finales
Op_H_a /= norme;
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n Op_H("<<a<<"): ";
Op_H_a.Affiche();
cout << "\n giB_tdt("<<a<<"): "; (*ex.giB_tdt)(a).Affiche();
cout << "\n giH_tdt("<<a<<"): "; (*ex.giH_tdt)(a).Affiche();
// affichage du repère tournée
CoordonneeB V(3);
for (int j=1;j<4;j++) V += (*ex.giB_tdt)(j) * Op_H_a(j);
cout << "\n op_H("<<a<<"): en absolu " << V ;
};
#endif
};
}
else
// transport de type covariant
{ // on calcule les coordonnées de la base O' dans la base duale
#ifdef MISE_AU_POINT
if (save_resul.O_B == NULL)
{cout << "\n *** erreur, le repere d'anisotropie n'est pas defini "
<<"on ne peut pas continuer " << flush ;
cout << " Projection_anisotrope_3D::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// BaseB& O_B = (*save_resul.O_B); // pour simplifier
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
for (int a=1; a<4;a++)
{ // tout d'abord la base non normalisée
CoordonneeB& Op_B_a = Op_B.CoordoB(a);
// calcul de la norme du vecteur Op_B_a
double norme = sqrt(Op_B_a * gijHH * Op_B_a);
// coordonnées finales
Op_B_a /= norme;
};
// maintenant on calcule la base Op_H correspondante
for (int i=1; i<4;i++)
Op_H.CoordoH(i) = Op_B.CoordoB(i) * gijHH;
};
// on calcul la matrice de passage de la base g_i vers la base O'_i
// O'_i = beta_i^{.j} * g_j
// et on a également \hat{\vec g}^i = {\beta}_{a}^{.i}~\hat{\vec O'}^a
// comme O'_i est déjà exprimé dans g_j, ses coordonnées sont directement béta
for (int i=1; i<4;i++)
for (int j=1; j<4;j++)
{ beta(i,j) = Op_H(i)(j);};
// puis on calcul les coordonnées de la base duale
beta_transpose = beta.Transpose();
gamma = beta_transpose.Inverse();
gamma_transpose = gamma.Transpose();
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n beta: ";beta.Affiche();
cout << "\n beta_transpose: ";beta_transpose.Affiche();
cout << "\n gamma: ";gamma.Affiche();
cout << "\n gamma_transpose: ";gamma_transpose.Affiche();
};
#endif
// changement de base (cf. théorie) : la matrice beta est telle que:
// gpB(i) = beta(i,j) * gB(j) <==> gp_i = beta_i^j * g_j
// calcul des coordonnées de la déformation dans le repère O'_i
#ifdef MISE_AU_POINT
if (affichage)
{// Tenseur3BB eps_p_BB(epsBB);
cout <<"\n eps_p_BB en g^i: "; epsBB.Ecriture(cout);
Tenseur3BB tiutiu;
epsBB.BaseAbsolue(tiutiu,*(ex.giH_tdt));
cout << "\n eps en absolu :";tiutiu.Ecriture(cout);
};
#endif
// il faut passer en 2 fois contravariant
Tenseur3HH eps_p_HH(epsHH); // init
#ifdef MISE_AU_POINT
if (affichage)
{cout <<"\n eps_p_HH en g_i: "; eps_p_HH.Ecriture(cout);
};
#endif
eps_p_HH.ChBase(gamma); // res = (gamma * res) * gamma.Transpose();
#ifdef MISE_AU_POINT
if (affichage)
{cout <<"\n eps_p_HH en O_p: "; eps_p_HH.Ecriture(cout);
cout <<"\n gijHH: "; gijHH.Ecriture(cout);
Tenseur3HH toto(eps_p_HH);
toto.ChBase(beta_transpose);
cout << "\n retour dans la base g^i de eps^{ij} :";toto.Ecriture(cout);
Tenseur3BB titi(gijBB * toto * gijBB);
cout << "\n def eps_ij :";titi.Ecriture(cout);
};
#endif
double untiers = 1./3.;
// dans le cas de fonctions nD récupération des valeurs
if (!Bll_fct_para)
//rappel: Tableau <Fonction_nD* > fct_para; // fonction nD éventuelle d'évolution des paramètres
{// un tableau de travail pour les valeurs sous forme indicée
Tableau < double * > coef(6);
coef(1) = &A1; coef(2) = & A2; coef(3) = &A3;
coef(4) = & B12; coef(5) = & B13; coef(6) = & B23;
// 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;
// on passe en revue les fonctions nD
for (int i=1;i<7;i++)
{if (fct_para(i) != NULL)
{ Fonction_nD* pt_fonct = fct_para(i); // pour simplifier
// on utilise la méthode générique de loi abstraite
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(pt_fonct,1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_expli
,NULL
,NULL
,NULL
);
/* // ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
);
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 1)
{ cout << "\nErreur : la fonction nD relative au parametre materiau " << i
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " Projection_anisotrope_3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
*/
// on récupère le premier élément du tableau uniquement
(*coef(i)) = tab_val(1);
};
};
};
// --- on calcule le tenseur de contrainte de référence et l'opérateur tangent ----
Tenseur3HH sig_ref_HH; // def et init
double compress_inter=0.; double cisaill_inter=0; // init
{ (save_resul.l_sigoHH)->Inita(0.);
// passage des informations spécifique à la loi liste_des_SaveResul
loi_interne->IndiqueSaveResult(save_resul.SaveResul_interne);
loi_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
loi_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
// calcul de la contrainte résultante,
d_sig_deps_3D_HHHH.Inita(0.); // peut-être inutile ***
loi_interne->Calcul_dsigma_deps(en_base_orthonormee,*(save_resul.l_sigoHH_t),DepsBB
,epsBB_tdt,delta_epsBB,jacobien_0,jacobien
,*(save_resul.l_sigoHH),d_sig_deps_3D_HHHH
,save_resul.l_energ,save_resul.l_energ_t
,compress_inter,cisaill_inter
,ex);
switch (completude_calcul)
{ case CONTRAINTE_ET_TANGENT: // cas normal
{if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{
double coef = save_resul.f_ponder;
compress_inter *= coef;
cisaill_inter *= coef;
energ = save_resul.l_energ * coef; // update des énergies totales
sig_ref_HH = coef*(*(save_resul.l_sigoHH)); // on coefficiente la contrainte
// update de l'opérateur tangent
d_sig_deps_3D_HHHH *= coef;
}
else
{energ = save_resul.l_energ; // update des énergies totales
sig_ref_HH = (*(save_resul.l_sigoHH)); // récup de la contrainte
};
break;
};
case CONTRAINTE_UNIQUEMENT:
{if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{
double coef = save_resul.f_ponder;
compress_inter *= coef;
cisaill_inter *= coef;
energ = save_resul.l_energ * coef; // update des énergies totales
sig_ref_HH = coef*(*(save_resul.l_sigoHH)); // on coefficiente la contrainte
}
else
{energ = save_resul.l_energ; // update des énergies totales
sig_ref_HH = (*(save_resul.l_sigoHH)); // récup de la contrainte
};
break;
};
case TANGENT_UNIQUEMENT:
{ // on remet à 0 l'énergie de ref car pas de prise en compte de la contrainte
save_resul.l_energ.Inita(0.);energ = save_resul.l_energ;
sig_ref_HH.Inita(0.);// on remet à 0 la contrainte de référence
(*(save_resul.l_sigoHH)) = sig_ref_HH;
// on ne tiens pas en compte du module de compressibilité ni de cisaillement
module_compressibilite = module_cisaillement = 0.;
if ((ponderation_nD_quelconque != NULL) && (save_resul.f_ponder != 1.))
{ double coef = save_resul.f_ponder;
// update de l'opérateur tangent
d_sig_deps_3D_HHHH *= coef;
};
// sinon rien, on utilisera l'opérateur tangent déjà calculé
break;
};
};
};
#ifdef MISE_AU_POINT
if (affichage)
{cout << "\n sig_refHH :";sig_ref_HH.Ecriture(cout);
};
if (affichage)
{ cout << "\n epsBB="<<epsBB << " gijBB_tdt= " << gijBB << " gijHH_tdt= " << gijHH<<" ";
Tenseur3BB titi;
epsBB.BaseAbsolue(titi, (*ex.giH_tdt));
cout << "\n epsBB= en absolue "; titi.Ecriture(cout);
Tenseur3HH toto;
sig_ref_HH.BaseAbsolue(toto,(*ex.giB_tdt));
cout << "\n sighh_ref= en absolue: ";toto.Ecriture(cout);
cout << "\n d_sigma_deps_ref_HHHH= ";
d_sig_deps_3D_HHHH.Affiche_bidim(cout);
};
#endif
// --- on calcule la projection du tenseur de contrainte dans le repère de travail
{ sigHH.Inita(0.);
if (completude_calcul != TANGENT_UNIQUEMENT)
{ for (int i=1;i<4;i++) for (int j=i;j<4;j++)
{for (int k=1;k<4;k++) for (int l=1;l<4;l++) for (int a=1;a<4;a++)for (int b=1;b<4;b++)
sigHH.Coor(i,j) += (*hij(a,b)) * beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)
* sig_ref_HH(k,l);
};
};
};
// dans le cas où on veut une sortie des grandeurs on sauvegarde
if (sortie_post)
{ if (save_resul.eps_loc_HH == NULL) {save_resul.eps_loc_HH = NevezTenseurHH(3);};
if (save_resul.sig_loc_HH == NULL) {save_resul.sig_loc_HH = NevezTenseurHH(3);};
if (save_resul.para_loi == NULL ) {save_resul.para_loi = new Vecteur (6);};
(*save_resul.eps_loc_HH) = eps_p_HH;
(*save_resul.sig_loc_HH) = sigHH; // le tenseur final en gi
// puis on le met dans le repère d'anisotropie
save_resul.sig_loc_HH->ChBase(gamma); // res = (gamma * res) * gamma.Transpose();
// paramètres de la loi
(*save_resul.para_loi)(1) = A1; (*save_resul.para_loi)(2) = A2;
(*save_resul.para_loi)(3) = A3;
(*save_resul.para_loi)(4) = B12; (*save_resul.para_loi)(5) = B13;
(*save_resul.para_loi)(6) = B23;
};
// passage dans la variance mixte
mise à jour 7.010
2023-05-03 17:23:49 +02:00
Tenseur3BH sigBH = gijBB * sigHH ;
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
double trace_sig = sigBH.Trace();
switch (cas_calcul)
{ case 0: // calcul normal (tous les termes)
{ // on ne fait rien de spécial
break;
}
case 1: // calcul de la partie déviatorique seule
{ sigBH -= (untiers * trace_sig) * IdBH3;
sigHH -= (untiers * trace_sig) * gijHH;
break;
}
case 2: // calcul de la partie sphérique seule
{ sigBH = (untiers * trace_sig) * IdBH3;
sigHH = (untiers * trace_sig) * gijHH;
break;
}
default:
{ cout << "\n erreur l'indicateur cas_calcul= " << cas_calcul << " n'a pas une valeur correcte !! "
<< "\n Projection_anisotrope_3D::Calcul_DsigmaHH_tdt (.... ";
Sortie(1);
}
};
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
{ Tenseur3HH toto;
sigHH.BaseAbsolue(toto,(*ex.giB_tdt));
cout << "\n sighh= en absolue ";toto.Ecriture(cout);
};
#endif
// traitement des énergies
{ // on commence par calculer l'incrément de l'énergie de référence
EnergieMeca delta_ener = save_resul.l_energ - save_resul.l_energ_t;
// on calcule le rapport des incréments linéarisés des énergies globales
double puiss_deltat = sigHH && delta_epsBB;
double puiss_ref_deltat = (*(save_resul.l_sigoHH)) && delta_epsBB;
double r = puiss_deltat / (puiss_ref_deltat + ConstMath::trespetit);
delta_ener *= r ; // l'incrément dilaté
energ = energ_t + delta_ener; // l'énergie finale
};
// -- calcul des modules
{
// on n'utilise plus la forme linéaire, mais à la place la variation relative de volume
// constaté, au sens d'une mesure logarithmique, sauf dans le cas où cette variation est trop petite
// calcul de la valeur de la variation relative de volume en log
double log_var_vol = log((*(ex.jacobien_tdt))/(*(ex.jacobien_0)));
// pour le module de compressibilité, choix entre les différents cas
if ((cas_calcul == 0) || (cas_calcul == 2))
{// calcul du module initiale : cas particulier
// celui d'une variation purement volumique (cf. théorie)
//K_{s} = \frac{1}{3} {g"}_{ab}~\left (\lambda(a,b) ~{g"}^{ab}\right ) K_{(ref)s}
// avec {g"}_{ab} et {g"}^{ab} les composantes de la métriques dans O'
// dans le cas ou la def est null
// le repère initial est orthonormé d'où un calcul très simple
double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
double module_compressibilite_initial = untiers * coeff * compress_inter;
double module_courant=module_compressibilite_initial; // init
if (log_var_vol > ConstMath::petit)
{module_courant = untiers * sigBH.Trace() / (log_var_vol);};
// maintenant on va essayer de faire un choix
if (module_courant >= module_compressibilite_initial * ratio_inf_module_compressibilite)
{ module_compressibilite = module_courant;}
else
// cas où le module n'est pas recevable (a priori)
{module_compressibilite = module_compressibilite_initial; };
// if (log_var_vol > 0.01 ) //ConstMath::petit)
// {module_compressibilite = untiers * sigBH.Trace() / (log_var_vol);}
// else // si la variation de volume est trop faible on passe par un cas particulier
// // celui d'une variation purement volumique (cf. théorie)
// {//K_{s} = \frac{1}{3} {g"}_{ab}~\left (\lambda(a,b) ~{g"}^{ab}\right ) K_{(ref)s}
// // avec {g"}_{ab} et {g"}^{ab} les composantes de la métriques dans O'
// // dans le cas ou la def est null
// // le repère initial est orthonormé d'où un calcul très simple
// double coeff = (*hij(1,1)) + (*hij(2,2)) + (*hij(3,3));
// module_compressibilite = untiers * coeff * compress_inter;
// };
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
cout << "\n module_compressibilite= " << module_compressibilite
<< ", module_loi_interne= " << compress_inter
<< " hij(1,1)= "<< (*hij(1,1))
<< " ,hij(2,2)= "<< (*hij(2,2))
<< " ,hij(3,3)= "<< (*hij(3,3))
<< " ,log_var_vol= "<< log_var_vol
<< " ,sigBH.Trace()= "<< sigBH.Trace()
<< flush;
if (Permet_affichage() > 5) cout << "\n cas_calcul= " << cas_calcul;
#endif
}
else
// dans le cas d'une loi purement déviatorique, le module de compressibilité = 0 (approximativement)
{module_compressibilite = 0.;
#ifdef MISE_AU_POINT
if (Permet_affichage() > 4)
cout << "\n module_compressibilite= " << module_compressibilite;
if (Permet_affichage() > 5) cout << "\n cas_calcul= " << cas_calcul;
#endif
};
// pour la partie cisaillement
if ((cas_calcul == 0) || (cas_calcul == 1))
{// on commence par calculer l'intensité du déviateur des déformations
Tenseur3BH eps_barre_BH(epsBH); // init
double I_eps = epsBH.Trace();
eps_barre_BH -= (untiers * I_eps) * IdBH3;
double II_eps = eps_barre_BH.II();
if (II_eps > ConstMath::petit)
{Tenseur3BH sig_barre_BH(sigBH);
sig_barre_BH -= (untiers * trace_sig) * IdBH3;
module_cisaillement = sig_barre_BH.II()/II_eps;
}
else // sinon c'est le cas à déformation nulle
{module_cisaillement = 0.5*untiers * (B12+B13+B23)
*cisaill_inter;
};
}
else
// en purement sphérique, le module est supposé nul
{module_cisaillement = 0.; };
};
// ----- calcul de l'opérateur tangent -----------
// cas de la variation des beta_a^{.i} par rapport aux eps_kl
// on va les stocker dans un tableau d_beta(a,i) de tenseur HH
Tableau2 <Tenseur3HH> d_beta_HH(3,3);
Tableau2 <Tenseur3HH> d_beta_transpose_HH(3,3);
if (type_transport == 0)
// transport de type contravariant
{ // on calcule les coordonnées de la base O' dans la base naturelle
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// BaseH& O_H = (*save_resul.O_H); // pour simplifier
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
for (int a=1;a<4;a++)
{// calcul de la norme du vecteur O'_a
CoordonneeH alpha_H_a = alpha_H(a);
double n_O_a = tab_norme(a);
for (int i=1;i<4;i++)
{d_beta_transpose_HH(i,a)
= d_beta_HH(a,i)
= -(alpha_H_a(i)
/(n_O_a*n_O_a*n_O_a))*Tenseur3HH::Prod_tensoriel(alpha_H_a,alpha_H_a);
};
};
}; // fin du transport 0
//**** à faire le second transport
// vérif d_beta OK
/* {// vérif en différence finie de d_beta
for (int a=1;a<4;a++)
{ CoordonneeH alpha_H_a = alpha_H(a);
//Tenseur3HH toto(Tenseur3HH::Prod_tensoriel(alpha_H_a,alpha_H_a));
//cout << "\n Prod_tensoriel(alpha_H_a,alpha_H_a)= ";
//toto.Ecriture(cout);
//Tenseur3HH titi;
//for (int j=1;j<4;j++)
// for (int i=1;i<4;i++)
// titi.Coor(i,j) = alpha_H_a(i) * alpha_H_a(j);
//cout << "\n Prod_tensoriel(alpha_H_a,alpha_H_a) numerique = ";
//titi.Ecriture(cout);
// calcul de la norme du vecteur O'_a
double d_beta_num=0.;
double n_O_a = tab_norme(a);
for (int i=1;i<4;i++)
{d_beta_num = d_beta_HH(a,i) && epsBB;
// avec les composantes
double inter=0.;
for (int k=1;k<4;k++) for (int l=1;l<4;l++)
inter += d_beta_HH(a,i)(k,l) * epsBB(k,l);
cout << "\n delta_beta("<<a<<","<<i<<")= "
<< (beta(a,i)-alpha_H_a(i))<< " num: "<<d_beta_num
<< " beta(a,i)= " << beta(a,i) << " alpha_H_a(i)= " << alpha_H_a(i)
<< " en_composante " << inter ;
}
}
};
Sortie(1);
*/
// on calcule maintenant la variation de gamma
Tableau2 <Tenseur3HH> d_gamma_HH(3,3);
for (int b=1;b<4;b++)
for (int j=1;j<4;j++)
{for (int a=1;a<4;a++)
for (int i=1;i<4;i++)
d_gamma_HH(b,j) -= (gamma(a,j) * gamma(b,i)) * d_beta_HH(a,i);
};
// vérification numérique de d_gamma
// a priori ok
/* // on a \hat{\vec g}_j = {\gamma}_{.j}^{b}~\hat{\vec O'}_b
{Mat_pleine gamma_0(3,3);
Mat_pleine O_p_0_ij(3,3);Mat_pleine O_p_tdt_ij(3,3);
for (int a=1;a<4;a++)
{CoordonneeH alpha_H_a = alpha_H(a);
const Tenseur3BB & gijBB_0 = *((Tenseur3BB*) ex.gijBB_0);
CoordonneeB alpha_B_a = alpha_H_a * gijBB_0;
double d_gamma_num=0.;
double n_O_a = tab_norme(a);
for (int i=1;i<4;i++)
{d_gamma_num = d_gamma_HH(a,i) && epsBB;
cout << "\n delta_gamma("<<a<<","<<i<<")= "
<< (gamma(a,i)-alpha_B_a(i))<< " num: "<<d_gamma_num
<< " gamma(a,i)= " << gamma(a,i) << " alpha_B_a(i)= " << alpha_B_a(i) ;
};
};
};
*/
// --- opérateur tangent final ----
// variation des contraintes (composantes dans le repère g_i ) / au eps_kl
// on calcul par composante car il y a peut-être des pb de symétrie qui sont difficiles à maîtriser
// en tensoriel avec un stockage symétrique (il faudra peut-être essayer d'améliorer via les tenseurs
// généraux) mais cela doit-être cependant assez efficace !
for (int m=1;m<4;m++) for (int n=1;n<4;n++)
{for (int i=1;i<4;i++) for (int j=i;j<4;j++)
{double inter = 0.;
for (int k=1;k<4;k++) for (int l=1;l<4;l++) for (int a=1;a<4;a++)for (int b=1;b<4;b++)
inter += (*hij(a,b)) * (
( d_beta_HH(a,i)(m,n) * beta(b,j) * gamma(a,m) * gamma(b,n)
+beta(a,i) * d_beta_HH(b,j)(m,n) * gamma(a,m) * gamma(b,n)
+beta(a,i) * beta(b,j) * d_gamma_HH(a,m)(m,n) * gamma(b,n)
+beta(a,i) * beta(b,j) * gamma(a,m) * d_gamma_HH(b,n)(m,n)
)* sig_ref_HH(k,l)
+ beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)
* d_sig_deps_3D_HHHH(k,l,m,n)
);
dsig_ijkl_HHHH.Change(i,j,m,n,inter);
};
};
#ifdef MISE_AU_POINT
if (affichage)
{ cout << "\n avant transfert de gene en symetrique: d_sigma_depsHHHH= ";
int e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++)
for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<dsig_ijkl_HHHH(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
};
#endif
// on passe au tenseur de retour
//**** pour voir
d_sigma_deps.TransfertDunTenseurGeneral(dsig_ijkl_HHHH.Symetrise1et2_3et4());
//d_sigma_deps.TransfertDunTenseurGeneral(dsig_ijkl_HHHH);
#ifdef MISE_AU_POINT
if (affichage)
{ cout << "\n au final: d_sigma_depsHHHH= ";
d_sigma_deps.Affiche_bidim(cout);
};
#endif
/*
//$$$$$$$$ vérif ************
{ cout << "\n $$$$$$$$ vérifs intermédiaires ************ ";
TenseurQ3geneHHHH dsig_ijkl_HHHH;
{for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++) for (int l=1;l<4;l++)
{ double inter= 0.;
for (int a=1;a<4;a++) for (int b=1;b<4;b++) for (int m=1;m<4;m++) for (int n=1;n<4;n++)
inter += (*hij(a,b)) * (
+ beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)* d_sig_deps_3D_HHHH(m,n,k,l)
);
dsig_ijkl_HHHH.Change(i,j,k,l,inter);
};
};
{ cout << "\n avant transfert de gene en symetrique: d_sigma_depsHHHH= ";
int e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++)
for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<dsig_ijkl_HHHH(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
};
// on passe au tenseur de retour
Tenseur3HHHH inter_HHHH;
inter_HHHH.TransfertDunTenseurGeneral(dsig_ijkl_HHHH.Symetrise1et2_3et4());
//d_sigma_deps.TransfertDunTenseurGeneral(dsig_ijkl_HHHH);
if (affichage)
{ cout << "\n au final: d_sigma_depsHHHH= ";
inter_HHHH.Affiche_bidim(cout);
};
// vérif en accroissement fini / aux eps en cours
Tenseur3HH delta_sig_num_HH = inter_HHHH && epsBB;
cout << "\n delta_sig_num_HH "; delta_sig_num_HH.Ecriture(cout);
cout << "\n sig_HH "; sigHH.Ecriture(cout);
{TenseurQ3geneHHHH toto;
for (int m=1;m<4;m++) for (int n=1;n<4;n++)
{for (int i=1;i<4;i++) for (int j=i;j<4;j++)
{double inter = 0.;
for (int k=1;k<4;k++) for (int l=1;l<4;l++) for (int a=1;a<4;a++)for (int b=1;b<4;b++)
inter += (*hij(a,b)) * (
( d_beta_HH(a,i)(m,n) * beta(b,j) * gamma(a,m) * gamma(b,n)
+beta(a,i) * d_beta_HH(b,j)(m,n) * gamma(a,m) * gamma(b,n)
+beta(a,i) * beta(b,j) * d_gamma_HH(a,m)(m,n) * gamma(b,n)
+beta(a,i) * beta(b,j) * gamma(a,m) * d_gamma_HH(b,n)(m,n)
)* sig_ref_HH(k,l)
+ beta(a,i) * beta(b,j) * gamma(a,k) * gamma(b,l)
* d_sig_deps_3D_HHHH(k,l,m,n)
);
toto.Change(i,j,m,n,inter);
};
};
cout << "\n autre calcul toto= ";
int e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++)
for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<toto(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
// on passe au tenseur de retour
// --- pour voir
d_sigma_deps.TransfertDunTenseurGeneral(toto.Symetrise1et2_3et4());
{ cout << "\n au final à partir de toto: d_sigma_depsHHHH= ";
d_sigma_deps.Affiche_bidim(cout);
};
// vérif en accroissement fini / aux eps en cours
Tenseur3HH delta_sig_num_HH = toto && epsBB;
cout << "\n autre_delta_sig_num_HH "; delta_sig_num_HH.Ecriture(cout);
cout << "\n sig_HH "; sigHH.Ecriture(cout);
};
}
//$$$$$$$$ fin vérif ************
*/
//cout << "\n d_sigma_deps= ";
//d_sigma_deps.Affiche_bidim(cout);
//
// d_sigma_deps = Tenseur3HHHH::Var_tenseur_dans_nouvelle_base
// (beta,var_sig_ab_kl_HHHH,d_beta_HH,sigHH);
//
//cout << "\n (2) d_sigma_deps= ";
//d_sigma_deps.Affiche_bidim(cout);
// Sortie(1);
LibereTenseur();
LibereTenseurQ();
};
// fonction interne utilisée par les classes dérivées de Loi_comp_abstraite
// pour répercuter les modifications de la température
// ici utiliser pour modifier la température des lois élémentaires
// l'Enum_dure: indique quel est la température courante : 0 t ou tdt
void Projection_anisotrope_3D::RepercuteChangeTemperature(Enum_dure temps)
{
loi_interne->temperature_0 = this->temperature_0;
loi_interne->temperature_t = this->temperature_t;
loi_interne->temperature_tdt = this->temperature_tdt;
loi_interne->dilatation=dilatation;
// on répercute également les déformations thermiques, qui ne sont utilisées
// telles quelles que pour certaines lois: ex: loi hyper-élastique
if (dilatation)
{// a- dimensionnement des tenseurs intermédiaires
int dim_tens = epsBB_therm->Dimension();
// -- cas de la déformation
if (loi_interne->epsBB_therm == NULL) { loi_interne->epsBB_therm = NevezTenseurBB(dim_tens);}
else if (loi_interne->epsBB_therm->Dimension() != dim_tens)
{ delete loi_interne->epsBB_therm;loi_interne->epsBB_therm = NevezTenseurBB(dim_tens);};
// -- cas de la vitesse de déformation
if (loi_interne->DepsBB_therm == NULL) { loi_interne->DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (loi_interne->DepsBB_therm->Dimension() != dim_tens)
{ delete loi_interne->DepsBB_therm;loi_interne->DepsBB_totale = NevezTenseurBB(dim_tens);};
// b- affectation des tenseurs
(*loi_interne->epsBB_therm)=(*epsBB_therm);
(*loi_interne->DepsBB_therm)=(*DepsBB_therm);
};
// on répercute sur les lois internes
loi_interne->RepercuteChangeTemperature(temps);
Loi_comp_abstraite * loi_ili = loi_interne; // pour le débug
switch (temps)
{ case TEMPS_0: {loi_interne->temperature = &loi_interne->temperature_0; break;}
case TEMPS_t: {loi_interne->temperature = &loi_interne->temperature_t; break;}
case TEMPS_tdt: {loi_interne->temperature = &loi_interne->temperature_tdt; break;}
default:
{ cout << "\n erreur, cas de temps non prevu !! "
<< "\n Projection_anisotrope_3D::RepercuteChangeTemperature(...";
Sortie(1);
};
};
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 7)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\n init temperatures:\n "
<< " loi_ili->temperature_0= " << loi_ili->temperature_0
<< " loi_ili->temperature_t= " << loi_ili->temperature_t
<< " loi_ili->temperature_tdt= " << loi_ili->temperature_tdt
<< " loi_ili->temperature= " << *(loi_ili->temperature)
<< " \n Projection_anisotrope_3D::RepercuteChangeTemperature(.."
<< endl;
};
#endif
};
// fonction surchargée dans les classes dérivée si besoin est
void Projection_anisotrope_3D::CalculGrandeurTravail
(const PtIntegMecaInterne& ptintmeca
,const Deformation & def,Enum_dure temps,const ThermoDonnee& dTP
,const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat
,const List_io<Ddl_etendu>* exclure_dd_etend
,const List_io<const TypeQuelconque *>* exclure_Q)
{
// récup et enregistrement dans les variables spécifiques au point calculé
SaveResulProjection_anisotrope_3D & save_resul = *((SaveResulProjection_anisotrope_3D*) saveResul);
// ******* prise en compte de la pondération: pour l'instant pas effective
if (ponderation_nD_quelconque != NULL)
{ // dans le cas où il y a une dépendance à des grandeurs locales
Ponderation_TypeQuelconque& pondnD = *ponderation_nD_quelconque;
list <SaveResul*> list_save; // liste inter
if (save_resul.SaveResul_interne != NULL)
list_save.push_back(save_resul.SaveResul_interne);
// on utilise la méthode générique de loi abstraite
Tableau <double> & tab_val = Loi_comp_abstraite::Loi_comp_Valeur_FnD_Evoluee
(pondnD.C_proport(),1 // une seule valeur attendue en retour
,ex_impli,ex_expli_tdt,ex_umat
,exclure_dd_etend
,exclure_Q
,&(list_save)
);
save_resul.f_ponder = tab_val(1); // on sauvegarde
};
// répercution sur la classe dérivée si besoin est
{// passage des informations spécifique
loi_interne->IndiqueSaveResult(save_resul.SaveResul_interne);
loi_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
loi_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
loi_interne->CalculGrandeurTravail(ptintmeca,def,temps,dTP
,ex_impli,ex_expli_tdt,ex_umat,exclure_dd_etend,exclure_Q);
};
};
// vérification et préparation de l'acces aux grandeurs locales
void Projection_anisotrope_3D::Verif_et_preparation_acces_grandeurs_locale()
{// 1) on demande aux lois de prendre en compte éventuellement les grandeurs locales nécessaires
list <EnumTypeQuelconque > listEnuQuelc; // def de la liste des grandeurs quelconque
// ---récup des grandeurs quelconques nécessaires
// tout d'abord au niveau de la pondération
if (ponderation_nD_quelconque != NULL)
{const Tableau <EnumTypeQuelconque>& tab_enu_quelc = ponderation_nD_quelconque->C_proport()->Tab_enu_quelconque();
int tail = tab_enu_quelc.Taille();
if (tail)
{for (int i = 1; i<= tail; i++)
listEnuQuelc.push_back(tab_enu_quelc(i));
};
};
// puis au niveau des fonctions paramètres si elles existent
for (int i=1;i<7;i++)
{if (fct_para(i) != NULL)
{// les grandeurs patentées quelconques
const Tableau <EnumTypeQuelconque> & ti_grandeur = fct_para(i)->Tab_enu_quelconque();
{int taille = ti_grandeur.Taille();
for (int j=1;j<=taille;j++)
listEnuQuelc.push_back(ti_grandeur(j));
};
// les énum étendues
const Tableau <Ddl_enum_etendu>& tab_enu = fct_para(i)->Tab_enu_etendu();
{int taille = tab_enu.Taille();
for (int j=1;j<=taille;j++)
{EnumTypeQuelconque enQ = Ddl_enum_etendu
::Equivalent_en_grandeur_quelconque(tab_enu(j));
if (enQ != RIEN_TYPEQUELCONQUE)
{listEnuQuelc.push_back(enQ);}
else
{ cout << "\n erreur, la grandeur " << tab_enu(j).NomGenerique()
<< " n'a pas d'equivalent en type quelconque "
<< " !! on ne peut pas continuer "
<< "\n Projection_anisotrope_3D::Verif_et_preparation_acces_grandeurs_locale(...";
Sortie(1);
};
};
};
};
};
// nettoyage de la liste
listEnuQuelc.sort();
listEnuQuelc.unique();
// -- on vérifie que la loi interne produira la grandeur quelconque
// ou this
list <EnumTypeQuelconque >::iterator jj,jjfin=listEnuQuelc.end();
for (jj=listEnuQuelc.begin();jj != jjfin; jj++)
{bool existe = false; // par défaut pb
{if (loi_interne->Existe_stockage_grandeurs_quelconques((*jj)))
{existe = true;
loi_interne->Activation_stockage_grandeurs_quelconques(listEnuQuelc);
};
if (this->Existe_stockage_grandeurs_quelconques((*jj)))
{existe = true;
this->Activation_stockage_grandeurs_quelconques(listEnuQuelc);
};
};
if (!existe)
{ cout << "\n *** erreur d'utilisation de grandeur locale "
<< " la grandeur "<<NomTypeQuelconque(*jj)
<< " n'est pas disponible au niveau de la loi interne "
<< " et au niveau de la loi de projection elle-meme "
<< " on ne peut pas continuer !"
<< "\n Projection_anisotrope_3D::Verif_et_preparation_acces_grandeurs_locale(... "
<< endl;
Sortie(1);
};
};
};
Reference in a new issue Copy permalink