3077 lines
158 KiB
C++
Executable file
3077 lines
158 KiB
C++
Executable file
// FICHIER : LoiContraintesPlanesDouble.cp
|
|
// CLASSE : LoiContraintesPlanesDouble
|
|
|
|
|
|
// This file is part of the Herezh++ application.
|
|
//
|
|
// The finite element software Herezh++ is dedicated to the field
|
|
// of mechanics for large transformations of solid structures.
|
|
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
|
|
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
|
|
//
|
|
// Herezh++ is distributed under GPL 3 license ou ultérieure.
|
|
//
|
|
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
|
|
// AUTHOR : Gérard Rio
|
|
// E-MAIL : gerardrio56@free.fr
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License,
|
|
// or (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty
|
|
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
// See the GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
//
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
|
|
|
|
|
//#include "Debug.h"
|
|
#include "LesLoisDeComp.h"
|
|
|
|
# include <iostream>
|
|
using namespace std; //introduces namespace std
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include "Sortie.h"
|
|
#include "TypeQuelconqueParticulier.h"
|
|
#include "TypeConsTens.h"
|
|
#include "NevezTenseurQ.h"
|
|
#include "Util.h"
|
|
#include "ExceptionsLoiComp.h"
|
|
#include "MotCle.h"
|
|
#include "MathUtil.h"
|
|
#include "MathUtil2.h"
|
|
|
|
#include "LoiContraintesPlanesDouble.h"
|
|
|
|
|
|
|
|
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
|
|
// ceci dans le cas où l'axe de traction simple est quelconque
|
|
// ici on suppose que:
|
|
// - tous les calculs sont en dim 3
|
|
// - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
|
|
//
|
|
// ViB et ViH : correspond à la base de traction: les vecteurs sont identiques
|
|
// par contre leurs composantes sont différentes car exprimés dans g^i et g_i
|
|
// telle que \vec ViB(1) = \vec ViH(1) = l'axe de traction
|
|
// La base Vi est supposée orthonormée
|
|
// en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormee
|
|
// le tenseur de déformation et son incrémentsont également en orthonormees
|
|
// si = false: les bases transmises sont utilisées
|
|
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a !!** non pas tout le temps a revoir
|
|
void LoiContraintesPlanesDouble::Calcul_dsigma_deps
|
|
(bool en_base_orthonormee,const BaseB * ViB,const BaseH * ViH
|
|
,const TenseurHH & sigHH_t,const TenseurBB& DepsBB
|
|
,const TenseurBB & epsBB_tdt,const TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
|
|
,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps
|
|
,EnergieMeca & energ,const EnergieMeca &
|
|
,double& module_compressibilite,double& module_cisaillement
|
|
,const Met_abstraite::Umat_cont& ex)
|
|
{ // récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
module_compressibilite=module_cisaillement=0.; // init
|
|
energ.Inita(0.); // initialisation des énergies mises en jeux
|
|
//récup de eps_mecaBB_t
|
|
// 1) dans le repère gi
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
|
|
// 2) dans le dernier repère V_a: a utiliser en sortie post, pas pour le calcul !!
|
|
Tenseur3BB& eps_P_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_P_mecaBB_t));
|
|
Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_resul.eps_P_mecaBB));
|
|
|
|
// on récupère la dimension de l'espace de travail qui doit-être de dimension 3
|
|
// ainsi que le nombre de vecteur: pour l'instant 3
|
|
#ifdef MISE_AU_POINT
|
|
{if (Permet_affichage() > 5)
|
|
cout << "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
|
|
if ((ViB->Dimension() != 3) || (ViH->Dimension() != 3))
|
|
{ cout << "\n***Erreur : la dimension de la base ViB = "<< ViB->Dimension()
|
|
<< " ou la dimension de la base ViH = " << ViH->Dimension()
|
|
<< " est different de 3, on ne peut pas continuer ! "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
|
|
Sortie(1);
|
|
};
|
|
if ((ViB->NbVecteur() != 3) || (ViH->NbVecteur() != 3))
|
|
{ cout << "\n***Erreur : le nombre de vecteur de la base ViB = "
|
|
<< ViB->NbVecteur()
|
|
<< " ou le nombre de vecteur de la base ViH = " << ViH->NbVecteur()
|
|
<< " est different de 3, on ne peut pas continuer ! "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
|
|
Sortie(1);
|
|
};
|
|
// idem pour l'Umat: là on ne test que les vecteurs à 0, pour alléger le test
|
|
if (ex.giB_0->Dimension() != 3)
|
|
{ cout << "\n***Erreur : la dimension de la base ex.giB_0 = "
|
|
<< ex.giB_0->Dimension()
|
|
<< " du conteneur Umat, est different de 3, on ne peut pas continuer ! "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
|
|
Sortie(1);
|
|
};
|
|
if (ex.giB_0->NbVecteur() != 3)
|
|
{ cout << "\n***Erreur : le nombre de vecteur de la base ex.giB_0 = "
|
|
<< ex.giB_0->NbVecteur()
|
|
<< " du conteneur Umat est different de 3, on ne peut pas continuer ! "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
#endif
|
|
|
|
// les coordonnées des vecteurs \vec V_e sont exprimées
|
|
// dans la base de travail g_i finale on a donc
|
|
// Vi_B(i) = {\beta'}_i^{.j} \hat{\vec g}_j = betaP(i,j) * giB_tdt(j)
|
|
// i indice de ligne, j indice de colonne
|
|
// \vec V_i = \vec V^i : les vecteurs sont égaux mais pas les composantes !
|
|
|
|
// 1 -- passage de la base \vec V_e à la base \hat{\vec g}_i
|
|
// et \vec V^e à la base \hat{\vec g}^i
|
|
// calcul des grandeurs gamma et beta:
|
|
// \hat{\vec g}^i = \gamma^i_{~.e}~\vec V^e
|
|
// \hat{\vec g}_i = \beta_i^{~.e}~\vec V_e
|
|
|
|
|
|
// 2 -- passage de la base \hat{\vec g}_i à la base \vec V_e
|
|
// et \hat{\vec g}^i à la base \vec V^e
|
|
//
|
|
// Vi_B(i) = {\beta'}_i^{.j} \hat{\vec g}_j
|
|
// Vi_H(i) = {\gamma'}^i_{.j} \hat{\vec g}^j
|
|
|
|
// gamma3D(i,e) = \gamma^i_{~.e} = \hat{\vec g}^i . \vec V_e
|
|
// ce qui donne que [gamma3D(i,e)] = [{\beta'}_i^{.j}]^T
|
|
for (int i=1;i<4;i++) // calcul de la matrice beta'
|
|
for (int e=1;e<4;e++) // au lieu de passer par un produit scalaire on récupère
|
|
{betaP_3D(i,e) = (*ViH)(i)(e); // directement les coordonnées
|
|
gammaP_3D(i,e) = (*ViB)(i)(e); // " "
|
|
};
|
|
|
|
// \beta_i^{~.e} = \hat{\vec g}_i . \vec V^e
|
|
// c'est aussi l'inverse de la transposé de gamma3D
|
|
// on choisit de calculer avec la formule
|
|
gamma3D = betaP_3D.Transpose(); // ok
|
|
beta3D = gammaP_3D.Transpose(); // betaP_3D.Inverse();
|
|
// gammaP_3D = beta3D.Transpose();
|
|
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 6)
|
|
{cout << "\n Vi_B(1)= "<<(*ViB)(1) << " Vi_B(2)= "<<(*ViB)(2) << " Vi_B(3)= "<<(*ViB)(3);
|
|
cout << "\n Vi_H(1)= "<<(*ViH)(1) << " Vi_H(2)= "<<(*ViH)(2) << " Vi_H(3)= "<<(*ViH)(3);
|
|
cout << "\n Vi_B(i) = {beta'}_i^{.j} hat{vec g}_j --> beta'= ";betaP_3D.Affiche();
|
|
cout << "\n Vi_H(i) = {gamma'}^i_{.j} hat{vec g}^j--> gamma'= ";gammaP_3D.Affiche();
|
|
cout << "\n hat{vec g}^i = gamma^i_{~.e}~vec V^e --> gamma= ";gamma3D.Affiche();
|
|
cout << "\n hat{vec g}_i = beta_i^{~.e}~vec V_e--> beta= ";beta3D.Affiche();
|
|
};
|
|
#endif
|
|
|
|
|
|
// on alimente les pointeurs de tenseurs d'entrée de la cinématique
|
|
epsBB_tdt_cin = ((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
|
|
DepsBB_cin = ((Tenseur3BB*) &DepsBB); // passage en dim 3
|
|
delta_epsBB_cin = ((Tenseur3BB*) &delta_epsBB); // passage en dim 3
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 6)
|
|
{cout << "\n epsBB_tdt_cin "<< (*epsBB_tdt_cin) << "\n DepsBB_cin "<< (*DepsBB_cin)
|
|
<< "\n delta_epsBB_cin "<< (*delta_epsBB_cin) << flush;
|
|
// en base orthonormee
|
|
Tenseur3BB epsBB_tdt_cin_inter;
|
|
epsBB_tdt_cin->BaseAbsolue(epsBB_tdt_cin_inter,(* ex.giH_tdt));
|
|
cout << "\n epsBB_tdt_cin= en absolue 3D : "<<epsBB_tdt_cin_inter;
|
|
};
|
|
#endif
|
|
|
|
// --- initialisation des déformations mécaniques dans le repère V_a
|
|
eps_P_BB = eps_mecaBB_t; // init: là on est dans le repère g^i
|
|
// on passe ensuite dans le repère V_a qui est donc le nouveau repère, pour ce passage
|
|
eps_P_BB.ChBase(betaP_3D);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 6)
|
|
{cout << "\n eps_mecaBB_t= "; eps_mecaBB_t.Ecriture(cout);
|
|
cout << "\n eps_P_BB apres change base donc en V_a "<< eps_P_BB;
|
|
};
|
|
#endif
|
|
|
|
// on calcul la déformation cinématique dans la direction V_1
|
|
// \epsilon'_{11}(t+\Delta t) = \epsilon_{kl}(t+\Delta t) ~\gamma^k_{~.1} ~\gamma^l_{~.1}
|
|
// idem pour l'accroissement de déformation et la vitesse de déformation
|
|
eps_P_BB.Coor(1,1) = 0.; Deps_P_BB.Coor(1,1) = 0.; delta_eps_P_BB.Coor(1,1) = 0.;
|
|
for (int k=1;k<4;k++) // calcul de la matrice gamma
|
|
for (int l=1;l<4;l++)
|
|
{eps_P_BB.Coor(1,1) += (*epsBB_tdt_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
|
|
Deps_P_BB.Coor(1,1) += (*DepsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
|
|
delta_eps_P_BB.Coor(1,1) += (*delta_epsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
|
|
};
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 6)
|
|
{cout << "\n eps_P_BB avec la cinematique "<< eps_P_BB;
|
|
cout << "\n Deps_P_BB avec la cinematique "<< Deps_P_BB;
|
|
cout << "\n delta_eps_P_BB avec la cinematique "<< delta_eps_P_BB << flush;
|
|
};
|
|
#endif
|
|
// on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
|
|
// sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
|
|
sig_HH_t_3D = sigHH_t;
|
|
// concernant la métrique relative au repère de travail, on utilise celle passée en paramètre
|
|
(*umat_cont_3D)= ex;
|
|
// concernant le jacobien, le jacobien du conteneur est supposé être celui 1D
|
|
sauve_jacobien1D_tdt = *(ex.jacobien_tdt);
|
|
|
|
// --- pour les contraintes passage en 3D
|
|
// bool plusZero = true;
|
|
// --- pour la cinématique
|
|
// passage des informations spécifique à la loi le_SaveResul
|
|
lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
|
|
lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
|
|
lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en coursVhVg
|
|
// la base de traction
|
|
ViB_3D = ViB; ViH_3D = ViH;
|
|
// calcul d'élements intermédiaires
|
|
V1V1_BB = Tenseur3BB::Prod_tensoriel((*ViB_3D)(1),(*ViB_3D)(1));
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
{ cout << "\n V1V1_BB: "<< V1V1_BB;
|
|
}
|
|
#endif
|
|
// pour (m,n) = (2,2); (3,3); (1,2), et pour (g,h) = (2,2); (3,3); et (1,2)
|
|
// on va utiliser des matrices intermédiaires tels que:
|
|
// in = 1 pour (2,2) , 2 pour (3,3) et 3 pour (1,2)
|
|
|
|
for (int gh=1;gh<4;gh++)
|
|
{//if (gh <3)
|
|
{VhVg_BB(gh) = Tenseur3BB::Prod_tensoriel((*ViB_3D)(jnd(gh)),(*ViB_3D)(ind(gh)));}
|
|
// else // donne des tenseurs non symétriques pour (1,2) -> on symétrise
|
|
// {VhVg_BB(gh) = 0.5* Tenseur3BB::Prod_tensoriel((*ViB_3D)(ind(gh)),(*ViB_3D)(jnd(gh)))
|
|
// + 0.5* Tenseur3BB::Prod_tensoriel((*ViB_3D)(jnd(gh)),(*ViB_3D)(ind(gh)));
|
|
// }
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
{ cout << "\n VhVg_BB("<<gh<<")= "<< VhVg_BB(gh);
|
|
}
|
|
#endif
|
|
};
|
|
|
|
// on sauvegarde les dimensions transverse, car elles vont être modifiées et si on doit
|
|
// repartir au départ, ce ne sera plus possible sans cette sauvegarde
|
|
double sauve_hsurh0=save_resul.hsurh0;
|
|
double sauve_bsurb0=save_resul.bsurb0;
|
|
// on initialise les jacobiens à 0 et t
|
|
jacobien_0_3D = jacobien_0; // à t=0 -> save_resul.hsurh0 * save_resul.bsurb0 == 1;
|
|
jacobien_t_3D = (*ex.jacobien_t) * save_resul.h_tsurh0 * save_resul.b_tsurb0;
|
|
|
|
// choix entre calcul avec une boucle locale de Newton on non.
|
|
// Dans le premier cas il nous faut la version d_sig/d_eps
|
|
// au lieu de d_sig/d_ddl
|
|
bool mauvaise_convergence = false; // init
|
|
bool newtonlocal = false;
|
|
switch (type_de_contrainte)
|
|
{case NEWTON_LOCAL: // cas de la version d_sig/d_eps
|
|
{ newtonlocal=true; // pour la gestion du cas par défaut
|
|
// initialisation du tenseur final de contrainte
|
|
sig_HH_3D = sigHH_tdt; // c-a-d de la solution: peut-être pas nécessaire ***
|
|
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
|
|
Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
|
|
|
|
// on appel la procédure de résolution de sig33(..)=0
|
|
if ((sortie_post)&&(save_resul.indicateurs_resolution.Taille()!= 2)) // dimensionnement éventuelle de la sortie d'indicateurs
|
|
save_resul.indicateurs_resolution.Change_taille(2);
|
|
// ----- pour ce faire on appelle une méthode de recherche de zero
|
|
// on initialise les inconnues à 0. On pourrait partir des def_P sauvegardés
|
|
// mais si on change d'orientation dans la pratique la convergence est moins bonne ??
|
|
val_initiale_3D.Zero(); // = save_resul.def_P_t;
|
|
// val_initiale_3D = save_resul.def_P_t;
|
|
// on impose que les grandeurs soient dans les limites admises
|
|
double& hsurh0 = save_resul.h_tsurh0;double& bsurb0 = save_resul.b_tsurb0; // init
|
|
// la méthode va mettre à jour les variations bsurb0 et hsurh0 et le jacobien
|
|
if (Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D))
|
|
{if (Permet_affichage() > 3)
|
|
cout << "\n pb limitation val initiales: LoiContraintesPlanesDouble::Calcul_dsigma_deps(... "
|
|
<< flush;
|
|
};
|
|
|
|
//// ---- debut
|
|
// // on vérifie que les anciennes transversales ne sont pas négatives
|
|
// if ((val_initiale(1) < 0.) || (val_initiale(2) < 0.) )
|
|
// { cout << "\n debug LoiContraintesPlanesDouble::Calcul_DsigmaHH_tdt : "
|
|
// << " initialement "
|
|
// << "\n save_resul.b_tsurb0 = " << save_resul.b_tsurb0
|
|
// << ", save_resul.h_tsurh0 " << save_resul.h_tsurh0 << endl;
|
|
// Sortie(1);
|
|
// };
|
|
////------ fin debug
|
|
|
|
residu_3D.Zero(); // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
|
|
derResidu_3D.Initialise(0.); // init de la matrice dérivée à 0.
|
|
try // on met le bloc sous surveillance
|
|
{
|
|
// résolution de l'équation constitutive d'avancement discrétisée en euler implicite
|
|
int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
|
|
|
|
// dans le cas où on utilise une précision qui est pilotée
|
|
{ // opération de transmission de la métrique: encapsulé ici
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Umat_cont* ex_expli = &ex;
|
|
|
|
if (fct_tolerance_residu != 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 = fct_tolerance_residu->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu->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 = fct_tolerance_residu->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 de pilotage de la precision "
|
|
<< " doit calculer un scalaire or le tableau de retour est de taille "
|
|
<< tab_val.Taille() << " ce n'est pas normal !\n";
|
|
cout << " LoiContraintesPlanesDouble::Calcul_dsigma_deps\n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// on récupère le premier élément du tableau uniquement
|
|
double tol = tab_val(1);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
cout << "\n Newton_prec: tol_= "<< tol;
|
|
#endif
|
|
alg_zero.Modif_prec_res_abs(tol);
|
|
};
|
|
|
|
if (fct_tolerance_residu_rel != 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 = fct_tolerance_residu_rel->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu_rel->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 = fct_tolerance_residu_rel->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 de pilotage de la precision relative "
|
|
<< " doit calculer un scalaire or le tableau de retour est de taille "
|
|
<< tab_val.Taille() << " ce n'est pas normal !\n";
|
|
cout << " LoiContraintesPlanesDouble::Calcul_dsigma_deps\n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// on récupère le premier élément du tableau uniquement
|
|
double tol_rel = tab_val(1);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 5)
|
|
cout << ", tol_relative= "<< tol_rel;
|
|
#endif
|
|
alg_zero.Modif_prec_res_rel(tol_rel);
|
|
};
|
|
};
|
|
|
|
bool conver=alg_zero.Newton_raphson
|
|
(*this,&LoiContraintesPlanesDouble::Residu_constitutif_3D
|
|
,&LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D
|
|
,val_initiale_3D,racine_3D,derResidu_3D,nb_incr_total,nb_iter_total
|
|
,maxi_delta_var_eps_sur_iter_pour_Newton);
|
|
if(sortie_post) // sauvegarde éventuelle des indicateurs
|
|
{save_resul.indicateurs_resolution(1)+=nb_incr_total;
|
|
save_resul.indicateurs_resolution(2)+=nb_iter_total;
|
|
};
|
|
// on vérifie qu'il n'y a pas de pb de convergence
|
|
double absracinemax=racine_3D.Max_val_abs();
|
|
if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
|
|
{ if (Permet_affichage() > 0)
|
|
{ cout << "\n**>> non convergence sur l'algo de la resolution de (sigma.V_2=0, sigma.V_3=0),\n";
|
|
if (ptintmeca_en_cours != NULL) ptintmeca_en_cours->Signature();
|
|
cout << ", nb_incr_total=" << nb_incr_total
|
|
<< ", nb_iter_total=" << nb_iter_total<<", val obtenues:"
|
|
<< "\n eps_P(2,2)= " << racine_3D(1) << ", "
|
|
<< " eps_P(3,3)= " << racine_3D(2) << ", "<< " eps_P(1,2)= " << racine_3D(3)
|
|
<< "\n precedentes: "<<eps_mecaBB_t(2,2)<<", "<<eps_mecaBB_t(3,3)<<", "<<eps_mecaBB_t(1,2);
|
|
cout << "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
|
|
};
|
|
// on garde en mémoire
|
|
mauvaise_convergence=true;
|
|
};
|
|
}
|
|
catch (ErrNonConvergence_Newton erreur)
|
|
{
|
|
if (Permet_affichage() > 0)
|
|
{ cout << "\n erreur de non convergence avec l'algo de la resolution de"
|
|
<< " (sig'22, sig'33, sig'12)=(0,0,0) "
|
|
<< " on obtient une valeur infinie ou NaN "
|
|
<< "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
|
|
};
|
|
// on garde en mémoire
|
|
mauvaise_convergence=true;
|
|
}
|
|
catch (ErrSortieFinale)
|
|
// cas d'une direction voulue vers la sortie
|
|
// on relance l'interuption pour le niveau supérieur
|
|
{ ErrSortieFinale toto;
|
|
throw (toto);
|
|
}
|
|
catch ( ... )
|
|
{ // dans le cas d'une erreur inconnue, on génère également une exception
|
|
if (Permet_affichage() > 0)
|
|
{ cout << "\n erreur non identifiee sur l'algo de la resolution de "
|
|
<< " (sig'22, sig'33, sig'12)=(0,0,0) "
|
|
<< "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
|
|
};
|
|
// on garde en mémoire
|
|
mauvaise_convergence=true;
|
|
};
|
|
if (!mauvaise_convergence)
|
|
{ // arrivée ici, cela veut dire que tout à bien fonctionné
|
|
save_resul.def_P = racine_3D; // récup de la solution dans le repère V_a
|
|
// on doit sauvegarder les def méca dans le repère g^i, ce sont directement
|
|
// les dernières composantes utilisées pour l'appel de la loi 3D
|
|
eps_mecaBB = eps_BB_3D; // dans g^i
|
|
eps_P_mecaBB = eps_P_BB; // dans V_a pour le post traitement
|
|
// on met à jour les modules
|
|
module_compressibilite= module_compressibilite_3D;
|
|
module_cisaillement= module_cisaillement_3D;
|
|
|
|
// retour du tenseurs résultat
|
|
sigHH_tdt = sig_HH_3D;
|
|
|
|
// calcul de l'opérateur tangent final: ici en 3D
|
|
// d sigma^{ij} / d epsilon_{kl}
|
|
Calcul2_dsigma_deps(d_sigma_deps);
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 6)
|
|
{ cout << "\nLoiContraintesPlanesDouble::Calcul_dsigma_deps: bonne convergence de Newton ";
|
|
TenseurHH* ptHH = NevezTenseurHH(sig_HH_3D);
|
|
sig_HH_3D.BaseAbsolue(*ptHH,giB_tdt_3D);
|
|
cout << "\n sigma apres CP2: ";
|
|
ptHH->Ecriture(cout);
|
|
delete ptHH;
|
|
};
|
|
if (Permet_affichage() > 5)
|
|
{ cout << "\n def_P(2,2) " << save_resul.def_P(1)
|
|
<< " def_P(3,3) " << save_resul.def_P(2)
|
|
<< " def_P(1,2) " << save_resul.def_P(3)
|
|
<< "\n en g^i: eps_mecaBB.Coor(2,2) " << eps_mecaBB.Coor(2,2)
|
|
<< " eps_mecaBB.Coor(3,3) " << eps_mecaBB.Coor(3,3)
|
|
<< " eps_mecaBB.Coor(1,2) " << eps_mecaBB.Coor(1,2)
|
|
<< "\n en V_a: eps_P_mecaBB.Coor(2,2) " << eps_P_mecaBB.Coor(2,2)
|
|
<< " eps_P_mecaBB.Coor(3,3) " << eps_P_mecaBB.Coor(3,3)
|
|
<< " eps_P_mecaBB.Coor(1,2) " << eps_P_mecaBB.Coor(1,2)
|
|
<< "\n module_compressibilite= " << module_compressibilite
|
|
<< " module_cisaillement= " << module_cisaillement;
|
|
// on affiche les épaisseur et largeur
|
|
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
|
|
cout << "\n hsurh0= "<< hsurh0 << " bsurb0= " << bsurb0 << flush;
|
|
};
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n dans le repere locale d_sigma_deps_3D= \n";
|
|
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<<")= "<<d_sigma_deps_3D(i,j,k,l) << " ; ";
|
|
if (e>6) {cout << "\n"; e=1;}
|
|
};
|
|
Tenseur3HHHH inter_HHHH;
|
|
d_sigma_deps_3D.ChangeBase(inter_HHHH,giB_tdt_3D);
|
|
cout << "\n dans le repere orthonormee d_sigma_deps_3D= \n";
|
|
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<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
|
|
if (e>6) {cout << "\n"; e=1;}
|
|
};
|
|
};
|
|
#endif
|
|
break;
|
|
};
|
|
// sinon on ne fait pas de break, donc on continue avec la méthode de perturbation
|
|
};
|
|
|
|
|
|
default: // cas de la version avec perturbation, on encore du cas où la méthode de Newton n'a pas convergé
|
|
{
|
|
// initialisation du tenseurs contrainte
|
|
sig_HH_3D.Inita(0.);
|
|
// dans le cas ou on est passé par la méthode de newton et que cela n'a pas
|
|
// convergé, on doit réinitialiser les variables de départ
|
|
if ((newtonlocal) && (Permet_affichage() > 3))
|
|
cout << "\n calcul direct sans Newton ";
|
|
save_resul.hsurh0 = sauve_hsurh0;
|
|
save_resul.bsurb0 = sauve_bsurb0;
|
|
|
|
// appel du calcul en perturbation
|
|
Cal_avec_perturbation2();
|
|
// on met à jour les modules
|
|
module_compressibilite= module_compressibilite_3D;
|
|
module_cisaillement= module_cisaillement_3D;
|
|
|
|
if (mauvaise_convergence)
|
|
{if (Permet_affichage() > 0)
|
|
{ cout << ": valeur final de h/h_0(t+dt)= " << save_resul.hsurh0
|
|
<< " et de b/b_0(t+dt)= " << save_resul.bsurb0;
|
|
};
|
|
// on annulle les dérivées des épaisseurs
|
|
save_resul.d_hsurh0.Zero();save_resul.d_bsurb0.Zero();
|
|
};
|
|
|
|
}
|
|
break;
|
|
};
|
|
|
|
// sig_HH_3D.Ecriture(cout);
|
|
|
|
(*save_resul.l_sigoHH) = sig_HH_3D; // sauvegarde en locale
|
|
// ---calcul des invariants de déformation et de vitesse de déformation
|
|
Calcul_invariants_et_def_cumul();
|
|
|
|
energ = save_resul.l_energ; // récup des énergies
|
|
|
|
// --- on remet à jour éventuellement l'épaisseur
|
|
// dans le cas d'une contrainte par perturbation
|
|
// if ((type_de_contrainte == PERTURBATION )||mauvaise_convergence)
|
|
if (type_de_contrainte == PERTURBATION )
|
|
// calcul de la déformation d'épaisseur correspondant à la condition de contraintes planes
|
|
{ //Tenseur1BH sigBH = gijBB_tdt * sigHH_tdt;
|
|
Tenseur3BH sigBH = gijBB_tdt_3D * sig_HH_3D;
|
|
Tenseur3BH sigBH_t = gijBB_t_3D * sig_HH_t_3D;
|
|
Calcul_eps_trans_parVarVolume2(sigBH_t,jacobien_0,module_compressibilite,jacobien,sigBH,*(ex.jacobien_t));
|
|
if ((Permet_affichage() > 0) && mauvaise_convergence)
|
|
{ cout << "\n valeur final de h/h_0(t+dt)= " << save_resul.hsurh0
|
|
<< " et de b/b_0(t+dt)= " << save_resul.bsurb0;
|
|
};
|
|
};
|
|
|
|
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 LoiContraintesPlanesDouble::RepercuteChangeTemperature(Enum_dure temps)
|
|
{ lois_interne->temperature_0 = this->temperature_0;
|
|
lois_interne->temperature_t = this->temperature_t;
|
|
lois_interne->temperature_tdt = this->temperature_tdt;
|
|
lois_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 (lois_interne->epsBB_therm == NULL) { lois_interne->epsBB_therm = NevezTenseurBB(dim_tens);}
|
|
else if (lois_interne->epsBB_therm->Dimension() != dim_tens)
|
|
{ delete lois_interne->epsBB_therm;lois_interne->epsBB_therm = NevezTenseurBB(dim_tens);};
|
|
// -- cas de la vitesse de déformation
|
|
if (lois_interne->DepsBB_therm == NULL) { lois_interne->DepsBB_therm = NevezTenseurBB(dim_tens);}
|
|
else if (lois_interne->DepsBB_therm->Dimension() != dim_tens)
|
|
{ delete lois_interne->DepsBB_therm;lois_interne->DepsBB_totale = NevezTenseurBB(dim_tens);};
|
|
// b- affectation des tenseurs
|
|
(*lois_interne->epsBB_therm)=(*epsBB_therm);
|
|
(*lois_interne->DepsBB_therm)=(*DepsBB_therm);
|
|
};
|
|
// puis en interne
|
|
lois_interne->RepercuteChangeTemperature(temps);
|
|
switch (temps)
|
|
{ case TEMPS_0:
|
|
{lois_interne->temperature = &lois_interne->temperature_0;
|
|
break;
|
|
}
|
|
case TEMPS_t:
|
|
{lois_interne->temperature = &lois_interne->temperature_t;
|
|
break;
|
|
}
|
|
case TEMPS_tdt:
|
|
{lois_interne->temperature = &lois_interne->temperature_tdt;
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n erreur, cas de temps non prevu !! "
|
|
<< "\n LoiContraintesPlanesDouble::RepercuteChangeTemperature(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
};
|
|
|
|
void LoiContraintesPlanesDouble::CalculGrandeurTravail
|
|
(const PtIntegMecaInterne& ptintmec
|
|
,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 des infos spécifiques
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
|
|
// par défaut on active les invariants et on intègre les infos def en 3D
|
|
if (!calcul_en_3D_via_direction_quelconque)
|
|
{ //dans le cas classique, on dispose de la variation d'épaisseur et de largeur pour déterminer les def
|
|
|
|
// ptintmec concerne des tenseurs d'ordres 1 donc il faut transformer en 3D
|
|
bool plusZero=true; // on ajoute des 0 pour les grandeurs manquantes dans une première étape
|
|
// on commence par transférer le tout
|
|
ptintmeca.Affectation_1D_a_3D(ptintmec,plusZero);
|
|
// ensuite on traite les points particuliers
|
|
// pour tout ce qui est contrainte: l'ajout de 0 est ok
|
|
|
|
// pour les déformations on met à jour en fonction des grandeurs sauvegardées
|
|
double eps33=0.; double eps33_t=0.; //init par défaut
|
|
double eps22=0.; double eps22_t=0.; //init par défaut
|
|
|
|
// -- maintenant on tient compte de l'élongation d'épaisseur
|
|
switch (type_de_deformation)
|
|
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
|
|
// cas d'une déformation d'Almansi
|
|
{ // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
|
|
eps33 = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
|
|
eps33_t = 0.5 * (save_resul.h_tsurh0 * save_resul.h_tsurh0 - 1.);
|
|
eps22 = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
|
|
eps22_t = 0.5 * (save_resul.b_tsurb0 * save_resul.b_tsurb0 - 1.);
|
|
};
|
|
break;
|
|
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
|
|
// cas d'une def logarithmique ou une approximation
|
|
{ eps33 = log(save_resul.hsurh0);
|
|
eps33_t = log(save_resul.h_tsurh0);
|
|
eps22 = log(save_resul.bsurb0);
|
|
eps22_t = log(save_resul.b_tsurb0);
|
|
};
|
|
break;
|
|
default :
|
|
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
|
|
<< Nom_type_deformation(type_de_deformation);
|
|
cout << "\n LoiContraintesPlanesDouble::CalculGrandeurTravail \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
(*ptintmeca.EpsBB()).Coor(3,3) = eps33;
|
|
(*ptintmeca.DeltaEpsBB()).Coor(3,3) = eps33 - eps33_t;
|
|
(*ptintmeca.EpsBB()).Coor(2,2) = eps22;
|
|
(*ptintmeca.DeltaEpsBB()).Coor(2,2) = eps22 - eps22_t;
|
|
|
|
(*ptintmeca.DepsBB()).Coor(3,3) = (eps33 - eps33_t) * unSurDeltat;
|
|
(*ptintmeca.DepsBB()).Coor(2,2) = (eps22 - eps22_t) * unSurDeltat;
|
|
}
|
|
else
|
|
// dans le cas où on travaille dans le repère gi en 3D on dispose directement des infos de déformation
|
|
{//récup de eps_mecaBB
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
|
|
// on sauve ce qu'il y a dans la direction 11 pour éviter les erreurs d'arrondis
|
|
double eps11= (*ptintmeca.EpsBB())(1,1);
|
|
double Deps11= (*ptintmeca.DepsBB())(1,1);
|
|
double delta_eps11= (*ptintmeca.DeltaEpsBB())(1,1);
|
|
// les manips entre tenseurs
|
|
(*ptintmeca.EpsBB()) = eps_mecaBB;
|
|
(*ptintmeca.DeltaEpsBB()) = eps_mecaBB - eps_mecaBB_t;
|
|
(*ptintmeca.DepsBB()) = (*ptintmeca.DeltaEpsBB()) * unSurDeltat;
|
|
// récup des valeurs sauvegardées (normalement ne devrait pas changer...)
|
|
(*ptintmeca.EpsBB()).Coor(1,1) = eps11;
|
|
(*ptintmeca.DeltaEpsBB()).Coor(1,1) = delta_eps11;
|
|
(*ptintmeca.DepsBB()).Coor(1,1) = Deps11;
|
|
};
|
|
|
|
// on met à jour les invariants 3D, on récupère ce qui a été sauvegardé
|
|
if (ptintmeca.Statut_Invariants_deformation()) // le conteneur des invariants a pu être effacé par Affectation_1D_a_3D
|
|
ptintmeca.EpsInvar() = save_resul.epsInvar;
|
|
if (ptintmeca.Statut_Invariants_vitesseDeformation()) // le conteneur des invariants a pu être effacé par Affectation_1D_a_3D
|
|
ptintmeca.DepsInvar() = save_resul.depsInvar;
|
|
// idem au niveau des déformations cumulées
|
|
ptintmeca.Deformation_equi() = save_resul.def_equi;
|
|
ptintmeca.Deformation_equi_t() = save_resul.def_equi_t;
|
|
|
|
|
|
// passage des informations spécifique à la loi le_SaveResul
|
|
lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
|
|
lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
|
|
lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
|
|
|
|
lois_interne->CalculGrandeurTravail(ptintmeca,def,temps,dTP
|
|
,ex_impli,ex_expli_tdt,ex_umat,exclure_dd_etend,exclure_Q); // passage à la loi 3D
|
|
};
|
|
|
|
|
|
// passage des grandeurs métriques de l'ordre 1 à 3, cas implicite
|
|
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Impli& ex)
|
|
{// on s'occupe du redimensionnement éventuel
|
|
// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
|
|
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
|
|
{if (ex.gradVmoyBB_t != NULL) {impli_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
|
|
if (ex.gradVmoyBB_tdt != NULL) {impli_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
|
|
if (ex.gradVBB_tdt != NULL) {impli_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
|
|
|
|
// on s'occupe tout d'abord des tableaux de vitesses qui peuvent ne pas exister
|
|
if (ex.d_gradVmoyBB_t != NULL) // cas où il existe
|
|
{ int tail=ex.d_gradVmoyBB_t->Taille();
|
|
if (d_gradVmoyBB_t_3D_P == NULL) {d_gradVmoyBB_t_3D_P = new Tableau<TenseurBB *> (tail);
|
|
impli_3D->d_gradVmoyBB_t = d_gradVmoyBB_t_3D_P;}
|
|
else {d_gradVmoyBB_t_3D_P->Change_taille(tail);};
|
|
d_gradVmoyBB_t_3D.Change_taille(tail);
|
|
for (int i=1;i<= tail;i++)
|
|
(*d_gradVmoyBB_t_3D_P)(i) = &(d_gradVmoyBB_t_3D(i));
|
|
};
|
|
if (ex.d_gradVmoyBB_tdt != NULL) // cas où il existe
|
|
{ int tail=ex.d_gradVmoyBB_tdt->Taille();
|
|
if (d_gradVmoyBB_tdt_3D_P == NULL) {d_gradVmoyBB_tdt_3D_P = new Tableau<TenseurBB *> (tail);
|
|
impli_3D->d_gradVmoyBB_tdt = d_gradVmoyBB_tdt_3D_P;}
|
|
else {d_gradVmoyBB_tdt_3D_P->Change_taille(tail);};
|
|
d_gradVmoyBB_tdt_3D.Change_taille(tail);
|
|
for (int i=1;i<= tail;i++)
|
|
(*d_gradVmoyBB_tdt_3D_P)(i) = &(d_gradVmoyBB_tdt_3D(i));
|
|
};
|
|
if (ex.d_gradVBB_t != NULL) // cas où il existe
|
|
{ int tail=ex.d_gradVBB_t->Taille();
|
|
if (d_gradVBB_t_3D_P == NULL) {d_gradVBB_t_3D_P = new Tableau<TenseurBB *> (tail);
|
|
impli_3D->d_gradVBB_t = d_gradVBB_t_3D_P;}
|
|
else {d_gradVBB_t_3D_P->Change_taille(tail);};
|
|
d_gradVBB_t_3D.Change_taille(tail);
|
|
for (int i=1;i<= tail;i++)
|
|
(*d_gradVBB_t_3D_P)(i) = &(d_gradVBB_t_3D(i));
|
|
};
|
|
if (ex.d_gradVBB_tdt != NULL) // cas où il existe
|
|
{ int tail=ex.d_gradVBB_tdt->Taille();
|
|
if (d_gradVBB_tdt_3D_P == NULL) {d_gradVBB_tdt_3D_P = new Tableau<TenseurBB *> (tail);
|
|
impli_3D->d_gradVBB_tdt = d_gradVBB_tdt_3D_P;}
|
|
else {d_gradVBB_tdt_3D_P->Change_taille(tail);};
|
|
d_gradVBB_tdt_3D.Change_taille(tail);
|
|
for (int i=1;i<= tail;i++)
|
|
(*d_gradVBB_tdt_3D_P)(i) = &(d_gradVBB_tdt_3D(i));
|
|
};
|
|
|
|
}; // fin de la partie dédiée à la vitesse
|
|
// maintenant on s'occupe uniquement du redimensionnement des tableaux restants
|
|
// -- on s'occupe des tableaux nécessaire à la métrique
|
|
int ta_ex_d_giB_tdt = ex.d_giB_tdt->Taille(); // la dimension
|
|
if (d_giB_tdt_3D.Taille() != ta_ex_d_giB_tdt)
|
|
{ d_giB_tdt_3D.Change_taille(ta_ex_d_giB_tdt);};
|
|
int ta_ex_d_giH_tdt = ex.d_giH_tdt->Taille(); // la dimension
|
|
if (d_giH_tdt_3D.Taille() != ta_ex_d_giH_tdt)
|
|
{ d_giH_tdt_3D.Change_taille(ta_ex_d_giH_tdt);};
|
|
int ta_ex_d_gijBB_tdt = ex.d_gijBB_tdt->Taille(); // la dimension
|
|
if (d_gijBB_tdt_3D.Taille() != ta_ex_d_gijBB_tdt)
|
|
{ d_gijBB_tdt_3D.Change_taille(ta_ex_d_gijBB_tdt); d_gijBB_tdt_3D_P.Change_taille(ta_ex_d_gijBB_tdt);
|
|
for (int i=1;i<=ta_ex_d_gijBB_tdt;i++) {d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));};
|
|
};
|
|
int ta_ex_d_gijHH_tdt = ex.d_gijHH_tdt->Taille();
|
|
if (d_gijHH_tdt_3D.Taille() != ta_ex_d_gijHH_tdt)
|
|
{ d_gijHH_tdt_3D.Change_taille(ta_ex_d_gijHH_tdt); d_gijHH_tdt_3D_P.Change_taille(ta_ex_d_gijHH_tdt);
|
|
for (int i=1;i<=ta_ex_d_gijHH_tdt;i++) {d_gijHH_tdt_3D_P(i) = &(d_gijHH_tdt_3D(i));};
|
|
};
|
|
if( d_jacobien_tdt_3D.Taille() != ex.d_jacobien_tdt->Taille())
|
|
{d_jacobien_tdt_3D.Change_taille(ex.d_jacobien_tdt->Taille());}
|
|
|
|
// le tableau à double dimension
|
|
if (ex.d2_gijBB_tdt != NULL)
|
|
{ int taille1=ex.d2_gijBB_tdt->Taille1(); int taille2=ex.d2_gijBB_tdt->Taille2();
|
|
if (d2_gijBB_tdt_3D_P == NULL)
|
|
{ d2_gijBB_tdt_3D_P = impli_3D->d2_gijBB_tdt = new Tableau2<TenseurBB *> (taille1,taille2);
|
|
d2_gijBB_tdt_3D.Change_taille(taille1,taille2);
|
|
for (int i=1;i<=taille1;i++) for (int j=1;j<= taille2;j++)
|
|
{ (*d2_gijBB_tdt_3D_P)(i,j) = &d2_gijBB_tdt_3D(i,j);};
|
|
}
|
|
else if ((taille1 != d2_gijBB_tdt_3D.Taille1()) || (taille2 != d2_gijBB_tdt_3D.Taille2()))
|
|
{ d2_gijBB_tdt_3D_P->Change_taille(taille1,taille2);
|
|
d2_gijBB_tdt_3D.Change_taille(taille1,taille2);
|
|
for (int i=1;i<=taille1;i++) for (int j=1;j<= taille2;j++)
|
|
{ (*d2_gijBB_tdt_3D_P)(i,j) = &d2_gijBB_tdt_3D(i,j);};
|
|
};
|
|
};
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
|
|
// si les grandeurs à 0 ont déjà été calculées, on les sauvegarde d'abord car elles vont
|
|
// être écrasées
|
|
// on commence par recopier les grandeurs de l'ordre 1 à 3
|
|
bool plusZero = true; // on complète avec des 0 dans un premier temps
|
|
int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
|
|
if ( save_resul.meti_ref_00.giB_ != NULL)
|
|
type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
|
|
impli_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3("
|
|
<< "\n 1) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
|
|
};
|
|
#endif
|
|
|
|
// maintenant on s'occupe de mettre à jour les grandeurs manquantes
|
|
// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
|
|
// et sont identiques pour les bases naturelles et duales
|
|
// Dans le cas d'un comportement isotrope, on choisit un repère par défaut
|
|
|
|
// si les repères associés ne sont pas définis on les définit
|
|
if ( save_resul.meti_ref_00.giB_ == NULL)
|
|
{ Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
|
|
save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
|
|
// on s'occupe de la création du premier repère
|
|
CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
|
|
giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
|
|
giB0_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
|
|
// deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
|
|
if (ex.giB_0->NbVecteur() > 1)
|
|
{ // a priori le second vecteur de la base est un vecteur d'un plan de ref
|
|
giB0_3 = (Util::ProdVec_coorB(giB_0_3D(1), (*(ex.giB_0))(2))).Normer();
|
|
giB0_2 = Util::ProdVec_coorB(giB0_3,giB0_1);
|
|
}
|
|
else
|
|
{ // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
|
|
MathUtil2::Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
|
|
};
|
|
////--- debug
|
|
// cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
|
|
// cout << "\n giB0_1: "; giB0_1.Affiche();
|
|
// cout << "\n giB0_2: "; giB0_2.Affiche();
|
|
// cout << "\n giB0_3: "; giB0_3.Affiche();
|
|
//// --- fin debug
|
|
giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
|
|
giH_0_3D.CoordoH(2) = giB_0_3D(2).Bas_haut();
|
|
giH_0_3D.CoordoH(3) = giB_0_3D(3).Bas_haut();
|
|
};
|
|
////--- debug
|
|
// cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
|
|
// cout << "\n giH_0(2): "; giH_0_3D(2).Affiche();
|
|
// cout << "\n giH_0(3): "; giH_0_3D(3).Affiche();
|
|
//// --- fin debug
|
|
// on met à jour les deux autres directions
|
|
// *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
|
|
// *** , mais il faudra surement améliorer
|
|
CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
|
|
giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
|
|
giBt_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
|
|
// deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
|
|
if (ex.giB_t->NbVecteur() > 1)
|
|
{ // a priori le second vecteur de la base est un vecteur d'un plan de ref
|
|
giBt_3 = (Util::ProdVec_coorB(giB_t_3D(1), (*(ex.giB_t))(2))).Normer();
|
|
giBt_2 = Util::ProdVec_coorB(giBt_3,giBt_1);
|
|
}
|
|
else
|
|
{ // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
|
|
MathUtil2::Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
|
|
};
|
|
giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
|
|
giH_t_3D.CoordoH(2) = giB_t_3D(2).Bas_haut();
|
|
giH_t_3D.CoordoH(3) = giB_t_3D(3).Bas_haut();
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n 2) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
|
|
};
|
|
#endif
|
|
|
|
|
|
CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
|
|
giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
|
|
giBtdt_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
|
|
// deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
|
|
if (ex.giB_tdt->NbVecteur() > 1)
|
|
{ // a priori le second vecteur de la base est un vecteur d'un plan de ref
|
|
giBtdt_3 = (Util::ProdVec_coorB(giB_tdt_3D(1), (*(ex.giB_tdt))(2))).Normer();
|
|
giBtdt_2 = Util::ProdVec_coorB(giBtdt_3,giBtdt_1);
|
|
}
|
|
else
|
|
{ // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
|
|
MathUtil2::Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
|
|
};
|
|
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
|
|
giH_tdt_3D.CoordoH(2) = giB_tdt_3D(2).Bas_haut();
|
|
giH_tdt_3D.CoordoH(3) = giB_tdt_3D(3).Bas_haut();
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n 3) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
|
|
};
|
|
#endif
|
|
// ici on calcul les vecteurs
|
|
// >>>> et la première partie de leurs variations, fonction du premier vecteur uniquement
|
|
|
|
// - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
|
|
gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
|
|
=gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=1.;
|
|
gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
|
|
=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
|
|
|
|
// lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0 et idem pour l'épaisseur
|
|
// ne sont pas affectés car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
|
|
// ensuite elles évolueront, idem pour la largeur
|
|
bool premier_passage = false;
|
|
if (save_resul.h_tsurh0 == 0.)
|
|
{ save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
|
|
save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
|
|
premier_passage = true;
|
|
};
|
|
if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
|
|
{ save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
|
|
if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
|
|
{ save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
|
|
// on calcul les vecteur 2 et 3 naturels et duals
|
|
giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
|
|
giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;
|
|
|
|
giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
|
|
giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n 4) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
|
|
};
|
|
#endif
|
|
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
|
|
// sont normées initialement, seule l'intensité
|
|
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
|
|
double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
|
|
gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
|
|
|
|
double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
|
|
gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
|
|
|
|
double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
|
|
gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;
|
|
|
|
double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
|
|
gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;
|
|
|
|
// -- cas des variations
|
|
|
|
if (save_resul.d_hsurh0.Taille() != ta_ex_d_giB_tdt)
|
|
{save_resul.d_hsurh0.Change_taille(ta_ex_d_giB_tdt); // donc à 0 la première fois
|
|
save_resul.d_bsurb0.Change_taille(ta_ex_d_giB_tdt);
|
|
};
|
|
|
|
// - les variations du vecteur normal
|
|
#ifdef MISE_AU_POINT
|
|
if (ex.d_giB_tdt->Taille() != ex.d_giH_tdt->Taille())
|
|
{ cout << "\n **** sans doute une erreur: (ex.d_giB_tdt->Taille() "<< ex.d_giB_tdt->Taille()
|
|
<< " != ex.d_giH_tdt->Taille()) " << ex.d_giH_tdt->Taille()
|
|
<< "\n LoiDeformationsPlanes::Passage_metrique_ordre1_vers_3(const Met_abstraite::Impli& ex)";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
|
|
// maintenant on s'occupe de la seconde partie des variations
|
|
// >>>> fonction du changement d'épaisseur et de largeur
|
|
|
|
if (!premier_passage)
|
|
{Vecteur& d_hsurh0 = save_resul.d_hsurh0;
|
|
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
|
|
for (int iddl=1;iddl<= ta_ex_d_giB_tdt;iddl++)
|
|
{ double d_hsurh0_iddl = d_hsurh0(iddl);
|
|
// gib(3) = h/h0 * N -> d_gib(3) = d (h/h0 * N) = d_(h/h0) * N + h/h0 * d_N
|
|
// actuellement, c'est la partie d_N qui est stockée dans d_gib(3), on rajoute donc la partie restante
|
|
d_giB_tdt_3D(iddl).CoordoB(3) *= save_resul.hsurh0;
|
|
d_giB_tdt_3D(iddl).CoordoB(3) += d_hsurh0_iddl * giB_normer_3_tdt_3D_sauve; // base naturelle
|
|
// gih(3) = 1/giB(3) -> d_gih(3)(3) = - d_giB(3)(3) / (giB(3)*giB(3))
|
|
d_giH_tdt_3D(iddl).CoordoH(3) = (d_giB_tdt_3D(iddl)(3) / (-hsurh0_carre)).Bas_haut() ; // base duale
|
|
// gijBB(3)= giB(3) * giB(3) -> d_gijBB(3)= 2. * d_giB(3) * giB(3)
|
|
d_gijBB_tdt_3D(iddl).Coor(3,3) = 2. * (d_giB_tdt_3D(iddl)(3).ScalBB(giB_tdt_3D(3)));
|
|
// gijHH(3)= giH(3) * giH(3) -> d_gijHH(3)= 2. * d_giH(3) * giH(3)
|
|
d_gijHH_tdt_3D (iddl).Coor(3,3) = 2. * (d_giH_tdt_3D(iddl)(3).ScalHH(giH_tdt_3D(3)));
|
|
|
|
double d_bsurb0_iddl = d_bsurb0(iddl);
|
|
// gib(2) = b/b0 * N -> d_gib(2) = d (b/b0 * N) = d_(b/b0) * N + b/b0 * d_N
|
|
// actuellement, c'est la partie d_N qui est stockée dans d_gib(3), on rajoute donc la partie restante
|
|
d_giB_tdt_3D(iddl).CoordoB(2) *= save_resul.bsurb0;
|
|
d_giB_tdt_3D(iddl).CoordoB(2) += d_bsurb0_iddl * giB_normer_3_tdt_3D_sauve; // base naturelle
|
|
// gih(2) = 1/giB(2) -> d_gih(2)(2) = - d_giB(3)(3) / (giB(3)*giB(3))
|
|
d_giH_tdt_3D(iddl).CoordoH(2) = (d_giB_tdt_3D(iddl)(2) / (-bsurb0_carre)).Bas_haut() ; // base duale
|
|
// gijBB(2)= giB(2) * giB(2) -> d_gijBB(2)= 2. * d_giB(2) * giB(2)
|
|
d_gijBB_tdt_3D(iddl).Coor(2,2) = 2. * (d_giB_tdt_3D(iddl)(2).ScalBB(giB_tdt_3D(2)));
|
|
// gijHH(2)= giH(2) * giH(2) -> d_gijHH(2)= 2. * d_giH(2) * giH(2)
|
|
d_gijHH_tdt_3D (iddl).Coor(2,2) = 2. * (d_giH_tdt_3D(iddl)(2).ScalHH(giH_tdt_3D(2)));
|
|
|
|
};
|
|
};
|
|
|
|
};
|
|
|
|
// passage des grandeurs métriques de l'ordre 1 à 3: cas explicite
|
|
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Expli_t_tdt& ex)
|
|
{ // on s'occupe du redimensionnement éventuel
|
|
// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
|
|
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
|
|
{if (ex.gradVmoyBB_t != NULL) {expli_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
|
|
if (ex.gradVmoyBB_tdt != NULL) {expli_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
|
|
if (ex.gradVBB_tdt != NULL) {expli_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
|
|
}; // fin de la partie dédiée à la vitesse
|
|
// la dimension des tableaux : on considère que tous les tableaux doivent avoir la même dimension
|
|
int ta_ex_d_gijBB_tdt = ex.d_gijBB_tdt->Taille(); // la dimension
|
|
// maintenant on s'occupe uniquement du redimensionnement des tableaux restants
|
|
int ta_d_gijBB_tdt_3D = d_gijBB_tdt_3D.Taille();
|
|
if (ta_d_gijBB_tdt_3D != ta_ex_d_gijBB_tdt)
|
|
{// cela veut dire que tous les tableaux sont mal dimensionnés
|
|
// -- on s'occupe des tableaux nécessaire à la métrique
|
|
d_gijBB_tdt_3D.Change_taille(ta_ex_d_gijBB_tdt); d_gijBB_tdt_3D_P.Change_taille(ta_ex_d_gijBB_tdt);
|
|
for (int i=1;i<=ta_ex_d_gijBB_tdt;i++) {d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));};
|
|
};
|
|
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
|
|
// on commence par recopier les grandeurs de l'ordre 1 à 3
|
|
bool plusZero = true; // on complète avec des 0 dans un premier temps
|
|
int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
|
|
if ( save_resul.meti_ref_00.giB_ != NULL)
|
|
type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
|
|
expli_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
|
|
|
|
// maintenant on s'occupe de mettre à jour les grandeurs manquantes
|
|
// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
|
|
// et sont identiques pour les bases naturelles et duales
|
|
// Dans le cas d'un comportement isotrope, on choisit un repère par défaut
|
|
|
|
// si les repères associés ne sont définis on les définit
|
|
if ( save_resul.meti_ref_00.giB_ == NULL)
|
|
{ Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
|
|
save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
|
|
// on s'occupe de la création du premier repère
|
|
CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
|
|
giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
|
|
giB0_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
|
|
// Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
|
|
giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
|
|
};
|
|
// on met à jour les deux autres directions
|
|
// *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
|
|
// *** , mais il faudra surement améliorer
|
|
CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
|
|
giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
|
|
giBt_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
|
|
// Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
|
|
giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
|
|
|
|
|
|
CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
|
|
giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
|
|
giBtdt_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
|
|
// Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
|
|
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
|
|
|
|
// - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
|
|
gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
|
|
=gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=1.;
|
|
gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
|
|
=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
|
|
|
|
// lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
|
|
// ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
|
|
// ensuite elles évolueront, idem pour la largeur
|
|
bool premier_passage = false;
|
|
if (save_resul.h_tsurh0 == 0.)
|
|
{ save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
|
|
save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
|
|
premier_passage = true;
|
|
};
|
|
if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
|
|
{ save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
|
|
if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
|
|
{ save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
|
|
// on calcul les vecteur 2 et 3 naturels et duals
|
|
giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
|
|
giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;
|
|
|
|
giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
|
|
giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
|
|
|
|
|
|
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
|
|
// sont normées initialement, seule l'intensité
|
|
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
|
|
double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
|
|
gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;
|
|
|
|
double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
|
|
gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;
|
|
|
|
double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
|
|
gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
|
|
|
|
double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
|
|
gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
|
|
|
|
};
|
|
|
|
|
|
// passage des grandeurs métriques de l'ordre 1 à 3: cas Umat
|
|
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Umat_cont& ex)
|
|
{ // on s'occupe du redimensionnement éventuel
|
|
// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
|
|
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
|
|
{if (ex.gradVmoyBB_t != NULL) {umat_cont_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
|
|
if (ex.gradVmoyBB_tdt != NULL) {umat_cont_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
|
|
if (ex.gradVBB_tdt != NULL) {umat_cont_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
|
|
}; // fin de la partie dédiée à la vitesse
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
|
|
// on commence par recopier les grandeurs de l'ordre 1 à 3
|
|
bool plusZero = true; // on complète avec des 0 dans un premier temps
|
|
int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
|
|
if ( save_resul.meti_ref_00.giB_ != NULL)
|
|
type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
|
|
umat_cont_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
|
|
|
|
// maintenant on s'occupe de mettre à jour les grandeurs manquantes
|
|
// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
|
|
// et sont identiques pour les bases naturelles et duales
|
|
// Dans le cas d'un comportement isotrope, on choisit un repère par défaut
|
|
|
|
// si les repères associés ne sont définis on les définit
|
|
if ( save_resul.meti_ref_00.giB_ == NULL)
|
|
{ Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
|
|
save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
|
|
// on s'occupe de la création du premier repère
|
|
CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
|
|
giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
|
|
giB0_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
|
|
MathUtil2::Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
|
|
// comme la base est orthonormée, les giH sont identiques aux giB
|
|
(*(save_resul.meti_ref_00.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_00.giB_)));
|
|
// on affecte les vecteurs qui seront transmis via l'umat
|
|
giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
|
|
BaseH& a0H = (*(save_resul.meti_ref_00.giH_)); // pour simplifier
|
|
giH_0_3D.CoordoH(1)=a0H(1);giH_0_3D.CoordoH(2)=a0H(2);giH_0_3D.CoordoH(3)=a0H(3);
|
|
};
|
|
////--- debug
|
|
// cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
|
|
// cout << "\n giH_0(2): "; giH_0_3D(2).Affiche();
|
|
// cout << "\n giH_0(3): "; giH_0_3D(3).Affiche();
|
|
//// --- fin debug
|
|
// on met à jour les deux autres directions
|
|
// *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
|
|
// *** , mais il faudra surement améliorer
|
|
CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
|
|
giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
|
|
giBt_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
|
|
MathUtil2::Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
|
|
// comme la base est orthonormée, les giH sont identiques aux giB
|
|
// ensuite plus bas, les normes seront changées pour tenir compte des variations d'épaisseur et de largeur
|
|
(*(save_resul.meti_ref_t.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_t.giB_)));
|
|
// on affecte les vecteurs qui seront transmis via l'umat
|
|
giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
|
|
BaseH& a_t_H = (*(save_resul.meti_ref_t.giH_)); // pour simplifier
|
|
giH_t_3D.CoordoH(1)=a_t_H(1);giH_t_3D.CoordoH(2)=a_t_H(2);giH_t_3D.CoordoH(3)=a_t_H(3);
|
|
|
|
//--- cas à tdt
|
|
CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
|
|
giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
|
|
giBtdt_1.Normer(); // on le norme
|
|
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
|
|
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
|
|
MathUtil2::Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
|
|
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
|
|
// comme la base est orthonormée ici, les giH sont identiques aux giB
|
|
// ensuite plus bas, les normes seront changées pour tenir compte des variations d'épaisseur et de largeur
|
|
(*(save_resul.meti_ref_tdt.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_tdt.giB_)));
|
|
// on affecte les vecteurs qui seront transmis via l'umat
|
|
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
|
|
BaseH& a_tdt_H = (*(save_resul.meti_ref_tdt.giH_)); // pour simplifier
|
|
giH_tdt_3D.CoordoH(1)=a_tdt_H(1);giH_tdt_3D.CoordoH(2)=a_tdt_H(2);giH_tdt_3D.CoordoH(3)=a_tdt_H(3);
|
|
|
|
// - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
|
|
gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
|
|
=gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=1.;
|
|
gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
|
|
=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
|
|
|
|
// lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
|
|
// ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
|
|
// ensuite elles évolueront, idem pour la largeur
|
|
bool premier_passage = false;
|
|
if (save_resul.h_tsurh0 == 0.)
|
|
{ save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
|
|
premier_passage = true;
|
|
};
|
|
if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
|
|
{ save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
|
|
// on peut avoir la même chose en b, indépendamment, lorsque l'on passe de contrainte plane
|
|
// à doublement plane
|
|
if (save_resul.b_tsurb0 == 0.)
|
|
{ save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
|
|
premier_passage = true;
|
|
};
|
|
if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
|
|
{ save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
|
|
// on calcul les vecteur 2 et 3 naturels et duals
|
|
giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
|
|
giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;
|
|
|
|
giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
|
|
giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
|
|
|
|
|
|
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
|
|
// sont normées initialement, seule l'intensité
|
|
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
|
|
double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
|
|
gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;
|
|
|
|
double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
|
|
gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;
|
|
|
|
double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
|
|
gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
|
|
|
|
double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
|
|
gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
|
|
|
|
};
|
|
|
|
|
|
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
|
|
// utilisée pour les méthodes 1 et 2
|
|
// entre autre: calcul du jacobien 3D en fonction de l'élongation 1D et des variations h/h0 et b/b0
|
|
|
|
void LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3
|
|
// (int cas, TenseurBB& DepsBB,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>* d_epsBB
|
|
(TenseurBB& DepsBB,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>* d_epsBB
|
|
,TenseurBB & delta_epsBB,const double& jacobien_0,double& jacobien
|
|
,Vecteur* d_jacobien_tdt,const double& jacobien_t)
|
|
|
|
{ // au début on complète avec des 0, puis dans un second temps on tient compte de la déformation d'épaisseur eps33
|
|
bool plusZero = true;
|
|
Deps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &DepsBB),plusZero);
|
|
eps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &epsBB_tdt),plusZero);
|
|
delta_eps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &delta_epsBB),plusZero);
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
if(d_epsBB != NULL)
|
|
{ int taille = d_epsBB->Taille();
|
|
// redimensionnement éventuel
|
|
int ta_d_eps_BB_3D = d_eps_BB_3D.Taille();
|
|
if (ta_d_eps_BB_3D != taille)
|
|
{ // cela veut dire que tous les tableaux sont mal dimensionnés
|
|
d_eps_BB_3D.Change_taille(taille); d_eps_BB_3D_P.Change_taille(taille);
|
|
for (int i=1;i<=taille;i++) {d_eps_BB_3D_P(i) = &(d_eps_BB_3D(i));}
|
|
};
|
|
// passage 1D 3D
|
|
for (int i=1;i<=taille;i++)
|
|
{ d_eps_BB_3D(i).Affectation_trans_dimension(*((Tenseur1BB*) (*d_epsBB)(i)),plusZero);}
|
|
if (save_resul.d_hsurh0.Taille() != taille)
|
|
save_resul.d_hsurh0.Change_taille(taille); // donc à 0 la première fois
|
|
if (save_resul.d_bsurb0.Taille() != taille)
|
|
save_resul.d_bsurb0.Change_taille(taille); // donc à 0 la première fois
|
|
|
|
};
|
|
// -- calcul du jacobien 3D
|
|
|
|
/* la partie cas = 2: ne fonctionne pas bien le choix est fait
|
|
d'utiliser les variations d'épaisseur et de largeur
|
|
// dans le cas d'un calcul dans une direction quelconque
|
|
// on utilise la déformation mécanique sauvegardée pour calculer les jacobiens
|
|
if (cas == 2)
|
|
{ //récup de eps_mecaBB_t
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
// --- on s'occupe des jacobiens
|
|
// à t le jacobien est calculée à partir de la déformation supposée d'Almansi
|
|
// ici il s'agit du jacobien relatif au comportement mécanique
|
|
// |J| = sqrt(det(2 * eps_ij_t + g_ij(0)))
|
|
// calcul de la matrice correspondante à (2 * eps_ij_t + g_ij(0))
|
|
gij_meca_BB_t = (2. * eps_mecaBB_t + (*(umat_cont_3D->gijBB_0)));
|
|
gij_meca_BB_t.Matrice_composante(mat_inter);
|
|
// calcul du jacobien
|
|
jacobien_t_3D = sqrt(Dabs(mat_inter.Determinant()));
|
|
|
|
// on fait de même pour le jacobien initial ce qui permet d'être cohérent
|
|
// avec la métrique initiale et la déformation au sens d'Almansi
|
|
(umat_cont_3D->gijBB_0)->Matrice_composante(mat_inter);
|
|
jacobien_0_3D = sqrt(Dabs(mat_inter.Determinant()));
|
|
// le jacobien à tdt sera calculé pendant la résolution
|
|
}
|
|
else if (cas == 1)
|
|
*/
|
|
{// première phase on se contente de recopier, ce qui suppose une épaisseur = 1.
|
|
jacobien_0_3D = jacobien_0; // -> partie 1D
|
|
jacobien_tdt_3D = sauve_jacobien1D_tdt = jacobien;
|
|
// ici, l'épaisseur initiale h0 est toujours 1., donc h/h0 = également h et idem pour b/b0
|
|
jacobien_tdt_3D *= save_resul.hsurh0 * save_resul.bsurb0;
|
|
// même travail pour t
|
|
jacobien_t_3D = jacobien_t;
|
|
jacobien_t_3D *= save_resul.h_tsurh0 * save_resul.b_tsurb0;
|
|
};
|
|
|
|
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
|
|
// cas des def, delta def et D
|
|
if (!calcul_en_3D_via_direction_quelconque)
|
|
{ //dans le cas classique, on dispose de la variation d'épaisseur et de largeur pour déterminer les def
|
|
switch (type_de_deformation)
|
|
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
|
|
// cas d'une déformation d'Almansi
|
|
{ // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
|
|
// dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
|
|
eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
|
|
delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
|
|
// idem pour la largeur
|
|
eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
|
|
delta_eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - save_resul.b_tsurb0 * save_resul.b_tsurb0);
|
|
Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3)* unSurDeltat;
|
|
Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2)* unSurDeltat;
|
|
|
|
|
|
// dans le cas où on calcule les variations, on intègre également les variations de eps33 et eps22
|
|
if (d_epsBB != NULL)
|
|
{ Vecteur& d_hsurh0 = save_resul.d_hsurh0;
|
|
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
|
|
int taille = d_epsBB->Taille();
|
|
for (int i=1;i<=taille;i++)
|
|
{d_eps_BB_3D(i).Coor(3,3)= save_resul.hsurh0 * d_hsurh0(i);
|
|
d_eps_BB_3D(i).Coor(2,2)= save_resul.bsurb0 * d_bsurb0(i);
|
|
}
|
|
};
|
|
};
|
|
break;
|
|
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
|
|
// cas d'une def logarithmique ou une approximation
|
|
{ eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
|
|
delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
|
|
// idem pour la largeur
|
|
eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0);
|
|
delta_eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0) - log(save_resul.b_tsurb0);
|
|
Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3)* unSurDeltat;
|
|
Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2)* unSurDeltat;
|
|
|
|
// dans le cas où on calcule les variations, on intègre également les variations de eps33 et eps22
|
|
if (d_epsBB != NULL)
|
|
{ Vecteur& d_hsurh0 = save_resul.d_hsurh0;
|
|
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
|
|
int taille = d_epsBB->Taille();
|
|
// d_eps_BB_3D(3,3) = d_hsurh0 * exp(save_resul.hsurh0)
|
|
for (int i=1;i<=taille;i++)
|
|
{d_eps_BB_3D(i).Coor(3,3)= d_hsurh0(i) * eps_BB_3D(3,3);
|
|
d_eps_BB_3D(i).Coor(2,2)= d_bsurb0(i) * eps_BB_3D(2,2);
|
|
}
|
|
};
|
|
};
|
|
break;
|
|
default :
|
|
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
|
|
<< Nom_type_deformation(type_de_deformation);
|
|
cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3 \n";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
// dans le cas où on travaille dans le repère gi en 3D on dispose directement des infos de déformation
|
|
{//récup de eps_mecaBB
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
|
|
// on sauve ce qu'il y a dans la direction 11 pour éviter les erreurs d'arrondis
|
|
double eps11= (*ptintmeca.EpsBB())(1,1);
|
|
double Deps11= (*ptintmeca.DepsBB())(1,1);
|
|
double delta_eps11= (*ptintmeca.DeltaEpsBB())(1,1);
|
|
// les manips entre tenseurs
|
|
(*ptintmeca.EpsBB()) = eps_mecaBB;
|
|
(*ptintmeca.DeltaEpsBB()) = eps_mecaBB - eps_mecaBB_t;
|
|
(*ptintmeca.DepsBB()) = (*ptintmeca.DeltaEpsBB()) * unSurDeltat;
|
|
// récup des valeurs sauvegardées (normalement ne devrait pas changer...)
|
|
(*ptintmeca.EpsBB()).Coor(1,1) = eps11;
|
|
(*ptintmeca.DeltaEpsBB()).Coor(1,1) = delta_eps11;
|
|
(*ptintmeca.DepsBB()).Coor(1,1) = Deps11;
|
|
};
|
|
|
|
if(d_jacobien_tdt != NULL)
|
|
{ Vecteur& d_hsurh0 = save_resul.d_hsurh0;
|
|
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
|
|
// jacobien_tdt_3D = hsurh0 * jacobien_1D
|
|
// --> d_jacobien_tdt_3D = d_hsurh0 * jacobien_1D + hsurh0 * d_jacobien_1D
|
|
d_jacobien_tdt_3D = d_hsurh0 * save_resul.bsurb0 * jacobien
|
|
+ save_resul.hsurh0 * d_bsurb0 * jacobien
|
|
+ save_resul.hsurh0 * save_resul.bsurb0 * (*d_jacobien_tdt);
|
|
};
|
|
|
|
//---debug
|
|
//cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3, eps_BB_3D(3,3) = "
|
|
// << eps_BB_3D(3,3) << " jacobien_tdt_3D= " << jacobien_tdt_3D;
|
|
// --- fin debug
|
|
|
|
|
|
};
|
|
|
|
|
|
// mise à jour des informations liées à la déformation de 1 vers 3
|
|
// uniquement méthode 1
|
|
void LoiContraintesPlanesDouble::Mise_a_jour_deformations_et_Jacobien_en_3D
|
|
()
|
|
|
|
{ // récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
// -- calcul du jacobien 3D
|
|
// première phase on se contente de recopier, ce qui suppose une épaisseur = 1. et une largeur = 1.
|
|
jacobien_tdt_3D = sauve_jacobien1D_tdt * save_resul.hsurh0 * save_resul.bsurb0;
|
|
// -- maintenant on tient compte de la déformation d'épaisseur
|
|
switch (type_de_deformation)
|
|
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
|
|
// cas d'une déformation d'Almansi
|
|
{ // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
|
|
// dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
|
|
eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
|
|
delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
|
|
// idem pour la direction 2
|
|
eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
|
|
delta_eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - save_resul.b_tsurb0 * save_resul.b_tsurb0);
|
|
};
|
|
break;
|
|
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
|
|
// cas d'une def logarithmique ou une approximation
|
|
{ eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
|
|
delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
|
|
// idem pour la direction 2
|
|
eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0);
|
|
delta_eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0) - log(save_resul.b_tsurb0);
|
|
};
|
|
break;
|
|
default :
|
|
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
|
|
<< Nom_type_deformation(type_de_deformation);
|
|
cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3 \n";
|
|
Sortie(1);
|
|
};
|
|
// -- cas de la vitesse de déformation
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2) * unSurDeltat;
|
|
Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3) * unSurDeltat;
|
|
|
|
};
|
|
|
|
|
|
// calcul de la déformation d'épaisseur et de largeur
|
|
// cas de la méthode 1
|
|
void LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume1(Tenseur3BH& sigBH_t,double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
|
|
,Tenseur3BH& sigBH,double& jaco_1D_t)
|
|
{ // récup des infos spécifiques
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
// on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
|
|
// avec trace(sigma) en 1D = idem en 3D car on est en contraintes doublement planes
|
|
// jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_0*h_0*b_0) = delta_V
|
|
// on considère de plus que l'on est isotrope donc h/h_0 = b/b_0
|
|
// en appelant B= trace(sigma)/3 /K_moy on a:
|
|
// (h/h_0)^2 = jaco1D_0 / jaco1D_tdt * 1./(1.-B) d'où la déformation logarithmique
|
|
if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
|
|
{ //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
|
|
double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
|
|
double exp_log_VsurV0 = exp(log_VsurV0);
|
|
double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
|
|
// on test les mini-maxi
|
|
if (hSurh_0 < mini_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
|
|
<< " on le ramene a la limite "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume(... " << endl;
|
|
};
|
|
hSurh_0 = mini_hsurh0;
|
|
};
|
|
if (hSurh_0 > maxi_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
|
|
<< " on le ramene a la limite "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume(... " << endl;
|
|
};
|
|
hSurh_0 = maxi_hsurh0;
|
|
};
|
|
// on met à jour épaisseur et largeur d'une même valeur
|
|
save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
|
|
save_resul.bsurb0 = hSurh_0; // élongation de largeur
|
|
|
|
////--- debug
|
|
// save_resul.hsurh0 = save_resul.h_tsurh0; // élongation d'épaisseur
|
|
// save_resul.bsurb0 = save_resul.b_tsurb0; // élongation d'épaisseur
|
|
////-----fin debug
|
|
|
|
|
|
}
|
|
else
|
|
{ save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
|
|
|
|
// ---- la suite est un essai, mais ne semble pas non plus fonctionner avec la loi critère
|
|
// car le jacobien à t n'est pas calculé ?? ce qui n'est pas normal
|
|
// ticket #120
|
|
|
|
// // a priori cela ne marche pas bien: on va tenter avec le delta jacobien
|
|
// // on utilise la relation: (V(tdt)-V(t))/V(tdt) = trace(delta_sigma)/3 /K_moy
|
|
// // avec trace(delta_sigma) en 1D = idem en 3D car on est en contraintes doublement planes
|
|
// // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_t*h_t*b_t) = delta_V
|
|
// // on considère de plus que l'on est isotrope donc delta_h/h = delta_b/b d'où h/h_t = b/b_t
|
|
// // en appelant B= trace(delta_sigma)/3 /K / (h_t*b_t) on a:
|
|
// // (h/h_t)^2 = jaco1D_t / (jaco1D_tdt * (1.-B))
|
|
// if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
|
|
// { //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
|
|
// double B = (sigBH(1,1)-sigBH_t(1,1)) / (3. * module_compressibilite );
|
|
// double hSurh_t = sqrt((jaco_1D_t / jaco_1D_tdt) / (1.-B)) ;
|
|
//
|
|
// save_resul.hsurh0 = hSurh_t * save_resul.h_tsurh0;
|
|
// save_resul.bsurb0 = hSurh_t * save_resul.b_tsurb0; // élongation d'épaisseur
|
|
// }
|
|
// else
|
|
// { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
|
|
|
|
|
|
};
|
|
|
|
// calcul de la déformation d'épaisseur et de largeur
|
|
// cas de la méthode 2
|
|
void LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume2(Tenseur3BH& sigBH_t,double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
|
|
,Tenseur3BH& sigBH,double& jaco_1D_t)
|
|
{ // récup des infos spécifiques
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
// on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
|
|
// avec trace(sigma) en 1D = idem en 3D car on est en contraintes doublement planes
|
|
// jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_0*h_0*b_0) = delta_V
|
|
// on considère de plus que l'on est isotrope donc h/h_0 = b/b_0
|
|
// en appelant B= trace(sigma)/3 /K_moy on a:
|
|
// (h/h_0)^2 = jaco1D_0 / jaco1D_tdt * 1./(1.-B) d'où la déformation logarithmique
|
|
if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
|
|
{ //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
|
|
double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
|
|
double exp_log_VsurV0 = exp(log_VsurV0);
|
|
double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
|
|
// on test les mini-maxi
|
|
if (hSurh_0 < mini_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
|
|
<< " on le ramene a la limite "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume2(... " << endl;
|
|
};
|
|
hSurh_0 = mini_hsurh0;
|
|
};
|
|
if (hSurh_0 > maxi_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
|
|
<< " on le ramene a la limite "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume1(... " << endl;
|
|
};
|
|
hSurh_0 = maxi_hsurh0;
|
|
};
|
|
// on met à jour épaisseur et largeur d'une même valeur
|
|
save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
|
|
save_resul.bsurb0 = hSurh_0; // élongation de largeur
|
|
|
|
////--- debug
|
|
// save_resul.hsurh0 = save_resul.h_tsurh0; // élongation d'épaisseur
|
|
// save_resul.bsurb0 = save_resul.b_tsurb0; // élongation d'épaisseur
|
|
////-----fin debug
|
|
|
|
|
|
}
|
|
else
|
|
{ save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
|
|
|
|
// ---- la suite est un essai, mais ne semble pas non plus fonctionner avec la loi critère
|
|
// car le jacobien à t n'est pas calculé ?? ce qui n'est pas normal
|
|
// ticket #120
|
|
|
|
// // a priori cela ne marche pas bien: on va tenter avec le delta jacobien
|
|
// // on utilise la relation: (V(tdt)-V(t))/V(tdt) = trace(delta_sigma)/3 /K_moy
|
|
// // avec trace(delta_sigma) en 1D = idem en 3D car on est en contraintes doublement planes
|
|
// // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_t*h_t*b_t) = delta_V
|
|
// // on considère de plus que l'on est isotrope donc delta_h/h = delta_b/b d'où h/h_t = b/b_t
|
|
// // en appelant B= trace(delta_sigma)/3 /K / (h_t*b_t) on a:
|
|
// // (h/h_t)^2 = jaco1D_t / (jaco1D_tdt * (1.-B))
|
|
// if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
|
|
// { //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
|
|
// double B = (sigBH(1,1)-sigBH_t(1,1)) / (3. * module_compressibilite );
|
|
// double hSurh_t = sqrt((jaco_1D_t / jaco_1D_tdt) / (1.-B)) ;
|
|
//
|
|
// save_resul.hsurh0 = hSurh_t * save_resul.h_tsurh0;
|
|
// save_resul.bsurb0 = hSurh_t * save_resul.b_tsurb0; // élongation d'épaisseur
|
|
// }
|
|
// else
|
|
// { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
|
|
|
|
|
|
};
|
|
|
|
// calcul des déformations d'épaisseur et de largeur et des variations
|
|
// méthode 1
|
|
void LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
|
|
,TenseurHH& sigHH_,Vecteur& d_jaco_1D
|
|
,Tableau <TenseurHH *>& d_sigHH,Tableau <TenseurBB *>& d_gijBB_tdt,TenseurBB & gijBB_)
|
|
{ // récup des infos spécifiques
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
// on utilise la relation: log(V/V0) = trace(sigma)/3 /K_moy
|
|
// avec trace(sigma) en 1D = idem en 3D car on est en contraintes planes double
|
|
// jacobien 1D = la longueur, donc le volume: V = jacobien 1D * h * b -> (jaco1D_tdt * h * b)/(jaco1D_0*h_0*b_0) = V/V0
|
|
// en appelant log_VsurV0= trace(sigma)/3 /K_moy on a:
|
|
// h/h_0 * b/b_0= jaco1D_0 / jaco1D_tdt * exp(log_VsurV0)
|
|
// ici sans info sup, on considère une loi isotrope d'où h/h_0 = b/b_0
|
|
// -> h/h_0 = (jaco1D_0 / jaco1D_tdt * exp(log_VsurV0))^(0.5)
|
|
if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
|
|
{ const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_); // pour simplifier
|
|
const Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_); // pour simplifier
|
|
Tenseur1BH sigBH = gijBB * sigHH;
|
|
double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
|
|
double exp_log_VsurV0 = exp(log_VsurV0);
|
|
double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
|
|
// on test les mini-maxi
|
|
if (hSurh_0 < mini_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
|
|
<< " on le ramene a la limite "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(... " << endl;
|
|
};
|
|
hSurh_0 = mini_hsurh0;
|
|
};
|
|
if (hSurh_0 > maxi_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
|
|
<< " on le ramene a la limite "
|
|
<< "\n LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(... " << endl;
|
|
};
|
|
hSurh_0 = maxi_hsurh0;
|
|
};
|
|
// on met à jour épaisseur et largeur d'une même valeur
|
|
save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
|
|
save_resul.bsurb0 = hSurh_0; // élongation de largeur
|
|
|
|
// maintenant le cas de la variation de h/h0
|
|
// on tient compte de la variation relative à la dépendance du jacobien et à la trace de sigma
|
|
// a) vérification et modif éventuelle de la taille
|
|
Vecteur& d_hsurh0 = save_resul.d_hsurh0;
|
|
// b) calcul de la variation
|
|
//2.* h/h_0 * d h/h0 = - djaco_1D * (jaco_1D_0 / ((jaco_1D_tdt)^2)) * exp(log_VsurV0)
|
|
// + d_(sigBH.Trace()) * (jaco_1D_0 / jaco_1D_tdt) * exp(log_VsurV0) / (3. * module_compressibilite )
|
|
int taille = d_hsurh0.Taille();
|
|
double fact1 = - jaco_1D_0 / (jaco_1D_tdt*jaco_1D_tdt) * exp_log_VsurV0 / (2. * hSurh_0);
|
|
double fact2 = (jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0 / (3. * module_compressibilite / (2. * hSurh_0));
|
|
////---debug
|
|
//cout << "\n Calcul_d_eps33_parVarVolume() module_compressibilite= " << module_compressibilite
|
|
// << " log_VsurV0= "<< log_VsurV0
|
|
// << " sigBH.Trace()= " << sigBH.Trace()
|
|
// << " hSurh_0 =" << hSurh_0
|
|
// << "\n fact1 =" << fact1 << " fact2 =" << fact2;
|
|
//// --- fin debug
|
|
|
|
for (int i=1;i<=taille;i++)
|
|
{
|
|
Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i))); // passage en dim 1
|
|
const Tenseur1BB & dgijBB = *((Tenseur1BB*)(d_gijBB_tdt(i))) ; // pour simplifier l'ecriture
|
|
Tenseur1BH d_sigBH = gijBB * dsigHH + dgijBB * sigHH;
|
|
d_hsurh0(i) = fact1 * d_jaco_1D(i) + fact2 * (d_sigBH.Trace());
|
|
//---debug
|
|
//cout << "\n Calcul_d_eps33_parVarVolume() "
|
|
// << " d_hsurh0("<< i<<")= "<< d_hsurh0(i);
|
|
// --- fin debug
|
|
};
|
|
// en isotrope les axes h et b sont équivalents
|
|
save_resul.d_bsurb0 = save_resul.d_hsurh0;
|
|
}
|
|
else
|
|
{ save_resul.hsurh0 = save_resul.bsurb0 = 1.;
|
|
save_resul.d_hsurh0.Zero();save_resul.d_bsurb0.Zero();
|
|
};
|
|
|
|
};
|
|
|
|
// calcul de la variation des déformations d'hauteur et de largeur, par rapport à la def_11, à cause
|
|
// des contraintes planes doubles: première méthode
|
|
void LoiContraintesPlanesDouble::Calcul_d_eps_eg_11()
|
|
{ // on prépare la résolution
|
|
// ,d_sig_ef_gh(2,2),d_eps_ef_11(2)
|
|
// la matrice
|
|
double a = d_sig_ef_gh(1,1) = d_sigma_deps_3D(2,2,2,2);
|
|
double b = d_sig_ef_gh(1,2) = d_sigma_deps_3D(2,2,3,3);
|
|
double c = d_sig_ef_gh(2,1) = d_sigma_deps_3D(3,3,2,2);
|
|
double d = d_sig_ef_gh(2,2) = d_sigma_deps_3D(3,3,3,3);
|
|
double determinant = a*d - b*c;
|
|
// //----- debug
|
|
// cout << "\n debug LoiContraintesPlanesDouble::Calcul_d_eps_eg_11() "<< "\n d_sig_ef_gh= ";
|
|
// d_sig_ef_gh.Affiche();
|
|
// //---fin debug
|
|
if (Dabs(determinant) > ConstMath::petit)
|
|
{ // le second membre
|
|
d_eps_ef_11(1)=-d_sigma_deps_3D(2,2,1,1);
|
|
d_eps_ef_11(2)=-d_sigma_deps_3D(3,3,1,1);
|
|
// //----- debug
|
|
// cout << "\n second membre= "<< d_eps_ef_11(1) << ", " << d_eps_ef_11(2);
|
|
// //---fin debug
|
|
// résolution
|
|
d_sig_ef_gh.Resol_systID(d_eps_ef_11);
|
|
// //----- debug
|
|
// cout << "\n solution= "<< d_eps_ef_11(1) << ", " << d_eps_ef_11(2);
|
|
// // vérif
|
|
// a = d_sig_ef_gh(1,1) = d_sigma_deps_3D(2,2,2,2);
|
|
// b = d_sig_ef_gh(1,2) = d_sigma_deps_3D(2,2,3,3);
|
|
// c = d_sig_ef_gh(2,1) = d_sigma_deps_3D(3,3,2,2);
|
|
// d = d_sig_ef_gh(2,2) = d_sigma_deps_3D(3,3,3,3);
|
|
// Vecteur truc(2);
|
|
// truc = d_sig_ef_gh * d_eps_ef_11;
|
|
// cout << "\n second membre correspondant "<< truc(1) << ", " << truc(2)<< flush;
|
|
// //---fin debug
|
|
}
|
|
else // sinon on met à 0
|
|
{d_eps_ef_11.Inita(0.);};
|
|
};
|
|
|
|
|
|
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
|
|
// correspondant aux cas 3D: valable pour les 2 méthodes
|
|
void LoiContraintesPlanesDouble::Calcul_invariants_et_def_cumul()
|
|
{ // récup des infos spécifiques
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat=0;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // non un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "; };
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
|
|
// ---calcul des invariants de déformation et de vitesse de déformation
|
|
Tenseur3BH espBH = eps_BB_3D * gijHH_tdt_3D;
|
|
Tenseur3BH delta_espBH = delta_eps_BB_3D * gijHH_tdt_3D;
|
|
save_resul.epsInvar(1) = espBH.Trace();
|
|
save_resul.epsInvar(2) = espBH.II();
|
|
save_resul.epsInvar(3) = espBH.Det();
|
|
save_resul.depsInvar(1) = unSurDeltat * delta_espBH.Trace();
|
|
save_resul.depsInvar(2) = unSurDeltat * delta_espBH.II();
|
|
save_resul.depsInvar(3) = unSurDeltat * delta_espBH.Det();
|
|
|
|
// cas des grandeurs cumulées (cf. la classe Deformation)
|
|
|
|
// tableau relatif aux différentes grandeurs de type def scalaires équivalentes
|
|
// def_equi(1) = deformation cumulée = somme des sqrt(2./3. * (delta_eps_barre_BH && delta_eps_barre_BH)) ;
|
|
// def_equi(2) = deformation duale de la contrainte de mises = sqrt(2./3. * (eps_barre_BH && eps_barre_BH)) ;
|
|
// def_equi(3) = niveau maxi atteind par def_equi(2)
|
|
// def_equi(4) = delta def cumulée = sqrt(2./3. * (delta_eps_barre_BH && delta_eps_barre_BH));
|
|
|
|
double delta_eps_equi = sqrt(2./3. * ( (delta_espBH && delta_espBH) - Sqr(delta_espBH.Trace()) /3. ));
|
|
save_resul.def_equi(1) = save_resul.def_equi_t(1) + delta_eps_equi ;
|
|
save_resul.def_equi(2) = sqrt(2./3. * (espBH && espBH)) ;
|
|
if (save_resul.def_equi(2) > save_resul.def_equi_t(3))
|
|
save_resul.def_equi(3) = save_resul.def_equi(2);
|
|
save_resul.def_equi(4) = delta_eps_equi;
|
|
|
|
};
|
|
|
|
// calcul de la fonction résidu de la résolution de l'équation constitutive: sig^ef(h/h0,b/b0) = 0
|
|
// avec "e,f" = 2 et 3
|
|
// h/h0, bb0 sont les inconnues du problème et sont l'élément d'entrée: <x(1),x(2)>
|
|
// l'argument test ramène
|
|
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
|
|
// fatal, qui invalide le calcul du résidu
|
|
Vecteur& LoiContraintesPlanesDouble::Residu_constitutif (const double & alpha,const Vecteur & x, int& test)
|
|
{
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
test = 1; // init par défaut
|
|
|
|
// %%% on met à jour les conteneurs locaux en fonction des paramètres d'entrée %%%
|
|
// -- prise en compte de la variation de la déformation d'épaisseur
|
|
// pour nous ici on travaille avec giB_03D qui est normé et = 1 c-a-d h0 = 1
|
|
// donc à tdt il faut multiplier par h/h0 à tdt
|
|
double bsurb0 = x(1);double hsurh0 = x(2);
|
|
// vérif des grandeurs
|
|
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
|
|
if (Limitation_h_b(hsurh0,bsurb0))
|
|
{if (Permet_affichage() > 3)
|
|
cout << "\n LoiContraintesPlanesDouble::Residu_constitutif(... " << endl;
|
|
erreur_sortie = true;
|
|
test = -1; // l'erreur n'est pas fatale mais elle signifie que ce n'est pas très normal !!
|
|
};
|
|
if (erreur_sortie && (Permet_affichage() > 4))
|
|
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
};
|
|
|
|
save_resul.hsurh0 = hsurh0; // sauvegarde de l'élongation d'épaisseur
|
|
save_resul.bsurb0 = bsurb0; // sauvegarde de l'élongation de largeur
|
|
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
|
|
giB_tdt_3D.CoordoB(3) = giBtdt_3 * hsurh0; giH_tdt_3D.CoordoH(3) = (giBtdt_3 / hsurh0).Bas_haut();
|
|
giB_tdt_3D.CoordoB(2) = giBtdt_2 * bsurb0; giH_tdt_3D.CoordoH(2) = (giBtdt_2 / bsurb0).Bas_haut();
|
|
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs
|
|
// et en largeurs sont normées initialement, seule l'intensité
|
|
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs et de largeur
|
|
double hsurh0_2 = hsurh0 * hsurh0;
|
|
gijBB_tdt_3D.Coor(3,3) = hsurh0_2; gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
|
|
double bsurb0_2 = bsurb0 * bsurb0;
|
|
gijBB_tdt_3D.Coor(2,2) = bsurb0_2; gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_2;
|
|
|
|
// mise à jour des informations liées à la déformation de 1 vers 3
|
|
Mise_a_jour_deformations_et_Jacobien_en_3D();
|
|
|
|
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
|
|
// correspondant aux cas 3D
|
|
Calcul_invariants_et_def_cumul();
|
|
|
|
// initialisation du comportement tangent / au def
|
|
d_sigma_deps_3D.Inita(0.);
|
|
// appel du calcul de sig et dsig
|
|
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas forcément en orthonormee
|
|
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
|
|
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D
|
|
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
|
|
,module_cisaillement_3D,*umat_cont_3D);
|
|
// retour des grandeurs voulues
|
|
residu(1) = sig_HH_3D(2,2);
|
|
residu(2) = sig_HH_3D(3,3);
|
|
return residu;
|
|
|
|
};
|
|
|
|
// calcul de la matrice tangente de la résolution de l'équation constitutive: sig^ef(h/h0,b/b0) = 0
|
|
// avec "e,f" = 2 et 3
|
|
// h/h0, bb0 sont les inconnues du problème et sont l'élément d'entrée: <x(1),x(2)>
|
|
// l'argument test ramène
|
|
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
|
|
// fatal, qui invalide le calcul du résidu et de la dérivée
|
|
Mat_abstraite& LoiContraintesPlanesDouble::Mat_tangente_constitutif
|
|
(const double & alphap,const Vecteur & x, Vecteur& residu, int& test)
|
|
{
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
test = 1; // init par défaut
|
|
|
|
// %%% on met à jour les conteneurs locaux en fonction des paramètres d'entrée %%%
|
|
// -- prise en compte de la variation de la déformation d'épaisseur et de largeur
|
|
// pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = b0 = 1
|
|
// donc à tdt il faut multiplier par h/h0 à tdt et idem pour b/b0
|
|
double bsurb0 = x(1);double hsurh0 = x(2);
|
|
// vérif des grandeurs
|
|
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
|
|
if (Limitation_h_b(hsurh0,bsurb0))
|
|
{if (Permet_affichage() > 3)
|
|
cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif(... " << endl;
|
|
erreur_sortie = true;
|
|
test = -1; // l'erreur n'est pas fatale mais elle signifie que ce n'est pas très normal !!
|
|
};
|
|
if (erreur_sortie && (Permet_affichage() > 4))
|
|
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
};
|
|
|
|
save_resul.hsurh0 = hsurh0; // sauvegarde de l'élongation d'épaisseur
|
|
save_resul.bsurb0 = bsurb0; // sauvegarde de l'élongation de largeur
|
|
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
|
|
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
|
|
giB_tdt_3D.CoordoB(3) = giBtdt_3 * hsurh0;
|
|
giH_tdt_3D.CoordoH(3) = (giBtdt_3 / hsurh0).Bas_haut();
|
|
giB_tdt_3D.CoordoB(2) = giBtdt_2 * bsurb0;
|
|
giH_tdt_3D.CoordoH(2) = (giBtdt_2 / bsurb0).Bas_haut();
|
|
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
|
|
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
|
|
double hsurh0_2 = hsurh0 * hsurh0;
|
|
gijBB_tdt_3D.Coor(3,3) = hsurh0_2;
|
|
gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
|
|
double bsurb0_2 = bsurb0 * bsurb0;
|
|
gijBB_tdt_3D.Coor(2,2) = bsurb0_2;
|
|
gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_2;
|
|
|
|
// mise à jour des informations liées à la déformation de 1 vers 3
|
|
Mise_a_jour_deformations_et_Jacobien_en_3D();
|
|
|
|
|
|
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
|
|
// correspondant aux cas 3D
|
|
Calcul_invariants_et_def_cumul();
|
|
|
|
// initialisation du comportement tangent / au def
|
|
d_sigma_deps_3D.Inita(0.);
|
|
// appel du calcul de sig et dsig
|
|
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas a priori en orthonormee
|
|
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
|
|
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D
|
|
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
|
|
,module_cisaillement_3D,*umat_cont_3D);
|
|
|
|
// retour des grandeurs voulues
|
|
residu(1) = sig_HH_3D(2,2);
|
|
residu(2) = sig_HH_3D(3,3);
|
|
// la matrice tangente dépend du type de mesure de déformation, avec la relation
|
|
// en appelant hb/hb0 les deux ddl: h/h0 et b/b0, i,j,k et l variant de 2 à 3
|
|
// d_sig^ij(h/h0,b/b0) / d_ddl = d_sig^ij(h/h0,b/b0)/d_eps_kl * d_eps_kl/d_ddl
|
|
// comme les repères restent orthogonales à tous moments il n'y a pas de couplage entre les directions 2 et 3
|
|
Vecteur d_eps_kl_sur_d_ddl(2,0.); // init
|
|
|
|
// -- choix en fonction du type de déformation
|
|
switch (type_de_deformation)
|
|
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
|
|
// cas d'une déformation d'Almansi
|
|
{ // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
|
|
// dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
|
|
// idem pour la direction 22
|
|
d_eps_kl_sur_d_ddl(1) = bsurb0;d_eps_kl_sur_d_ddl(2) = hsurh0;
|
|
};
|
|
break;
|
|
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
|
|
// cas d'une def logarithmique ou une approximation: eps_BB_3D(3,3) = log(hsurh0);
|
|
{ d_eps_kl_sur_d_ddl(1) = log(bsurb0);d_eps_kl_sur_d_ddl(2) = log(hsurh0);
|
|
};
|
|
break;
|
|
default :
|
|
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
|
|
<< Nom_type_deformation(type_de_deformation);
|
|
cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
// d'où la dérivée du résidu, en considérant qu'il n'y a pas de couplage entre les directions 2 et 3
|
|
// c-a-d d_eps33_sur_d_bSurb0 = 0 et d_eps22_sur_d_hSurh0 = 0
|
|
// non c'est faux, du à la partie sphérique des contraintes par exemple, il faut tout considérer
|
|
derResidu(1,1) = d_sigma_deps_3D(2,2,2,2) * d_eps_kl_sur_d_ddl(1);
|
|
derResidu(1,2) = d_sigma_deps_3D(2,2,3,3) * d_eps_kl_sur_d_ddl(2);
|
|
derResidu(2,1) = d_sigma_deps_3D(3,3,2,2) * d_eps_kl_sur_d_ddl(1);
|
|
derResidu(2,2) = d_sigma_deps_3D(3,3,3,3) * d_eps_kl_sur_d_ddl(2);
|
|
return derResidu;
|
|
};
|
|
|
|
//--- cas de la résolution de l'équation sigij(eps_meca)*V_j{.l}=0, l=2 et 3
|
|
// calcul de la fonction résidu de la résolution de l'équation constitutive
|
|
// l'argument test ramène
|
|
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
|
|
// fatal, qui invalide le calcul du résidu
|
|
Vecteur& LoiContraintesPlanesDouble::Residu_constitutif_3D(const double & alpha,const Vecteur & x, int& test)
|
|
{
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
test = 1; // init par défaut
|
|
|
|
//récup de eps_mecaBB_t
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
|
|
// les inconnues sont les déformations mécaniques {epsilon'}_kl;
|
|
// (k,l) = (22) (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
|
|
// c-a-d les def qui vont servir dans l'appel de la loi 3D
|
|
eps_P_BB.Coor(1,2) = x(3);eps_P_BB.Coor(2,2) = x(1);eps_P_BB.Coor(3,3) = x(2);
|
|
eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
|
|
if (Permet_affichage() > 4)
|
|
{cout << "\n eps_P_BB (dans Va) avant limitation : " ;eps_P_BB.Ecriture(cout);}
|
|
|
|
// vérif des grandeurs et calcul des variations d'épaisseur et de largeur
|
|
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
|
|
// la méthode va mettre à jour les variations bsurb0 et hsurh0 et le jacobien
|
|
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
|
|
int info_limit = Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D);
|
|
if (info_limit)
|
|
{if (Permet_affichage() > 3)
|
|
cout << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(... " << endl;
|
|
erreur_sortie = true;test = info_limit - 2;
|
|
};
|
|
if (erreur_sortie && (Permet_affichage() > 4))
|
|
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n eps_BB_3D: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) << flush;
|
|
};
|
|
|
|
// dans la direction V_1 c'est la déformation issue de la cinématique
|
|
// elle a déjà été calculée, et les autres composantes demeurent nulle
|
|
// donc eps_P_BB est au complet,
|
|
eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
|
|
eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
|
|
// on s'occupe maintenant de l'incrément de déformation
|
|
delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
|
|
// au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
|
|
delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
#ifdef MISE_AU_POINT
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "
|
|
<< "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(...";
|
|
};
|
|
#endif
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
// ici on considère que la déformation est de type Almansi
|
|
Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
|
|
|
|
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
|
|
// correspondant aux cas 3D
|
|
Calcul_invariants_et_def_cumul();
|
|
|
|
/* // --- on s'occupe des jacobiens
|
|
// à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
|
|
// ici il s'agit du jacobien relatif au comportement mécanique
|
|
// |J| = sqrt(det(2 * eps_ij + g_ij(0)))
|
|
// calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
|
|
gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
|
|
gij_meca_BB.Matrice_composante(mat_inter);
|
|
// calcul du jacobien
|
|
jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
|
|
*/
|
|
// initialisation du comportement tangent / au def
|
|
d_sigma_deps_3D.Inita(0.);
|
|
// appel du calcul de sig et dsig
|
|
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas en orthonormee
|
|
// car on travaille dans la base naturelle
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(..."
|
|
<< "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) ;
|
|
// en base orthonormee
|
|
Tenseur3BB epsBB_3D_inter;
|
|
eps_BB_3D.BaseAbsolue(epsBB_3D_inter,*(umat_cont_3D->giH_tdt));
|
|
cout << "\n eps_BB_3D en absolue 3D (dans Va): "<<epsBB_3D_inter;
|
|
cout << "\n jacobien_tdt_3D= "<< jacobien_tdt_3D
|
|
<< "\n jacobien_0= " << jacobien_0_3D
|
|
<< flush;
|
|
};
|
|
#endif
|
|
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
|
|
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D
|
|
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
|
|
,module_cisaillement_3D,*umat_cont_3D);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n sig_HH_3D: " << sig_HH_3D;
|
|
cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
|
|
};
|
|
#endif
|
|
|
|
//-------- calcul du residu -----------
|
|
// residu(mi) = V_{mj} * \sigma^{ij}(\epsilon'_{kl}) : m=2 et 3
|
|
|
|
CoordonneeH V_sig_H(sig_HH_3D * (*ViB_3D)(2));
|
|
// on considère que les densités d'effort interne dans le plan, la composante 23 devant
|
|
// être nulle compte tenu du tenseur des déformations
|
|
// residu_3D(1) = V_sig_H(1);residu_3D(3) = V_sig_H(2);
|
|
residu_3D(1) = V_sig_H * (*ViB_3D)(2); // composante 22
|
|
residu_3D(3) = (*ViB_3D)(1) * V_sig_H ; // composante 12
|
|
// essai avec le mixte
|
|
// residu_3D(3) = 0.5 * ((*ViB_3D)(1) * V_sig_H + V_sig_H * (*ViB_3D)(1)); // composante 12
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n en g_i: sig_HH * ViB(2)= " << V_sig_H
|
|
<< " df_ds_22= " << residu_3D(1) << " df_ds_21= " << residu_3D(3)
|
|
<< flush;
|
|
};
|
|
#endif
|
|
|
|
V_sig_H = sig_HH_3D * (*ViB_3D)(3);
|
|
// residu_3D(2) = V_sig_H(3); // seule la composante dans la direction 3 est utilisée
|
|
residu_3D(2) = V_sig_H * (*ViB_3D)(3); // composante 33
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n en g_i: sig_HH * ViB(3)= " << V_sig_H << " df_ds_33= " << residu_3D(2);
|
|
// pour info on calcul la projection symétrique
|
|
double f_12 = (*ViB_3D)(1)* sig_HH_3D * (*ViB_3D)(2);
|
|
cout << " df_ds_12= " << f_12 << flush;
|
|
};
|
|
#endif
|
|
|
|
return residu_3D;
|
|
|
|
};
|
|
//--- cas de la résolution de l'équation sigij(eps_meca)*V_j{.l}=0, l=2 et 3
|
|
// calcul de la matrice tangente de la résolution de l'équation constitutive
|
|
// l'argument test ramène
|
|
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
|
|
// fatal, qui invalide le calcul du résidu et de la dérivée
|
|
Mat_abstraite& LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D
|
|
(const double & alphap,const Vecteur & x, Vecteur& residu_local, int& test)
|
|
{
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
test = 1; // init par défaut
|
|
|
|
//récup de eps_mecaBB_t
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
|
|
// les inconnues sont les déformations mécaniques {epsilon'}_kl;
|
|
// (k,l) = (22) (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
|
|
// c-a-d les def qui vont servir dans l'appel de la loi 3D
|
|
eps_P_BB.Coor(1,2) = x(3);eps_P_BB.Coor(2,2) = x(1);eps_P_BB.Coor(3,3) = x(2);
|
|
eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
|
|
|
|
// vérif des grandeurs et calcul des variations d'épaisseur et de largeur
|
|
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
|
|
// la méthode va mettre à jour les valeurs de bsurb0 et hsurh0, et le jacobien
|
|
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
|
|
if (Permet_affichage() > 4)
|
|
{cout << "\n eps_P_BB (dans Va) avant limitation : " ;eps_P_BB.Ecriture(cout);};
|
|
int info_limit = Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D);
|
|
if (info_limit)
|
|
{if (Permet_affichage() > 3)
|
|
cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(... " << endl;
|
|
erreur_sortie = true;test = info_limit - 2;
|
|
};
|
|
if (erreur_sortie && (Permet_affichage() > 4))
|
|
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) << flush;
|
|
};
|
|
|
|
// dans la direction V_1 c'est la déformation issue de la cinématique
|
|
// elle a déjà été calculée, et les autres composantes demeurent nulle
|
|
// donc eps_P_BB est au complet,
|
|
eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
|
|
eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
|
|
// on s'occupe maintenant de l'incrément de déformation
|
|
delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
|
|
// au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
|
|
delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
#ifdef MISE_AU_POINT
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "
|
|
<< "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(...";
|
|
};
|
|
#endif
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
// ici on considère que la déformation est de type Almansi
|
|
Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
|
|
|
|
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
|
|
// correspondant aux cas 3D
|
|
Calcul_invariants_et_def_cumul();
|
|
|
|
/* // --- on s'occupe des jacobiens
|
|
// à t=0 on considère que c'est le même jacobien que transmis
|
|
jacobien_0_3D = *(umat_cont_3D->jacobien_0);
|
|
// à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
|
|
// ici il s'agit du jacobien relatif au comportement mécanique
|
|
// |J| = sqrt(det(2 * eps_ij + g_ij(0)))
|
|
// calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
|
|
gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
|
|
gij_meca_BB.Matrice_composante(mat_inter);
|
|
// calcul du jacobien
|
|
jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
|
|
*/
|
|
// initialisation du comportement tangent / au def
|
|
d_sigma_deps_3D.Inita(0.);
|
|
// appel du calcul de sig et dsig
|
|
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas a priori en orthonormee
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(..."
|
|
<< "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) ;
|
|
// en base orthonormee
|
|
Tenseur3BB epsBB_3D_inter;
|
|
eps_BB_3D.BaseAbsolue(epsBB_3D_inter,*(umat_cont_3D->giH_tdt));
|
|
cout << "\n eps_BB_3D en absolue 3D (dans Va): "<<epsBB_3D_inter;
|
|
cout << "\n jacobien_tdt_3D= "<< jacobien_tdt_3D
|
|
<< "\n jacobien_0= " << jacobien_0_3D
|
|
<< flush;
|
|
};
|
|
#endif
|
|
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
|
|
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D
|
|
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
|
|
,module_cisaillement_3D,*umat_cont_3D);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n sig_HH_3D: " << sig_HH_3D;
|
|
cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
|
|
};
|
|
#endif
|
|
|
|
//-------- calcul du residu -----------
|
|
// residu(mi) = V_{mj} * \sigma^{ij}(\epsilon'_{kl}) : m=2 et 3
|
|
|
|
//residu(1) = V_2 . \sigma . V_2 = V_{2i} ~\sigma^{ij} ~V_{2j} = ( V_2 \otimes V_2) : \sigma
|
|
//residu(2) = V_3 . \sigma . V_3 = V_{3i} ~\sigma^{ij} ~V_{3j} = ( V_3 \otimes V_3) : \sigma
|
|
//residu(3) = V_1 . \sigma . V_2 = V_{1i} ~\sigma^{ij} ~V_{2j} = ( V_1 \otimes V_2) : \sigma
|
|
|
|
|
|
CoordonneeH V_sig_H(sig_HH_3D * (*ViB_3D)(2)); // = \sigma . V_2
|
|
// ici on utilise le vecteur passé en argument
|
|
residu_local(1) = V_sig_H * (*ViB_3D)(2); // composante 22 = V_2 . \sigma . V_2
|
|
// residu(3) = V_sig_H * (*ViB_3D)(1); // composante 21
|
|
residu_local(3) = (*ViB_3D)(1) * V_sig_H ; // composante 12 = V_1 . \sigma . V_2
|
|
// essai avec le mixte
|
|
// residu_local(3) = 0.5 * ((*ViB_3D)(1) * V_sig_H + V_sig_H * (*ViB_3D)(1)); // composante 12
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n en g_i: sig_HH * ViB(2)= " << V_sig_H
|
|
<< " df_ds_22= " << residu_3D(1) << " df_ds_21= " << residu_3D(3)
|
|
<< flush;
|
|
};
|
|
#endif
|
|
|
|
V_sig_H = sig_HH_3D * (*ViB_3D)(3); // = \sigma . V_3
|
|
residu_local(2) = V_sig_H * (*ViB_3D)(3); // composante 33 = V_3 . \sigma . V_3
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n en g_i: sig_HH * ViB(3)= " << V_sig_H << " df_ds_33= " << residu_3D(2);
|
|
// pour info on calcul la projection symétrique
|
|
double f_12 = (*ViB_3D)(1)* sig_HH_3D * (*ViB_3D)(2);
|
|
cout << " df_ds_12= " << f_12 << flush;
|
|
};
|
|
#endif
|
|
|
|
// calcul de la matrice tangente
|
|
// formule générique
|
|
//\frac{\partial \left( \sigma^{ij}(\epsilon'_{kl}) ~ V_{mj} \right )}{\partial \epsilon'_{rs}}
|
|
// = V_{mj}~\frac{\partial \sigma^{ij}(\epsilon_{kl}) }{\partial \epsilon_{op}}~V^r_{~.o}~ V^s_{~.p}
|
|
// ici on doit utiliser uniquement:
|
|
//\frac{\partial \left( \sigma^{3j}(\epsilon'_{kl}) ~ V_{3j} \right ) }{\partial \epsilon'_{rs}}
|
|
//\frac{\partial \left( \sigma^{\alpha j}(\epsilon'_{kl}) ~ V_{2j} \right ) }{\partial \epsilon'_{rs}}
|
|
// avec \alpha =1,2 ; (k,l) et (r,s) } = (2,2) ; (3,3) ; (1,2)
|
|
|
|
// d (V_{mi} ~\sigma^{ij} ~V_{nj} \right ) / d \epsilon'_{gh}
|
|
// = V_{mi}~V_{nj}~ d \sigma^{ij} / d \epsilon_{op} ~V^r_{~.o}~ V^s_{~.p}
|
|
// pour (m,n) et (g,h) = (2,2); (3,3); et (1,2)
|
|
|
|
|
|
derResidu_3D.Initialise(0.); // init
|
|
// for (int j=1;j<4;j++) // j représente la composante de l'axe de projection
|
|
// { for (int i=1;i<4;i++) // i représente le classement de (\epsilon'_{kl})
|
|
// // i= 1 pour (k,l) = (2,2), i=2 pour (k,l) = (3,3), i=3 pour (k,l) = (1,2)
|
|
// {// d (sig(2,j) V(2)(j)) / d (eps')(k,l)
|
|
// derResidu_3D(1,i) += (*ViB_3D)(2)(j) *
|
|
// (d_sigma_deps_3D(2,j,1,1) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(1)
|
|
// + d_sigma_deps_3D(2,j,1,2) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(2)
|
|
// + d_sigma_deps_3D(2,j,2,1) * (*ViH_3D)(indx(i))(2) * (*ViH_3D)(indy(i))(1)
|
|
// + d_sigma_deps_3D(2,j,3,3) * (*ViH_3D)(indx(i))(3) * (*ViH_3D)(indy(i))(3)
|
|
// );
|
|
// // d (sig(1,j) V(2)(j)) / d (eps')(k,l)
|
|
// derResidu_3D(3,i) += (*ViB_3D)(2)(j) *
|
|
// (d_sigma_deps_3D(1,j,1,1) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(1)
|
|
// + d_sigma_deps_3D(1,j,1,2) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(2)
|
|
// + d_sigma_deps_3D(1,j,2,1) * (*ViH_3D)(indx(i))(2) * (*ViH_3D)(indy(i))(1)
|
|
// + d_sigma_deps_3D(1,j,3,3) * (*ViH_3D)(indx(i))(3) * (*ViH_3D)(indy(i))(3)
|
|
// );
|
|
// // d (sig(3,j) V(2)(j)) / d (eps')(k,l)
|
|
// derResidu_3D(2,i) += (*ViB_3D)(3)(j) *
|
|
// (d_sigma_deps_3D(3,j,1,1) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(1)
|
|
// + d_sigma_deps_3D(3,j,1,2) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(2)
|
|
// + d_sigma_deps_3D(3,j,2,1) * (*ViH_3D)(indx(i))(2) * (*ViH_3D)(indy(i))(1)
|
|
// + d_sigma_deps_3D(3,j,3,3) * (*ViH_3D)(indx(i))(3) * (*ViH_3D)(indy(i))(3)
|
|
// );
|
|
//
|
|
// };
|
|
// };
|
|
|
|
// pour (m,n) et (g,h) = (2,2); (3,3); et (1,2)
|
|
// d (V_{mi} ~\sigma^{ij} ~V_{nj} ) / d\epsilon'_{gh} =
|
|
// (V_{mi} V_{nj}) (d \sigma^{ij} / d \epsilon_{op}) (V^g_{~.o} V^h_{~.p})
|
|
// pour 2,2 et 3,3 pas de pb de symétrie
|
|
// pour 1,2
|
|
|
|
|
|
for (int mn=1;mn<4;mn++) // mn représente la composante de l'axe de projection
|
|
{// d'abord les termes diagonaux
|
|
for (int gh=1;gh<3;gh++) // gh représente le classement de (\epsilon'_{gh})
|
|
derResidu_3D(mn,gh) = VhVg_BB(mn) && (d_sigma_deps_3D && VhVg_BB(gh));
|
|
// maintenant les termes non diagonaux à doubler
|
|
derResidu_3D(mn,3) = 2.* VhVg_BB(mn) && (d_sigma_deps_3D && VhVg_BB(3));
|
|
};
|
|
|
|
return derResidu_3D;
|
|
};
|
|
|
|
|
|
// limitation des variations d'épaisseurs et largeurs
|
|
// ramène true s'il y a eu une modif: sert pour la méthode 1
|
|
bool LoiContraintesPlanesDouble::Limitation_h_b(double& hsurh0, double& bsurb0)
|
|
{ bool retour = false;
|
|
// il ne faut pas que la nouvelle épaisseur soit inf au mini
|
|
if (hsurh0 < mini_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est inferieur a la limite fixee "<<mini_hsurh0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
hsurh0 = mini_hsurh0;retour = true;
|
|
};
|
|
// il ne faut pas que la nouvelle largeur soit inf au mini
|
|
if (bsurb0 < mini_bsurb0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est inferieur a la limite fixee "<<mini_bsurb0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
bsurb0 = mini_bsurb0;retour = true;
|
|
};
|
|
// on regarde également pour les maxi
|
|
if (hsurh0 > maxi_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est superieur a la limite fixee "<<maxi_hsurh0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
hsurh0 = maxi_hsurh0;retour = true;
|
|
};
|
|
if (bsurb0 > maxi_bsurb0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est superieur a la limite fixee "<<maxi_bsurb0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
bsurb0 = maxi_bsurb0;retour = true;
|
|
};
|
|
// retour
|
|
return retour;
|
|
};
|
|
|
|
|
|
// calcul de l'opérateur tangent (2ième méthode)
|
|
// d sigma^{ij} / d epsilon_{kl}
|
|
void LoiContraintesPlanesDouble::Calcul2_dsigma_deps(TenseurHHHH& d_sigma_deps)
|
|
{
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n dans le repere locale d_sigma_deps_3D= \n";
|
|
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<<")= "<<d_sigma_deps_3D(i,j,k,l) << " ; ";
|
|
if (e>6) {cout << "\n"; e=1;}
|
|
};
|
|
Tenseur3HHHH inter_HHHH;
|
|
d_sigma_deps_3D.ChangeBase(inter_HHHH,giB_tdt_3D);
|
|
cout << "\n dans le repere orthonormee d_sigma_deps_3D= \n";
|
|
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<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
|
|
if (e>6) {cout << "\n"; e=1;}
|
|
};
|
|
};
|
|
#endif
|
|
// on calcule deux grandeurs intermédiaires
|
|
// a) d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}
|
|
// on récupère que le pointeur, ce qui évite d'utiliser l'opérateur de copie du tenseur
|
|
// !!!!**** a vérifier -> oui c'est ok
|
|
// const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
|
|
|
|
Tenseur3HH dSigDeps11P_HH(d_sigma_deps_3D && V1V1_BB);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
{ cout << "\n (*ViB_3D)(1)= "<<(*ViB_3D)(1);
|
|
cout << "\n V1V1_BB: "<< V1V1_BB
|
|
<< "\n dSigDeps11P_HH: "<< dSigDeps11P_HH
|
|
<< flush;
|
|
};
|
|
#endif
|
|
|
|
// b) d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^g_{~.o'}~ V^h_{~.p'}
|
|
Tableau <Tenseur3HH* > dsigdeps_gh_P(3);
|
|
for (int gh=1;gh<4;gh++)
|
|
{dsigdeps_gh_P(gh) = ((Tenseur3HH*) &( d_sigma_deps_3D && VhVg_BB(gh)));
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
{ cout << "\n VhVg_BB("<<gh<<")= "<< VhVg_BB(gh);
|
|
cout << "\n dsigdeps_gh_P(gh): "<< (*dsigdeps_gh_P(gh)) << flush;
|
|
}
|
|
#endif
|
|
};
|
|
|
|
// c) mat1(mn,gh) = [ V_{mi}~(d sigma^{ij} / d epsilon_{op}(3D)~~V^g_{~.o}~ V^h_{~.p}) ~~V_{nj}]
|
|
// SM(mn) = ( V_{mi}~~ d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}~~V_{nj} )
|
|
Mat_pleine mat1(3,3);
|
|
Vecteur SM(3);
|
|
for (int mn=1;mn<4;mn++)
|
|
{SM(mn) = VhVg_BB(mn) && dSigDeps11P_HH;
|
|
for (int gh=1;gh<4;gh++)
|
|
mat1(mn,gh) = VhVg_BB(mn) && (*dsigdeps_gh_P(gh));
|
|
};
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
{ cout << "\n mat1 "<< mat1;
|
|
cout << "\n SM: "<< SM << flush;
|
|
}
|
|
#endif
|
|
|
|
// d) calcul des variations des déformations transversales
|
|
// en fonction de la déformation longitudinale
|
|
// d epsilon'_gh / d epsilon'_11 =
|
|
// - [ V_{mi}~(d sigma^{ij} / d epsilon_{op}(3D)~~V^g_{~.o}~ V^h_{~.p}) ~~V_{nj}]^{-1}~
|
|
// * ( V_{mi}~~ d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}~~V_{nj} )
|
|
Vecteur depsP_gh_depsP_11_BB(3);
|
|
depsP_gh_depsP_11_BB = - mat1.Inverse() * SM;
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
{ cout << "\n depsP_gh_depsP_11_BB "<< depsP_gh_depsP_11_BB << flush;
|
|
}
|
|
#endif
|
|
|
|
// e) d sigma'^{11} / d epsilon'_11 (1D~CP) =
|
|
// = V_{1i}~V_{1j}~d sigma^{ij} / d epsilon_{op}(3D)~~V^1_{~.o}~ V^1_{~.p}
|
|
// + V_{1i}~V_{1j}~d sigma^{ij} / d epsilon_{op}(3D)~~~V^g_{~.o}~ V^h_{~.p}
|
|
// * d epsilon'_gh / d epsilon'_11
|
|
// ici on se sert du fait que l'on a déjà calculé:
|
|
// dSigDeps11P_HH = d sigma^{ij} / d epsilon_{op}(3D)~~V^1_{~.o}~ V^1_{~.p}
|
|
|
|
Tenseur1HHHH d_sigP11_depsP11_HHHH
|
|
( (V1V1_BB && dSigDeps11P_HH)
|
|
+ (V1V1_BB &&
|
|
( (*dsigdeps_gh_P(1))*depsP_gh_depsP_11_BB(1)
|
|
+ (*dsigdeps_gh_P(2))*depsP_gh_depsP_11_BB(2)
|
|
+ (*dsigdeps_gh_P(3))*depsP_gh_depsP_11_BB(3)
|
|
)
|
|
)
|
|
);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
cout << "\n d_sigP11_depsP11_HHHH "<< d_sigP11_depsP11_HHHH << flush;
|
|
#endif
|
|
// f) d_sigma^{ij} / d_epsilon_{kl} \approx
|
|
// V_1^{~.i}~V_1^{~.j}~d sigma'^{11} / d epsilon'_11 (1D~CP)~\gamma^{k}_{~.1} ~\gamma^{l}_{~.1}
|
|
//fonctions static définissant le produit tensoriel de deux vecteurs
|
|
// si les vecteurs sont égaux le tenseur est symétrique sinon il est non symétrique
|
|
TenseurHH & inter1_HH = Tenseur3HH::Prod_tensoriel((*ViH_3D)(1),(*ViH_3D)(1));
|
|
Coordonnee3H u_H(gamma3D.Colonne(1)); // qui correspond aux \gamma^{k}_{~.1}
|
|
TenseurHH & inter2_HH = Tenseur3HH::Prod_tensoriel(u_H,u_H);
|
|
d_sigma_deps = Tenseur3HHHH::Prod_tensoriel(inter1_HH,inter2_HH);
|
|
d_sigma_deps *= d_sigP11_depsP11_HHHH(1,1,1,1);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ Tenseur3HHHH inter_HHHH;
|
|
cout << "\n dans le repere locale d_sigma_deps= \n";
|
|
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<<")= "<<d_sigma_deps(i,j,k,l) << " ; ";
|
|
if (e>6) {cout << "\n"; e=1;}
|
|
};
|
|
Tenseur3HHHH inter1_HHHH(d_sigma_deps);
|
|
inter1_HHHH.ChangeBase(inter_HHHH,giB_tdt_3D);
|
|
cout << "\n dans le repere orthonormee d_sigma_deps= \n";
|
|
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<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
|
|
if (e>6) {cout << "\n"; e=1;}
|
|
};
|
|
};
|
|
#endif
|
|
|
|
};
|
|
|
|
|
|
// vérif des grandeurs et calcul des valeurs de l'épaisseur et de la largeur
|
|
// ainsi que calcul du jacobien final
|
|
// = 0 pas de modif
|
|
// = 1 modif due aux maxi permis de h et b, = 2 si def d'entrée incorrecte
|
|
// cas de la 2ième méthode
|
|
int LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b
|
|
(double& hsurh0,Tenseur3BB& eps_P_BB_3D, double& bsurb0,double& jacobien_tdt_3DD, const double & jacobien_0_3DD)
|
|
{ int retour = 0; // init par défaut
|
|
bool modif_hsurh0 = false; bool modif_bsurb0 = false; // init
|
|
|
|
// à partir de la valeur des déformations en cours on regarde les variations transversales correspondantes
|
|
// ici les déformations eps_BB_3D sont exprimées dans le repère V_a donc orthonormée
|
|
switch (type_de_deformation)
|
|
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
|
|
// cas d'une déformation d'Almansi
|
|
{ // on est dans le repère V_a, V_3 est direction principale du tenseur de def -> formule classique
|
|
// epsBB33 = 1/2 * (1. -(h_0/h_tdt)**2) avec un maxi de o.5
|
|
if (eps_P_BB_3D(3,3) > 0.5)
|
|
{if (Permet_affichage() > 0)
|
|
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
|
|
<< "\n **erreur: la deformation d'epaisseur dans V_a "<<eps_P_BB_3D(3,3)
|
|
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
|
|
<< " la nouvelle epaisseur, on garde donc l'ancienne " << flush;
|
|
retour = 2; // là le pb est sérieux
|
|
}
|
|
else // sinon cas normal
|
|
{hsurh0 = 1./sqrt(1.-2.*eps_P_BB_3D(3,3));
|
|
// il ne faut pas que la nouvelle épaisseur soit inf au mini
|
|
if (hsurh0 < mini_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est inferieur a la limite fixee "<<mini_hsurh0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
hsurh0 = mini_hsurh0; modif_hsurh0 = true;
|
|
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
|
|
};
|
|
// on regarde également pour les maxi
|
|
if (hsurh0 > maxi_hsurh0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est superieur a la limite fixee "<<maxi_hsurh0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
hsurh0 = maxi_hsurh0;modif_hsurh0 = true;
|
|
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
|
|
};
|
|
// si on a modifié via les extrèmes il faut mettre à jour la def
|
|
if (modif_hsurh0)
|
|
eps_BB_3D.Coor(3,3) = 0.5 * (1.-hsurh0 * hsurh0);
|
|
};
|
|
// -- suite seulement si retour est différent de 2 ---
|
|
// c-a-d que l'on n'a pas eu de pb sérieux sur h
|
|
if (retour < 2)
|
|
{// pour les deux autres def on n'est pas en directions principales mais c'est un plan principal
|
|
// on calcule les def principales via la technique de Mohr (possible car on est en repère ortho)
|
|
double epsP_c = 0.5 * (eps_P_BB_3D(1,1) + eps_P_BB_3D(2,2)); // le centre
|
|
double R = sqrt(Sqr(eps_P_BB_3D(2,2) - epsP_c) + Sqr(eps_P_BB_3D(1,2)));
|
|
double epsP_I = epsP_c + R; // la plus grande
|
|
double epsP_II = epsP_c - R; // la deuxième mais qui en toute rigueur pourrait être intercalé avec eps_P_BB_3D(3,3)
|
|
// d'où la variation de jacobien à condition de ne pas avoir de pb
|
|
bool pb = false;
|
|
// if ((epsP_I > 0.5) || (epsP_II > 0.5))
|
|
if (epsP_I > 0.5)
|
|
{if (Permet_affichage() > 0)
|
|
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
|
|
<< "\n **erreur: la deformation eps'_I "<<epsP_I
|
|
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
|
|
<< " la variation de volume " << flush;
|
|
epsP_I = 0.5-ConstMath::trespetit; // on le met au maxi
|
|
retour = 2; // là le pb est sérieux
|
|
pb = true;
|
|
};
|
|
if (epsP_II > 0.5)
|
|
{if (Permet_affichage() > 0)
|
|
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
|
|
<< "\n **erreur: la deformation eps'_II "<<epsP_I
|
|
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
|
|
<< " la variation de volume " << flush;
|
|
epsP_II = 0.5-ConstMath::trespetit; // on le met au maxi
|
|
retour = 2; // là le pb est sérieux
|
|
pb = true;
|
|
};
|
|
// on considère que la direction tranversale est celle correspondante à celle différente de
|
|
// eps_P_BB_3D(1,1) donc on en déduit une variation de b mais qui ne sera celle de V_2 que
|
|
// dans un cas isotrope, sinon c'est celle de la direction principale de def dans le plan !
|
|
if (retour < 2) // là encore on ne modifie que si retour est < 2
|
|
// car si retour == 2 on va obtenir des bsurb0 infinie
|
|
|
|
// eps_P_BB_3D(1,1), c'est la déformation cinématique imposée selon V1
|
|
// si elle est > eps_P_BB_3D(2,2), cela veut dire que epsP_I est de son côté
|
|
// par contre si c'est < , cela veut dire que c'est epsP_II qui est de son côté
|
|
// d'où le choix qui suit
|
|
{ if (eps_P_BB_3D(1,1) > eps_P_BB_3D(2,2))
|
|
{bsurb0 = 1./sqrt(1.-2.*epsP_II);}
|
|
else // cas d'une compression
|
|
{bsurb0 = 1./sqrt(1.-2.*epsP_I);};
|
|
};
|
|
|
|
// il ne faut pas que la nouvelle largeur soit inf au mini
|
|
if (bsurb0 < mini_bsurb0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est inferieur a la limite fixee "<<mini_bsurb0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
bsurb0 = mini_bsurb0; modif_bsurb0 = true;
|
|
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
|
|
};
|
|
if (bsurb0 > maxi_bsurb0)
|
|
{if (Permet_affichage() > 3)
|
|
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est superieur a la limite fixee "<<maxi_bsurb0
|
|
<< " on le ramene a la limite ";
|
|
};
|
|
bsurb0 = maxi_bsurb0;modif_bsurb0 = true;
|
|
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
|
|
};
|
|
// si on a modifié via les extrèmes il faut mettre à jour la def
|
|
// ici il s'agit uniquement de limiter eps_BB_3D(2,2) à des valeurs non absurdes donc on
|
|
// fait l'approximation (fausse en ortho mais cela ne devrait pas être un gros pb)
|
|
if (modif_bsurb0)
|
|
{eps_BB_3D.Coor(2,2) = 0.5 * (1.-bsurb0 * bsurb0);
|
|
// on est obligé de refaire une passe sur les def principales qui vont ^tre utilisées
|
|
// pour la variation de jacobien
|
|
epsP_c = 0.5 * (eps_P_BB_3D(1,1) + eps_P_BB_3D(2,2)); // le centre
|
|
R = sqrt(Sqr(eps_P_BB_3D(2,2) - epsP_c) + Sqr(eps_P_BB_3D(1,2)));
|
|
epsP_I = epsP_c + R; // la plus grande
|
|
epsP_II = epsP_c - R; // la deuxième mais qui en toute rigueur pourrait être intercalé avec eps_P_BB_3D(3,3)
|
|
// d'où la variation de jacobien à condition de ne pas avoir de pb
|
|
if (epsP_I > 0.5)
|
|
{if (Permet_affichage() > 0)
|
|
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
|
|
<< "\n **erreur: la deformation eps'_I (deuxieme passe!!) "<<epsP_I
|
|
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
|
|
<< " la variation de volume " << flush;
|
|
epsP_I = 0.5-ConstMath::trespetit; // on le met au maxi
|
|
retour = 2; // là le pb est sérieux
|
|
pb = true;
|
|
};
|
|
if (epsP_II > 0.5)
|
|
{if (Permet_affichage() > 0)
|
|
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
|
|
<< "\n **erreur: la deformation eps'_II (deuxieme passe!!) "<<epsP_I
|
|
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
|
|
<< " la variation de volume " << flush;
|
|
epsP_II = 0.5-ConstMath::trespetit; // on le met au maxi
|
|
retour = 2; // là le pb est sérieux
|
|
pb = true;
|
|
};
|
|
};
|
|
|
|
if (!pb) // cas normal
|
|
{double var_vol = 1./sqrt((1-2.*eps_P_BB_3D(3,3))*(1-2.*epsP_I)*(1-2.*epsP_II));
|
|
// jacobien_0_3DD contient uniquement le jacobien_0_1D et on part par défaut
|
|
// avec ici h0 = b0 = 1, c-a-d ici dans la loi de comp: on ne calcule en fait que h/h0 et b/b0
|
|
// la véritable épaisseur et largeur est obtenue au niveau de l'élément fini
|
|
jacobien_tdt_3DD = jacobien_0_3DD * var_vol;
|
|
};
|
|
// sinon on laise en l'état
|
|
}; // fin du cas retour < 2
|
|
};
|
|
break;
|
|
// -- pour l'instant je laisse le log en attente
|
|
// case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
|
|
// // cas d'une def logarithmique ou une approximation
|
|
// { hsurh0 = exp(eps_BB_3D(3,3));
|
|
// bsurb0 = exp(eps_BB_3D(2,2));
|
|
// };
|
|
// break;
|
|
default :
|
|
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
|
|
<< Nom_type_deformation(type_de_deformation);
|
|
cout << "\n LoiContraintesPlanesDouble::Limitation2_h_b(.. \n";
|
|
Sortie(1);
|
|
};
|
|
//---- debug
|
|
//jacobien_tdt_3DD = jacobien_0_3DD; // pour voir si c'est un pb de jacob ien
|
|
// retour
|
|
return retour;
|
|
};
|
|
|
|
|
|
// calcul de l'état final, dans le cas de la méthode de perturbation: version 2
|
|
void LoiContraintesPlanesDouble::Cal_avec_perturbation2()
|
|
{
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
|
|
//récup de eps_mecaBB_t
|
|
// 1) dans le repère gi
|
|
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
|
|
// 2) dans le dernier repère V_a: a utiliser en sortie post, pas pour le calcul !!
|
|
Tenseur3BB& eps_P_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_P_mecaBB_t));
|
|
Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_resul.eps_P_mecaBB));
|
|
|
|
// les inconnues sont les déformations mécaniques {epsilon'}_kl;
|
|
// (k,l) = (22) (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
|
|
// c-a-d les def qui vont servir dans l'appel de la loi 3D
|
|
// contrairement au cas de la méthode de Newton, ici on utilise les grandeurs stockées à t
|
|
// sauf pour la cinématique
|
|
// --- initialisation des déformations mécaniques dans le repère V_a
|
|
eps_P_BB = eps_mecaBB_t; // là on est dans le repère g^i
|
|
// on passe dans le repère V_a qui est donc le nouveau repère, pour ce passage
|
|
eps_P_BB.ChBase(betaP_3D);
|
|
|
|
// on calcul la déformation cinématique dans la direction V_1
|
|
// \epsilon'_{11}(t+\Delta t) = \epsilon_{kl}(t+\Delta t) ~\gamma^k_{~.1} ~\gamma^l_{~.1}
|
|
// idem pour l'accroissement de déformation et la vitesse de déformation
|
|
eps_P_BB.Coor(1,1) = 0.; Deps_P_BB.Coor(1,1) = 0.; delta_eps_P_BB.Coor(1,1) = 0.;
|
|
eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
|
|
for (int k=1;k<4;k++) // calcul de la matrice gamma
|
|
for (int l=1;l<4;l++)
|
|
{eps_P_BB.Coor(1,1) += (*epsBB_tdt_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
|
|
Deps_P_BB.Coor(1,1) += (*DepsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
|
|
delta_eps_P_BB.Coor(1,1) += (*delta_epsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
|
|
};
|
|
|
|
// vérif des grandeurs et calcul de la nouvelle épaisseur et largeur
|
|
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
|
|
// la méthode va mettre à jour les valeurs bsurb0 et hsurh0 et du jacobien
|
|
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
|
|
if (Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D))
|
|
{if (Permet_affichage() > 3)
|
|
cout << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(... " << endl;
|
|
erreur_sortie = true;
|
|
};
|
|
if (erreur_sortie && (Permet_affichage() > 4))
|
|
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
};
|
|
|
|
// dans la direction V_1 c'est la déformation issue de la cinématique
|
|
// elle a déjà été calculée, et les autres composantes demeurent nulle
|
|
// donc eps_P_BB est au complet,
|
|
eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
|
|
eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
|
|
// on s'occupe maintenant de l'incrément de déformation
|
|
delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
|
|
// au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
|
|
delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
|
|
// recup de l'incrément de temps
|
|
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
|
|
double unSurDeltat;
|
|
if (Abs(deltat) >= ConstMath::trespetit)
|
|
{unSurDeltat = 1./deltat;}
|
|
else
|
|
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
|
|
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
|
|
#ifdef MISE_AU_POINT
|
|
if (unSurDeltat < 0)
|
|
{ cout << "\n le pas de temps est négatif !! "
|
|
<< "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(...";
|
|
};
|
|
#endif
|
|
unSurDeltat = ConstMath::tresgrand;
|
|
};
|
|
// ici on considère que la déformation est de type Almansi
|
|
Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
|
|
|
|
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
|
|
// correspondant aux cas 3D
|
|
Calcul_invariants_et_def_cumul();
|
|
|
|
// --- on s'occupe des jacobiens
|
|
// à t=0 on considère que c'est le même jacobien que transmis
|
|
jacobien_0_3D = *(umat_cont_3D->jacobien_0);
|
|
// à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
|
|
// ici il s'agit du jacobien relatif au comportement mécanique
|
|
// |J| = sqrt(det(2 * eps_ij + g_ij(0)))
|
|
// calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
|
|
gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
|
|
gij_meca_BB.Matrice_composante(mat_inter);
|
|
// calcul du jacobien
|
|
jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
|
|
|
|
// initialisation du comportement tangent / au def
|
|
d_sigma_deps_3D.Inita(0.);
|
|
// appel du calcul de sig et dsig
|
|
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas en orthonormee
|
|
// car on travaille dans la base naturelle
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(..."
|
|
<< "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
|
|
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
|
|
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
|
|
};
|
|
#endif
|
|
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
|
|
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D
|
|
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
|
|
,module_cisaillement_3D,*umat_cont_3D);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{cout << "\n sig_HH_3D: " << sig_HH_3D;
|
|
cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
|
|
};
|
|
#endif
|
|
// sauvegarde des résultats
|
|
save_resul.def_P(1) = eps_P_BB(2,2); // récup de la solution dans le repère V_a
|
|
save_resul.def_P(2) = eps_P_BB(3,3); // "
|
|
save_resul.def_P(3) = eps_P_BB(1,2); // "
|
|
// on doit sauvegarder les def méca dans le repère g^i, ce sont directement
|
|
// les dernières composantes utilisées pour l'appel de la loi 3D
|
|
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
|
|
eps_mecaBB = eps_BB_3D; // en g^i
|
|
eps_P_mecaBB = eps_P_BB; // dans V_a pour le post traitement
|
|
};
|
|
|
|
|
|
// passage entre le calcul classique en 1D de Calcul_dsigma_deps et le calcul 3D
|
|
// méthode 2
|
|
void LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_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 &
|
|
,double& module_compressibilite,double& module_cisaillement
|
|
,const Met_abstraite::Umat_cont& ex
|
|
)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_deps";
|
|
#endif
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
module_compressibilite=module_cisaillement=0.; // init
|
|
energ.Inita(0.); // initialisation des énergies mises en jeux
|
|
|
|
// --- pour les contraintes passage en 3D
|
|
bool plusZero = true;
|
|
// on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
|
|
// sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
|
|
sig_HH_t_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_t),plusZero);
|
|
// --- pour la cinématique
|
|
// passage des informations spécifique à la loi le_SaveResul
|
|
// initialisation du tenseurs contrainte
|
|
sig_HH_3D.Inita(0.);
|
|
sig_HH_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_tdt),plusZero);
|
|
// passage des métriques de l'ordre 1 vers 3
|
|
Passage_metrique_ordre1_vers_3(ex);
|
|
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
|
|
Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
|
|
Vecteur* titi = NULL; // " " "
|
|
Passage_deformation_volume_ordre1_vers_3(DepsBB,epsBB_tdt,toto,delta_epsBB
|
|
,jacobien_0,jacobien,titi,*(ex.jacobien_t));
|
|
|
|
// -- construction de la base ortho associée à la direction de traction/compression
|
|
// ici il s'agit pour le premier vecteur de la direction 1,
|
|
// et compte tenu de Passage_metrique_ordre1_vers_3, le repère gi est orthogonale
|
|
// on va donc le normer
|
|
BaseB Vi_B(3); BaseH Vi_H(3);
|
|
BaseB * ViB = &Vi_B;BaseH * ViH = &Vi_H;
|
|
// tout d'abord on récupère la base gi actuelle et on la norme:
|
|
// on se sert de Vi_B comme stockage transitoire
|
|
// cout << "\n base giB(1) "<< (* ex.giB_tdt)(1)
|
|
// << "\n base giB(2) "<< (* ex.giB_tdt)(2);
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n base giB_tdt_3D(1) "<< giB_tdt_3D(1)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
|
|
};
|
|
#endif
|
|
|
|
|
|
Vi_B.CoordoB(1) = giB_tdt_3D(1);Vi_B.CoordoB(1).Normer();
|
|
Vi_B.CoordoB(2) = giB_tdt_3D(2);Vi_B.CoordoB(2).Normer();
|
|
Vi_B.CoordoB(3) = giB_tdt_3D(3);Vi_B.CoordoB(3).Normer();
|
|
// calcul des coordonnées locales de V_a en contravariants:
|
|
// V_a^{~i} = \vec V_a . \vec g^j
|
|
for (int a=1;a<4;a++) for (int i=1;i<4;i++)
|
|
Vi_H.CoordoH(a)(i) = Vi_B.CoordoB(a) * giH_tdt_3D(i);
|
|
// calcul des coordonnées contra
|
|
// Vi_B(a) = g_ij * Vi_H(a)
|
|
for (int a=1;a<4;a++)
|
|
Vi_B.CoordoB(a) = gijBB_t_3D * Vi_H(a);
|
|
|
|
// appel de la résolution 3D
|
|
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
|
|
// ceci dans le cas où l'axe de traction simple est quelconque
|
|
// ici on suppose que:
|
|
// - tous les calculs sont en dim 3
|
|
// - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
|
|
//
|
|
// ViB et ViH : correspond à la base de traction telle que ViB(1) = ViH(1) = l'axe de traction
|
|
// La base Vi est supposée orthonormée
|
|
// en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormee
|
|
// le tenseur de déformation et son incrémentsont également en orthonormees
|
|
// si = false: les bases transmises sont utilisées
|
|
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
|
|
Calcul_dsigma_deps (en_base_orthonormee,ViB,ViH,sig_HH_t_3D,Deps_BB_3D
|
|
,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D,save_resul.l_energ,save_resul.l_energ_t
|
|
,module_compressibilite,module_cisaillement,*umat_cont_3D);
|
|
// passage des tenseurs résultats à l'ordre 1
|
|
((Tenseur1HH*) &sigHH_tdt)->Affectation_trans_dimension(sig_HH_3D,false);
|
|
// maintenant on renseigne le tenseur de sortie
|
|
bool pluszero = true;
|
|
d_sigma_deps.Affectation_trans_dimension(d_sigma_deps_3D, pluszero);
|
|
|
|
};
|
|
|
|
// passage entre le calcul classique en 1D de Calcul_DsigmaHH_tdt et le calcul 3D
|
|
// méthode 2
|
|
void LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_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 &
|
|
,double& module_compressibilite,double& module_cisaillement
|
|
,const Met_abstraite::Impli& ex)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 7)
|
|
cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_td";
|
|
#endif
|
|
// récup du conteneur spécifique
|
|
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
|
|
module_compressibilite=module_cisaillement=0.; // init
|
|
energ.Inita(0.); // initialisation des énergies mises en jeux
|
|
|
|
// --- pour les contraintes passage en 3D
|
|
bool plusZero = true;
|
|
// on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
|
|
// sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
|
|
sig_HH_t_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_t),plusZero);
|
|
// --- pour la cinématique
|
|
// passage des informations spécifique à la loi le_SaveResul
|
|
// initialisation du tenseurs contrainte
|
|
sig_HH_3D.Inita(0.);
|
|
sig_HH_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_tdt),plusZero);
|
|
// passage des métriques de l'ordre 1 vers 3
|
|
Passage_metrique_ordre1_vers_3(ex);
|
|
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
|
|
Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
|
|
Vecteur* titi = NULL; // " " "
|
|
Passage_deformation_volume_ordre1_vers_3(DepsBB,epsBB_tdt,toto,delta_epsBB
|
|
,jacobien_0,jacobien,titi,*(ex.jacobien_t));
|
|
|
|
// -- construction de la base ortho associée à la direction de traction/compression
|
|
// ici il s'agit pour le premier vecteur de la direction 1,
|
|
// et compte tenu de Passage_metrique_ordre1_vers_3, le repère gi est orthogonale
|
|
// on va donc le normer
|
|
BaseB Vi_B(3); BaseH Vi_H(3);
|
|
BaseB * ViB = &Vi_B;BaseH * ViH = &Vi_H;
|
|
// tout d'abord on récupère la base gi actuelle et on la norme:
|
|
// on se sert de Vi_B comme stockage transitoire
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 8)
|
|
{ cout << "\n base giB_tdt_3D(1) "<< giB_tdt_3D(1)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
|
|
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
|
|
};
|
|
#endif
|
|
|
|
Vi_B.CoordoB(1) = giB_tdt_3D(1);Vi_B.CoordoB(1).Normer();
|
|
Vi_B.CoordoB(2) = giB_tdt_3D(2);Vi_B.CoordoB(2).Normer();
|
|
Vi_B.CoordoB(3) = giB_tdt_3D(3);Vi_B.CoordoB(3).Normer();
|
|
// calcul des coordonnées locales de V_a en contravariants:
|
|
// V_a^{~i} = \vec V_a . \vec g^j
|
|
for (int a=1;a<4;a++) for (int i=1;i<4;i++)
|
|
Vi_H.CoordoH(a)(i) = Vi_B.CoordoB(a) * giH_tdt_3D(i);
|
|
// calcul des coordonnées contra
|
|
// Vi_B(a) = g_ij * Vi_H(a)
|
|
for (int a=1;a<4;a++)
|
|
Vi_B.CoordoB(a) = gijBB_t_3D * Vi_H(a);
|
|
|
|
#ifdef MISE_AU_POINT
|
|
if (Permet_affichage() > 4)
|
|
{cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_td( ";
|
|
cout << "\n base Vi_B(1) "<< Vi_B(1) << "\n base Vi_B(2) "<< Vi_B(2)
|
|
<< "\n base Vi_B(3) "<< Vi_B(3)
|
|
<< "\n base Vi_H(1) "<< Vi_H(1) << "\n base Vi_H(2) "<< Vi_H(2)
|
|
<< "\n base Vi_H(3) "<< Vi_H(3) << flush;
|
|
};
|
|
#endif
|
|
|
|
|
|
// appel de la résolution 3D
|
|
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
|
|
// ceci dans le cas où l'axe de traction simple est quelconque
|
|
// ici on suppose que:
|
|
// - tous les calculs sont en dim 3
|
|
// - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
|
|
//
|
|
// ViB et ViH : correspond à la base de traction telle que ViB(1) = ViH(1) = l'axe de traction
|
|
// La base Vi est supposée orthonormée
|
|
// en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormee
|
|
// le tenseur de déformation et son incrémentsont également en orthonormees
|
|
// si = false: les bases transmises sont utilisées
|
|
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
|
|
bool en_base_orthonormee = false;
|
|
Calcul_dsigma_deps (en_base_orthonormee,ViB,ViH,sig_HH_t_3D,Deps_BB_3D
|
|
,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
|
|
,sig_HH_3D,d_sigma_deps_3D,save_resul.l_energ,save_resul.l_energ_t
|
|
,module_compressibilite,module_cisaillement,*umat_cont_3D);
|
|
// passage des tenseurs résultats à l'ordre 1
|
|
((Tenseur1HH*) &sigHH_tdt)->Affectation_trans_dimension(sig_HH_3D,false);
|
|
// maintenant on renseigne le tenseur de sortie
|
|
bool pluszero = true;
|
|
|
|
// d_sigma / d_ddl = d_sigma / d_eps * d_eps / d_ddl
|
|
// ici l'opérateur tangent correspond à la sensibilité de d sigma'^11 / d eps'_11
|
|
// et la direction V_1 c'est la direction de g_1 (à la norme près donc
|
|
// d_sigma / d_eps normalement c'est == à d sigma'^11 / d eps'_11
|
|
// du coup ici on n'a pas à itérer sur toutes les déformations
|
|
int nb_ddlPlus1 = d_epsBB.Taille()+1;
|
|
for (int i_ddl=1;i_ddl<nb_ddlPlus1;i_ddl++)
|
|
d_sigHH(i_ddl)->Coor(1,1) = (d_sigma_deps_3D(1,1,1,1) * (*d_epsBB(i_ddl)))(1,1);
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|