1285 lines
58 KiB
C++
1285 lines
58 KiB
C++
// FICHIER : ElemPoint.cc
|
|
// CLASSE : ElemPoint
|
|
|
|
// This file is part of the Herezh++ application.
|
|
//
|
|
// The finite element software Herezh++ is dedicated to the field
|
|
// of mechanics for large transformations of solid structures.
|
|
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
|
|
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
|
|
//
|
|
// Herezh++ is distributed under GPL 3 license ou ultérieure.
|
|
//
|
|
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
|
|
// AUTHOR : Gérard Rio
|
|
// E-MAIL : gerardrio56@free.fr
|
|
//
|
|
// This program is free software: you can redistribute it and/or modify
|
|
// it under the terms of the GNU General Public License as published by
|
|
// the Free Software Foundation, either version 3 of the License,
|
|
// or (at your option) any later version.
|
|
//
|
|
// This program is distributed in the hope that it will be useful,
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty
|
|
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
// See the GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU General Public License
|
|
// along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
//
|
|
// For more information, please consult: <https://herezh.irdl.fr/>.
|
|
|
|
//#include "Debug.h"
|
|
|
|
|
|
# include <iostream>
|
|
using namespace std; //introduces namespace std
|
|
#include <stdlib.h>
|
|
#include "Sortie.h"
|
|
|
|
|
|
#include "ElemPoint.h"
|
|
#include "FrontPointF.h"
|
|
#include "TypeConsTens.h"
|
|
#include "Def_Umat.h"
|
|
#include "ExceptionsLoiComp.h"
|
|
|
|
//----------------------------------------------------------------
|
|
// def des donnees commune a tous les elements
|
|
// la taille n'est pas defini ici car elle depend de la lecture
|
|
//----------------------------------------------------------------
|
|
ElemPoint::DonneeCommune * ElemPoint::doCo_Point = NULL;
|
|
ElemPoint::UneFois ElemPoint::unefois_Point;
|
|
ElemPoint::NombresConstruire ElemPoint::nombre_V;
|
|
ElemPoint::ConstrucElemPoint ElemPoint::construcElemPoint;
|
|
|
|
// constructeur définissant les nombres (de noeud, de point d'integ ..)
|
|
// utilisé dans la construction des éléments
|
|
ElemPoint::NombresConstruire::NombresConstruire()
|
|
{ nbne = 1; // le nombre de noeud de l'élément
|
|
nbiEr = 1;// le nombre de point d'intégration pour le calcul d'erreur
|
|
nbiMas = 1; // le nombre de point d'intégration pour le calcul de la matrice masse
|
|
};
|
|
|
|
// ---- definition du constructeur de la classe conteneur de donnees communes ------------
|
|
ElemPoint::DonneeCommune::DonneeCommune (GeomPoint& pteg,DdlElement& tab,DdlElement& tabErr,DdlElement& tab_Err1Sig,
|
|
Met_ElemPoint& met_point,
|
|
Tableau <Vecteur *> & resEr,Mat_pleine& raidEr,
|
|
GeomPoint& pteEr,Vecteur& residu_int,Mat_pleine& raideur_int,
|
|
Tableau <Vecteur* > & residus_extN,Tableau <Mat_pleine* >& raideurs_extN,
|
|
Mat_pleine& mat_masse,GeomPoint& pteMa,UmatAbaqus& umatAbaqusqus,int dim_tenseur,int nbi) :
|
|
ptpoint(pteg),tab_ddl(tab),tab_ddlErr(tabErr),tab_Err1Sig11(tab_Err1Sig)
|
|
,met_ElemPoint(met_point)
|
|
,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)
|
|
,pteEr(pteEr)
|
|
,residu_interne(residu_int),raideur_interne(raideur_int)
|
|
,residus_externeN(residus_extN),raideurs_externeN(raideurs_extN)
|
|
,matrice_masse(mat_masse)
|
|
,pteMas(pteMa)
|
|
,umatAbaqus(umatAbaqusqus),inne()
|
|
{
|
|
int nbddl = tab.NbDdl();
|
|
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_tenseur);
|
|
};
|
|
int tailledeps = d_epsBB.Taille();
|
|
for (int i=1;i<= tailledeps; i++)
|
|
d_epsBB(i) = NevezTenseurBB (dim_tenseur);
|
|
int tailledsig = d_sigHH.Taille();
|
|
for (int j=1;j<= tailledsig; j++)
|
|
d_sigHH(j) = NevezTenseurHH (dim_tenseur);
|
|
// liaisons entre inne et umatAbaqus
|
|
inne.incre = &umatAbaqus.nb_increment;
|
|
inne.step = &umatAbaqus.nb_step;
|
|
inne.nbe = &umatAbaqus.nb_elem;
|
|
inne.nbpti = &umatAbaqus.nb_pt_integ;
|
|
inne.temps_tdt = &umatAbaqus.temps_tdt;
|
|
inne.delta_t = &umatAbaqus.delta_t;
|
|
inne.nom_loi = &umatAbaqus.nom_materiau;
|
|
};
|
|
|
|
ElemPoint::DonneeCommune::DonneeCommune(DonneeCommune& a) :
|
|
ptpoint(a.ptpoint),tab_ddl(a.tab_ddl),tab_ddlErr(a.tab_ddlErr),tab_Err1Sig11(a.tab_Err1Sig11)
|
|
,met_ElemPoint(a.met_ElemPoint),matGeom(a.matGeom),matInit(a.matInit)
|
|
,d2_epsBB(a.d2_epsBB),resErr(a.resErr),raidErr(a.raidErr),pteEr(a.pteEr)
|
|
,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)
|
|
,matrice_masse(a.matrice_masse),pteMas(a.pteMas)
|
|
,umatAbaqus(a.umatAbaqus),inne()
|
|
{ int nbddl = d_sigHH.Taille();
|
|
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++)
|
|
d_epsBB(i) = NevezTenseurBB (*(a.d_epsBB(i)));
|
|
int tailledsig = d_sigHH.Taille();
|
|
for (int j=1;j<= tailledsig; j++)
|
|
d_sigHH(j) = NevezTenseurHH (*(d_sigHH(j)));
|
|
// liaisons entre inne et umatAbaqus
|
|
inne.incre = &umatAbaqus.nb_increment;
|
|
inne.step = &umatAbaqus.nb_step;
|
|
inne.nbe = &umatAbaqus.nb_elem;
|
|
inne.nbpti = &umatAbaqus.nb_pt_integ;
|
|
inne.temps_tdt = &umatAbaqus.temps_tdt;
|
|
inne.delta_t = &umatAbaqus.delta_t;
|
|
inne.nom_loi = &umatAbaqus.nom_materiau;
|
|
};
|
|
ElemPoint::DonneeCommune::~DonneeCommune()
|
|
{ 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);
|
|
};
|
|
|
|
// ---------- fin definition de la classe conteneur de donnees communes ------------
|
|
// -+-+ definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
|
|
ElemPoint::UneFois::UneFois () : // constructeur par défaut
|
|
doCoMemb(NULL),CalResPrem_t(0),CalResPrem_tdt(0),CalimpPrem(0),dualSortbiel(0)
|
|
,CalSMlin_t(0),CalSMlin_tdt(0),CalSMRlin(0)
|
|
,CalDynamique(0),CalPt_0_t_tdt(0)
|
|
,nbelem_in_Prog(0)
|
|
{};
|
|
ElemPoint::UneFois::~UneFois ()
|
|
{ delete doCoMemb;
|
|
};
|
|
|
|
// -+-+ fin definition de la classe contenant tous les indicateurs qui sont modifiés une seule fois -+-+-+
|
|
|
|
|
|
// ---------- fin definition de la classe conteneur de donnees communes ------------
|
|
|
|
// Constructeur par defaut
|
|
// les tenseurs on par défaut (dimension==-1) la dimension de l'espace
|
|
// sinon il ont la dimension donnée
|
|
ElemPoint::ElemPoint (int dimension) :
|
|
ElemMeca(),nbi(1),nb_appelsCalculUmat(0),lesPtMecaInt()
|
|
|
|
{ unefois=&unefois_Point; // le cas d'un vrai ElemPoint
|
|
if (unefois->nbelem_in_Prog == 0)
|
|
{ unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
|
|
}
|
|
else // sinon on construit
|
|
{id_interpol=CONSTANT; // donnees de la classe mere
|
|
id_geom=POINT;
|
|
Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE); //mise à jour du type de pb
|
|
tab_noeud.Change_taille(nombre_V.nbne);
|
|
// initialisation par défaut
|
|
int dim_tenseur = ParaGlob::Dimension();
|
|
if (dimension != -1) dim_tenseur = dimension;
|
|
Init (dim_tenseur);
|
|
unefois->nbelem_in_Prog++;
|
|
};
|
|
};
|
|
|
|
// Constructeur fonction d'un numero de maillage et d'identification
|
|
ElemPoint::ElemPoint (int num_mail,int num_id,int dimension) :
|
|
ElemMeca(num_mail,num_id,CONSTANT,POINT),nbi(1),nb_appelsCalculUmat(0)
|
|
,lesPtMecaInt()
|
|
{unefois=&unefois_Point; // le cas d'un vrai ElemPoint
|
|
if (unefois->nbelem_in_Prog == 0)
|
|
{ unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
|
|
}
|
|
else // sinon on construit
|
|
{tab_noeud.Change_taille(nombre_V.nbne);
|
|
// initialisation par défaut
|
|
int dim_tenseur = ParaGlob::Dimension();
|
|
if (dimension != -1) dim_tenseur = dimension;
|
|
unefois->nbelem_in_Prog++;
|
|
Init (dim_tenseur);
|
|
Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE); //mise à jour du type de pb
|
|
};
|
|
};
|
|
|
|
//------- pour des classes dérivées ------
|
|
// CONSTRUCTEURS :
|
|
// les tenseurs on par défaut (dimension==-1) la dimension de l'espace
|
|
// sinon il ont la dimension donnée
|
|
ElemPoint::ElemPoint (ElemPoint::UneFois& unefoisfois, Enum_geom nouveau_id, int dimension) :
|
|
ElemMeca(),nbi(1),nb_appelsCalculUmat(0),lesPtMecaInt()
|
|
|
|
{ unefois=&unefoisfois; // le cas d'un vrai ElemPoint
|
|
if (unefois->nbelem_in_Prog == 0)
|
|
{ unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
|
|
}
|
|
else // sinon on construit
|
|
{id_interpol=CONSTANT; // donnees de la classe mere
|
|
id_geom=nouveau_id;
|
|
Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE); //mise à jour du type de pb
|
|
tab_noeud.Change_taille(nombre_V.nbne);
|
|
// initialisation par défaut
|
|
int dim_tenseur = ParaGlob::Dimension();
|
|
if (dimension != -1) dim_tenseur = dimension;
|
|
unefois->nbelem_in_Prog++;
|
|
Init (dim_tenseur);
|
|
};
|
|
};
|
|
// Constructeur fonction d'un numero de maillage et d'identification
|
|
ElemPoint::ElemPoint (ElemPoint::UneFois& unefoisfois, Enum_geom nouveau_id, int num_mail,int num_id,int dimension) :
|
|
ElemMeca(num_mail,num_id,CONSTANT,nouveau_id),nbi(1),nb_appelsCalculUmat(0)
|
|
,lesPtMecaInt()
|
|
{unefois=&unefoisfois; // le cas d'un vrai ElemPoint
|
|
if (unefois->nbelem_in_Prog == 0)
|
|
{ unefois->nbelem_in_Prog++; // au premier passage on se contente d'incrémenter
|
|
}
|
|
else // sinon on construit
|
|
{tab_noeud.Change_taille(nombre_V.nbne);
|
|
// initialisation par défaut
|
|
int dim_tenseur = ParaGlob::Dimension();
|
|
if (dimension != -1) dim_tenseur = dimension;
|
|
unefois=&unefoisfois; // le cas d'un vrai ElemPoint
|
|
unefois->nbelem_in_Prog++;
|
|
Init (dim_tenseur);
|
|
Element::Change_TypeProblem(MECA_SOLIDE_DEFORMABLE); //mise à jour du type de pb
|
|
};
|
|
};
|
|
|
|
|
|
ElemPoint::ElemPoint (const ElemPoint& a) :
|
|
ElemMeca (a),nbi(a.nbi),nb_appelsCalculUmat(0)
|
|
,lesPtMecaInt(a.lesPtMecaInt)
|
|
// Constructeur de copie
|
|
// a priori si on utilise le constructeur de copie, donc il y a déjà un élément
|
|
// par contre a priori on ne doit pas faire une copie du premier élément
|
|
{unefois= a.unefois; // le cas d'un vrai ElemPoint
|
|
if (unefois->nbelem_in_Prog == 1)
|
|
{ cout << "\n **** erreur pour l'element ElemPoint, le constructeur de copie ne doit pas etre utilise"
|
|
<< " pour le premier element !! " << endl;
|
|
Sortie (1);
|
|
}
|
|
else
|
|
{// initialisation
|
|
unefois = a.unefois;
|
|
int dim_tenseur=lesPtMecaInt(1).EpsBB()->Dimension();
|
|
unefois->nbelem_in_Prog++;
|
|
Init (dim_tenseur);
|
|
};
|
|
};
|
|
|
|
ElemPoint::~ElemPoint ()
|
|
// Destructeur
|
|
{ LibereTenseur();
|
|
unefois->nbelem_in_Prog--;
|
|
Destruction();
|
|
};
|
|
|
|
// Lecture des donnees de la classe sur fichier
|
|
void
|
|
ElemPoint::LectureDonneesParticulieres
|
|
(UtilLecture * entreePrinc,Tableau<Noeud *> * tabMaillageNoeud)
|
|
{ int nb;
|
|
tab_noeud.Change_taille(1);
|
|
int i=1;
|
|
// on ne continue que si le tableau de noeud passé en paramètre a une taille non nulle
|
|
if ( (*tabMaillageNoeud).Taille() != 0)
|
|
{*(entreePrinc->entree) >> nb;
|
|
if ( ((entreePrinc->entree)->rdstate() == 0) || ((entreePrinc->entree)->eof()) )
|
|
// ici il n'y a qu'une seule information à lire, donc dès le début on a
|
|
// io_state == 1
|
|
|
|
// 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 != 2) 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 << "ElemPoint::LectureDonneesParticulieres";
|
|
Affiche();
|
|
Sortie (1);
|
|
};
|
|
// construction du tableau de ddl des noeuds de ElemPoint
|
|
ConstTabDdl();
|
|
};
|
|
};
|
|
|
|
// on associe un noeud à l'élément, remplace le noeud existant, s'il existe déjà
|
|
void ElemPoint::Associer_noeud (Noeud * noeu)
|
|
{ // on dimensionne le tableau de noeud éventuellement
|
|
tab_noeud.Change_taille(1);
|
|
// on réaffecte
|
|
tab_noeud(1)=noeu;
|
|
// construction du tableau de ddl des noeuds de ElemPoint
|
|
ConstTabDdl();
|
|
};
|
|
|
|
// 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 & ElemPoint::Point_physique(const Coordonnee& c_int,Coordonnee & co,Enum_dure temps)
|
|
{ ElemPoint::DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier l'écriture
|
|
// a) on commence par définir les bonnes grandeurs dans la métrique
|
|
if( !(unefois->CalPt_0_t_tdt ))
|
|
{ unefois->CalPt_0_t_tdt += 1;
|
|
Tableau<Enum_variable_metrique> tab(3);
|
|
tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
// b) calcul de l'interpolation
|
|
const Vecteur& phi = doCo->ptpoint.Phi_point(c_int);
|
|
// c) calcul du point
|
|
switch (temps)
|
|
{ case TEMPS_0 : co = doCo->met_ElemPoint.PointM_0(tab_noeud,phi); break;
|
|
case TEMPS_t : co = doCo->met_ElemPoint.PointM_t(tab_noeud,phi); break;
|
|
case TEMPS_tdt : co = doCo->met_ElemPoint.PointM_tdt(tab_noeud,phi); break;
|
|
};
|
|
// 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 ElemPoint::Point_physique(const Coordonnee& c_int,Tableau <Coordonnee> & t_co)
|
|
{ ElemPoint::DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier l'écriture
|
|
// a) on commence par définir les bonnes grandeurs dans la métrique
|
|
if( !(unefois->CalPt_0_t_tdt ))
|
|
{ unefois->CalPt_0_t_tdt += 1;
|
|
Tableau<Enum_variable_metrique> tab(3);
|
|
tab(1)=iM0;tab(2)=iMt;tab(3)=iMtdt;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
// b) calcul de l'interpolation
|
|
const Vecteur& phi = doCo->ptpoint.Phi_point(c_int);
|
|
// c) calcul des point
|
|
switch (t_co.Taille())
|
|
{ case 3 : t_co(3) = doCo->met_ElemPoint.PointM_tdt(tab_noeud,phi);
|
|
case 2 : t_co(2) = doCo->met_ElemPoint.PointM_t(tab_noeud,phi);
|
|
case 1 : t_co(1) = doCo->met_ElemPoint.PointM_0(tab_noeud,phi);
|
|
};
|
|
};
|
|
|
|
// Calcul du residu local à t ou tdt en fonction du booleen atdt
|
|
Vecteur* ElemPoint::CalculResidu (bool atdt,const ParaAlgoControle & pa)
|
|
{ ElemPoint::DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier l'écriture
|
|
Tableau <TenseurBB *>& d_epsBB = (doCo->d_epsBB);// "
|
|
// 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;
|
|
doCo->met_ElemPoint.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;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
};
|
|
// initialisation du résidu
|
|
residu->Zero();
|
|
Vecteur poids =(doCo->ptpoint).TaWi(); // poids d'interpolation = 2
|
|
ElemMeca::Cal_explicit ( doCo->tab_ddl,d_epsBB,1,poids,pa,atdt);
|
|
return residu;
|
|
};
|
|
|
|
// Calcul du residu local et de la raideur locale,
|
|
// pour le schema implicite
|
|
Element::ResRaid ElemPoint::Calcul_implicit (const ParaAlgoControle & pa)
|
|
{ ElemPoint::DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier l'écriture
|
|
Tableau <TenseurBB *>& d_epsBB = (doCo->d_epsBB);// "
|
|
Tableau <TenseurHH *>& d_sigHH = (doCo->d_sigHH);// "
|
|
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;
|
|
doCo->met_ElemPoint.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;
|
|
};
|
|
|
|
// initialisation du résidu
|
|
residu->Zero();
|
|
// initialisation de la raideur
|
|
raideur->Zero();
|
|
|
|
Vecteur poids =(doCo->ptpoint).TaWi();
|
|
// --- pour l'instant il n'y a pas de calcul, on retourne les matrices vides
|
|
// ElemMeca::Cal_implicit
|
|
// (doCo->tab_ddl, d_epsBB,(doCo->d2_epsBB),d_sigHH,nbi,poids,pa,cald_Dvirtuelle);
|
|
Element::ResRaid el;
|
|
el.res = residu;
|
|
el.raid = raideur;
|
|
return el;
|
|
};
|
|
|
|
// calcul de l'UMat pour abaqus
|
|
void ElemPoint::CalculUmatAbaqus(ParaAlgoControle & pa)
|
|
{ // -- récupération des variables Umat abaqus
|
|
UmatAbaqus& umatA = unefois->doCoMemb->umatAbaqus; // pour simplifier
|
|
// cas du premier passage
|
|
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;
|
|
unefois->doCoMemb->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
// -- appel de la loi de comportement
|
|
int ni=umatA.nb_pt_integ; // num du pt d'integ (pour simplifier)
|
|
// on change les coordonnées du point associé au premier pti
|
|
// (a priori ce n'est pas utilisé ici mais on le fait quand même pour le futur)
|
|
tab_noeud(1)->Change_coord2(umatA.coor_pt);
|
|
// on change la température au noeud en fonction de ce qui existe à l'umat
|
|
// car ce sera utilisé pour les grandeurs interpolées aux noeuds
|
|
if (tab_noeud(1)->Existe_ici(TEMP))
|
|
{tab_noeud(1)->Change_val_t(TEMP,umatA.temper_t);
|
|
tab_noeud(1)->Change_val_tdt(TEMP,umatA.temper_t + umatA.delta_temper);
|
|
};
|
|
|
|
// on change le numéro d'intégration pour les calculs d'interpolation c-a-d des acces
|
|
// aux tableaux phi et dphi : pour le point il n'y a qu'un pti !!
|
|
def->ChangeNumInteg(1); //ni);
|
|
|
|
// on associe la déformation à la bonne sauvegarde
|
|
def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
|
|
// mise à jour de la métrique avec les infos de umat
|
|
((Def_Umat*)def)->Mise_a_jourUmat(umatA);
|
|
// valeurs aux pt d'integ
|
|
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
|
|
|
|
// on regarde si c'est le premier incrément ou pas (d'où le calcul du premier pas)
|
|
premier_calcul_meca_impli_expli=false;
|
|
if (umatA.nb_increment <= 1) premier_calcul_meca_impli_expli = true;
|
|
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
|
|
if (loiTP != NULL) {sTP = tabSaveTP(ni);}; // au pt d'integ si elles existes
|
|
#ifdef MISE_AU_POINT
|
|
if ( ParaGlob::NiveauImpression() > 5)
|
|
{cout << "\n ElemPoint::CalculUmatAbaqus(... "
|
|
<< " dilatation= "<< dilatation
|
|
<< flush;
|
|
};
|
|
#endif
|
|
// on gère les exceptions éventuelles en mettant le bloc sous surveillance
|
|
try
|
|
{ // appel du comportement Umat
|
|
loiComp->ComportementUmat(tabSaveDon(ni),*def,ptIntegMeca,pa,sTP,loiTP,dilatation,umatA
|
|
,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
|
|
{ // l'idée est que dans le cas où il y a un pb avec la loi on signale qu'il faut diminuer le pas de temps ??
|
|
// c'est surtout pour une utilisation avec Abaqus
|
|
umatA.pnewdt = sqrt(2.)/2.;
|
|
if (ParaGlob::NiveauImpression() >= 1)
|
|
cout << "\n warning: CalculUmatAbaqus: modification du pas de temps -> du parametre pnewt";
|
|
}
|
|
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() >= 1)
|
|
{cout << "\n warning: exception generee par la loi de comportement sous Umat mais dont la prise en compte "
|
|
<< " n'est pas prevu !, on demande une diminution du pas de temps de 0.5 ";
|
|
// << " n'est pas prevu !, on ne fait rien et on continue le calcul";
|
|
umatA.pnewdt = 0.5;
|
|
if (ParaGlob::NiveauImpression() >= 4) cout << "\n ElemPoint::CalculUmatAbaqus(..";
|
|
};
|
|
};
|
|
|
|
// comptabilisation des appels si c'est le premier point d'intégration
|
|
if (ni == 1)
|
|
nb_appelsCalculUmat++;
|
|
};
|
|
|
|
// initialisation éventuelle: ajout de point d'intégration si nécessaire
|
|
// les tenseurs on par défaut (dim_tens==-1) la dimension de l'espace
|
|
void ElemPoint::InitialisationUmatAbaqus_interne(int dim_tens)
|
|
{ // dans le cas où le nombre de point d'intégration lue est supérieur à celui de l'élément
|
|
// on agrandi les tableaux (qui sont des tableaux de pointeurs
|
|
if (unefois->doCoMemb->umatAbaqus.nb_pt_integ > nbi)
|
|
{ int dim_tenseur = ParaGlob::Dimension();
|
|
if (dim_tens != -1) dim_tenseur = dim_tens;
|
|
ChangeNombrePtinteg(unefois->doCoMemb->umatAbaqus.nb_pt_integ,dim_tenseur);
|
|
};
|
|
};
|
|
|
|
// Calcul de la matrice masse pour l'élément
|
|
Mat_pleine * ElemPoint::CalculMatriceMasse (Enum_calcul_masse type_calcul_masse)
|
|
{ ElemPoint::DonneeCommune* doCo = 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;
|
|
doCo->met_ElemPoint.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 = doCo->tab_ddl.NbDdl();
|
|
(doCo->matrice_masse).Initialise (nbddl,nbddl,0.);
|
|
}
|
|
};
|
|
Vecteur poids =(doCo->pteMas).TaWi(); // poids d'intégration = 2
|
|
// appel de la routine générale
|
|
ElemMeca::Cal_Mat_masse (doCo->tab_ddl,type_calcul_masse,
|
|
nombre_V.nbiMas,(doCo->pteMas).TaPhi(),nombre_V.nbne
|
|
,poids);
|
|
return mat_masse;
|
|
};
|
|
|
|
//------- calcul d'erreur, remontée des contraintes -------------------
|
|
// 1) calcul du résidu et de la matrice de raideur pour le calcul d'erreur
|
|
Element::Er_ResRaid ElemPoint::ContrainteAuNoeud_ResRaid()
|
|
{ DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
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;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
// appel du programme général
|
|
ElemMeca:: SigmaAuNoeud_ResRaid(tab_noeud.Taille()
|
|
,(doCo->ptpoint).TaPhi()
|
|
,(doCo->ptpoint).TaWi()
|
|
,doCo-> resErr,doCo->raidErr
|
|
,(doCo->pteEr).TaPhi()
|
|
,(doCo->pteEr).TaWi());
|
|
return (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr)));
|
|
};
|
|
|
|
// 2) remontée aux erreurs aux noeuds
|
|
Element::Er_ResRaid ElemPoint::ErreurAuNoeud_ResRaid()
|
|
{ DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
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;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
// appel du programme général
|
|
ElemMeca::Cal_ErrAuxNoeuds(tab_noeud.Taille(), (doCo->ptpoint).TaPhi(),
|
|
(doCo->ptpoint).TaWi(),doCo-> resErr );
|
|
return (Element::Er_ResRaid( &(doCo-> resErr),&(doCo->raidErr)));
|
|
};
|
|
|
|
// actualisation des ddl et des grandeurs actives de t+dt vers t
|
|
void ElemPoint::TdtversT()
|
|
{ for (int ni=1;ni<= nbi; ni++)
|
|
{ lesPtMecaInt.TdtversT();
|
|
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
|
|
// normalement ça correspond à la convergence, on met à -1 le nombre d'appel
|
|
// ainsi au premier appel on aura l'itération 0
|
|
nb_appelsCalculUmat = -1;
|
|
};
|
|
// actualisation des ddl et des grandeurs actives de t vers tdt
|
|
void ElemPoint::TversTdt()
|
|
{ for (int ni=1;ni<= nbi; ni++)
|
|
{ lesPtMecaInt.TversTdt();
|
|
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. Ce calcul n'est disponible
|
|
// qu'une fois la remontée aux contraintes effectuées sinon aucune
|
|
// action. En retour la valeur de l'erreur sur l'élément
|
|
// type indique le type de calcul d'erreur :
|
|
void ElemPoint::ErreurElement(int type,double& errElemRelative
|
|
,double& numerateur, double& denominateur)
|
|
{ DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
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;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
// appel du programme général
|
|
ElemMeca::Cal_ErrElem(type,errElemRelative,numerateur,denominateur,
|
|
tab_noeud.Taille(),(doCo->ptpoint).TaPhi(),
|
|
(doCo->ptpoint).TaWi(),
|
|
(doCo->pteEr).TaPhi(),(doCo->pteEr).TaWi());
|
|
};
|
|
|
|
//============= 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 ElemPoint::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 la ElemPoint
|
|
// car on n'a pas dimensionné les tableaux qui suivent !!
|
|
switch (cas)
|
|
{ case 1 : // ------- on récupère tout -------------------------
|
|
{ string toto;
|
|
int nbi_old = nbi;
|
|
ent >> toto >> nbi ; // récup du nombre de point d'intégration
|
|
// construction du tableau de ddl des noeuds de ElemPoint
|
|
ConstTabDdl();
|
|
// des déformations et contraintes éventuelles
|
|
int dim_tens;
|
|
ent >> dim_tens; // lecture de l'indicateur pour savoir s'il y a des tenseurs ou pas
|
|
if (nbi != nbi_old)
|
|
ChangeNombrePtinteg(nbi,dim_tens);
|
|
// récupération
|
|
lesPtMecaInt.Lecture_base_info(ent,cas);
|
|
break;
|
|
}
|
|
case 2 : // ----------- lecture uniquement de se qui varie --------------------
|
|
{ // récup info
|
|
// des déformations et contraintes éventuelles
|
|
int dim_tens;string toto;
|
|
ent >> toto >> dim_tens; // lecture de l'indicateur pour savoir s'il y a des tenseurs ou pas
|
|
// récupération
|
|
lesPtMecaInt.Lecture_base_info(ent, cas);
|
|
break;
|
|
}
|
|
default :
|
|
{ cout << "\nErreur : valeur incorrecte du type de lecture !\n";
|
|
cout << "ElemPoint::Lecture_base_info(ofstream& sort,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 ElemPoint::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 la ElemPoint
|
|
// sinon on essaie de sauvegarder
|
|
switch (cas)
|
|
{ case 1 : // ------- on sauvegarde tout -------------------------
|
|
{sort << "\n nbpti= " << nbi << " ";
|
|
sort << lesPtMecaInt.DimTens() << " ";
|
|
// des tenseurs déformation et contrainte,
|
|
lesPtMecaInt.Ecriture_base_info(sort,cas);
|
|
break;
|
|
}
|
|
case 2 : // ----------- sauvegarde uniquement de se qui varie --------------------
|
|
{ // la déformation et la contrainte
|
|
sort << "\n dim_tenseur: " << lesPtMecaInt.DimTens() << " ";
|
|
// des tenseurs déformation et contrainte,
|
|
lesPtMecaInt.Ecriture_base_info(sort,cas);
|
|
break;
|
|
}
|
|
default :
|
|
{ cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
|
|
cout << "ElemPoint::Ecriture_base_info(ofstream& sort,int cas)"
|
|
<< " cas= " << cas << endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
};
|
|
|
|
|
|
// Calcul de la matrice géométrique et initiale
|
|
ElemMeca::MatGeomInit ElemPoint::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & pa)
|
|
{ ElemPoint::DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier l'écriture
|
|
Tableau <TenseurBB *>& d_epsBB = (doCo->d_epsBB);// "
|
|
Tableau <TenseurHH *>& d_sigHH = (doCo->d_sigHH);// "
|
|
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;
|
|
doCo->met_ElemPoint.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;
|
|
};
|
|
// Par simplicité
|
|
Mat_pleine & matGeom = doCo->matGeom;
|
|
Mat_pleine & matInit = doCo->matInit;
|
|
// mise à zéro de la matrice géométrique
|
|
matGeom.Initialise();
|
|
Vecteur poids =(doCo->ptpoint).TaWi(); // poids d'interpolation = 2
|
|
ElemMeca::Cal_matGeom_Init (matGeom,matInit,doCo->tab_ddl, d_epsBB,
|
|
doCo->d2_epsBB,d_sigHH,1,poids,pa,cald_Dvirtuelle);
|
|
return MatGeomInit(&matGeom,&matInit);
|
|
} ;
|
|
|
|
// retourne les tableaux de ddl associés aux noeuds, gere par l'element
|
|
// ce tableau et specifique a l'element
|
|
const DdlElement & ElemPoint::TableauDdl() const
|
|
{ DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
return doCo->tab_ddl; };
|
|
// liberation de la place pointee
|
|
void ElemPoint::Libere ()
|
|
{Element::Libere (); // liberation de residu et raideur
|
|
LibereTenseur() ; // liberation des tenseur intermediaires
|
|
};
|
|
|
|
// acquisition ou modification d'une loi de comportement
|
|
void ElemPoint::DefLoi (LoiAbstraiteGeneral * NouvelleLoi)
|
|
{ // verification du type de loi
|
|
int dim=ParaGlob::Dimension(); // pour l'instant par défaut
|
|
|
|
// non avec les CP on passe en 2D au niveau de la loi, donc a priori on peut avoir plusieurs dimensions
|
|
// if (NouvelleLoi->Dimension_loi() != dim)
|
|
// { cout << "\n Erreur, la loi de comportement a utiliser avec des ElemPoints";
|
|
// cout << " est de type = " << (NouvelleLoi->Dimension_loi()) << " au lieu de " << dim << "!!! " << endl;
|
|
// Sortie(1);
|
|
// };
|
|
|
|
// cas d'une loi mécanique
|
|
if (GroupeMecanique(NouvelleLoi->Id_categorie()))
|
|
{loiComp = (Loi_comp_abstraite *) NouvelleLoi;
|
|
// initialisation du stockage particulier, ici 1 pt d'integ
|
|
tabSaveDon(1) = loiComp->New_et_Initialise();
|
|
// idem pour le type de déformation mécanique associé
|
|
tabSaveDefDon(1) = def->New_et_Initialise();
|
|
// définition du type de déformation associé à la loi
|
|
loiComp->Def_type_deformation(*def);
|
|
// on active les données particulières nécessaires au fonctionnement de la loi de comp
|
|
loiComp->Activation_donnees(tab_noeud,dilatation,lesPtMecaInt);
|
|
};
|
|
// cas d'une loi thermo physique
|
|
if (GroupeThermique(NouvelleLoi->Id_categorie()))
|
|
{loiTP = (CompThermoPhysiqueAbstraite *) NouvelleLoi;
|
|
// initialisation du stockage particulier, ici 1 pt d'integ
|
|
tabSaveTP(1) = loiTP->New_et_Initialise();
|
|
// on active les données particulières nécessaires au fonctionnement de la loi de comp
|
|
loiTP->Activation_donnees(tab_noeud);
|
|
};
|
|
// cas d'une loi de frottement
|
|
if (GroupeFrottement(NouvelleLoi->Id_categorie()))
|
|
loiFrot = (CompFrotAbstraite *) NouvelleLoi;
|
|
};
|
|
|
|
// test si l'element est complet
|
|
int ElemPoint::TestComplet()
|
|
{ int res = ElemMeca::TestComplet(); // test dans la fonction mere
|
|
return res;
|
|
};
|
|
|
|
// ajout du tableau de ddl des noeuds de ElemPoint
|
|
void ElemPoint::ConstTabDdl()
|
|
{
|
|
Tableau <Ddl> ta(ParaGlob::Dimension());
|
|
int posi = Id_nom_ddl("X1") -1;
|
|
int dim = ParaGlob::Dimension();
|
|
for (int i =1; i<= dim; i++)
|
|
{Ddl inter((Enum_ddl(i+posi)),0.,LIBRE);
|
|
ta(i) = inter;
|
|
};
|
|
// attribution des ddls aux noeuds
|
|
tab_noeud(1)->PlusTabDdl(ta);
|
|
};
|
|
|
|
// procesure permettant de completer l'element apres
|
|
// sa creation avec les donnees du bloc transmis
|
|
// peut etre appeler plusieurs fois
|
|
Element* ElemPoint::Complete(BlocGen & bloc,LesFonctions_nD* lesFonctionsnD)
|
|
{ // complétion avec bloc
|
|
return ElemMeca::Complete_ElemMeca(bloc,lesFonctionsnD);
|
|
};
|
|
|
|
// lecture de données diverses sur le flot d'entrée
|
|
void ElemPoint::LectureContraintes(UtilLecture * entreePrinc)
|
|
{ DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
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;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
ElemMeca::LectureDesContraintes (true,entreePrinc,lesPtMecaInt.TabSigHH_t());
|
|
};
|
|
|
|
// retour des contraintes en absolu retour true si elle existe sinon false
|
|
bool ElemPoint::ContraintesAbsolues(Tableau <Vecteur>& tabSig)
|
|
{ DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
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;
|
|
doCo->met_ElemPoint.PlusInitVariables(tab) ;
|
|
};
|
|
ElemMeca::ContraintesEnAbsolues(true,lesPtMecaInt.TabSigHH_t(),tabSig);
|
|
return true;
|
|
};
|
|
|
|
// affichage dans la sortie transmise, des variables duales "nom"
|
|
// dans le cas ou nom est vide, affichage de "toute" les variables
|
|
void ElemPoint::AfficheVarDual(ofstream& sort, Tableau<string>& nom)
|
|
{// affichage de l'entête de l'element
|
|
sort << "\n******************************************************************";
|
|
sort << "\n Element bielette (2 noeuds 1 point d'integration) ";
|
|
sort << "\n******************************************************************";
|
|
// appel de la procedure de elem meca
|
|
if (!(unefois->dualSortbiel ) && (unefois->CalimpPrem ))
|
|
{ VarDualSort(sort,nom,1,1);
|
|
unefois->dualSortbiel += 1;
|
|
}
|
|
else if ((unefois->dualSortbiel) && (unefois->CalimpPrem ))
|
|
VarDualSort(sort,nom,1,11);
|
|
else if (!(unefois->dualSortbiel ) && (unefois->CalResPrem_tdt ))
|
|
{ VarDualSort(sort,nom,1,2);
|
|
unefois->dualSortbiel += 1;
|
|
}
|
|
else if ((unefois->dualSortbiel ) && (unefois->CalResPrem_tdt ))
|
|
VarDualSort(sort,nom,1,12);
|
|
// sinon on ne fait rien
|
|
};
|
|
|
|
// 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 ElemPoint::CoordPtInteg(Enum_dure temps,Enum_ddl enu,int iteg,bool& erreur)
|
|
{ // ici les coordonnées du point d'intégration et celles du noeud sont identique
|
|
// donc on ramène les coordonnées du noeud associé
|
|
switch (temps) {
|
|
case TEMPS_0: {return (tab_noeud(1)->Coord0()); break;};
|
|
case TEMPS_t: {return (tab_noeud(1)->Coord1()); break;};
|
|
case TEMPS_tdt:{ return (tab_noeud(1)->Coord2()); break;};
|
|
};
|
|
};
|
|
|
|
// récupération des valeurs au numéro d'ordre = iteg pour
|
|
// les grandeur enu
|
|
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
|
|
Tableau <double> ElemPoint::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->dualSortbiel) && (unefois->CalimpPrem))
|
|
{ cas=1;unefois->dualSortbiel += 1;
|
|
}
|
|
else if ((unefois->dualSortbiel) && (unefois->CalimpPrem))
|
|
{ cas = 11;}
|
|
else if (!(unefois->dualSortbiel) && (unefois->CalResPrem_tdt))
|
|
{ cas=2;unefois->dualSortbiel += 1;
|
|
}
|
|
else if ((unefois->dualSortbiel) && (unefois->CalResPrem_tdt))
|
|
{ cas = 12;}
|
|
// sinon pour l'instant pb, car il faut définir des variable dans la métrique
|
|
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 prevu, unefois->dualSortbiel= " << unefois->dualSortbiel
|
|
<< " unefois->CalimpPrem= " << unefois->CalimpPrem
|
|
<< "\n ElemPoint::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
|
|
void ElemPoint::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->dualSortbiel) && (unefois->CalimpPrem))
|
|
{ cas=1;unefois->dualSortbiel += 1;
|
|
}
|
|
else if ((unefois->dualSortbiel) && (unefois->CalimpPrem))
|
|
{ cas = 11;}
|
|
else if (!(unefois->dualSortbiel) && (unefois->CalResPrem_tdt))
|
|
{ cas=2;unefois->dualSortbiel = 1;
|
|
}
|
|
else if ((unefois->dualSortbiel) && (unefois->CalResPrem_tdt))
|
|
{ cas = 12;}
|
|
// sinon pour l'instant pb, car il faut définir des variable dans la métrique
|
|
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->dualSortbiel= " << unefois->dualSortbiel
|
|
<< " unefois->CalimpPrem= " << unefois->CalimpPrem
|
|
<< "\n ElemPoint::ValTensorielle_a_diff_temps(Enum_dure enu_t...";
|
|
Sortie(1);
|
|
};
|
|
ElemMeca::Valeurs_Tensorielles(absolue,enu_t,enu,iteg,cas);
|
|
};
|
|
|
|
// Calcul des frontieres de l'element
|
|
// creation des elements frontieres et retour du tableau de ces elements
|
|
// la création n'a lieu qu'au premier appel
|
|
// ou lorsque l'on force le paramètre force a true
|
|
// dans ce dernier cas seul les frontière effacées sont recréée
|
|
Tableau <ElFrontiere*> const & ElemPoint::Frontiere(bool force)
|
|
{ int cas = 3; // on veut des points
|
|
return Frontiere_elemeca(cas,force);
|
|
// { // le calcul et la création ne sont effectués qu'au premier appel
|
|
// // ou lorsque l'on veut forcer une recréation
|
|
// if (((ind_front_lin == 0) && (ind_front_surf == 0) && (ind_front_point == 0))
|
|
// || force )
|
|
//// if ((ind_front_point == 0) || force || (ind_front_point == 2))
|
|
// { // dimensionnement des tableaux intermediaires
|
|
// Tableau <Noeud *> tab(1); // les noeuds des points frontieres
|
|
// DdlElement ddelem(1); // les ddlelements des points frontieres
|
|
// int tail=1;
|
|
// tabb.Change_taille(tail); // le tableau total de frontières
|
|
//
|
|
// // le point
|
|
//// tab(1) = tab_noeud(1);
|
|
//// ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(1));
|
|
//// if (tabb(1+posi_tab_front_point) == NULL)
|
|
//// tabb(1+posi_tab_front_point) = new FrontPointF (tab,ddelem);
|
|
// // mise à jour des indicateurs
|
|
//// ind_front_point = 1;
|
|
// }
|
|
//
|
|
// return tabb;
|
|
};
|
|
|
|
//// ramène la frontière point
|
|
//// éventuellement création des frontieres points de l'element et stockage dans l'element
|
|
//// si c'est la première fois sinon il y a seulement retour de l'elements
|
|
//// a moins que le paramètre force est mis a true
|
|
//// dans ce dernier cas la frontière effacéee est recréée
|
|
//// num indique le numéro du point à créer (numérotation EF)
|
|
//ElFrontiere* const ElemPoint::Frontiere_points(int num,bool force)
|
|
// { // le calcul et la création ne sont effectués qu'au premier appel
|
|
// // ou lorsque l'on veut forcer une recréation
|
|
// #ifdef MISE_AU_POINT
|
|
// if (num != 1)
|
|
// { cout << "\n *** erreur, pour les ElemPoints il n'y a qu'une frontieres point ! "
|
|
// << "\n Frontiere_points(int num,bool force)";
|
|
// Sortie(1);
|
|
// }
|
|
// #endif
|
|
//
|
|
// if ((ind_front_point == 0) || force || (ind_front_point == 2))
|
|
// {Tableau <Noeud *> tab(1); // les noeuds des points frontieres
|
|
// DdlElement ddelem(1); // les ddlelements des points frontieres
|
|
// // on regarde si la frontière point existe sinon on la crée
|
|
// if (ind_front_point == 1)
|
|
// return (ElFrontiere*)tabb(posi_tab_front_point+num);
|
|
// else if ( ind_front_point == 2)
|
|
// // cas où certaines frontières existent
|
|
// if (tabb(posi_tab_front_point+num) != NULL)
|
|
// return (ElFrontiere*)tabb(posi_tab_front_point+num);
|
|
// // dans tous les autres cas on construit la frontière point
|
|
// // on commence par dimensionner le tableau de frontière
|
|
// if (ind_front_point==0)
|
|
// // cas où aucune frontière existe, on crée pour les points
|
|
// { tabb.Change_taille(1); posi_tab_front_point = 0;};
|
|
// // dans les deux autres cas ((ind_front_lin == 0) && (ind_front_point>0)) et
|
|
// // ((ind_front_lin > 0) && (ind_front_point>0)), les points existant déjà on n'a rien n'a faire
|
|
// // on définit le point
|
|
// // premier point
|
|
//// tab(1) = tab_noeud(1);
|
|
//// ddelem.Change_un_ddlNoeudElement(1,doCo->tab_ddl(1));
|
|
//// if (tabb(1+posi_tab_front_point) == NULL)
|
|
//// tabb(1+posi_tab_front_point) = new FrontPointF (tab,ddelem);
|
|
// ind_front_point=1; // mise à jour de l'indicateur
|
|
// };
|
|
// return (ElFrontiere*)tabb(num+posi_tab_front_point);
|
|
// };
|
|
//
|
|
//// ramène la frontière linéique
|
|
//// éventuellement création des frontieres linéique de l'element et stockage dans l'element
|
|
//// si c'est la première fois et en 3D sinon il y a seulement retour de l'elements
|
|
//// a moins que le paramètre force est mis a true
|
|
//// dans ce dernier cas la frontière effacéee est recréée
|
|
//// num indique le numéro de l'arête à créer (numérotation EF)
|
|
//// ici normalement la fonction ne doit pas être appelée
|
|
//ElFrontiere* const ElemPoint::Frontiere_lineique(int ,bool )
|
|
// { cout << "\n *** erreur, pour les ElemPoints il n'y a pas de frontiere lineique ! "
|
|
// << "\n Frontiere_lineique(int ,bool force = false)";
|
|
// Sortie(1);
|
|
// return NULL;
|
|
// };
|
|
//
|
|
//// ramène la frontière surfacique
|
|
//// éventuellement création des frontieres surfacique de l'element et stockage dans l'element
|
|
//// si c'est la première fois sinon il y a seulement retour de l'elements
|
|
//// a moins que le paramètre force est mis a true
|
|
//// dans ce dernier cas la frontière effacéee est recréée
|
|
//// num indique le numéro de la surface à créer (numérotation EF)
|
|
//// ici normalement la fonction ne doit pas être appelée
|
|
//ElFrontiere* const ElemPoint::Frontiere_surfacique(int ,bool )
|
|
// { cout << "\n *** erreur, pour les ElemPoints il n'y a pas de frontiere surface ! "
|
|
// << "\n Frontiere_surfacique(int ,bool force = false)";
|
|
// Sortie(1);
|
|
// return NULL;
|
|
// };
|
|
|
|
// =====>>>> methodes privées appelees par les classes dérivees <<<<=====
|
|
|
|
// fonction d'initialisation servant au niveau du constructeur
|
|
// dim_tenseur: = dimension des tenseurs
|
|
void ElemPoint::Init(int dim_tenseur)
|
|
{ // le fait de mettre les pointeurs a null permet
|
|
// de savoir que l'element n'est pas complet
|
|
tab_noeud.Change_taille(nombre_V.nbne);
|
|
for (int i =1;i<= nombre_V.nbne;i++) tab_noeud(i) = NULL;
|
|
// definition des donnees communes aux ElemPointxxx
|
|
// a la premiere definition d'une instance
|
|
if (unefois->doCoMemb == NULL)
|
|
unefois->doCoMemb = ElemPoint::Def_DonneeCommune(dim_tenseur);
|
|
// unefois->doCoMemb = doCo ;
|
|
met = &(unefois->doCoMemb->met_ElemPoint); // met est defini dans elemeca
|
|
DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier
|
|
// def pointe sur la deformation specifique a l'element pour le calcul mecanique
|
|
def = new Def_Umat(*met,tab_noeud,(doCo->ptpoint).TaDphi(),(doCo->ptpoint).TaPhi());
|
|
// idem pour la remontee aux contraintes et le calcul d'erreur
|
|
defEr = new Def_Umat(*met,tab_noeud,(doCo->pteEr).TaDphi(),(doCo->pteEr).TaPhi());
|
|
// idem pour la remontee aux contraintes et le calcul d'erreur
|
|
defMas = new Def_Umat(*met,tab_noeud,(doCo->pteMas).TaDphi(),(doCo->pteMas).TaPhi());
|
|
// idem pour le calcul de second membre
|
|
defArete.Change_taille(1); // 1 arrête utilisée pour le second membre
|
|
// la déformation sera construite si nécessaire au moment du calcul de second membre
|
|
defArete(1) = NULL;
|
|
|
|
//dimensionnement des deformations et contraintes
|
|
lesPtMecaInt.Change_taille_PtIntegMeca(nbi,dim_tenseur);
|
|
// attribution des numéros de référencement dans le conteneur
|
|
for (int ni = 1; ni<= nbi; ni++)
|
|
{lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage);
|
|
lesPtMecaInt(ni).Change_Nb_ele(this->num_elt);
|
|
lesPtMecaInt(ni).Change_Nb_pti(ni);
|
|
};
|
|
|
|
lesPtIntegMecaInterne = &lesPtMecaInt; // association avec le pointeur d'ElemMeca
|
|
|
|
// stockage des donnees particulieres de la loi de comportement au point d'integ
|
|
tabSaveDon.Change_taille(nbi);
|
|
tabSaveTP.Change_taille(nbi);
|
|
tabSaveDefDon.Change_taille(nbi);
|
|
tab_energ.Change_taille(nbi);
|
|
tab_energ_t.Change_taille(nbi);
|
|
// initialisation des pointeurs définis dans la classe Element concernant les résidus et
|
|
// raideur
|
|
// --- cas de la puissance interne ---
|
|
residu = &(doCo->residu_interne); // residu local
|
|
raideur = &(doCo->raideur_interne); // raideur locale
|
|
// --- cas de la dynamique -----
|
|
mat_masse = &(doCo->matrice_masse);
|
|
// --- cas des efforts externes concernant les noeuds ------
|
|
res_extN = &(doCo->residus_externeN); // pour les résidus et second membres
|
|
raid_extN= &(doCo->raideurs_externeN);// pour les raideurs
|
|
};
|
|
|
|
|
|
// fonction privee
|
|
// dans cette fonction il ne doit y avoir que les données communes !!
|
|
ElemPoint::DonneeCommune* ElemPoint::Def_DonneeCommune(int dim_tenseur)
|
|
{ int nbn = nombre_V.nbne;
|
|
// interpollation : element geometrique de base correspondant: 1 pt integ, 1 noeuds
|
|
GeomPoint pteg ;
|
|
// degre de liberte
|
|
int dim = ParaGlob::Dimension();
|
|
DdlElement tab_ddl(nbn,dim);
|
|
int posi = Id_nom_ddl("X1") -1;
|
|
for (int i =1; i<= dim; i++)
|
|
for (int j=1; j<= nbn; j++)
|
|
// tab_ddl (j,i) = Enum_ddl(i+posi);
|
|
tab_ddl.Change_Enum(j,i,Enum_ddl(i+posi));
|
|
// cas des ddl éléments secondaires pour le calcul d'erreur
|
|
// def du nombre de composantes du tenseur de contrainte en absolu
|
|
int nbcomposante = 0; // init ParaGlob::NbCompTens ();
|
|
int nb_vect=0; // init: nb de vecteur de la base
|
|
switch (dim_tenseur)
|
|
{ case 3: nbcomposante=6;nb_vect=3; break;
|
|
case 2: nbcomposante=3;nb_vect=2; break;
|
|
case 1: nbcomposante=1;nb_vect=1; break;
|
|
default: break;
|
|
};
|
|
DdlElement tab_ddlErr(nbn,nbcomposante);
|
|
posi = Id_nom_ddl("SIG11") -1;
|
|
for (int j=1; j<= nbn; j++)
|
|
{// on definit le nombre de composante de sigma en absolu
|
|
if (nbcomposante == 3)
|
|
// dans l'énumération les composantes ne sont pas a suivre
|
|
{//tab_ddlErr (j,1) = Enum_ddl(1+posi); // cas de SIG11
|
|
//tab_ddlErr (j,2) = Enum_ddl(2+posi); // cas de SIG22
|
|
//tab_ddlErr (j,3) = Enum_ddl(4+posi); // cas de SIG12
|
|
tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi)); // cas de SIG11
|
|
tab_ddlErr.Change_Enum(j,2,Enum_ddl(2+posi)); // cas de SIG22
|
|
tab_ddlErr.Change_Enum(j,3,Enum_ddl(4+posi)); // cas de SIG12
|
|
}
|
|
else if (nbcomposante == 6)
|
|
// les composantes sont a suivre dans l'enumération
|
|
{for (int i= 1;i<= nbcomposante; i++)
|
|
// tab_ddlErr (j,i) = Enum_ddl(i+posi);
|
|
tab_ddlErr.Change_Enum(j,i,Enum_ddl(i+posi));
|
|
}
|
|
else
|
|
// tab_ddlErr (j,1) = Enum_ddl(1+posi); // uniquement SIG11
|
|
{tab_ddlErr.Change_Enum(j,1,Enum_ddl(1+posi));}; // uniquement SIG11
|
|
};
|
|
// egalement pour tab_Err1Sig11, def d'un tableau de un ddl : enum SIG11
|
|
// par noeud
|
|
DdlElement tab_Err1Sig11(nbn,DdlNoeudElement(SIG11));
|
|
// toujours pour le calcul d'erreur definition des fonctions d'interpolation
|
|
// pour le calcul du hession de la fonctionnelle :
|
|
// 1 points d'integration et 1 noeuds
|
|
GeomPoint ptegEr ;
|
|
// pour le calcul de la matrice masse definition des fonctions d'interpolation
|
|
// 1 points d'integration et 1 noeuds, en particulier pour le calcul de la masse consistante
|
|
GeomPoint ptegMa ;
|
|
// def metrique
|
|
// on definit les variables a priori toujours utiles
|
|
Tableau<Enum_variable_metrique> tab(24);
|
|
tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt ;
|
|
tab(4) = igiB_0; tab(5) = igiB_t; tab(6) = igiB_tdt;
|
|
tab(7) = igiH_0; tab(8) = igiH_t; tab(9) = igiH_tdt ;
|
|
tab(10)= igijBB_0; tab(11)= igijBB_t; tab(12)= igijBB_tdt;
|
|
tab(13)= igijHH_0; tab(14)= igijHH_t; tab(15)= igijHH_tdt ;
|
|
tab(16)= id_gijBB_tdt; tab(17)= id_giH_tdt; tab(18)= id_gijHH_tdt;
|
|
tab(19)= idMtdt ; tab(20)= id_jacobien_tdt;tab(21)= id2_gijBB_tdt;
|
|
tab(22)= igradVBB_tdt; tab(23) = iVtdt; tab(24)= idVtdt;
|
|
// dim du pb , nb de vecteur de la base , tableau de ddl et la def de variables
|
|
Met_ElemPoint metri(ParaGlob::Dimension(),nb_vect,tab_ddl,tab,nbn) ;
|
|
// ---- cas du calcul d'erreur sur sigma ou epsilon
|
|
// les tenseurs sont exprimees en absolu donc nombre de composante fonction
|
|
// de la dimension absolue
|
|
Tableau <Vecteur *> resEr(nbcomposante);
|
|
for (int i = 1;i<= nbcomposante; i++)
|
|
resEr(i)=new Vecteur (nbn); // une composante par noeud
|
|
Mat_pleine raidEr(nbn,nbn); // la raideur pour l'erreur
|
|
// dimensionnement des différents résidus et raideurs pour le calcul mécanique
|
|
int nbddl = tab_ddl.NbDdl();
|
|
Vecteur residu_int(nbddl); Mat_pleine raideur_int(nbddl,nbddl);
|
|
// cas de la dynamique
|
|
Mat_pleine matmasse(1,nbddl); // a priori on dimensionne en diagonale
|
|
// il y a une extrémité a priori
|
|
Tableau <Vecteur* > residus_extN(1); residus_extN(1) = new Vecteur(dim);
|
|
// pas d'arêtes, un noeud
|
|
Tableau <Mat_pleine* > raideurs_extN(1);raideurs_extN(1) = new Mat_pleine(dim,dim);
|
|
// -- particularités pour la routine Umat pour Abaqus
|
|
UmatAbaqus umatAbaqus(ParaGlob::Dimension(),dim_tenseur);
|
|
// definition de la classe static contenant toute les variables communes aux ElemPoints
|
|
int nbi=1; // par défaut on met un pt d'integ, ensuite il faudra remettre à jour
|
|
// les données communes si le nb de pt augmente
|
|
DonneeCommune* doCoco = new DonneeCommune(pteg,tab_ddl,tab_ddlErr,tab_Err1Sig11,metri,resEr,raidEr,ptegEr
|
|
,residu_int,raideur_int,residus_extN,raideurs_extN
|
|
,matmasse,ptegMa,umatAbaqus,dim_tenseur,nbi);
|
|
return doCoco;
|
|
};
|
|
|
|
// destructions de certaines grandeurs pointées, créées au niveau de l'initialisation
|
|
void ElemPoint::Destruction()
|
|
{ // tout d'abord l'idée est de détruire certaines grandeurs pointées que pour le dernier élément
|
|
if ((unefois->nbelem_in_Prog == 0)&& (unefois->doCoMemb != NULL))
|
|
// if (unefois->unefois->nbelem_in_Prog == 0)
|
|
// cas de la destruction du dernier élément
|
|
{ ElemPoint::DonneeCommune* doCo = unefois->doCoMemb; // pour simplifier l'écriture
|
|
int resErrTaille = doCo->resErr.Taille();
|
|
for (int i=1;i<= resErrTaille;i++)
|
|
delete doCo->resErr(i);
|
|
delete doCo->residus_externeN(1);
|
|
delete doCo->raideurs_externeN(1);
|
|
};
|
|
};
|
|
|
|
// changement du nombre de point d'intégration
|
|
void ElemPoint::ChangeNombrePtinteg(int nevez_nbi, int dim_tens)
|
|
{ // les numéroq
|
|
int old_nbi = nbi;
|
|
nbi = nevez_nbi;
|
|
// stockage des donnees particulieres de la loi de comportement au point d'integ
|
|
tabSaveDon.Change_taille(nbi);
|
|
tabSaveTP.Change_taille(nbi);
|
|
tabSaveDefDon.Change_taille(nbi);
|
|
// puis les tenseurs internes
|
|
lesPtMecaInt.Change_taille_PtIntegMeca(nbi,dim_tens);
|
|
// attribution des numéros de référencement dans le conteneur
|
|
for (int ni = 1; ni<= nbi; ni++)
|
|
{lesPtMecaInt(ni).Change_Nb_mail(this->num_maillage);
|
|
lesPtMecaInt(ni).Change_Nb_ele(this->num_elt);
|
|
lesPtMecaInt(ni).Change_Nb_pti(ni);
|
|
};
|
|
|
|
for (int ni = old_nbi+1; ni <= nbi; ni++)
|
|
{ // cas d'une loi mécanique
|
|
// initialisation du stockage particulier pour le nouveau pt d'integ
|
|
tabSaveDon(ni) = loiComp->New_et_Initialise();
|
|
// idem pour le type de déformation mécanique associé
|
|
tabSaveDefDon(ni) = def->New_et_Initialise();
|
|
// cas éventuelle d'une loi thermo physique
|
|
if (loiTP != NULL)
|
|
{// initialisation du stockage particulier, ici 1 pt d'integ
|
|
tabSaveTP(ni) = loiTP->New_et_Initialise();
|
|
};
|
|
};
|
|
};
|