1489 lines
76 KiB
C++
Executable file
1489 lines
76 KiB
C++
Executable file
// FICHIER : SfeMembT.cc
|
|
// CLASSE : SfeMembT
|
|
|
|
|
|
// 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-2021 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 <iostream>
|
|
#include <cmath>
|
|
using namespace std; //introduces namespace std
|
|
#include <stdlib.h>
|
|
#include "Sortie.h"
|
|
#include "MathUtil.h"
|
|
|
|
#include "SfeMembT.h"
|
|
#include "TypeConsTens.h"
|
|
#include "FrontPointF.h"
|
|
#include "Loi_Umat.h"
|
|
#include "ExceptionsElemMeca.h"
|
|
#include "ExceptionsLoiComp.h"
|
|
#include "CharUtil.h"
|
|
|
|
// =========================== variables communes ===============================
|
|
// a priori les pointeurs de fonction pointent sur rien au debut
|
|
TenseurHH* SfeMembT::sig_bulk_pourSfe_HH = NULL; // variable de travail pour le bulk
|
|
|
|
// =========================== fin variables communes ===========================
|
|
// ---- definition du constructeur de la classe conteneur de donnees communes ------------
|
|
|
|
SfeMembT::DonnComSfe::DonnComSfe (ElemGeomC0* elCentre,const GeomSeg& seg,const DdlElement& tab
|
|
,DdlElement& tabErr,DdlElement& tab_Err1Sig,DdlElement& t_ddlXi_C
|
|
,const Met_Sfe1& met_gene,
|
|
Tableau <Vecteur *> & resEr,Mat_pleine& raidEr,
|
|
ElemGeomC0* elEr,GeomSeg& segEr,ElemGeomC0* elS,GeomSeg& segmS
|
|
,Vecteur& residu_int,Mat_pleine& raideur_int,
|
|
Tableau <Vecteur* > & residus_extN,Tableau <Mat_pleine* >& raideurs_extN,
|
|
Tableau <Vecteur* > & residus_extA,Tableau <Mat_pleine* >& raideurs_extA,
|
|
Tableau <Vecteur* >& residus_extS,Tableau <Mat_pleine* >& raideurs_extS,
|
|
Mat_pleine& mat_masse,ElemGeomC0* elMas,const GeomSeg& segMas,int nbi,
|
|
Tableau <EnuTypeCL> const & tabTypeCL,Tableau <Coordonnee3> const & vplan
|
|
,Tableau <int>& nMetddl,const Met_abstraite& met_cent ) :
|
|
eleCentre(elCentre),segment(seg),tab_ddl(tab),tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig)
|
|
,tab_ddlXi_C(t_ddlXi_C),met_SfeMembT(met_gene)
|
|
,matGeom(tab.NbDdl(),tab.NbDdl()),matInit(tab.NbDdl(),tab.NbDdl())
|
|
,d_epsBB(tab.NbDdl()),d_sigHH(tab.NbDdl()),d2_epsBB(nbi)
|
|
,resErr(resEr),raidErr(raidEr),eleEr(elEr),segmentEr(segEr),eleS(elS),segS(segmS)
|
|
,residu_interne(residu_int),raideur_interne(raideur_int)
|
|
,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN)
|
|
,residus_externeA(residus_extA),raideurs_externeA(raideurs_extA)
|
|
,residus_externeS(residus_extS),raideurs_externeS(raideurs_extS)
|
|
,matrice_masse(mat_masse),eleMas(elMas),segmentMas(segMas)
|
|
,tabType_rienCL(tabTypeCL),vplan_rien(vplan)
|
|
,nMetVTab_ddl(nMetddl),met_surf_cent(met_cent)
|
|
,sfematD(NULL),sferesD(NULL),noeud_a_prendre_en_compte(NULL)
|
|
{
|
|
int nbddl = tab.NbDdl();
|
|
int dim_tens = met_gene.Nbvec_des_bases();
|
|
for (int ni=1;ni<=nbi;ni++)
|
|
{d2_epsBB(ni).Change_taille(nbddl);
|
|
for (int i1=1; i1<= nbddl; i1++)
|
|
for (int i2=1; i2<= nbddl; i2++)
|
|
d2_epsBB(ni)(i1,i2) = NevezTenseurBB (dim_tens);
|
|
};
|
|
// int tailledeps = d_epsBB.Taille();
|
|
// for (int i=1;i<= tailledeps; i++)
|
|
for (int i=1;i<= nbddl; i++)
|
|
d_epsBB(i) = NevezTenseurBB (dim_tens);
|
|
// int tailledsig = d_sigHH.Taille();
|
|
// for (int j=1;j<= tailledsig; j++)
|
|
for (int j=1;j<= nbddl; j++)
|
|
d_sigHH(j) = NevezTenseurHH (dim_tens);
|
|
};
|
|
// Il ne doit exister qu'un exemplaire de la classe, donc au niveau du constructeur de copie les tableaux
|
|
// de pointeurs pointes sur les mêmes entités que la valeur passée en paramètre
|
|
SfeMembT::DonnComSfe::DonnComSfe(DonnComSfe& a) :
|
|
eleCentre(a.eleCentre),segment(a.segment)
|
|
,tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr),tab_Err1Sig11(a.tab_Err1Sig11)
|
|
,tab_ddlXi_C(a.tab_ddlXi_C),met_SfeMembT(a.met_SfeMembT),matGeom(a.matGeom)
|
|
,matInit(a.matInit),d2_epsBB(a.d2_epsBB),resErr(a.resErr),raidErr(a.raidErr)
|
|
,eleEr(a.eleEr),segmentEr(a.segmentEr),eleS(a.eleS),segS(a.segS)
|
|
,d_epsBB(a.d_epsBB),d_sigHH(a.d_sigHH)
|
|
,residu_interne(a.residu_interne),raideur_interne(a.raideur_interne)
|
|
,residus_externeN(a.residus_externeN),raideurs_externeN(a.raideurs_externeN)
|
|
,residus_externeA(a.residus_externeA),raideurs_externeA(a.raideurs_externeA)
|
|
,residus_externeS(a.residus_externeS),raideurs_externeS(a.raideurs_externeS)
|
|
,matrice_masse(a.matrice_masse),eleMas(a.eleMas),segmentMas(a.segmentMas)
|
|
,tabType_rienCL(a.tabType_rienCL),vplan_rien(a.vplan_rien)
|
|
,nMetVTab_ddl(a.nMetVTab_ddl),met_surf_cent(a.met_surf_cent)
|
|
,sfematD(NULL),sferesD(NULL),noeud_a_prendre_en_compte(NULL)
|
|
|
|
{ int nbddl = tab_ddl.NbDdl();
|
|
int nbi=d2_epsBB.Taille();
|
|
for (int ni=1;ni<=nbi;ni++)
|
|
for (int i1=1; i1<= nbddl; i1++)
|
|
for (int i2=1; i2<= nbddl; i2++)
|
|
d2_epsBB(ni)(i1,i2) = NevezTenseurBB (*(a.d2_epsBB(ni)(i1,i2)));
|
|
// int tailledeps = d_epsBB.Taille();
|
|
// for (int i=1;i<= tailledeps; i++)
|
|
for (int i=1;i<= nbddl; i++)
|
|
d_epsBB(i) = NevezTenseurBB (*(a.d_epsBB(i)));
|
|
// int tailledsig = d_sigHH.Taille();
|
|
// for (int j=1;j<= tailledsig; j++)
|
|
for (int j=1;j<= nbddl; j++)
|
|
d_sigHH(j) = NevezTenseurHH (*(d_sigHH(j)));
|
|
|
|
// cas d'une stabilisation éventuelle
|
|
if (a.sfematD != NULL)
|
|
sfematD = new Mat_pleine(*a.sfematD);
|
|
if (a.sferesD != NULL)
|
|
sferesD = new Vecteur (*a.sferesD);
|
|
if (a.noeud_a_prendre_en_compte != NULL)
|
|
noeud_a_prendre_en_compte = new Tableau<int> (*a.noeud_a_prendre_en_compte);
|
|
};
|
|
|
|
// certaines grandeurs pointées qui n'exitent qu'en un seul exemplaire
|
|
// ne seront supprimées que par l'ordre destruction,
|
|
SfeMembT::DonnComSfe::DonnComSfe::~DonnComSfe()
|
|
{ int nbddl = tab_ddl.NbDdl();
|
|
int nbi=d2_epsBB.Taille();
|
|
for (int ni=1;ni<=nbi;ni++)
|
|
for (int i1=1; i1<= nbddl; i1++)
|
|
for (int i2=1; i2<= nbddl; i2++)
|
|
delete d2_epsBB(ni)(i1,i2);
|
|
int tailledeps = d_epsBB.Taille();
|
|
for (int i=1;i<= tailledeps; i++)
|
|
delete d_epsBB(i);
|
|
int tailledsig = d_sigHH.Taille();
|
|
for (int j=1;j<= tailledsig; j++)
|
|
delete d_sigHH(j);
|
|
// cas d'une stabilisation éventuelle
|
|
if (sfematD != NULL)
|
|
delete sfematD;
|
|
if (sferesD != NULL)
|
|
delete sferesD;
|
|
if (noeud_a_prendre_en_compte != NULL)
|
|
delete noeud_a_prendre_en_compte;
|
|
|
|
};
|
|
|
|
|
|
// ---------- fin definition de la classe conteneur de donnees communes ------------
|
|
// -+-+ definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
|
|
SfeMembT::UneFois::UneFois () : // constructeur par défaut
|
|
doCoMemb(NULL),CalResPrem_t(0),CalResPrem_tdt(0),CalimpPrem(0),dualSortSfe(0)
|
|
,CalSMlin_t(0),CalSMlin_tdt(0),CalSMRlin(0)
|
|
,CalSMsurf_t(0),CalSMsurf_tdt(0),CalSMRsurf(0)
|
|
,CalSMvol_t(0),CalSMvol_tdt(0),CalSMvol(0)
|
|
,CalDynamique(0),CalPt_0_t_tdt(0)
|
|
,nbelem_in_Prog(0)
|
|
{};
|
|
SfeMembT::UneFois::~UneFois ()
|
|
{ delete doCoMemb;
|
|
};
|
|
|
|
// -+-+ fin definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
|
|
|
|
// CONSTRUCTEURS :
|
|
// Constructeur par defaut
|
|
SfeMembT::SfeMembT () :
|
|
ElemMeca(),lesPtMecaInt(),donnee_specif(),areteTypeCL(NULL),vplan(NULL)
|
|
,t_N_centre(),unefois(NULL),nombre(NULL)
|
|
{ lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca
|
|
};
|
|
// Constructeur fonction d'un numero
|
|
// d'identification , d'identificateur d'interpolation et de geometrie
|
|
SfeMembT::SfeMembT (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,string info) :
|
|
ElemMeca(num_mail,num_id,id_interp_elt,id_geom_elt,info),lesPtMecaInt(),donnee_specif()
|
|
,areteTypeCL(NULL),vplan(NULL),t_N_centre(),unefois(NULL),nombre(NULL)
|
|
{ lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca
|
|
};
|
|
|
|
// Constructeur fonction d'un numero de maillage et d'identification,
|
|
// du tableau de connexite des noeuds, d'identificateur d'interpolation et de geometrie
|
|
SfeMembT::SfeMembT (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,
|
|
const Tableau<Noeud *>& tab,string info) :
|
|
ElemMeca(num_mail,num_id,tab,id_interp_elt,id_geom_elt,info),lesPtMecaInt(),donnee_specif()
|
|
,areteTypeCL(NULL),vplan(NULL),t_N_centre(),unefois(NULL),nombre(NULL)
|
|
{ lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca
|
|
};
|
|
// Constructeur de copie
|
|
SfeMembT::SfeMembT (const SfeMembT& SfeM) :
|
|
ElemMeca (SfeM)
|
|
,lesPtMecaInt(SfeM.lesPtMecaInt)
|
|
,donnee_specif(SfeM.donnee_specif),unefois(SfeM.unefois),nombre(SfeM.nombre)
|
|
,areteTypeCL(SfeM.areteTypeCL),vplan(SfeM.vplan),t_N_centre(SfeM.t_N_centre)
|
|
{ lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca
|
|
// --- cas des différentes puissances et ... ---
|
|
SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
residu = &(CoSfe->residu_interne); // residu local
|
|
raideur = &(CoSfe->raideur_interne); // raideur locale
|
|
// --- cas de la dynamique -----
|
|
mat_masse = &(CoSfe->matrice_masse);
|
|
// --- cas des efforts externes concernant noeuds ------
|
|
res_extN = &(CoSfe->residus_externeN); // pour les résidus et second membres
|
|
raid_extN= &(CoSfe->raideurs_externeN);// pour les raideurs
|
|
// --- cas des efforts externes concernant les aretes ------
|
|
res_extA = &(CoSfe->residus_externeA); // pour les résidus et second membres
|
|
raid_extA= &(CoSfe->raideurs_externeA);// pour les raideurs
|
|
// --- cas des efforts externes concernant les faces ------
|
|
res_extS= &(CoSfe->residus_externeS); // pour les résidus et second membres
|
|
raid_extS= &(CoSfe->raideurs_externeS); // pour les raideurs
|
|
// pour les conditions limites on les crée si besoin
|
|
if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL)))
|
|
{ // si ça ne pointe pas sur la valeur par défaut, alors c'est qu'il y a
|
|
// des CL particulières, on les crée donc à l'identique
|
|
areteTypeCL = new Tableau <EnuTypeCL> (*(SfeM.areteTypeCL));
|
|
vplan = new Tableau <Coordonnee3> (*(SfeM.vplan));
|
|
};
|
|
};
|
|
|
|
// DESTRUCTEUR :
|
|
SfeMembT::~SfeMembT ()
|
|
{ /* // on va libérer les déformations d'arrête et de surface qui sont un peu particulières
|
|
if (defArete.Taille() != 0)
|
|
if (defArete(1) != NULL)
|
|
{ delete defArete(1);
|
|
// on remet à null les pointeurs pour un appel correcte du destructeur d'elemMeca
|
|
for (int ia=1;ia<= 12; ia++)
|
|
defArete(ia) = NULL;
|
|
}
|
|
if (defSurf.Taille() != 0)
|
|
if (defSurf(1) != NULL)
|
|
{ delete defSurf(1);
|
|
// on remet à null les pointeurs pour un appel correcte du destructeur d'elemMeca
|
|
for (int ia=1;ia<= 12; ia++)
|
|
defSurf(ia) = NULL;
|
|
} */
|
|
// concernant tous les éléments qui sont définit qu'une seule fois, leur destruction est effectuée avec la méthode destruction
|
|
// appelée par les classes descendantes
|
|
LibereTenseur();
|
|
};
|
|
|
|
|
|
// Lecture des donnees de la classe sur fichier
|
|
void SfeMembT::LectureDonneesParticulieres
|
|
(UtilLecture * entreePrinc,Tableau<Noeud *> * tabMaillageNoeud)
|
|
{ int nb;int nbnte = nombre->nbnte;
|
|
tab_noeud.Change_taille(nbnte);
|
|
for (int i=1; i<= nbnte; i++)
|
|
{ *(entreePrinc->entree) >> nb;
|
|
if ((entreePrinc->entree)->rdstate() == 0)
|
|
// pour mémoire ici on a
|
|
/* enum io_state
|
|
{ badbit = 1<<0, // -> 1 dans rdstate()
|
|
eofbit = 1<<1, // -> 2
|
|
failbit = 1<<2, // -> 4
|
|
goodbit = 0 // -> O
|
|
};*/
|
|
tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
|
|
#ifdef ENLINUX
|
|
else if ((entreePrinc->entree)->fail())
|
|
// on a atteind la fin de la ligne et on appelle un nouvel enregistrement
|
|
{ entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement
|
|
*(entreePrinc->entree) >> nb;
|
|
tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
|
|
}
|
|
#else
|
|
/* #ifdef SYSTEM_MAC_OS_X_unix
|
|
else if ((entreePrinc->entree)->fail())
|
|
// on a atteind la fin de la ligne et on appelle un nouvel enregistrement
|
|
{ entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement
|
|
*(entreePrinc->entree) >> nb;
|
|
tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture normale
|
|
}
|
|
#else */
|
|
else if ((entreePrinc->entree)->eof())
|
|
// la lecture est bonne mais on a atteind la fin de la ligne
|
|
{ tab_noeud(i) = (*tabMaillageNoeud)(nb); // lecture
|
|
// si ce n'est pas la fin de la lecture on appelle un nouvel enregistrement
|
|
if (i != nbnte) entreePrinc->NouvelleDonnee(); // lecture d'un nouvelle enregistrement
|
|
}
|
|
// #endif
|
|
#endif
|
|
else // cas d'une erreur de lecture
|
|
{ cout << "\n erreur de lecture inconnue ";
|
|
entreePrinc->MessageBuffer("** lecture des données particulières **");
|
|
cout << "SfeMembT::LecDonPart";
|
|
Affiche(2);
|
|
Sortie (1);
|
|
}
|
|
};
|
|
// construction du tableau de ddl des noeuds de SfeMembT
|
|
ConstTabDdl();
|
|
// construction du tableau des noeuds du centre
|
|
t_N_centre.Change_taille(nombre->nbnce);
|
|
for (int i=1;i<= nombre->nbnce;i++)
|
|
t_N_centre(i)=tab_noeud(i);
|
|
};
|
|
|
|
// calcul d'un point dans l'élément réel en fonction des coordonnées dans l'élément de référence associé
|
|
// temps: indique si l'on veut les coordonnées à t = 0, ou t ou tdt
|
|
// 1) cas où l'on utilise la place passée en argument
|
|
Coordonnee & SfeMembT::Point_physique(const Coordonnee& ,Coordonnee & co,Enum_dure )
|
|
{ cout << "\n methode non encore implantee !!"
|
|
<< "\n Coordonnee & SfeMembT::Point_physique(...";
|
|
Sortie(1);
|
|
// d) retour
|
|
return co;
|
|
};
|
|
|
|
// 3) cas où l'on veut les coordonnées aux 1, 2 ou trois temps selon la taille du tableau t_co
|
|
void SfeMembT::Point_physique(const Coordonnee& ,Tableau <Coordonnee> & )
|
|
{ cout << "\n methode non encore implantee !!"
|
|
<< "\n void SfeMembT::Point_physique(...";
|
|
Sortie(1);
|
|
};
|
|
|
|
// vérification que la courbure ne soit pas anormale au temps enu
|
|
// inf_normale : indique en entrée le det mini pour la courbure en locale
|
|
// retour 1: si tout est ok,
|
|
// 0: une courbure trop grande a été détecté
|
|
// retour = -1 : il y a un pb inconnue dans le calcul
|
|
int SfeMembT::Test_courbure_anormale3(Enum_dure enu,double inf_normale)
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
if (unefois->CalimpPrem == 0)
|
|
{ unefois->CalimpPrem = 1;
|
|
Tableau<Enum_variable_metrique> tab(20);
|
|
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
|
|
tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
|
|
tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt;
|
|
tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt;
|
|
tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
|
|
tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
};
|
|
// la variable de retour initialisée à 1 par défaut
|
|
int retour = 1;
|
|
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
|
|
try
|
|
{ // puis les indices de boucle
|
|
int indice = 1;
|
|
int nbtotsurf = CoSfe->eleCentre->Nbi();
|
|
int nbtotepais =(CoSfe->segment).Nbi();
|
|
for (int nisur =1; nisur<= nbtotsurf; nisur++) // boucle sur les pt d'integ de surface
|
|
{for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur
|
|
{((DeformationSfe1*) def)->ChangeNumIntegSfe1(nisur,niepais);
|
|
((DeformationSfe1*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
|
|
retour *= ((DeformationSfe1*) def)->Test_courbure_anormale(enu,inf_normale);
|
|
if (!retour)
|
|
break; // on arrête dès qu'un pb est repéré
|
|
}; // fin boucle sur les points d'intégration
|
|
if (!retour)
|
|
break; // on arrête dès qu'un pb est repéré
|
|
};
|
|
}
|
|
catch (ErrSortieFinale)
|
|
// cas d'une direction voulue vers la sortie
|
|
// on relance l'interuption pour le niveau supérieur
|
|
{ ErrSortieFinale toto;
|
|
throw (toto);
|
|
}
|
|
catch ( ... )
|
|
{ if (ParaGlob::NiveauImpression() > 3)
|
|
{cout << "\n warning: erreur generee par un element lors du test de verification de la courbure ";
|
|
};
|
|
retour = -1;
|
|
}
|
|
return retour;
|
|
};
|
|
|
|
// récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu
|
|
// ici il s'agit de grandeurs scalaire,
|
|
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
|
|
Tableau <double> SfeMembT::Valeur_a_diff_temps(bool absolue,Enum_dure enu_t,const List_io<Ddl_enum_etendu>& enu,int iteg)
|
|
{ // appel de la procedure de elem meca
|
|
int cas;
|
|
if (!(unefois->dualSortSfe ) && (unefois->CalimpPrem ))
|
|
{ cas=1;unefois->dualSortSfe += 1;
|
|
}
|
|
else if ((unefois->dualSortSfe ) && (unefois->CalimpPrem ))
|
|
{ cas = 11;}
|
|
else if (!(unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))
|
|
{ cas=2;unefois->dualSortSfe += 1;
|
|
}
|
|
else if ((unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))
|
|
{ cas = 12;}
|
|
// sinon problème
|
|
else
|
|
{ cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
|
|
<< " pour inialiser les conteneurs des tenseurs resultats ";
|
|
if (ParaGlob::NiveauImpression() >= 4)
|
|
cout << "\n cas non prévu, unefois->dualSortSfe= " << unefois->dualSortSfe
|
|
<< " unefois->CalimpPrem= " << unefois->CalimpPrem
|
|
<< "\n SfeMembT::Valeur_a_diff_temps(Enum_dure enu_t...";
|
|
Sortie(1);
|
|
};
|
|
return ElemMeca::Valeur_multi(absolue,enu_t,enu,iteg,cas);
|
|
};
|
|
|
|
// récupération des valeurs au numéro d'ordre = iteg pour les grandeurs enu
|
|
// ici il s'agit de grandeurs tensorielles, le retour s'effectue dans la liste
|
|
// de conteneurs quelconque associée
|
|
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
|
|
void SfeMembT::ValTensorielle_a_diff_temps(bool absolue,Enum_dure enu_t,List_io<TypeQuelconque>& enu,int iteg)
|
|
{ // appel de la procedure de elem meca
|
|
int cas;
|
|
if (!(unefois->dualSortSfe ) && (unefois->CalimpPrem ))
|
|
{ cas=1;unefois->dualSortSfe += 1;
|
|
}
|
|
else if ((unefois->dualSortSfe ) && (unefois->CalimpPrem ))
|
|
{ cas = 11;}
|
|
else if (!(unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))
|
|
{ cas=2;unefois->dualSortSfe += 1;
|
|
}
|
|
else if ((unefois->dualSortSfe ) && (unefois->CalResPrem_tdt ))
|
|
{ cas = 12;}
|
|
// sinon problème
|
|
else
|
|
{ cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
|
|
<< " pour inialiser les conteneurs des tenseurs resultats ";
|
|
if (ParaGlob::NiveauImpression() >= 4)
|
|
cout << "\n cas non prévu, unefois->dualSortSfe= " << unefois->dualSortSfe
|
|
<< " unefois->CalimpPrem= " << unefois->CalimpPrem
|
|
<< "\n SfeMembT::ValTensorielle_a_diff_temps(Enum_dure enu_t...";
|
|
Sortie(1);
|
|
};
|
|
ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas);
|
|
};
|
|
|
|
// Calcul du residu local à t ou tdt en fonction du booleen atdt
|
|
Vecteur* SfeMembT::CalculResidu (bool atdt,const ParaAlgoControle & pa)
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
Tableau <TenseurBB *>& d_epsBB = (CoSfe->d_epsBB);// "
|
|
energie_totale.Initialisation_differenciee(energie_totale_t); // init de l'énergie totale sur l'élément
|
|
TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
|
|
if (ElemMeca::Bulk_visco_actif())
|
|
{ if (sig_bulk_pourSfe_HH == NULL) {sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);}
|
|
else if (sig_bulk_pourSfe_HH->Dimension() != sigHH_t_1->Dimension())
|
|
{delete sig_bulk_pourSfe_HH; sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);};
|
|
};
|
|
// dimensionnement de la metrique
|
|
if (!atdt)
|
|
{if( !(unefois->CalResPrem_t ))
|
|
{ unefois->CalResPrem_t += 1;
|
|
Tableau<Enum_variable_metrique> tab(10);
|
|
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
|
|
tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
|
|
tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
};}
|
|
else
|
|
{if( !(unefois->CalResPrem_tdt ))
|
|
{unefois->CalResPrem_tdt += 1;
|
|
Tableau<Enum_variable_metrique> tab(11);
|
|
tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;tab(4) = igijBB_tdt;
|
|
tab(5) = igijHH_tdt; tab(6) = id_giB_tdt; tab(7) = id_gijBB_tdt ;
|
|
tab(8) = id_gijBB_t ;tab(9) = id_gijBB_tdt ;
|
|
tab(10) = igradVBB_tdt;tab(11) = iVtdt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
};};
|
|
// dimensionnement du residu
|
|
int nbddl = CoSfe->tab_ddl.NbDdl();
|
|
if ( residu == NULL)
|
|
residu = new Vecteur(nbddl); // cas du premier passage
|
|
else
|
|
for (int i =1; i<= nbddl; i++) // cas des autres passages
|
|
(*residu)(i) = 0.; // on initialise a zero
|
|
double jacobien; // def du jacobien
|
|
int indice = 1;
|
|
int nbtotsurf =CoSfe->eleCentre->Nbi();
|
|
int nbtotepais =(CoSfe->segment).Nbi();
|
|
// préparation du calcul de l'épaisseur finale (qui dépend de tous les points d'intégration)
|
|
double somme_epaisseurs_finales = 0;
|
|
// on utilise une épaisseur intermédiaire pour la métrique
|
|
// de manière à utiliser l'épaisseur de l'incrément précédent
|
|
Epai epais_pour_metrique;
|
|
|
|
// on doit boucler à l'extérieur sur les pt de surface et à l'intérieur sur les pt d'épaisseurs
|
|
// pour profiter du fait que une fois un pt calculer, on n'a pas à recalculer la courbure !!!!!
|
|
for (int nisur =1; nisur<= nbtotsurf; nisur++) // boucle sur les pt d'integ de surface
|
|
for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur
|
|
|
|
|
|
{ if ((donnee_specif.epais != NULL) // on change les épaisseurs si elles sont stockées
|
|
&& (!Loi_rien(loiComp->Id_comport()))
|
|
)
|
|
{ epais_pour_metrique = *(donnee_specif.epais);
|
|
epais_pour_metrique.epaisseur_tdt = epais_pour_metrique.epaisseur_t;
|
|
// ((DeformationSfe1*) def)->Change_epaisseur(*donnee_specif.epais); // au niveau de l'élément
|
|
((DeformationSfe1*) def)->Change_epaisseur(epais_pour_metrique);
|
|
};
|
|
((DeformationSfe1*) def)->ChangeNumIntegSfe1(nisur,niepais);
|
|
((DeformationSfe1*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
|
|
|
|
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
|
|
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
|
|
TenseurBB & DepsBB_ = *(ptIntegMeca.DepsBB());
|
|
EnergieMeca& energ = tab_energ(indice);
|
|
EnergieMeca& energ_t = tab_energ_t(indice);
|
|
double poids = CoSfe->eleCentre->Wi(nisur) * (CoSfe->segment).Wi(niepais) ;
|
|
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
|
|
if (loiTP != NULL) {sTP = tabSaveTP(indice);}; // au pt d'integ si elles existes
|
|
if ((loiComp->Id_comport()==LOI_VIA_UMAT) // cas de l'utilisation d'une loi umat
|
|
|| (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
|
|
((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);
|
|
double epaisseur_tdt = 0.; // ne sert que pour les éléments 2D avec épaisseurs stockées à l'élément
|
|
double epaisseur_pris_en_compte = 0.; // celle qui sera réellement prise en compte
|
|
|
|
// -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
|
|
try // ce qui permet de donner les numéros d'éléments et de pti
|
|
{if (atdt)
|
|
{ const Met_abstraite::Expli_t_tdt& ex = loiComp->Cal_explicit_tdt
|
|
(tabSaveDon(indice), *def,CoSfe->tab_ddl,ptIntegMeca,d_epsBB,jacobien
|
|
,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
|
|
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
|
|
if (Bulk_visco_actif())
|
|
ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_ ,(*sig_bulk_pourSfe_HH));
|
|
// on tient compte de la variation d'épaisseur si l'épaisseur est stockée dans l'élément
|
|
// et que l'on n'a pas une loi rien
|
|
if ((donnee_specif.epais != NULL) && (!Loi_rien(loiComp->Id_comport())))
|
|
{ switch (loiComp->Comportement_3D_CP_DP_1D())
|
|
{ case COMP_CONTRAINTES_PLANES :
|
|
// en contrainte plane on recalcule l'épaisseur
|
|
{ CalEpaisseurAtdtExp_et_vol_pti(donnee_specif.epais->epaisseur0,pa,donnee_specif.epais->epaisseur_t
|
|
,ptIntegMeca,epaisseur_tdt,ex);
|
|
//------ le 7 nov 2016: changement car je ne veux pas utiliser la variation d'épaisseur
|
|
// qui est calculée dans les lois. C'est sans doute possible, bien qu'avec les lois additives et mélange
|
|
// ce n'est pas sûr que cela veuille dire quelque chose. Donc on essaie de sans passer
|
|
// De plus c'est cohérent avec les éléments membranes ...
|
|
|
|
// ---- remodif le 1 avril 2017 car la méthode précédente fait que l'élément converge très mal !!
|
|
// je ne sais pas trop pourquoi (avec une loi de déformation plane, cela converge immédiatement !)
|
|
// // en fait, on récupère la variation d'épaisseur calculée via la condition de contrainte plane
|
|
// // on utilise uniquement la variation d'épaisseur calculée à la surface médiane ??? pour l'instant
|
|
// epaisseur_tdt = loiComp->HsurH0(tabSaveDon(nisur)) * donnee_specif.epais->epaisseur0;
|
|
|
|
// il faudrait peut-être utiliser finalement la variation d'épaisseur calculée par la loi, comme pour les éléments
|
|
// membranes qui ont été modifiées pour tenir compte de la différence de def entre géométrique et la def méca
|
|
// Par exemple lors de critère de rupture, on aura le même pb qu'avec les plis
|
|
//*** à faire
|
|
epaisseur_pris_en_compte = donnee_specif.epais->epaisseur_tdt;
|
|
break;
|
|
}
|
|
case COMP_DEFORMATIONS_PLANES :
|
|
// en deformation plane, l'épaisseur ne change pas
|
|
{ epaisseur_tdt = donnee_specif.epais->epaisseur0;
|
|
ptIntegMeca.Volume_pti() *= donnee_specif.epais->epaisseur0;
|
|
epaisseur_pris_en_compte = donnee_specif.epais->epaisseur0;
|
|
break;
|
|
}
|
|
default :
|
|
cout << "\nErreur : cas de loi de comportement : " << Nom_comp((loiComp->Id_comport()))
|
|
<< " non prevu pour le calcul le l'epaisseur avec l'element SFE ";
|
|
cout << "\n SfeMembT::Calcul_implicit(.... \n";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
somme_epaisseurs_finales += epaisseur_tdt; //pour le calcul final
|
|
}
|
|
else
|
|
{ const Met_abstraite::Expli& ex = loiComp->Cal_explicit_t
|
|
(tabSaveDon(indice), *def,CoSfe->tab_ddl,ptIntegMeca,d_epsBB,jacobien
|
|
,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
|
|
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
|
|
if (Bulk_visco_actif())
|
|
ModifContrainteBulk(TEMPS_t,*ex.gijHH_t,sigHH,DepsBB_ ,(*sig_bulk_pourSfe_HH));
|
|
// pour l'instant on expédie un message, car il y a des chances que cette partie disparaisse !!
|
|
cout << "\n attention *** calcul non operationnelle en explicite t, on ne tient pas compte de la"
|
|
<< " variation d'epaisseur !! (a modifier ) ";
|
|
};
|
|
}
|
|
catch (ErrNonConvergence_loiDeComportement excep)
|
|
// cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
|
|
// d'une loi de comportement incrémentale
|
|
{ cout << "\n erreur de loi de comportement element= " << this->Num_elt()
|
|
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
|
|
throw excep;
|
|
}
|
|
catch (ErrSortieFinale)
|
|
// cas d'une direction voulue vers la sortie
|
|
// on relance l'interuption pour le niveau supérieur
|
|
{ ErrSortieFinale toto;
|
|
throw (toto);
|
|
}
|
|
catch ( ... )
|
|
{ cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
|
|
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
|
|
throw (Err_inconnue_ElemMeca());
|
|
};
|
|
|
|
if ((jacobien <= 0.)|| (std::isinf(jacobien))||( std::isnan(jacobien))) // vérif du jacobien
|
|
{ if ((std::isinf(jacobien))||( std::isnan(jacobien)))
|
|
// on met un message quelque soit le niveau d'impression
|
|
{ cout << "\n ********** attention on a trouve un jacobien infini ou nan !! = ("<<jacobien<<")********* " << endl; };
|
|
if (ParaGlob::NiveauImpression() >= 1)
|
|
{ cout << "\n ********** attention on a trouve un jacobien negatif!! ("<<jacobien<<")********* " << endl; };
|
|
};
|
|
if ((jacobien <= 0.) && (pa.JabobienNegatif() != 0)) // vérif du jacobien et traitement adéquate si besoin
|
|
{ switch (pa.JabobienNegatif())
|
|
{ case 1:
|
|
{ cout << "\n on annulle sa contribution. Element nb: " << Num_elt() << " ni_surf: " << nisur
|
|
<< " ni_ep: " << niepais;
|
|
break;
|
|
}
|
|
// dans les autres cas on ne fait rien
|
|
};
|
|
}
|
|
else
|
|
{ // on prend en compte toutes les variations
|
|
// dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
|
|
// épaisseur_tdt ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
|
|
// pour éviter des tests et une recopie du code on laisse epaisseur_tdt, avec une valeur 1
|
|
// ce qui fait qu'elle n'intervient pas
|
|
if (donnee_specif.epais == NULL) epaisseur_pris_en_compte = 1.;
|
|
// il ne faut pas interpoler l'énergie à t, car elle était associée à un volume qui a changé
|
|
// donc il ne faut considérer que l'accroissement d'énergie
|
|
energie_totale.Ajout_differenciee(energ,energ_t,(poids * jacobien * epaisseur_pris_en_compte));
|
|
// energie_totale += (energ-energ_t) * (poids * jacobien * epaisseur_tdt);
|
|
// energie_totale += energ * (poids * jacobien * epaisseur_tdt);
|
|
|
|
//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
|
|
//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
|
|
Delta_Jacobien_anormal(tabSaveDefDon(indice),nisur,niepais);
|
|
|
|
for (int j =1; j<= nbddl; j++)
|
|
{int jloc = CoSfe->nMetVTab_ddl(j); // adressage indirecte
|
|
(*residu)(jloc) -= (sigHH && (*(d_epsBB(j)))) * jacobien * poids * epaisseur_pris_en_compte;
|
|
volume += jacobien * poids * epaisseur_tdt; // là on prend la vrai épaisseur
|
|
};
|
|
};
|
|
if (ElemMeca::Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
|
|
{ sigHH -= (*sig_bulk_pourSfe_HH); } // ceci pour la sauvegarde ou autre utilisation
|
|
|
|
}; // fin boucle sur les points d'intégration
|
|
// on calcul éventuellement l'épaisseur moyenne finale si c'est un élément 2D
|
|
if ((donnee_specif.epais != NULL)&&(!Loi_rien(loiComp->Id_comport()))) // cas 2D
|
|
donnee_specif.epais->epaisseur_tdt = somme_epaisseurs_finales / (indice - 1);
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
return residu;
|
|
};
|
|
|
|
// Calcul du residu local et de la raideur locale,
|
|
// pour le schema implicite
|
|
Element::ResRaid SfeMembT::Calcul_implicit (const ParaAlgoControle & pa)
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
Tableau <TenseurBB *>& d_epsBB = (CoSfe->d_epsBB);// "
|
|
Tableau <TenseurHH *>& d_sigHH = (CoSfe->d_sigHH);// "
|
|
energie_totale.Initialisation_differenciee(energie_totale_t); // init de l'énergie totale sur l'élément
|
|
E_elem_bulk_tdt = E_elem_bulk_t; P_elem_bulk = 0.; // init pour l'énergie et la puissance associées au bulk
|
|
// init éventuel de la contrainte de bulk viscosity
|
|
TenseurHH* sigHH_t_1 = (*lesPtIntegMecaInterne)(1).SigHH_t(); // simplification d'écriture
|
|
if (ElemMeca::Bulk_visco_actif())
|
|
{ if (sig_bulk_pourSfe_HH == NULL) {sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);}
|
|
else if (sig_bulk_pourSfe_HH->Dimension() != sigHH_t_1->Dimension())
|
|
{delete sig_bulk_pourSfe_HH; sig_bulk_pourSfe_HH = NevezTenseurHH(*sigHH_t_1);};
|
|
};
|
|
|
|
// init éventuel des intégrales de volume et temps
|
|
Init_Integ_vol_et_temps();
|
|
// --- init pour l'init de la stabilisation de membrane-biel éventuelle
|
|
const Met_abstraite::Impli* ex_final=NULL; // contiendra après les boucles la dernière métrique
|
|
Mise_a_jour_A_calculer_force_stab(); // permettra ensuite de savoir où le calcul doit être fait
|
|
|
|
bool cald_Dvirtuelle = false;
|
|
if (unefois->CalimpPrem == 0)
|
|
{ unefois->CalimpPrem = 1;
|
|
Tableau<Enum_variable_metrique> tab(20);
|
|
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igiB_tdt; tab(4) = igijBB_0;
|
|
tab(5) = igijBB_t;tab(6) = igijBB_tdt; tab(7) = igijHH_tdt; tab(8) = id_giB_tdt;
|
|
tab(9) = id_gijBB_tdt ;tab(10) = igiH_tdt;tab(11) = id_giH_tdt;
|
|
tab(12) = id_gijHH_tdt;tab(13) = id_jacobien_tdt;tab(14) = id2_gijBB_tdt;
|
|
tab(15) = id_gijBB_t ;tab(16) = id_gijBB_tdt ;tab(17) = idMtdt ;
|
|
tab(18) = igradVBB_tdt; tab(19) = iVtdt; tab(20) = idVtdt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
// on ne calcul la dérivée de la déformation virtuelle qu'une fois
|
|
// car elle est constante dans le temps et indépendante des coordonnées
|
|
cald_Dvirtuelle=true;
|
|
};
|
|
// dimensionnement du residu
|
|
int nbddl = CoSfe->tab_ddl.NbDdl();
|
|
if ( residu == NULL)
|
|
residu = new Vecteur(nbddl); // cas du premier passage
|
|
else
|
|
for (int i =1; i<= nbddl; i++) // cas des autres passages
|
|
(*residu)(i) = 0.; // on initialise a zero
|
|
// dimensionnement de la raideur
|
|
if ( raideur == NULL)
|
|
raideur = new Mat_pleine(nbddl,nbddl); // cas du premier passage
|
|
else
|
|
for (int i =1; i<= nbddl; i++) // cas des autres passages
|
|
for (int j=1; j<= nbddl; j++) //
|
|
(*raideur)(i,j) = 0.; // on initialise a zero
|
|
double jacobien; // def du jacobien
|
|
Vecteur d_jacobien_tdt;
|
|
int indice = 1;
|
|
int nbtotsurf = CoSfe->eleCentre->Nbi();
|
|
int nbtotepais =(CoSfe->segment).Nbi();
|
|
// préparation du calcul de l'épaisseur finale (qui dépend de tous les points d'intégration)
|
|
double somme_epaisseurs_finales = 0;
|
|
// on utilise une épaisseur intermédiaire pour la métrique
|
|
// de manière à utiliser l'épaisseur de l'incrément précédent
|
|
Epai epais_pour_metrique;
|
|
//--debug
|
|
//if (num_elt == 89)
|
|
// { cout << "\n SfeMembT::Calcul_implicit debug ";
|
|
// cout << endl;
|
|
// };
|
|
|
|
// on doit boucler à l'extérieur sur les pt de surface et à l'intérieur sur les pt d'épaisseurs
|
|
// pour profiter du fait que une fois un pt calculer, on n'a pas à recalculer la courbure !!!!!
|
|
// ne sert plus** int nbint = nbtotsurf*nbtotepais; // nb total de pti
|
|
for (int nisur =1; nisur<= nbtotsurf; nisur++) // boucle sur les pt d'integ de surface
|
|
for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur
|
|
{if ((donnee_specif.epais != NULL) // on change les épaisseurs si elles sont stockées
|
|
&& (!Loi_rien(loiComp->Id_comport()))
|
|
)
|
|
{ epais_pour_metrique = *(donnee_specif.epais);
|
|
epais_pour_metrique.epaisseur_tdt = epais_pour_metrique.epaisseur_t;
|
|
// ((DeformationSfe1*) def)->Change_epaisseur(*donnee_specif.epais); // au niveau de l'élément
|
|
((DeformationSfe1*) def)->Change_epaisseur(epais_pour_metrique);
|
|
};
|
|
((DeformationSfe1*) def)->ChangeNumIntegSfe1(nisur,niepais);
|
|
((DeformationSfe1*) def)->Mise_a_jour_data_specif(tabSaveDefDon(indice));
|
|
|
|
|
|
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(indice);
|
|
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
|
|
TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
|
|
EnergieMeca& energ = tab_energ(indice);
|
|
EnergieMeca& energ_t = tab_energ_t(indice);
|
|
double poids = (CoSfe->eleCentre->Wi(nisur) * (CoSfe->segment).Wi(niepais)) * 0.5 ;
|
|
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
|
|
if (loiTP != NULL) {sTP = tabSaveTP(indice);}; // au pt d'integ si elles existes
|
|
if ((loiComp->Id_comport()==LOI_VIA_UMAT) // cas de l'utilisation d'une loi umat
|
|
|| (loiComp->Id_comport()==LOI_VIA_UMAT_CP))
|
|
((Loi_Umat*) loiComp)->Mise_a_jour_nbe_nbptinteg(this->Num_elt(),indice);
|
|
|
|
// -- on gère les exceptions éventuelles en mettant le bloc sous surveillance
|
|
const Met_abstraite::Impli* ex_pointe = NULL; // c'est pour construire ex
|
|
try // ce qui permet de donner les numéros d'éléments et de pti
|
|
{ ex_pointe = &loiComp->Cal_implicit(tabSaveDon(indice), *def,CoSfe->tab_ddl
|
|
,ptIntegMeca,d_epsBB,jacobien,d_jacobien_tdt
|
|
,d_sigHH,pa,sTP,loiTP,dilatation,energ,energ_t
|
|
,premier_calcul_meca_impli_expli);
|
|
}
|
|
catch (ErrNonConvergence_loiDeComportement excep)
|
|
// cas d'une d'une erreur survenue à cause d'une non convergence pour la résolution
|
|
// d'une loi de comportement incrémentale
|
|
{ cout << "\n erreur de loi de comportement element= " << this->Num_elt()
|
|
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
|
|
throw excep;
|
|
}
|
|
catch (ErrSortieFinale)
|
|
// cas d'une direction voulue vers la sortie
|
|
// on relance l'interuption pour le niveau supérieur
|
|
{ ErrSortieFinale toto;
|
|
throw (toto);
|
|
}
|
|
catch ( ... )
|
|
{ cout << "\n erreur inconnue de loi de comportement, element= " << this->Num_elt()
|
|
<< " point d'integration de surface= "<<nisur << " d'epaisseur= "<< niepais << endl;
|
|
throw (Err_inconnue_ElemMeca());
|
|
};
|
|
const Met_abstraite::Impli& ex=(*ex_pointe);
|
|
ex_final = ex_pointe; // stockage pour la stabilisation éventuelle de membraneBiel
|
|
|
|
// on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
|
|
Tableau2 <TenseurBB *>& d2_epsBB = (CoSfe->d2_epsBB)(indice); // pour simplifier
|
|
|
|
if (cald_Dvirtuelle)
|
|
def->Cal_var_def_virtuelle(false,d2_epsBB);
|
|
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
|
|
if (Bulk_visco_actif())
|
|
ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_tdt,(*sig_bulk_pourSfe_HH) );
|
|
double rap=1.; // pour le calcul éventuel du rapport de jacobien actuel
|
|
if (pa.MaxiVarJacobien() > 1.) // cas de la demande du calcul du rapport de jacobien
|
|
{if (Dabs(jacobien) > (*(ex.jacobien_0))) {rap = Dabs(jacobien) / (*(ex.jacobien_0));}
|
|
else {rap = (*(ex.jacobien_0)) / Dabs(jacobien) ;}
|
|
};
|
|
// on tient compte de la variation d'épaisseur pour ce point d'integ
|
|
double epaisseur_tdt = 0.; // ne sert que pour les éléments 2D avec épaisseurs stockées à l'élément
|
|
double epaisseur_pris_en_compte = 0.; // celle qui sera réellement prise en compte
|
|
if ((donnee_specif.epais != NULL) && (!Loi_rien(loiComp->Id_comport())) )
|
|
{ switch (loiComp->Comportement_3D_CP_DP_1D())
|
|
{ case COMP_CONTRAINTES_PLANES :
|
|
{
|
|
// en contrainte plane on recalcule l'épaisseur
|
|
CalEpaisseurAtdtImp_et_vol_pti(donnee_specif.epais->epaisseur0
|
|
,pa,donnee_specif.epais->epaisseur_t
|
|
,ptIntegMeca,epaisseur_tdt,ex);
|
|
//------ le 7 nov 2016: changement car je ne veux pas utiliser la variation d'épaisseur
|
|
// qui est calculée dans les lois. C'est sans doute possible, bien qu'avec les lois additives et mélange
|
|
// ce n'est pas sûr que cela veuille dire quelque chose. Donc on essaie de sans passer
|
|
// De plus c'est cohérent avec les éléments membranes ...
|
|
|
|
// ---- remodif le 1 avril 2017 car la méthode précédente fait que l'élément converge très mal !!
|
|
// je ne sais pas trop pourquoi (avec une loi de déformation plane, cela converge immédiatement !)
|
|
// // en fait, on récupère la variation d'épaisseur calculée via la condition de contrainte plane
|
|
// // on utilise uniquement la variation d'épaisseur calculée à la surface médiane ??? pour l'instant
|
|
// epaisseur_tdt = loiComp->HsurH0(tabSaveDon(nisur)) * donnee_specif.epais->epaisseur0;
|
|
// ************$la ligne qui suit est uniquement pour faire passer qq calcul ---> mais c'est faux !!!
|
|
// ************$la ligne qui suit est uniquement pour faire passer qq calcul ---> mais c'est faux !!!
|
|
// ************$la ligne qui suit est uniquement pour faire passer qq calcul ---> mais c'est faux !!!
|
|
|
|
// epaisseur_tdt = donnee_specif.epais->epaisseur0;
|
|
// on utilise l'épaisseur de l'incrément précédent
|
|
epaisseur_pris_en_compte = donnee_specif.epais->epaisseur_tdt;
|
|
break;
|
|
}
|
|
case COMP_DEFORMATIONS_PLANES :
|
|
// en deformation plane, l'épaisseur ne change pas
|
|
{ epaisseur_tdt = donnee_specif.epais->epaisseur0;
|
|
ptIntegMeca.Volume_pti() *= donnee_specif.epais->epaisseur0;
|
|
epaisseur_pris_en_compte = epaisseur_tdt;
|
|
break;
|
|
}
|
|
default :
|
|
cout << "\nErreur : cas de loi de comportement : " << Nom_comp((loiComp->Id_comport()))
|
|
<< " non prevu pour le calcul le l'epaisseur avec l'element SFE ";
|
|
cout << "\n SfeMembT::Calcul_implicit(.... \n";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
somme_epaisseurs_finales += epaisseur_tdt; //pour le calcul final
|
|
|
|
//--------- pour le débug --------
|
|
// Deformation & defS = *defSurf(1); // pour simplifier l'écriture
|
|
// int nbddls = ddlS.NbDdl();
|
|
|
|
|
|
//------- pour le débug fin -------
|
|
|
|
|
|
if ((jacobien <= 0.)|| (std::isinf(jacobien))||( std::isnan(jacobien))) // vérif du jacobien
|
|
{ if ((std::isinf(jacobien))||( std::isnan(jacobien)))
|
|
// on met un message quelque soit le niveau d'impression
|
|
{ cout << "\n ********** attention on a trouve un jacobien infini ou nan !! = ("<<jacobien<<")********* " << endl; };
|
|
if (ParaGlob::NiveauImpression() >= 1)
|
|
{ cout << "\n ********** attention on a trouve un jacobien negatif!! ("<<jacobien<<")********* " << endl; };
|
|
};
|
|
if ((jacobien <= 0.) && (pa.JabobienNegatif()!= 0)) // vérif du jacobien et traitement adéquate si besoin
|
|
{ cout << "\n ParaGlob::NiveauImpression() = "<<ParaGlob::NiveauImpression() << " jacobien= " << jacobien << endl;
|
|
switch (pa.JabobienNegatif())
|
|
{ case 1:
|
|
{ cout << "\n on annulle sa contribution. Element nb: " << Num_elt() << " ni_surf: " << nisur
|
|
<< " ni_ep: " << niepais;
|
|
break;
|
|
}
|
|
case 2:
|
|
{ // on génère une exception
|
|
throw ErrJacobienNegatif_ElemMeca();break;
|
|
}
|
|
// dans les autres cas on ne fait rien
|
|
};
|
|
}
|
|
else if ((pa.MaxiVarJacobien() > 1.) && ( rap > pa.MaxiVarJacobien()))
|
|
{ if (ParaGlob::NiveauImpression() >= 3)
|
|
{ cout << "\n *** attention la variation maximale du jacobien est atteinte !! *** "
|
|
<< "\n ele= "<< Num_elt() << " ni_surf: " << nisur
|
|
<< " ni_ep: " << niepais << " jacobien= " << jacobien << " jacobien_0= "
|
|
<< (*(ex.jacobien_0)) << endl;
|
|
};
|
|
// on génère une exception
|
|
throw ErrVarJacobienMini_ElemMeca();break;
|
|
}
|
|
else
|
|
{ // on prend en compte toutes les variations
|
|
|
|
// dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
|
|
// épaisseur_tdt ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
|
|
// pour éviter des tests et une recopie du code on laisse epaisseur_tdt, avec une valeur 1
|
|
// ce qui fait qu'elle n'intervient pas
|
|
if (donnee_specif.epais == NULL) epaisseur_pris_en_compte = 1.;
|
|
|
|
// on continue que s'il s'agit d'une loi différente de rien
|
|
if (!Loi_rien(loiComp->Id_comport()))
|
|
{// il ne faut pas interpoler l'énergie à t, car elle était associée à un volume qui a changé
|
|
// donc il ne faut considérer que l'accroissement d'énergie
|
|
energie_totale.Ajout_differenciee(energ,energ_t,(poids * jacobien * epaisseur_pris_en_compte));
|
|
// energie_totale += (energ-energ_t) * (poids * jacobien * epaisseur_tdt);
|
|
// energie_totale += energ * (poids * jacobien * epaisseur_tdt);
|
|
// l'ordre des ddl vient de la métrique: tous le xi, puis les épaisseurs, par contre dans le stockage
|
|
// local on doit avoir noeud1: ddlxi puis epais, noeud2: idem etc pour les noeuds centraux, puis ddl xi
|
|
// pour les noeuds externes. Si pas de ddl epais, tous les noeuds sont traités de la même manière
|
|
// on utilise un adressage indirecte
|
|
////--- debug
|
|
//if (Dabs(epaisseur_pris_en_compte - donnee_specif.epais->epaisseur0) > 0.01)
|
|
// {cout << "\n debug SfeMembT::Calcul_implicit: epai= "
|
|
// << epaisseur_pris_en_compte << " epai_0= "<< donnee_specif.epais->epaisseur0 << endl;
|
|
// };
|
|
//// --- fin debug
|
|
|
|
for (int j =1; j<= nbddl; j++) // 1ere boucle sur les ddl
|
|
{ int jloc = CoSfe->nMetVTab_ddl(j); // adressage indirecte
|
|
(*residu)(jloc) -= ((*(d_epsBB(j))) && sigHH ) * poids * jacobien* epaisseur_pris_en_compte;
|
|
/*//--debug
|
|
if ((*residu)(jloc) > 10.)
|
|
{ cout << "\n SfeMembT::Calcul_implicit debug ";
|
|
cout << "\n (*ex.giB_0)(1): " ;(*ex.giB_0)(1).Affiche();
|
|
cout << " (*ex.giB_0)(2): " ;(*ex.giB_0)(2).Affiche();
|
|
cout << " (*ex.giH_0)(1): " ;(*ex.giH_0)(1).Affiche();
|
|
cout << " (*ex.giH_0)(2): " ;(*ex.giH_0)(2).Affiche();
|
|
cout << "\n (*ex.giBB_0): " ;(*ex.gijBB_0).Ecriture(cout);
|
|
cout << "\n (*ex.giHH_0): " ;(*ex.gijHH_0).Ecriture(cout);
|
|
|
|
cout << "\n (*(d_epsBB(j))): " ;(*(d_epsBB(j))).Ecriture(cout);
|
|
cout << "\n sigHH: " ;sigHH.Ecriture(cout);
|
|
cout << endl;
|
|
|
|
ex_pointe = &loiComp->Cal_implicit(tabSaveDon(indice), *def,CoSfe->tab_ddl
|
|
,ptIntegMeca,d_epsBB,jacobien,d_jacobien_tdt
|
|
,d_sigHH,pa,sTP,loiTP,dilatation,energ,energ_t
|
|
,premier_calcul_meca_impli_expli);
|
|
|
|
Sortie(1);
|
|
};
|
|
//--fin debug */
|
|
//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
|
|
//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
|
|
Delta_Jacobien_anormal(tabSaveDefDon(indice),nisur,niepais);
|
|
|
|
volume += jacobien * poids * epaisseur_tdt; // là on prend la vraie épaisseur
|
|
for (int i =1; i<= nbddl; i++) // 2ere boucle sur les ddl
|
|
{ int iloc = CoSfe->nMetVTab_ddl(i); // adressage indirecte
|
|
(*raideur)(jloc,iloc) += ( ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien
|
|
+ (sigHH && (*((CoSfe->d2_epsBB)(indice)(j,i)))) * jacobien
|
|
+ (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i)
|
|
) * poids * epaisseur_pris_en_compte;
|
|
|
|
// on utilise toujour la variation de D, donc c'est maintenant directement intégré dans l'expression précédente, et la suite est commentée
|
|
// // on intègre la variation seconde si c'est demandé
|
|
// if (pa.Var_D())
|
|
// (*raideur)(jloc,iloc) += (sigHH && (*((CoSfe->d2_epsBB)(indice)(j,i))))
|
|
// * jacobien * poids * epaisseur_tdt;
|
|
// ---- pour test débug
|
|
// if (i < 10)
|
|
// (*raideur)(j,i) += (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i) * poids * donnee_specif.epaisseur_tdt;
|
|
};
|
|
};
|
|
};
|
|
// --- cas éventuel d'une stabilisation membrane-biel ---
|
|
// ici il s'agit de la contribution précise à chaque pti
|
|
if (pt_StabMembBiel != NULL)
|
|
if (!(pt_StabMembBiel->Aa_calculer()))
|
|
{ double poids_vol = poids * epaisseur_pris_en_compte;
|
|
Cal_implicit_StabMembBiel
|
|
(nisur,*ex_final,nbtotsurf,poids_vol,unefois->doCoMemb->noeud_a_prendre_en_compte);
|
|
};
|
|
};// fin du if sur la valeur non négative du jacobien
|
|
////--debug
|
|
//for (int i = 1; i<= 18; i++)
|
|
// { if ((*residu)(i) != 0. )
|
|
// { cout << "\n SfeMembT::Calcul_implicit debug ";
|
|
// cout << endl;
|
|
// };
|
|
////if (num_elt == 176)
|
|
//// { cout << "\n SfeMembT::Calcul_implicit debug ";
|
|
//// cout << endl;
|
|
// };
|
|
////--fin debug
|
|
if (ElemMeca::Bulk_visco_actif()) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
|
|
{ sigHH -= (*sig_bulk_pourSfe_HH); } // ceci pour la sauvegarde ou autre utilisation
|
|
|
|
}; // fin boucle sur les points d'intégration
|
|
#ifdef MISE_AU_POINT
|
|
if (ParaGlob::NiveauImpression() > 9)
|
|
{ cout << "\n Raideur et second membre locaux: ? (o/n) "; string rep(" ");
|
|
// procédure de lecture avec prise en charge d'un retour chariot
|
|
rep = lect_return_defaut(false,"n");
|
|
|
|
if ((rep == "0")||(rep == "o"))
|
|
{cout << "\n Raideur et second membre locaux: elem= " << Num_elt_const()
|
|
<< ", maillage= " << Num_maillage();
|
|
cout << "\n raideur: ";
|
|
raideur->Affiche();
|
|
cout << "\n second membre: " << (*residu);
|
|
};
|
|
};
|
|
#endif
|
|
// --- intervention dans le cas d'une stabilisation d'hourglass
|
|
// stabilisation pour un calcul implicit
|
|
if (type_stabHourglass)
|
|
Cal_implicit_hourglass();
|
|
// --- cas éventuel d'une stabilisation membrane-biel ---
|
|
// ici il s'agit soit du calcul approché d'initialisation, soit de la fin du calcul après la boucle
|
|
// modif: 10 janvier 2019 non c'est le calcul correct une fois la raideur calculée
|
|
if (pt_StabMembBiel != NULL)
|
|
{if (pt_StabMembBiel->Aa_calculer())
|
|
{Cal_implicit_StabMembBiel
|
|
(0,*ex_final, nbtotsurf,volume,unefois->doCoMemb->noeud_a_prendre_en_compte);}
|
|
else {Cal_implicit_StabMembBiel
|
|
(-1,*ex_final, nbtotsurf,volume,unefois->doCoMemb->noeud_a_prendre_en_compte);} ;
|
|
};
|
|
#ifdef MISE_AU_POINT
|
|
if (ParaGlob::NiveauImpression() > 9)
|
|
{ if ((type_stabHourglass)||(pt_StabMembBiel != NULL))
|
|
{cout << "\n apres stabilisation: Raideur et second membre locaux: ? (o/n) "; string rep(" ");
|
|
// procédure de lecture avec prise en charge d'un retour chariot
|
|
rep = lect_return_defaut(false,"n");
|
|
|
|
if ((rep == "0")||(rep == "o"))
|
|
{cout << "\n Raideur et second membre locaux: elem= " << Num_elt_const()
|
|
<< ", maillage= " << Num_maillage();
|
|
cout << "\n raideur: ";
|
|
raideur->Affiche();
|
|
cout << "\n second membre: " << (*residu);
|
|
};
|
|
};
|
|
};
|
|
#endif
|
|
|
|
////--debug
|
|
//for (int i = 1; i<= 18; i++)
|
|
// { if ((*residu)(i) != 0. )
|
|
// { cout << "\n SfeMembT::Calcul_implicit debug ";
|
|
// cout << endl;
|
|
// };
|
|
////if (num_elt == 176)
|
|
//// { cout << "\n SfeMembT::Calcul_implicit debug ";
|
|
//// cout << endl;
|
|
// };
|
|
////--fin debug
|
|
|
|
|
|
// on calcul éventuellement l'épaisseur moyenne finale si c'est un élément 2D
|
|
if ((donnee_specif.epais != NULL) && (!Loi_rien(loiComp->Id_comport())) )// cas 2D
|
|
donnee_specif.epais->epaisseur_tdt = somme_epaisseurs_finales / (indice - 1);
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = residu;
|
|
el.raid = raideur;
|
|
return el;
|
|
};
|
|
|
|
// Calcul de la matrice masse pour l'élément
|
|
// celle-ci est calculée uniquement pour l'élément centrale et en considérant
|
|
// que la masse volumique est constante dans l'épaisseur
|
|
Mat_pleine * SfeMembT::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse)
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
// dimensionement de la métrique si nécessaire
|
|
if (!(unefois->CalDynamique ))
|
|
{ unefois->CalDynamique += 1;
|
|
Tableau<Enum_variable_metrique> tab(5);
|
|
tab(1) = igiB_0; tab(2) = igiB_tdt; tab(3) = igijBB_0;
|
|
tab(4) = igijBB_tdt; tab(5) = igradVmoyBB_t;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
// on vérifie le bon dimensionnement de la matrice
|
|
if (type_calcul_masse == MASSE_CONSISTANTE)
|
|
// dans le cas où la masse est consistante il faut la redimensionner
|
|
{ int nbddl = CoSfe->tab_ddl.NbDdl();
|
|
(CoSfe->matrice_masse).Initialise (nbddl,nbddl,0.);
|
|
}
|
|
};
|
|
if ((donnee_specif.epais != NULL) // on change les épaisseurs si elles sont stockées
|
|
&& (!Loi_rien(loiComp->Id_comport()))
|
|
)
|
|
((DeformationSfe1*) def)->Change_epaisseur(*donnee_specif.epais); // au niveau de l'élément
|
|
// appel de la routine générale
|
|
////************ non, ce n'est pas correcte pour les éléments 3D, il faut revoir le calcul de la masse
|
|
//// en utilisant d(OM)/d_ddl au lieu de phi^r
|
|
ElemMeca::Cal_Mat_masse (CoSfe->tab_ddl,type_calcul_masse,
|
|
nombre->nbisMas,CoSfe->eleMas->TaPhi(),nombre->nbnce
|
|
,CoSfe->eleMas->TaWi());
|
|
if ((donnee_specif.epais != NULL) // prise en compte de l'épaisseur finale, si elle est stockée
|
|
&& (!Loi_rien(loiComp->Id_comport()))
|
|
)
|
|
(*mat_masse) *= donnee_specif.epais->epaisseur_tdt; // au niveau de l'élément, sinon on est en 3D, donc c'est
|
|
// intégré dans le jacobien
|
|
return mat_masse;
|
|
};
|
|
|
|
//============= lecture écriture dans base info ==========
|
|
|
|
// cas donne le niveau de la récupération
|
|
// = 1 : on récupère tout
|
|
// = 2 : on récupère uniquement les données variables (supposées comme telles)
|
|
void SfeMembT::Lecture_base_info
|
|
(ifstream& ent,const Tableau<Noeud *> * tabMaillageNoeud,const int cas)
|
|
{// tout d'abord appel de la lecture de la classe elem_meca
|
|
ElemMeca::Lecture_bas_inf(ent,tabMaillageNoeud,cas);
|
|
// traitement du cas particulier de l'élément centrale
|
|
string nom;
|
|
switch (cas)
|
|
{ case 1 : // ------- on récupère tout -------------------------
|
|
{ // construction du tableau de ddl des noeuds de l'élément centrale
|
|
ConstTabDdl();
|
|
// récup contraintes et déformation
|
|
lesPtMecaInt.Lecture_base_info(ent,cas);
|
|
|
|
// les données spécifiques
|
|
ent >> nom;
|
|
if (nom == "epaisStockeDansElement")
|
|
{ // cas où les épaisseurs sont stockées dans l'élément
|
|
if (donnee_specif.epais == NULL) donnee_specif.epais = new Epai;
|
|
ent >> nom >> donnee_specif.epais->epaisseur0
|
|
>> nom >> donnee_specif.epais->epaisseur_t
|
|
>> nom >> donnee_specif.epais->epaisseur_tdt;
|
|
}
|
|
else
|
|
// cas d'un stockage au niveau des noeuds
|
|
{ if (donnee_specif.epais != NULL) delete donnee_specif.epais;};
|
|
// CL éventuelles :
|
|
ent >> nom ;
|
|
bool existe_CL; ent >> existe_CL;
|
|
if (existe_CL)
|
|
{if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL)))
|
|
{ // si les tableaux n'existent pas on les crée
|
|
if (areteTypeCL == (&(unefois->doCoMemb->tabType_rienCL)))
|
|
areteTypeCL = new Tableau <EnuTypeCL> (3);
|
|
if (vplan == (&(unefois->doCoMemb->vplan_rien)))
|
|
vplan = new Tableau <Coordonnee3> (3);
|
|
// et on récupère les infos
|
|
ent >> nom; (*areteTypeCL)(1) = Id_nomTypeCL(nom); ent >> (*vplan)(1);
|
|
ent >> nom; (*areteTypeCL)(2) = Id_nomTypeCL(nom); ent >> (*vplan)(2);
|
|
ent >> nom; (*areteTypeCL)(3) = Id_nomTypeCL(nom); ent >> (*vplan)(3);
|
|
}
|
|
}
|
|
else
|
|
{// on supprime les tableaux locaux éventuels et on pointe vers les valeurs par défaut
|
|
if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL))) delete areteTypeCL;
|
|
areteTypeCL = (&(unefois->doCoMemb->tabType_rienCL));
|
|
if (vplan != (&(unefois->doCoMemb->vplan_rien))) delete vplan;
|
|
vplan = (&(unefois->doCoMemb->vplan_rien));
|
|
};
|
|
break;
|
|
}
|
|
case 2 : // ----------- lecture uniquement de se qui varie --------------------
|
|
{ // récup contraintes et déformation
|
|
lesPtMecaInt.Lecture_base_info(ent,cas);
|
|
// les données spécifiques
|
|
bool avecEpaisDansElement = true;
|
|
ent >> avecEpaisDansElement;
|
|
if (avecEpaisDansElement)
|
|
{ if (donnee_specif.epais == NULL) donnee_specif.epais = new Epai;
|
|
ent >> nom >> donnee_specif.epais->epaisseur_tdt;
|
|
donnee_specif.epais->epaisseur_t = donnee_specif.epais->epaisseur_tdt;
|
|
};
|
|
|
|
// CL éventuelles : on ne les lit pas car elles sont généré éventuellement
|
|
// à chaque fois que l'on met les CL, sinon elles ne changent pas
|
|
break;
|
|
}
|
|
default :
|
|
{ cout << "\nErreur : valeur incorrecte du type de lecture !\n";
|
|
cout << "SfeMembT::::Lecture_base_info(ifstream& ent,const "
|
|
<< "Tableau<Noeud *> * tabMaillageNoeud,const int cas)"
|
|
<< " cas= " << cas << endl;
|
|
Sortie(1);
|
|
}
|
|
}
|
|
};
|
|
|
|
// cas donne le niveau de sauvegarde
|
|
// = 1 : on sauvegarde tout
|
|
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
|
|
void SfeMembT::Ecriture_base_info(ofstream& sort,const int cas)
|
|
{// tout d'abord appel de l'écriture de la classe elem_meca
|
|
ElemMeca::Ecriture_bas_inf(sort,cas);
|
|
// traitement du cas particulier de l'élément central
|
|
switch (cas)
|
|
{ case 1 : // ------- on sauvegarde tout -------------------------
|
|
{
|
|
// des tenseurs déformation et contrainte,
|
|
lesPtMecaInt.Ecriture_base_info(sort,cas);
|
|
// les données spécifiques
|
|
if (donnee_specif.epais != NULL)
|
|
{ sort << "\n epaisStockeDansElement " << " epaisseur0= " << donnee_specif.epais->epaisseur0
|
|
<< " epaisseur_t= " << donnee_specif.epais->epaisseur_t
|
|
<< " epaisseur_tdt= " << donnee_specif.epais->epaisseur_tdt;
|
|
}
|
|
else { sort << "\n epaisStockeAuxNoeuds "; };
|
|
// CL éventuelles :
|
|
sort << " CL= ";
|
|
if (areteTypeCL != (&(unefois->doCoMemb->tabType_rienCL)))
|
|
{ // si ça ne pointe pas sur la valeur par défaut, alors c'est qu'il y a
|
|
// des CL particulières,
|
|
sort << " 1 ";
|
|
sort << NomTypeCL((*areteTypeCL)(1)) << " " << (*vplan)(1)
|
|
<< NomTypeCL((*areteTypeCL)(2)) << " " << (*vplan)(2)
|
|
<< NomTypeCL((*areteTypeCL)(3)) << " " << (*vplan)(3);
|
|
}
|
|
else
|
|
{sort << " 0 ";
|
|
};
|
|
|
|
break;
|
|
}
|
|
case 2 : // ----------- sauvegarde uniquement de se qui varie --------------------
|
|
{ // des tenseurs déformation et contrainte,
|
|
lesPtMecaInt.Ecriture_base_info(sort,cas);
|
|
// les données spécifiques
|
|
if (donnee_specif.epais != NULL)
|
|
{ sort << "\n 1 epaisseur_tdt= " << donnee_specif.epais->epaisseur_tdt;}
|
|
else
|
|
{ sort << "\n 0 ";};
|
|
// CL éventuelles : on ne les sort pas car elles sont généré éventuellement
|
|
// à chaque fois que l'on met les CL, sinon elles ne changent pas
|
|
break;
|
|
}
|
|
default :
|
|
{ cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
|
|
cout << "SfeMembT::::Ecriture_base_info(ofstream& sort,const int cas)"
|
|
<< " cas= " << cas << endl;
|
|
Sortie(1);
|
|
}
|
|
};
|
|
};
|
|
|
|
//------- calcul d'erreur, remontée des contraintes -------------------
|
|
|
|
// calcul du résidu et de la matrice de raideur pour le calcul d'erreur
|
|
// en fait on ne calculera l'erreur que sur la contrainte sur toute la surface
|
|
// mais sur le premier pt d'integ d'épaisseur
|
|
Element::Er_ResRaid SfeMembT::ContrainteAuNoeud_ResRaid()
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
// dimensionnement de la metrique
|
|
if(!( unefois->CalResPrem_t ))
|
|
{unefois->CalResPrem_t += 1;
|
|
Tableau<Enum_variable_metrique> tab(10);
|
|
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
|
|
tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
|
|
tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
};
|
|
// appel du programme général
|
|
int tabn_taille = tab_noeud.Taille();
|
|
ElemMeca:: SigmaAuNoeud_ResRaid(tabn_taille,
|
|
CoSfe->eleCentre->TaPhi(),
|
|
CoSfe->eleCentre->TaWi(),
|
|
CoSfe-> resErr,CoSfe->raidErr,
|
|
CoSfe->eleEr->TaPhi(),
|
|
CoSfe->eleEr->TaWi());
|
|
|
|
// --- on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
|
|
// dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
|
|
// épaisseur_t ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
|
|
// pour éviter des tests et une recopie du code on laisse epaisseur_t, avec une valeur 1
|
|
// ce qui fait qu'elle n'intervient pas
|
|
double epaisseur = 0.; // init
|
|
if (donnee_specif.epais == NULL) {epaisseur = 1.;}
|
|
else {epaisseur = donnee_specif.epais->epaisseur_t;};
|
|
// on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
|
|
int taille_resErr=CoSfe->resErr.Taille();
|
|
for (int i=1;i<= taille_resErr;i++)
|
|
(*CoSfe-> resErr(i)) *= epaisseur;
|
|
CoSfe->raidErr *= epaisseur;
|
|
|
|
return (Element::Er_ResRaid( &(CoSfe-> resErr),&(CoSfe->raidErr)));
|
|
} ;
|
|
|
|
// calcul du résidu pour le calcul d'erreur ramenée aux noeuds
|
|
Element::Er_ResRaid SfeMembT::ErreurAuNoeud_ResRaid()
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
// dimensionnement de la metrique
|
|
if(!( unefois->CalResPrem_t ))
|
|
{ unefois->CalResPrem_t += 1;
|
|
Tableau<Enum_variable_metrique> tab(10);
|
|
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
|
|
tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
|
|
tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
};
|
|
// appel du programme général
|
|
int tabn_taille = tab_noeud.Taille();
|
|
ElemMeca::Cal_ErrAuxNoeuds(tabn_taille,
|
|
CoSfe->eleCentre->TaPhi(),CoSfe->eleCentre->TaWi()
|
|
,CoSfe-> resErr );
|
|
|
|
// --- on tient compte de l'épaisseur à t, supposée déjà calculée ou récupérée
|
|
// dans le cas d'un élément 3D, épaisseur est stockée aux noeuds, donc ici la variable
|
|
// épaisseur_t ne devrait pas exister, en particulier le jacobien intègre déjà l'épaisseur
|
|
// pour éviter des tests et une recopie du code on laisse epaisseur_t, avec une valeur 1
|
|
// ce qui fait qu'elle n'intervient pas
|
|
double epaisseur = 0.; // init
|
|
if (donnee_specif.epais == NULL) {epaisseur = 1.;}
|
|
else {epaisseur = donnee_specif.epais->epaisseur_t;};
|
|
int taille_resErr=CoSfe->resErr.Taille();
|
|
for (int i=1;i<= taille_resErr;i++)
|
|
(*CoSfe-> resErr(i)) *= epaisseur;
|
|
CoSfe->raidErr *= epaisseur;
|
|
|
|
return (Element::Er_ResRaid( &(CoSfe-> resErr),&(CoSfe->raidErr)));
|
|
} ;
|
|
|
|
// actualisation des ddl et des grandeurs actives de t+dt vers t
|
|
void SfeMembT::TdtversT()
|
|
{ lesPtMecaInt.TdtversT(); // contrainte
|
|
if (donnee_specif.epais != NULL)
|
|
donnee_specif.epais->epaisseur_t = donnee_specif.epais->epaisseur_tdt;
|
|
int taille = lesPtMecaInt.NbPti();
|
|
for (int ni=1;ni<= taille; ni++)
|
|
{ if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->TdtversT();
|
|
if (tabSaveTP(ni) != NULL) tabSaveTP(ni)->TdtversT();
|
|
if (tabSaveDefDon(ni) != NULL) tabSaveDefDon(ni)->TdtversT();
|
|
}
|
|
ElemMeca::TdtversT_(); // appel de la procédure mère
|
|
};
|
|
// actualisation des ddl et des grandeurs actives de t vers tdt
|
|
void SfeMembT::TversTdt()
|
|
{ lesPtMecaInt.TversTdt(); // contrainte
|
|
if (donnee_specif.epais != NULL)
|
|
donnee_specif.epais->epaisseur_tdt = donnee_specif.epais->epaisseur_t;
|
|
int taille = lesPtMecaInt.NbPti();
|
|
for (int ni=1;ni<= taille; ni++)
|
|
{ if (tabSaveDon(ni) != NULL) tabSaveDon(ni)->TversTdt();
|
|
if (tabSaveTP(ni) != NULL) tabSaveTP(ni)->TversTdt();
|
|
if (tabSaveDefDon(ni) != NULL) tabSaveDefDon(ni)->TversTdt();
|
|
}
|
|
ElemMeca::TversTdt_(); // appel de la procédure mère
|
|
};
|
|
|
|
// calcul de l'erreur sur l'élément.
|
|
// en fait on ne calculera l'erreur que sur la contrainte sur toute la surface
|
|
// mais sur le premier pt d'integ d'épaisseur
|
|
void SfeMembT::ErreurElement(int type,double& errElemRelative
|
|
,double& numerateur, double& denominateur)
|
|
{ SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
|
|
// dimensionnement de la metrique
|
|
if(!( unefois->CalResPrem_t ))
|
|
{ unefois->CalResPrem_t += 1;
|
|
Tableau<Enum_variable_metrique> tab(10);
|
|
tab(1) = igiB_0; tab(2) = igiB_t; tab(3) = igijBB_0;tab(4) = igijBB_t;
|
|
tab(5) = igijHH_t; tab(6) = id_giB_t; tab(7) = id_gijBB_t ;
|
|
tab(8) = id_gijBB_t ;tab(9) = igradVBB_t; tab(10) = iVt;
|
|
CoSfe->met_SfeMembT.PlusInitVariables(tab) ;
|
|
};
|
|
// appel du programme général
|
|
ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur,
|
|
tab_noeud.Taille(),CoSfe->eleCentre->TaPhi(),
|
|
CoSfe->eleCentre->TaWi(),
|
|
CoSfe->eleEr->TaPhi(),CoSfe->eleEr->TaWi());
|
|
} ;
|
|
|
|
// mise à jour de la boite d'encombrement de l'élément, suivant les axes I_a globales
|
|
// en retour coordonnées du point mini dans retour.Premier() et du point maxi dans .Second()
|
|
// la méthode est différente de la méthode générale car il faut prendre en compte l'épaisseur de l'élément
|
|
const DeuxCoordonnees& SfeMembT::Boite_encombre_element(Enum_dure temps)
|
|
{ int tab_taille= tab_noeud.Taille();
|
|
// on commence par calculer la boite d'encombrement pour l'élément médian
|
|
Element::Boite_encombre_element( temps);
|
|
// ensuite on augment sytématiquement dans toutes directions d'une valeur h/2 majorée
|
|
// h étant l'épaisseur maxi constatée
|
|
double epais_maxi=0.; // init
|
|
if (donnee_specif.epais != NULL)
|
|
{ // cas où les épaisseurs sont stockées sur l'élément
|
|
epais_maxi = donnee_specif.epais->epaisseur_tdt;
|
|
}
|
|
else
|
|
{ // cas où les épaisseurs sont stockées aux noeuds
|
|
for (int ine=1;ine<=nombre->nbnce;ine++) // balayage des noeuds
|
|
epais_maxi = MaX(epais_maxi,tab_noeud(ine)->Valeur_tdt(EPAIS));
|
|
};
|
|
|
|
double hSur2maj = epais_maxi * 0.5
|
|
* (ParaGlob::param->ParaAlgoControleActifs().Extra_boite_prelocalisation()-1.);
|
|
// ajout d'un extra dans toutes les directions
|
|
hSur2maj += ParaGlob::param->ParaAlgoControleActifs().Ajout_extra_boite_prelocalisation();
|
|
// mise à jour
|
|
boite_encombre.Premier().Ajout_meme_valeur(-hSur2maj); // le min
|
|
boite_encombre.Second().Ajout_meme_valeur(hSur2maj); // le max
|
|
// retour
|
|
return boite_encombre;
|
|
};
|
|
|
|
// ramène le nombre de points d'intégration de surface correspondant à un type énuméré
|
|
// ramène 0 si l'élément n'est pas une plaque ou coque
|
|
int SfeMembT::NbPtIntegSurface(Enum_ddl ddl ) const
|
|
{ switch (ddl)
|
|
{ case ERREUR : return nombre->nbisEr; break;
|
|
default : return nombre->nbis; break;
|
|
};
|
|
};
|
|
// ramène le nombre de points d'intégration en épaisseur correspondant à un type énuméré
|
|
// ramène 0 si l'élément n'est pas une plaque ou coque
|
|
int SfeMembT::NbPtIntegEpaiss(Enum_ddl ddl) const
|
|
{ switch (ddl)
|
|
{ case ERREUR : return nombre->nbieEr; break;
|
|
default : return nombre->nbie; break;
|
|
};
|
|
};
|
|
// ramene l'element geometrique de surface correspondant au ddl passé en paramètre
|
|
// ou null si ce n'est pas définie, dans ce cas si l'élément géométrique de surface est 2D
|
|
// cela signifie qu'il faut se référer à l'élément générique: ElementGeometrie(...
|
|
ElemGeomC0* SfeMembT::ElementGeometrieSurface(Enum_ddl ddl ) const
|
|
{ switch (ddl)
|
|
{ case ERREUR : return (unefois->doCoMemb->eleEr); break;
|
|
default : return (unefois->doCoMemb->eleCentre); break;
|
|
};
|
|
};
|
|
// ramene l'element geometrique d'épaisseur correspondant au ddl passé en paramètre
|
|
// ou null si ce n'est pas définie, dans ce cas si l'élément géométrique de surface est 2D
|
|
// cela signifie que tout est constant (identique) dans l'épaisseur
|
|
ElemGeomC0* SfeMembT::ElementGeometrieEpaiss(Enum_ddl ddl ) const
|
|
{ switch (ddl)
|
|
{ case ERREUR : return &(unefois->doCoMemb->segmentEr); break;
|
|
default : return &(unefois->doCoMemb->segment); break;
|
|
};
|
|
};
|
|
|
|
// retourne un numero d'ordre d'un point le plus près ou est exprimé la grandeur enum
|
|
// par exemple un point d'intégration, mais n'est utilisable qu'avec des méthodes particulières
|
|
// par exemple CoordPtInteg, ou Valeur_a_diff_temps
|
|
// car le numéro d'ordre peut-être différent du numéro d'intégration au sens classique
|
|
// temps: dit si c'est à 0 ou t ou tdt
|
|
int SfeMembT::PointLePlusPres(Enum_dure temps,Enum_ddl enu, const Coordonnee& M)
|
|
{ int iret;
|
|
// récupération des éléments géométriques correspondant à Enu
|
|
ElemGeomC0* elSur = this->ElementGeometrieSurface(enu);
|
|
ElemGeomC0* elEpa = this->ElementGeometrieEpaiss(enu);
|
|
// récupération des tableaux des fonctions d'interpolations
|
|
const Tableau <Vecteur> & tabphiS = elSur->TaPhi();
|
|
const Tableau < Mat_pleine > & tabDphiS = elSur->TaDphi();
|
|
const Tableau <Vecteur> & tabphiH = elEpa->TaPhi();
|
|
|
|
// on boucle sur les pt d'integ pour avoir le point le plus près
|
|
int nbptSur = tabphiS.Taille();
|
|
int nbptEpa = tabphiH.Taille();
|
|
|
|
Coordonnee P; iret=1; double dist= ConstMath::tresgrand;
|
|
|
|
int ipt=1;
|
|
for (int iptH = 1;iptH <= nbptEpa;iptH++)
|
|
for (int iptS = 1; iptS <= nbptSur;iptS++,ipt++)
|
|
{ switch (temps)
|
|
{ case TEMPS_0 : P = ((Met_Sfe1*) met)->PointSfe1M_0
|
|
(tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
|
|
case TEMPS_t : P = ((Met_Sfe1*) met)->PointSfe1M_t
|
|
(tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
|
|
case TEMPS_tdt : P = ((Met_Sfe1*) met)->PointSfe1M_tdt
|
|
(tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
|
|
}
|
|
double di=(M-P).Norme();
|
|
if (di < dist) { dist = di; iret = ipt;};
|
|
}
|
|
return iret; // retour
|
|
};
|
|
|
|
// recuperation des coordonnées du point de numéro d'ordre iteg pour
|
|
// la grandeur enu
|
|
// temps: dit si c'est à 0 ou t ou tdt
|
|
// si erreur retourne erreur à true
|
|
Coordonnee SfeMembT::CoordPtInteg(Enum_dure temps,Enum_ddl enu,int iteg,bool& erreur)
|
|
{ Coordonnee ptret; // le retour
|
|
// récupération des éléments géométriques correspondant à Enu
|
|
ElemGeomC0* elSur = this->ElementGeometrieSurface(enu);
|
|
ElemGeomC0* elEpa = this->ElementGeometrieEpaiss(enu);
|
|
// récupération des tableaux des fonctions d'interpolations
|
|
const Tableau <Vecteur> & tabphiS = elSur->TaPhi();
|
|
const Tableau < Mat_pleine > & tabDphiS = elSur->TaDphi();
|
|
const Tableau <Vecteur> & tabphiH = elEpa->TaPhi();
|
|
// def du numéro d'integ en surface et épaisseur
|
|
int nbptSur = tabphiS.Taille();
|
|
int nbptEpa = tabphiH.Taille();
|
|
int iptS = iteg / nbptEpa +1; // division entière
|
|
int iptH = iteg - (iptS-1) * nbptEpa ;
|
|
// calcul du point
|
|
if ((iteg < 0) || (iteg > (nbptSur*nbptEpa)))
|
|
{ erreur = false;}
|
|
else
|
|
{ switch (temps)
|
|
{ case TEMPS_0 : ptret = ((Met_Sfe1*) met)->PointSfe1M_0
|
|
(tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
|
|
case TEMPS_t : ptret = ((Met_Sfe1*) met)->PointSfe1M_t
|
|
(tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
|
|
case TEMPS_tdt : ptret = ((Met_Sfe1*) met)->PointSfe1M_tdt
|
|
(tab_noeud,tabphiS(iptS),tabDphiS(iptS),donnee_specif.epais,tabphiH(iptH)); break;
|
|
}
|
|
}
|
|
return ptret; // retour
|
|
};
|