Gérard Rio
ece68c42db
Introduction du calcul de l'accélération interpolée, au niveau de la métrique Intro de nouveaux affichages au niveau des fonctions nD
3583 lines
182 KiB
C++
3583 lines
182 KiB
C++
|
|
// 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 <cmath>
|
|
#include "ElemMeca.h"
|
|
#include <iomanip>
|
|
#include "ConstMath.h"
|
|
#include "Util.h"
|
|
#include "Coordonnee2.h"
|
|
#include "Coordonnee3.h"
|
|
#include "ParaAlgoControle.h"
|
|
#include "ExceptionsElemMeca.h"
|
|
#include "Loi_Umat.h"
|
|
#include "ExceptionsLoiComp.h"
|
|
#include "CharUtil.h"
|
|
|
|
// =========================== variables communes ===============================
|
|
// a priori les pointeurs de fonction pointent sur rien au debut
|
|
const Coordonnee & (Met_abstraite::*(ElemMeca::PointM))
|
|
(const Tableau<Noeud *>& tab_noeud,const Vecteur& Phi) ;
|
|
void (Met_abstraite::*(ElemMeca::BaseND))
|
|
(const Tableau<Noeud *>& tab_noeud,const Mat_pleine& dphi,
|
|
const Vecteur& phi,BaseB& bB,BaseH& bH);
|
|
int ElemMeca::bulk_viscosity = 0; // indique si oui ou non on utilise le bulk viscosity
|
|
double ElemMeca::c_traceBulk = 0.; // coeff de la trace de D pour le bulk
|
|
double ElemMeca::c_carreBulk = 0.; // coeff du carré de la trace de D pour le bulk
|
|
TenseurHH* ElemMeca::sig_bulkHH = NULL; // variable de travail pour le bulk
|
|
// ---------- pour le calcul de la masse pour la méthode de relaxation dynamique
|
|
Ddl_enum_etendu ElemMeca::masse_relax_dyn; // au début rien, def à la première utilisation
|
|
// pour la gestion d'hourglass: dimensionnement dans ElemMeca::Init_hourglass_comp
|
|
Tableau< Vecteur> ElemMeca::vec_hourglass; // tableau de vecteurs de travail, éventuellement défini
|
|
// pour la stabilisation de membrane: dimensionnement lors de l'utilisation
|
|
//Tableau< Vecteur> ElemMeca::vec_StabMembBiel; // tableau de vecteurs de travail, éventuellement défini
|
|
|
|
// =========================== fin variables communes ===========================
|
|
// -------------------- stabilisation transversale éventuelle de membrane ou de biel -------------------
|
|
// utilisée et alimentées par les classes dérivées: ElemMeca sert de conteneur ce qui permet
|
|
// d'éviter de redéfinir des conteneurs locaux aux éléments membranes et biel
|
|
// par défaut
|
|
ElemMeca::StabMembBiel::StabMembBiel() :
|
|
type_stabMembrane(STABMEMBRANE_BIEL_NON_DEFINIE)
|
|
,gamma(0.),pt_fct_gamma(NULL),stab_ref(1.),activite_en_explicite(true)
|
|
,F_StabMembBiel(),F_StabMembBiel_t()
|
|
,E_StabMembBiel(),E_StabMembBiel_t()
|
|
,nom_fctnD(NULL),nom_fctnD_2(NULL)
|
|
,affichage_stabilisation(0),gestion_maxi_mini(false)
|
|
,beta(1.) // valeur arbitraire pour l'instant
|
|
,f_mini(ConstMath::petit) // valeur arbitraire pour l'instant
|
|
,d_maxi(1000.) // déplacement arbitraire de 1000
|
|
,pt_fct_ponder_Ft(NULL)
|
|
{};
|
|
// fonction du nombre de noeud ou de point d'interpolation
|
|
ElemMeca::StabMembBiel::StabMembBiel(int nbnoe) :
|
|
type_stabMembrane(STABMEMBRANE_BIEL_NON_DEFINIE)
|
|
,gamma(0.),pt_fct_gamma(NULL),stab_ref(1.),activite_en_explicite(true)
|
|
,F_StabMembBiel(nbnoe),F_StabMembBiel_t(nbnoe)
|
|
,E_StabMembBiel(nbnoe),E_StabMembBiel_t(nbnoe)
|
|
,nom_fctnD(NULL),nom_fctnD_2(NULL)
|
|
,affichage_stabilisation(0),gestion_maxi_mini(false)
|
|
,beta(1.) // valeur arbitraire pour l'instant
|
|
,f_mini(ConstMath::petit) // valeur arbitraire pour l'instant
|
|
,d_maxi(1000.) // déplacement arbitraire de 100
|
|
,pt_fct_ponder_Ft(NULL)
|
|
{};
|
|
// fonction d'une valeur numérique et des noms éventuels de fonction
|
|
ElemMeca::StabMembBiel::StabMembBiel(double v,string * sgamma,string* sponder)
|
|
: type_stabMembrane(STABMEMBRANE_BIEL_NON_DEFINIE)
|
|
,gamma(v),pt_fct_gamma(NULL),F_StabMembBiel()
|
|
,E_StabMembBiel(),E_StabMembBiel_t()
|
|
,nom_fctnD(NULL),nom_fctnD_2(NULL)
|
|
,stab_ref(1.),activite_en_explicite(true)
|
|
,affichage_stabilisation(0),gestion_maxi_mini(false)
|
|
,beta(1.) // valeur arbitraire pour l'instant
|
|
,f_mini(ConstMath::petit) // valeur arbitraire pour l'instant
|
|
,d_maxi(1000.) // déplacement arbitraire de 100
|
|
,pt_fct_ponder_Ft(NULL)
|
|
{if (sgamma != NULL) nom_fctnD= new string (*sgamma);
|
|
if (sponder != NULL) nom_fctnD_2= new string (*sponder);
|
|
};
|
|
// de copie
|
|
ElemMeca::StabMembBiel::StabMembBiel(const StabMembBiel& a):
|
|
type_stabMembrane(a.type_stabMembrane)
|
|
,gamma(a.gamma),pt_fct_gamma(a.pt_fct_gamma),stab_ref(a.stab_ref)
|
|
,activite_en_explicite(a.activite_en_explicite)
|
|
,F_StabMembBiel(a.F_StabMembBiel),F_StabMembBiel_t(a.F_StabMembBiel_t)
|
|
,E_StabMembBiel(a.E_StabMembBiel),E_StabMembBiel_t(a.E_StabMembBiel_t)
|
|
,nom_fctnD(NULL),nom_fctnD_2(NULL)
|
|
,affichage_stabilisation(a.affichage_stabilisation)
|
|
,gestion_maxi_mini(a.gestion_maxi_mini)
|
|
,beta(a.beta)
|
|
,f_mini(a.f_mini)
|
|
,d_maxi(a.d_maxi),pt_fct_ponder_Ft(a.pt_fct_ponder_Ft)
|
|
{if (a.nom_fctnD != NULL) nom_fctnD = new string (*(a.nom_fctnD));
|
|
if (a.nom_fctnD_2 != NULL) nom_fctnD_2 = new string (*(a.nom_fctnD_2));
|
|
};
|
|
ElemMeca::StabMembBiel::~StabMembBiel()
|
|
{pt_fct_gamma=NULL; // la fonction est forcément globale, donc on ne la détruie pas
|
|
if (nom_fctnD != NULL) delete nom_fctnD;
|
|
// idem pour la pondération
|
|
if (nom_fctnD_2 != NULL) delete nom_fctnD_2;
|
|
};
|
|
ElemMeca::StabMembBiel& ElemMeca::StabMembBiel::operator= (const ElemMeca::StabMembBiel& a)
|
|
{type_stabMembrane = a.type_stabMembrane;
|
|
gamma=a.gamma;pt_fct_gamma=a.pt_fct_gamma;
|
|
stab_ref=a.stab_ref;activite_en_explicite = a.activite_en_explicite;
|
|
F_StabMembBiel=a.F_StabMembBiel;F_StabMembBiel_t=a.F_StabMembBiel_t;
|
|
E_StabMembBiel=a.E_StabMembBiel;E_StabMembBiel_t = a.E_StabMembBiel_t;
|
|
nom_fctnD=NULL;pt_fct_ponder_Ft=a.pt_fct_ponder_Ft;
|
|
affichage_stabilisation = a.affichage_stabilisation;
|
|
gestion_maxi_mini=a.gestion_maxi_mini;
|
|
beta=a.beta;f_mini=a.f_mini;d_maxi=a.d_maxi;
|
|
if (a.nom_fctnD != NULL) nom_fctnD = new string (*(a.nom_fctnD));
|
|
if (a.nom_fctnD_2 != NULL) nom_fctnD_2 = new string (*(a.nom_fctnD_2));
|
|
return (*this);
|
|
};
|
|
|
|
// énergie totale à t+dt développée sur l'élément pour la stabilisation
|
|
double ElemMeca::StabMembBiel::EE_total_StabMembBiel() const
|
|
{ double retour = 0; // init
|
|
int taille = E_StabMembBiel.Taille();
|
|
for (int i=1; i<=taille;i++)
|
|
retour += E_StabMembBiel(i);
|
|
return retour;
|
|
};
|
|
|
|
// énergie totale à t développée sur l'élément pour la stabilisation
|
|
double ElemMeca::StabMembBiel::EE_total_StabMembBiel_t() const
|
|
{ double retour = 0; // init
|
|
int taille = E_StabMembBiel_t.Taille();
|
|
for (int i=1; i<=taille;i++)
|
|
retour += E_StabMembBiel_t(i);
|
|
return retour;
|
|
};
|
|
|
|
// actualisation de la force de stabilisation de t+dt vers t
|
|
void ElemMeca::StabMembBiel::TdtversT()
|
|
{F_StabMembBiel_t = F_StabMembBiel;
|
|
E_StabMembBiel_t = E_StabMembBiel;
|
|
};
|
|
// actualisation de la force de stabilisation de t vers tdt
|
|
void ElemMeca::StabMembBiel::TversTdt()
|
|
{F_StabMembBiel = F_StabMembBiel_t;
|
|
E_StabMembBiel = E_StabMembBiel_t;
|
|
};
|
|
|
|
// 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 ElemMeca::StabMembBiel::Lecture_bas_inf(ifstream& ent,const int cas)
|
|
{string toto,nom;
|
|
switch (cas)
|
|
{ case 1 : // ------- on récupère tout -------------------------
|
|
{ // récup du type de stabilisation
|
|
ent >> toto;
|
|
type_stabMembrane = Id_Enum_StabMembraneBiel(toto);
|
|
// récup pour alpha
|
|
ent >> toto;
|
|
if (toto == "par_fonction_nD")
|
|
{ ent >> nom;
|
|
if (nom_fctnD == NULL)
|
|
nom_fctnD = new string(nom);
|
|
else *nom_fctnD = nom;
|
|
gamma = 0.; // init au cas ou
|
|
pt_fct_gamma = NULL; // même si diff de NULL
|
|
// il n'y a pas delete à faire car il s'agit de fonction globale
|
|
// ensuite il faudra définir la fonction de stabilisation
|
|
// au moment de compléter l'élément
|
|
}
|
|
else if (toto == "par_valeur")
|
|
{ent >> gamma;}
|
|
else
|
|
{ cout << "\n *** erreur en lecture, on attendait"
|
|
<< " par_fonction_nD ou par_valeur et on a lue "
|
|
<< toto ;
|
|
Sortie(1);
|
|
};
|
|
// activite_en_explicite
|
|
ent >> toto >> activite_en_explicite;
|
|
// pondération éventuelle
|
|
ent >> toto;
|
|
if (toto == "ponderation_F_t")
|
|
{ ent >> nom;
|
|
if (nom_fctnD_2 == NULL)
|
|
nom_fctnD_2 = new string(nom);
|
|
else *nom_fctnD_2 = nom;
|
|
pt_fct_ponder_Ft = NULL;// même si diff de NULL
|
|
// il n'y a pas delete à faire car il s'agit de fonction globale
|
|
// ensuite il faudra définir la fonction de pondération
|
|
// au moment de compléter l'élément
|
|
};
|
|
// affichage et les limitations
|
|
ent >> toto >> affichage_stabilisation
|
|
>> toto >> gestion_maxi_mini
|
|
>> toto >> beta >> toto >> f_mini >> toto >> d_maxi;
|
|
break;
|
|
}
|
|
case 2 : // ----------- lecture uniquement de se qui varie --------------------
|
|
{ ent >> toto >> F_StabMembBiel;
|
|
ent >> toto >> E_StabMembBiel;
|
|
break;
|
|
}
|
|
default :
|
|
{ cout << "\nErreur : valeur incorrecte du type de lecture !\n";
|
|
cout << "ElemMeca::StabMembBiel::Lecture_bas_inf(ifstream& ent,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 ElemMeca::StabMembBiel::Ecriture_bas_inf(ofstream& sort,const int cas)
|
|
{switch (cas)
|
|
{ case 1 : // ------- on sauvegarde tout -------------------------
|
|
{ sort << " "<<Nom_StabMembraneBiel(type_stabMembrane);
|
|
if (pt_fct_gamma != NULL)
|
|
{sort << " par_fonction_nD " << pt_fct_gamma->NomFonction()<<" ";}
|
|
else {sort << " par_valeur " << gamma<< " ";};
|
|
sort << " activite_en_explicite "<< activite_en_explicite;
|
|
if (pt_fct_ponder_Ft != NULL)
|
|
{sort << " ponderation_F_t " << pt_fct_ponder_Ft->NomFonction()<<" ";}
|
|
else {sort << " pas_de_ponderation_F_t " ;};
|
|
sort << "\n affichage_stabilisation= "<<affichage_stabilisation
|
|
<< " gestion_maxi_mini= "<< gestion_maxi_mini
|
|
<< " beta= "<< beta << " f_mini= "
|
|
<< f_mini << " d_maxi= " << d_maxi << " ";
|
|
break;
|
|
}
|
|
case 2 : // ----------- sauvegarde uniquement de se qui varie --------------------
|
|
{ sort << " F_StabMembBiel "<< F_StabMembBiel;
|
|
sort << " E_StabMembBiel "<< E_StabMembBiel
|
|
<< " ";
|
|
break;
|
|
}
|
|
default :
|
|
{ cout << "\nErreur : valeur incorrecte du type d'écriture !\n";
|
|
cout << "ElemMeca::StabMembBiel::Ecriture_bas_inf(ofstream& sort,const int cas)"
|
|
<< " cas= " << cas << endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
};
|
|
|
|
// complète en attribuant les fct nD si besoin, ceci après une lecture
|
|
// initiale base_info
|
|
void ElemMeca::StabMembBiel::Complete_StabMembBiel(LesFonctions_nD* lesFonctionsnD)
|
|
{ // on complète si les pointeurs de fonction sont nulle et qu'il y a des noms de fonction
|
|
if (nom_fctnD != NULL)
|
|
if (pt_fct_gamma == NULL)
|
|
pt_fct_gamma = lesFonctionsnD->Trouve(*nom_fctnD);
|
|
if (nom_fctnD_2 != NULL)
|
|
if (pt_fct_ponder_Ft == NULL)
|
|
pt_fct_ponder_Ft = lesFonctionsnD->Trouve(*nom_fctnD_2);
|
|
};
|
|
|
|
// changement du nombre de noeuds ou de pti
|
|
void ElemMeca::StabMembBiel::Change_nb_pti(int nbnoe)
|
|
{F_StabMembBiel.Change_taille(nbnoe);
|
|
F_StabMembBiel_t.Change_taille(nbnoe);
|
|
E_StabMembBiel.Change_taille(nbnoe);
|
|
E_StabMembBiel_t.Change_taille(nbnoe);
|
|
};
|
|
// -------------------- fin conteneur stabilisation transversale éventuelle de membrane ou de biel -------------------
|
|
|
|
|
|
// Constructeur par defaut
|
|
ElemMeca::ElemMeca () :
|
|
Element(),tabSaveDon(),tabSaveTP(),tabSaveDefDon(),defSurf(),defArete()
|
|
,masse_volumique(masse_volumique_defaut),dilatation(0),tab_energ(),tab_energ_t()
|
|
,energie_totale(),energie_totale_t()
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,premier_calcul_meca_impli_expli(true),lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
// hourglass
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
// stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL;
|
|
sigErreur_relative = NULL;
|
|
};
|
|
|
|
// Constructeur utile quand le numero d'identification de l'element est connu
|
|
ElemMeca::ElemMeca (int num_mail,int num_id) :
|
|
Element(num_mail,num_id),tabSaveDon(),tabSaveTP(),tabSaveDefDon(),defSurf(),defArete()
|
|
,masse_volumique(masse_volumique_defaut),dilatation(0),tab_energ(),tab_energ_t()
|
|
,energie_totale(),energie_totale_t()
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,premier_calcul_meca_impli_expli(true)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
// hourglass
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
// stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL; sigErreur_relative= NULL;
|
|
};
|
|
|
|
// Constructeur utile quand le numero et le tableau des noeuds
|
|
// de l'element sont connu
|
|
ElemMeca::ElemMeca (int num_mail,int num_id,const Tableau<Noeud *>& tab):
|
|
Element(num_mail,num_id,tab)
|
|
,tabSaveDon(),tabSaveTP(),tabSaveDefDon(),defSurf(),defArete()
|
|
,masse_volumique(masse_volumique_defaut),dilatation(0),tab_energ(),tab_energ_t()
|
|
,energie_totale(),energie_totale_t()
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,premier_calcul_meca_impli_expli(true)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
// hourglass
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
// stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL; sigErreur_relative = NULL;
|
|
};
|
|
|
|
|
|
// Constructeur utile quand le numero d'identification est connu,
|
|
// ainsi que la geometrie et le type d'interpolation de l'element
|
|
ElemMeca::ElemMeca (int num_mail,int num_id,Enum_interpol id_interp_elt,Enum_geom id_geom_elt,string info):
|
|
Element (num_mail,num_id,id_interp_elt,id_geom_elt,MECA_SOLIDE_DEFORMABLE,info)
|
|
,tabSaveDon(),tabSaveTP(),tabSaveDefDon()
|
|
,defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(0)
|
|
,tab_energ(),tab_energ_t(),energie_totale(),energie_totale_t()
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,premier_calcul_meca_impli_expli(true)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
//--- stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL; sigErreur_relative = NULL;
|
|
};
|
|
|
|
ElemMeca::ElemMeca (int num_mail,int num_id,char* nom_interpol,char* nom_geom,string info):
|
|
// Constructeur utile quand le numero d'identification est connu,
|
|
// ainsi que la geometrie et le type d'interpolation de l'element
|
|
Element(num_mail,num_id,nom_interpol,nom_geom,NomElemTypeProblem(MECA_SOLIDE_DEFORMABLE),info)
|
|
,tabSaveDon()
|
|
,tabSaveTP(),tabSaveDefDon(),defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(0)
|
|
,tab_energ(),tab_energ_t(),energie_totale(),energie_totale_t(),premier_calcul_meca_impli_expli(true)
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
//--- stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL; sigErreur_relative = NULL;
|
|
};
|
|
|
|
ElemMeca::ElemMeca (int num_mail,int num_id,const Tableau<Noeud *>& tab,Enum_interpol id_interp_elt,
|
|
Enum_geom id_geom_elt,string info):
|
|
// Constructeur utile quand toutes les donnees sont connues
|
|
Element (num_mail,num_id,tab,id_interp_elt,id_geom_elt,MECA_SOLIDE_DEFORMABLE,info)
|
|
,tabSaveDon(),tabSaveTP()
|
|
,tabSaveDefDon(),defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(0)
|
|
,tab_energ(),tab_energ_t(),energie_totale(),energie_totale_t(),premier_calcul_meca_impli_expli(true)
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
//--- stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL; sigErreur_relative = NULL;
|
|
};
|
|
|
|
// Constructeur utile quand toutes les donnees sont connues
|
|
ElemMeca::ElemMeca (int num_mail,int num_id,const Tableau<Noeud *>& tab,char* nom_interpol,
|
|
char* nom_geom,string info):
|
|
Element (num_mail,num_id,tab,nom_interpol,nom_geom,NomElemTypeProblem(MECA_SOLIDE_DEFORMABLE),info)
|
|
,tabSaveDon()
|
|
,tabSaveTP(),tabSaveDefDon(),defSurf(),defArete(),masse_volumique(masse_volumique_defaut),dilatation(0)
|
|
,tab_energ(),tab_energ_t(),energie_totale(),energie_totale_t(),premier_calcul_meca_impli_expli(true)
|
|
,E_elem_bulk_t(0.),E_elem_bulk_tdt(0.),P_elem_bulk(0.),E_Hourglass(0.)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
,type_stabHourglass(STABHOURGLASS_NON_DEFINIE),tab_elHourglass()
|
|
,coefStabHourglass(0.),raid_hourglass_transitoire(NULL)
|
|
//--- stab transversale
|
|
,pt_StabMembBiel(NULL),matD(NULL),resD(NULL),maxi_F_t(0.)
|
|
{ loiComp = NULL;loiTP=NULL;loiFrot=NULL;
|
|
met = NULL;
|
|
def = NULL;
|
|
defEr = NULL;
|
|
defMas = NULL;
|
|
sigErreur = NULL; sigErreur_relative = NULL;
|
|
};
|
|
|
|
// Constructeur de copie
|
|
ElemMeca::ElemMeca (const ElemMeca& elt):
|
|
Element (elt),defSurf(elt.defSurf),defArete(elt.defArete)
|
|
,masse_volumique(elt.masse_volumique),dilatation(elt.dilatation)
|
|
,tab_energ(elt.tab_energ),tab_energ_t(elt.tab_energ_t)
|
|
,energie_totale(elt.energie_totale),energie_totale_t(elt.energie_totale_t)
|
|
,E_elem_bulk_t(elt.E_elem_bulk_t),E_elem_bulk_tdt(elt.E_elem_bulk_tdt),P_elem_bulk(elt.P_elem_bulk)
|
|
//--- hourglass
|
|
,E_Hourglass(elt.E_Hourglass)
|
|
//--- stab transversale
|
|
,pt_StabMembBiel(NULL),matD(elt.matD),resD(elt.resD),maxi_F_t(0.)
|
|
,premier_calcul_meca_impli_expli(true)
|
|
,lesPtIntegMecaInterne(NULL)
|
|
,lesChargeExtSurEle(NULL)
|
|
,type_stabHourglass(elt.type_stabHourglass),tab_elHourglass()
|
|
,coefStabHourglass(elt.coefStabHourglass),raid_hourglass_transitoire(NULL)
|
|
{ // pour la loi on pointe sur la même que l'élément de référence, car de toute
|
|
// la loi ne contient pas les données spécifique
|
|
loiComp = (elt.loiComp);loiTP = (elt.loiTP);loiFrot = (elt.loiFrot);
|
|
// idem pour la métrique qui est a priori une métrique générique
|
|
met = (elt.met);
|
|
// la déformation est spécifique à l'élément
|
|
if (elt.def!=NULL) def = elt.def->Nevez_deformation(tab_noeud); // *def = *(elt.def);
|
|
if (elt.defEr!=NULL) defEr = elt.defEr->Nevez_deformation(tab_noeud); //*defEr = *(elt.defEr);
|
|
if (elt.defMas!=NULL) defMas = elt.defMas->Nevez_deformation(tab_noeud); //*defMas = *(elt.defMas);
|
|
// les données spécifiques à la déformation mécanique
|
|
int idefdtr = (elt.tabSaveDefDon).Taille();
|
|
tabSaveDefDon.Change_taille(idefdtr);
|
|
for (int i=1;i<= idefdtr; i++)
|
|
{ tabSaveDefDon(i) = def->New_et_Initialise();
|
|
if (tabSaveDefDon(i)!= NULL) *tabSaveDefDon(i) = *(elt.tabSaveDefDon)(i);
|
|
};
|
|
// def: cas des arrêtes et des faces si besoin
|
|
int tailledefSurf = defSurf.Taille();
|
|
int itab = 1; // indice pour parcourir tabb
|
|
for (int idsu =1;idsu<=tailledefSurf;idsu++,itab++)
|
|
{if ((defSurf(idsu) != NULL) && (tabb.Taille() != 0)) // dans le cas des coques
|
|
// on peut avoir defSurf qui est défini, mais pas la frontière donc dans ce cas
|
|
// on ne peut pas utiliser cette copie
|
|
{ defSurf(idsu) = elt.defSurf(idsu)->Nevez_deformation(tabb(itab)->TabNoeud());}
|
|
else // sinon on supprime, il sera reconstruit si besoin est
|
|
{ // dans le cas où on a qu'une seule surface a priori et que tabb n'est pas définit
|
|
// on considère que le tableau de noeud associé est celui de l'élément
|
|
delete(defSurf(idsu));defSurf(idsu)=NULL;
|
|
};
|
|
};
|
|
int tailledefArete = defArete.Taille();
|
|
for (int idA =1;idA<=tailledefArete;idA++)
|
|
{if ((defArete(idA) != NULL) && (tabb.Taille() != 0))
|
|
defArete(idA) = elt.defArete(idA)->Nevez_deformation(tabb(itab)->TabNoeud());
|
|
else {delete(defArete(idA));defArete(idA)=NULL; }
|
|
};
|
|
// les donnees de la loi de comportement mécanique
|
|
int dtr = (elt.tabSaveDon).Taille();
|
|
tabSaveDon.Change_taille(dtr);
|
|
for (int i=1;i<= dtr; i++)
|
|
{ tabSaveDon(i) = loiComp->New_et_Initialise();
|
|
if (tabSaveDon(i)!= NULL) *tabSaveDon(i) = *(elt.tabSaveDon)(i);
|
|
};
|
|
// les donnees de la loi de comportement thermo physique
|
|
int dtrTP = (elt.tabSaveTP).Taille();
|
|
tabSaveTP.Change_taille(dtrTP);
|
|
for (int i=1;i<= dtrTP; i++)
|
|
{ tabSaveTP(i) = NULL; // init
|
|
if (loiTP != NULL ) tabSaveTP(i) = loiTP->New_et_Initialise();
|
|
if (tabSaveTP(i)!= NULL) *tabSaveTP(i) = *(elt.tabSaveTP)(i);
|
|
};
|
|
// et autre
|
|
if (elt.sigErreur != NULL)
|
|
{ sigErreur = new double;
|
|
(*sigErreur) = *(elt.sigErreur);
|
|
sigErreur_relative = new double;
|
|
*sigErreur_relative = *(elt.sigErreur_relative);
|
|
}
|
|
else
|
|
{sigErreur = sigErreur_relative = NULL;};
|
|
//---- stabilisation transversale éventuelle
|
|
if (elt.pt_StabMembBiel != NULL)
|
|
pt_StabMembBiel = new StabMembBiel(*(elt.pt_StabMembBiel));
|
|
};
|
|
|
|
// Destructeur
|
|
ElemMeca::~ElemMeca ()
|
|
{ // les donnees pour la loi de comportement mécanique
|
|
int imaxi = tabSaveDon.Taille();
|
|
for (int i=1;i<=imaxi; i++)
|
|
if (tabSaveDon(i) != NULL) {delete tabSaveDon(i) ; tabSaveDon(i)=NULL;}
|
|
// les donnees pour la loi de comportement thermo physique
|
|
int imaxTP = tabSaveTP.Taille();
|
|
for (int i=1;i<=imaxTP; i++)
|
|
if (tabSaveTP(i) != NULL) {delete tabSaveTP(i) ; tabSaveTP(i)=NULL;}
|
|
// les donnees à la déformation mécanique
|
|
int idefmaxi = tabSaveDefDon.Taille();
|
|
for (int i=1;i<=idefmaxi; i++)
|
|
if (tabSaveDefDon(i) != NULL) {delete tabSaveDefDon(i) ; tabSaveDefDon(i)=NULL;}
|
|
// les déformations sont spécifiquent aux elements, on les detruits si necessaire
|
|
if (def != NULL) { delete def; def = NULL;}
|
|
if (defEr != NULL) {delete defEr; defEr = NULL;}
|
|
if (defMas != NULL) {delete defMas; defMas = NULL;}
|
|
int nbdefSurf = defSurf.Taille();
|
|
for (int i = 1; i<= nbdefSurf; i++)
|
|
if (defSurf(i) != NULL) {delete defSurf(i); defSurf(i) = NULL;}
|
|
int nbdefArete = defArete.Taille();
|
|
for (int i = 1; i<= nbdefArete; i++)
|
|
if (defArete(i) != NULL) {delete defArete(i);defArete(i) = NULL;}
|
|
// par contre la metrique et la loi de comportement sont communes a une classe d'element
|
|
// donc ici met et loicomp ne sont que des pointeurs dont il ne faut pas supprimer les objets
|
|
// pointes
|
|
if (sigErreur!= NULL) {delete sigErreur; sigErreur = NULL;
|
|
delete sigErreur_relative; sigErreur_relative=NULL;
|
|
};
|
|
// --- cas des éléments de stabilisation d'hourglass ---
|
|
if (tab_elHourglass.Taille() != 0)
|
|
{ int taillhour = tab_elHourglass.Taille();
|
|
for (int i=1;i<= taillhour;i++)
|
|
{ if (tab_elHourglass(i) != NULL) delete tab_elHourglass(i);};
|
|
};
|
|
if (raid_hourglass_transitoire != NULL) delete raid_hourglass_transitoire;
|
|
// --- cas de la stabilisation transversale ---
|
|
if (pt_StabMembBiel != NULL)
|
|
delete pt_StabMembBiel;
|
|
//---- suppression du dimensionnement des charges externes éventuellement ----
|
|
if (lesChargeExtSurEle != NULL) delete lesChargeExtSurEle;
|
|
};
|
|
|
|
// =========================== methodes publiques ===============================
|
|
|
|
// calcul si un point est a l'interieur de l'element ou non
|
|
// il faut que M est la dimension globale
|
|
// les trois fonctions sont pour l'etude a t=0, t et tdt
|
|
// retour : =0 le point est externe, =1 le point est interne ,
|
|
// = 2 le point est sur la frontière à la précision près
|
|
// coor_locales : s'il est différent de NULL, est affecté des coordonnées locales calculées,
|
|
// uniquement précises si le point est interne
|
|
// dans le cas où l'on veut un fonctionnement différent, c'est surchargé dans les
|
|
// éléments fils
|
|
int ElemMeca::Interne_0(const Coordonnee& M,Coordonnee* coor_locales)
|
|
{ PointM = &Met_abstraite::PointM_0;
|
|
BaseND = &Met_abstraite::BaseND_0;
|
|
return Interne(TEMPS_0,M,coor_locales);
|
|
};
|
|
int ElemMeca::Interne_t(const Coordonnee& M,Coordonnee* coor_locales)
|
|
{ PointM = &Met_abstraite::PointM_t;
|
|
BaseND = &Met_abstraite::BaseND_t;
|
|
return Interne(TEMPS_t,M,coor_locales);
|
|
};
|
|
int ElemMeca::Interne_tdt(const Coordonnee& M,Coordonnee* coor_locales)
|
|
{ PointM = &Met_abstraite::PointM_tdt;
|
|
BaseND = &Met_abstraite::BaseND_tdt;
|
|
return Interne(TEMPS_tdt,M,coor_locales);
|
|
};
|
|
|
|
// test si l'element est complet
|
|
int ElemMeca::TestComplet()
|
|
{ int ret = Element::TestComplet(); // teste de la class mere
|
|
if (loiComp == NULL)
|
|
{ cout << " \n la loi n'est pas defini dans l element " << num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
else if (loiComp->TestComplet() != 1)
|
|
{ cout << " \n la loi n'est pas complete dans l element " << num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
// dans le cas d'une loi thermo physique on vérifie qu'elle est complète
|
|
if (loiTP != NULL)
|
|
{ if (loiTP->TestComplet() != 1)
|
|
{ cout << " \n la loi thermo physique n'est pas complete dans l element " << num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
};
|
|
// dans le cas où l'on prend en compte de la dilatation, on vérifie qu'une loi thermo plysique
|
|
// a été défini
|
|
if (dilatation)
|
|
{if (loiTP == NULL)
|
|
{cout << " \n L element " << num_elt
|
|
<< " prend en compte la dilatation, alors qu'aucune loi thermo physique n'est defini, "
|
|
<< " il faut donc en definir une ! \n ";
|
|
ret = 0;
|
|
}
|
|
}
|
|
// dans le cas d'une loi de frottement on vérifie qu'elle est complète
|
|
if (loiFrot != NULL)
|
|
{ if (loiFrot->TestComplet() != 1)
|
|
{ cout << " \n la loi de frottement n'est pas complete dans l element " << num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
};
|
|
// suite des datas
|
|
if ( masse_volumique == masse_volumique_defaut)
|
|
{ cout << "\n la masse volumique n'est pas defini dans l element " << num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
if (met == NULL)
|
|
{ cout << " \n la metrique n'est pas defini dans l element " << num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
if (def == NULL)
|
|
{ cout << " \n le type de deformation n'est pas defini dans l element "
|
|
<< num_elt << '\n';
|
|
ret = 0;
|
|
}
|
|
return ret;
|
|
};
|
|
|
|
// test pour savoir si le calcul de contrainte en absolu est possible
|
|
bool ElemMeca::ContrainteAbsoluePossible()
|
|
{ // il faut qu'il soit complet et il faut que les coordonnées a t =0 et t existe
|
|
bool retour = true;
|
|
if (TestComplet() !=0)
|
|
{// on regarde la définition des coordonnées
|
|
int nbmaxnoeud=tab_noeud.Taille();
|
|
for (int i=1;i<=nbmaxnoeud;i++)
|
|
if((!tab_noeud(i)->ExisteCoord0())||(!tab_noeud(i)->ExisteCoord1()))
|
|
{ retour = false; break;}
|
|
}
|
|
else
|
|
retour = false;
|
|
return retour;
|
|
};
|
|
|
|
|
|
// --------- calculs utils dans le cadre de la recherche du flambement linéaire
|
|
// dans un premier temps uniquement virtuelles, ensuite se sera virtuelle pure pour éviter
|
|
// les oubli de définition ----> AMODIFIER !!!!!!!!
|
|
// Calcul de la matrice géométrique et initiale
|
|
ElemMeca::MatGeomInit ElemMeca::MatricesGeometrique_Et_Initiale (const ParaAlgoControle & )
|
|
{ cout << "\n attention le calcul de la matrice g\'eom\'etrique et de la matrice "
|
|
<< " initiale ne sont pas implant\'ees \n" ;
|
|
this->Affiche(1);
|
|
cout << "\nvoid ElemMeca::MatricesGeometrique_Et_Initiale (Mat_pleine &,Mat_pleine &)"
|
|
<< endl;
|
|
Sortie(1);
|
|
return MatGeomInit(NULL,NULL); // pour supprimer le warning
|
|
};
|
|
|
|
// 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
|
|
// = 1 : erreur = (int (delta sigma):(delta sigma) dv)/(int sigma:sigma dv)
|
|
// le numerateur et le denominateur sont tel que :
|
|
// errElemRelative = numerateur / denominateur , si denominateur different de 0
|
|
// sinon denominateur = numerateur si numerateur est different de 0, sinon
|
|
// tous sont nuls mais on n'effectue pas la division
|
|
//!!!!!!!!!!!!! pour l'instant en virtuelle il faudra après en
|
|
// virtuelle pure !!!!!!!!!!!!!!!!!!!!!!!!!!! donc le code est à virer
|
|
void ElemMeca::ErreurElement(int ,double& ,double& , double& )
|
|
{ cout << "\n attention le calcul de l'erreur sur l'element n'est pas implant\'e "
|
|
<< " pour cet element \n" ;
|
|
this->Affiche(1);
|
|
cout << "\nElemMeca::ErreurElement(int....."
|
|
<< endl;
|
|
Sortie(1);
|
|
};
|
|
|
|
// ajout des ddl d'erreur pour les noeuds de l'élément
|
|
void ElemMeca::Plus_ddl_Erreur()
|
|
{int nbne = tab_noeud.Taille();
|
|
for (int ne=1; ne<= nbne; ne++) // pour chaque noeud
|
|
{ Ddl ddl(ERREUR,0.,LIBRE);
|
|
tab_noeud(ne)->PlusDdl(ddl);
|
|
}
|
|
};
|
|
// inactive les ddl d'erreur
|
|
void ElemMeca::Inactive_ddl_Erreur()
|
|
{ int nbne = tab_noeud.Taille();
|
|
for (int ne=1; ne<= nbne; ne++)
|
|
tab_noeud(ne)->Met_hors_service(ERREUR);
|
|
};
|
|
// active les ddl d'erreur
|
|
void ElemMeca::Active_ddl_Erreur()
|
|
{ int nbne = tab_noeud.Taille();
|
|
for (int ne=1; ne<= nbne; ne++)
|
|
tab_noeud(ne)->Met_en_service(ERREUR);
|
|
};
|
|
// retourne un tableau de ddl element, correspondant à la
|
|
// composante d'erreur -> ERREUR, pour chaque noeud
|
|
// -> utilisé pour l'assemblage de la raideur d'erreur
|
|
//dans cette version tous les noeuds sont supposés avoi un ddl erreur
|
|
// dans le cas contraire il faut redéfinir la fonction dans l'élément terminal
|
|
DdlElement ElemMeca::Tableau_de_ERREUR() const
|
|
{ return DdlElement(tab_noeud.Taille(), DdlNoeudElement(ERREUR) );
|
|
};
|
|
|
|
// =========================== methodes protegees ===============================
|
|
|
|
|
|
// calcul si un point est a l'interieur de l'element ou non
|
|
// il faut que M est la dimension globale
|
|
// retour : = 0 le point est externe, = 1 le point est interne ,
|
|
// = 2 le point est sur la frontière à la précision près
|
|
// coor_locales : s'il est différent de NULL, est affecté des coordonnées locales calculées,
|
|
|
|
int ElemMeca::Interne(Enum_dure temps,const Coordonnee& M,Coordonnee* coor_locales)
|
|
{ int dim_glob = M.Dimension(); // dim de M = dim de l'espace géométrique
|
|
ElemGeomC0& elgeom = this->ElementGeometrique(); // la geometrie
|
|
int dim_loc = elgeom.Dimension(); // dim de l'élément de référence
|
|
#ifdef MISE_AU_POINT
|
|
if (ParaGlob::Dimension() != dim_glob)
|
|
{ cout << "\n *** erreur, la dimension du point = (" << dim_glob
|
|
<< ") et de l'espace geometrique (" << ParaGlob::Dimension() << ") sont differentes! ";
|
|
cout << "\nElemMeca::Interne(Coordonnee& M) " << endl;
|
|
this->Affiche();
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// la methode consiste tout d'abord a calculer les coordonnees locales correspondants
|
|
// a M. Pour cela on utilise un processus iteratif
|
|
// ---- 1 :init de la recherche
|
|
const Coordonnee refL(dim_loc) ; // coordonnees locales du point de reference, initialisees a zero
|
|
// - construction du point local en fonction de l'existance ou non de coor_locales
|
|
Coordonnee* pt_ret=NULL; Coordonnee theta_local(refL);
|
|
if (coor_locales != NULL) {pt_ret = coor_locales;} else {pt_ret = &theta_local;};
|
|
Coordonnee& theta = *pt_ret; // init de AL qui representera les coordonnes locales de M
|
|
// - fin construction du point local
|
|
Vecteur ph_i = elgeom.Phi_point(refL); // fonctions d'interpolation en refL
|
|
Mat_pleine dph_i = elgeom.Dphi_point(refL); // derivees des fonctions d'interpolation en refL
|
|
// l'image de theta par l'interpolation, de dimension : espace géométrique
|
|
Coordonnee A = (met->*PointM) (tab_noeud,ph_i); //
|
|
// calcul de la base naturelle et duale en refL, en fonction des coord a t
|
|
BaseB bB(dim_glob,dim_loc);BaseH bH(dim_glob,dim_loc);
|
|
(met->*BaseND)(tab_noeud,dph_i,ph_i,bB,bH);
|
|
// calcul des coordonnees locales du point M dans le repere en refl
|
|
Coordonnee AM = M - A; // en global
|
|
for (int i=1;i<=dim_loc;i++)
|
|
theta(i) = AM * bH(i).Coor(); // ok même si dim_loc < dim_glob -> projection sur le plan tangent
|
|
Coordonnee theta_repere(theta); // le point où l'on calcul le repère
|
|
Coordonnee theta_initial(theta); Coordonnee delta_theta(theta);
|
|
// dans le cas où le point est externe à l'élément, on limite le repère de calcul au point externe
|
|
// de l'élément dans la direction de theta
|
|
if (!(elgeom.Interieur(theta)))
|
|
{ // calcul d'un point extreme de l'élément dans le sens de M
|
|
theta_repere = elgeom.Maxi_Coor_dans_directionGM(theta);
|
|
};
|
|
// calcul du point correspondant à theta dans l'element, qui remplace A
|
|
ph_i = elgeom.Phi_point(theta_repere); // fonctions d'interpolation en A
|
|
Coordonnee A_sauve= A;
|
|
A = (met->*PointM) (tab_noeud,ph_i);
|
|
AM = M-A; // nouveau gap
|
|
|
|
// ---- 2 :iteration
|
|
// on boucle avec un test de stabilité des coordonnées locales
|
|
int compteur = 1; // contrôle du maxi d'itération
|
|
int nb_exterieur=0; // contrôle du maxi de fois que l'on trouve consécutivement le point à l'extérieur
|
|
// - récup des paramètres de controle
|
|
double delta_thetai_maxi = ParaGlob::param->ParaAlgoControleActifs().PointInterneDeltaThetaiMaxi();
|
|
int nb_boucle_sur_delta_thetai= ParaGlob::param->ParaAlgoControleActifs().PointInterneNbBoucleSurDeltaThetai();
|
|
int nb_max_a_l_exterieur= ParaGlob::param->ParaAlgoControleActifs().PointInterneNbExterne();
|
|
double prec_tehetai_interne= ParaGlob::param->ParaAlgoControleActifs().PointInternePrecThetaiInterne();
|
|
// - boucle de recherche des thetai
|
|
while (delta_theta.Norme() >= delta_thetai_maxi)
|
|
{dph_i = elgeom.Dphi_point(theta_repere); // derivees des fonctions d'interpolation en theta_repere
|
|
ph_i = elgeom.Phi_point(theta_repere); // fonctions d'interpolation en theta_repere (déjà calculé la première fois)
|
|
// calcul de la base naturelle et duale en fonction des coord a t
|
|
(met->*BaseND)(tab_noeud,dph_i,ph_i,bB,bH);
|
|
// amelioration des coordonnees locales du point M dans le repere en theta
|
|
for (int i=1;i<=dim_loc;i++)
|
|
delta_theta(i)=AM * bH(i).Coor();
|
|
// theta += delta_theta; erreur je pense: le 11 juin 2007
|
|
// le nouveau theta = la position du repère + le deltat theta
|
|
theta = theta_repere + delta_theta;
|
|
|
|
theta_repere=theta;
|
|
if (!(elgeom.Interieur(theta))) // si le point est à l'extérieur
|
|
{ // calcul d'un point extreme de l'élément dans le sens de M
|
|
theta_repere = elgeom.Maxi_Coor_dans_directionGM(theta);
|
|
nb_exterieur++;
|
|
}
|
|
else // si nb_exterieur est diff de zero, on le remet à zéro, car on est de nouveau à l'intérieur
|
|
{ if (nb_exterieur != 0) nb_exterieur=0; };
|
|
|
|
// si cela fait nb_max_a_l_exterieur fois que l'on est à l'extérieur on arrête
|
|
if ( nb_exterieur >= nb_max_a_l_exterieur) break;
|
|
// calcul du point correspondant dans l'element, qui remplace theta
|
|
ph_i = elgeom.Phi_point(theta_repere); // fonctions d'interpolation en theta
|
|
A = (met->*PointM) (tab_noeud,ph_i);
|
|
AM = M-A; // nouveau gap
|
|
compteur++;
|
|
if (compteur > nb_boucle_sur_delta_thetai)
|
|
{ // cela signifie que l'on n'arrive pas à converger, petit message si besoin
|
|
// et on choisit le point initial pour le teste
|
|
if (ParaGlob::NiveauImpression() >= 4)
|
|
cout << " \n ** warning, l'algorithme de recherche d'un point a l'interieur d'un "
|
|
<< "element ne converge pas, utilisation du repere a l'origine"
|
|
<< "\nElemMeca::Interne(Coordonnee& M)" << endl;
|
|
if (elgeom.Interieur(theta_initial)) {theta=theta_initial; }
|
|
else {theta=theta_initial; };
|
|
break;
|
|
};
|
|
};
|
|
// // à la sortie de la boucle on calcul les thétas correspondants à M final
|
|
// dph_i = elgeom.Dphi(theta_repere); // derivees des fonctions d'interpolation en theta
|
|
// ph_i = elgeom.Phi(theta_repere); // fonctions d'interpolation en theta
|
|
// // calcul de la base naturelle et duale en fonction des coord a t
|
|
// (met->*BaseND)(tab_noeud,dph_i,ph_i,bB,bH);
|
|
// // amelioration des coordonnees locales du point M dans le repere en theta
|
|
// for (int i=1;i<=dim_glob;i++)
|
|
// theta(i) += AM * bH(i).Vect();
|
|
//
|
|
// -- maintenant on regarde si le point est interne ou pas a l'element
|
|
int retour=0;
|
|
if (elgeom.Interieur(theta))
|
|
{ // on regarde si le point est intérieur d'une toute petite précision
|
|
// calcul d'un point extreme de l'élément dans le sens de M
|
|
theta_repere = elgeom.Maxi_Coor_dans_directionGM(theta);
|
|
double diff =(theta_repere-theta).Norme();
|
|
if (diff <= prec_tehetai_interne)
|
|
{retour = 2;} // !-! à l'erreur près, le point est sur la frontière
|
|
else {retour = 1;}; // !-! point à l'intérieur
|
|
}
|
|
else {return 0;}; // !-! point à l'extérieur -> retour directe
|
|
// --- cas particuliers ---
|
|
// -- dans le cas où "retour" est différent de 0, on regarde les cas particuliers
|
|
switch (dim_glob)
|
|
{case 3:
|
|
{switch (dim_loc)
|
|
{ case 3: break; // cas pas de pb
|
|
case 2: // cas des "coques ou plaques" on regarde si le point est réellement à l'intérieur (+- 1/2 épaisseur)
|
|
{ double epaisseur = Epaisseurs(temps,theta);
|
|
A = (met->*PointM) (tab_noeud,ph_i); // le point projeté sur la surface de l'élément
|
|
AM = M-A; double delta_hauteur = 0.5 * epaisseur - AM.Norme();
|
|
if (Dabs(delta_hauteur) <= prec_tehetai_interne)
|
|
{retour = 2;} // !-! à l'erreur près, le point est sur la frontière
|
|
else if (delta_hauteur < 0.)
|
|
{retour = 0;} // !-! point externe
|
|
else {retour = 1;}; // !-! point à l'intérieur
|
|
break;
|
|
}
|
|
default:
|
|
cout << " \n *** erreur : cas non implemente pour l'instant "
|
|
<< " dim_glob= " << dim_glob << " dim_loc= " << dim_loc
|
|
<< "\n ElemMeca::Interne(... ";
|
|
Sortie(1);
|
|
};// -- fin du switch dim_loc pour dim_glob=3
|
|
break;
|
|
} // fin dim_glob=3
|
|
case 2:
|
|
{switch (dim_loc)
|
|
{case 2: break; // pas de pb
|
|
case 1: case 3: // cas des poutres et volume (pb non consistance) pour l'instant non traité
|
|
// donc on laisse aller à la valeur par défaut
|
|
default:
|
|
cout << " \n *** erreur : cas non implemente pour l'instant "
|
|
<< " dim_glob= " << dim_glob << " dim_loc= " << dim_loc
|
|
<< "\n ElemMeca::Interne(... ";
|
|
Sortie(1);
|
|
};// -- fin du switch dim_loc pour dim_glob=3
|
|
break;
|
|
} // fin dim_glob=2
|
|
case 1: // la dimension locale est forcément 1, et pas de pb
|
|
{break;
|
|
}
|
|
default: // default pour le switch sur di_glob
|
|
cout << " \n *** erreur : cas non implemente pour l'instant "
|
|
<< " dim_glob= " << dim_glob << " dim_loc= " << dim_loc
|
|
<< "\n ElemMeca::Interne(... ";
|
|
Sortie(1);
|
|
break;
|
|
}; // -- fin du switch sur dim_glob
|
|
// retour
|
|
return retour;
|
|
};
|
|
|
|
// Calcul du residu local et de la raideur locale,
|
|
// pour le schema implicite
|
|
// cald_Dvirtuelle = indique si l'on doit calculer la dérivée de la vitesse de déformation virtuelle
|
|
void ElemMeca::Cal_implicit (DdlElement & tab_ddl
|
|
,Tableau <TenseurBB *> & d_epsBB,Tableau < Tableau2 <TenseurBB *> > d2_epsBB_
|
|
,Tableau <TenseurHH *>& d_sigHH,int nbint
|
|
,const Vecteur& poids,const ParaAlgoControle & pa,bool cald_Dvirtuelle)
|
|
{
|
|
int nbddl = tab_ddl.NbDdl();
|
|
double jacobien; // def du jacobien
|
|
volume = 0. ; // init
|
|
Vecteur d_jacobien_tdt;
|
|
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 (bulk_viscosity)
|
|
{ if (sig_bulkHH == NULL) {sig_bulkHH = NevezTenseurHH(*sigHH_t_1);}
|
|
else if (sig_bulkHH->Dimension() != sigHH_t_1->Dimension())
|
|
{delete sig_bulkHH; sig_bulkHH = NevezTenseurHH(*sigHH_t_1);};
|
|
};
|
|
// init éventuel des intégrales de volume et temps
|
|
Init_Integ_vol_et_temps();
|
|
/* 19 nov: en fait il ne faut pas prendre en compte l'épaisseur méca: cf. doc théorique
|
|
// --- init éventuel pour les épaisseurs et largeurs
|
|
double epaisseur_0 = EpaisseurMoyenne(TEMPS_0 );
|
|
double epaisseur = epaisseur_0;double epaisseur_moyenne = 0.; // init
|
|
*/
|
|
// --- 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
|
|
|
|
for (int ni =1; ni<= nbint; ni++) // boucle sur les pt d'integ
|
|
{def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
|
|
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
|
|
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
|
|
TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
|
|
EnergieMeca& energ = tab_energ(ni);
|
|
EnergieMeca& energ_t = tab_energ_t(ni);
|
|
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
|
|
if (loiTP != NULL) {sTP = tabSaveTP(ni);}; // 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(),ni);
|
|
|
|
// -- 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(ni)
|
|
, *def,tab_ddl,ptIntegMeca,d_epsBB
|
|
,jacobien,d_jacobien_tdt,d_sigHH
|
|
,pa,sTP,loiTP,dilatation,energ,energ_t
|
|
,premier_calcul_meca_impli_expli);
|
|
// //---debug
|
|
// {double maxsig = sigHH.MaxiComposante();
|
|
// if ((!isfinite(maxsig)) || (isnan(maxsig)) )
|
|
// cout << "\n debug ElemMeca::Cal_implicit "
|
|
// << " maxsig= "
|
|
// << (maxsig)
|
|
// << flush;
|
|
//
|
|
//
|
|
// }
|
|
// //--- fin debug
|
|
|
|
|
|
}
|
|
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= "<<ni << 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= "<<ni << 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 = d2_epsBB_(ni); // 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
|
|
DeuxDoubles delta_ener_bulk_vol; // init
|
|
if (bulk_viscosity)
|
|
delta_ener_bulk_vol=ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_tdt,(*sig_bulkHH) );
|
|
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) ;}
|
|
};
|
|
//// -- pour le debug
|
|
// if ((pa.Variables_de_temps().TempsCourant()==0.021)&&(Num_elt()==55)&&(ni==7))
|
|
// {// on affiche les infos des noeuds
|
|
// for (int i=1;i<= tab_noeud.Taille();i++) tab_noeud(i)->Affiche();
|
|
// // on affiche les vecteurs de bases à 0 et t
|
|
// cout << "\n gi_0: " << *(ex.giB_0)<< "\n gi_tdt: " << *(ex.giB_tdt);
|
|
// cout << "\n debug : ElemMeca::Cal_implicit ";
|
|
// };
|
|
////--fin debug
|
|
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;
|
|
// là on ne peut plus rien faire donc on génère une exception
|
|
switch (pa.JabobienNegatif())
|
|
{ case 1:
|
|
{ cout << "\n on annulle sa contribution. Element nb: " << Num_elt() << " nb d'integ : " << ni;
|
|
break;
|
|
}
|
|
case 2:
|
|
{ // on génère une exception
|
|
throw ErrJacobienNegatif_ElemMeca();break;
|
|
}
|
|
// dans les autres cas on ne fait rien
|
|
};
|
|
};
|
|
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() << " nb d'integ : " << ni;
|
|
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() >= 6)
|
|
{// on affiche les infos des noeuds
|
|
for (int i=1;i<= tab_noeud.Taille();i++) tab_noeud(i)->Affiche();
|
|
// on affiche les vecteurs de bases à 0 et t
|
|
cout << "\n gi_0: " << *(ex.giB_0)<< "\n gi_tdt: " << *(ex.giB_tdt);
|
|
};
|
|
if (ParaGlob::NiveauImpression() >= 3)
|
|
{ cout << "\n *** attention la variation maximal du jacobien est atteinte !! *** "
|
|
<< "\n ele= "<< Num_elt() << " nbi= " << ni << " 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
|
|
double poid_jacobien= poids(ni) * jacobien;
|
|
// calcul des intégrales de volume et de volume + temps
|
|
// cas = 1 : pour un appel après un calcul implicit
|
|
// cas = 2 : pour un appel après un calcul explicit
|
|
Calcul_Integ_vol_et_temps(1,poid_jacobien,ni);
|
|
|
|
/* 19 nov: en fait il ne faut pas prendre en compte l'épaisseur méca: cf. doc théorique
|
|
//---cas éventuel de la prise en compte de la variation d'une épaisseur, ou section
|
|
switch (loiComp->Comportement_3D_CP_DP_1D())
|
|
{ case COMP_CONTRAINTES_PLANES :
|
|
// en contrainte plane on recalcule l'épaisseur
|
|
{epaisseur = epaisseur_0 * loiComp->HsurH0(tabSaveDon(ni));
|
|
epaisseur_moyenne += epaisseur;
|
|
poid_jacobien *= epaisseur;
|
|
break;
|
|
}
|
|
case COMP_DEFORMATIONS_PLANES :
|
|
// en deformation plane, l'épaisseur ne change pas
|
|
{epaisseur = epaisseur_0 ;
|
|
epaisseur_moyenne += epaisseur;
|
|
poid_jacobien *= epaisseur;
|
|
break;
|
|
}
|
|
default : // on n'a rien à faire
|
|
break;
|
|
};
|
|
*/
|
|
volume += poid_jacobien;ptIntegMeca.Volume_pti()=poid_jacobien;
|
|
// 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,poid_jacobien);
|
|
// energie_totale += (energ-energ_t) * poid_jacobien;
|
|
E_elem_bulk_tdt += delta_ener_bulk_vol.un * poid_jacobien;
|
|
P_elem_bulk += delta_ener_bulk_vol.deux * poid_jacobien;
|
|
// calcul du résidu et variation
|
|
for (int j =1; j<= nbddl; j++) // 1ere boucle sur les ddl
|
|
{(*residu)(j) -= ((*(d_epsBB(j))) && sigHH ) * poid_jacobien;
|
|
for (int i =1; i<= nbddl; i++) // 2ere boucle sur les ddl
|
|
(*raideur)(j,i) += ( ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien
|
|
+ (sigHH && (*(d2_epsBB(j,i)))) * jacobien
|
|
+ (sigHH && (*(d_epsBB(j)))) * d_jacobien_tdt(i)
|
|
) * poids(ni);
|
|
};
|
|
};
|
|
// --- 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()))
|
|
Cal_implicit_StabMembBiel(ni,*ex_final,nbint,poids(ni),NULL);
|
|
////debug
|
|
//cout << "\n ElemMeca::Cal_implicit debug : ";
|
|
//raideur->Affiche();
|
|
//// fin debug
|
|
|
|
};// fin du if sur la valeur non négative du jacobien
|
|
if (bulk_viscosity) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
|
|
{ sigHH -= (*sig_bulkHH); } // 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, nbint,volume,NULL);}
|
|
else {Cal_implicit_StabMembBiel(-1,*ex_final, nbint,volume,NULL);} ;
|
|
};
|
|
#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 ---
|
|
//if ((num_elt == 3)&&((integ_vol_typeQuel != NULL) || (integ_vol_t_typeQuel != NULL)))
|
|
// {cout << "\n debug ElemMeca::Cal_implicit ( ";
|
|
// if (integ_vol_typeQuel != NULL)
|
|
// {int taille = integ_vol_typeQuel->Taille();
|
|
// Tableau <TypeQuelconque> & tab_integ = *integ_vol_typeQuel; // pour simplifier
|
|
// for (int il =1;il <= taille ;il++)
|
|
// cout << "\n tab_integ("<<il<<")= "<< *(tab_integ(il).Grandeur_pointee());
|
|
// };
|
|
// // 2) intégration de volume et en temps
|
|
// if (integ_vol_t_typeQuel != NULL)
|
|
// {int taille = integ_vol_t_typeQuel->Taille();
|
|
// Tableau <TypeQuelconque> & tab_integ = *integ_vol_t_typeQuel; // pour simplifier
|
|
// Tableau <TypeQuelconque> & tab_integ_t = *integ_vol_t_typeQuel_t; // ""
|
|
// for (int il =1;il <= taille ;il++)
|
|
// cout << "\n tab_integ_et_temps("<<il<<")= "<< *(tab_integ(il).Grandeur_pointee());
|
|
// };
|
|
// cout << endl;
|
|
// };
|
|
//
|
|
//
|
|
////---- fin debug ----
|
|
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
};
|
|
|
|
// Calcul uniquement du residu local a l'instant t ou pas suivant la variable atdt
|
|
// pour le schema explicit par exemple
|
|
void ElemMeca::Cal_explicit (DdlElement & tab_ddl,Tableau <TenseurBB *> & d_epsBB
|
|
,int nbint,const Vecteur& poids,const ParaAlgoControle & pa,bool atdt)
|
|
|
|
{ double jacobien; // def du jacobien
|
|
int nbddl = tab_ddl.NbDdl();
|
|
energie_totale.Initialisation_differenciee(energie_totale_t); // init de l'énergie totale sur l'élément
|
|
volume = 0. ; // init
|
|
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 (bulk_viscosity)
|
|
{ if (sig_bulkHH == NULL) {sig_bulkHH = NevezTenseurHH(*sigHH_t_1);}
|
|
else if (sig_bulkHH->Dimension() != sigHH_t_1->Dimension())
|
|
{delete sig_bulkHH; sig_bulkHH = NevezTenseurHH(*sigHH_t_1);};
|
|
};
|
|
// --- init éventuel des intégrales de volume et temps
|
|
Init_Integ_vol_et_temps();
|
|
/* 19 nov: en fait il ne faut pas prendre en compte l'épaisseur méca: cf. doc théorique
|
|
// --- init éventuel pour les épaisseurs et largeurs
|
|
double epaisseur_0 = EpaisseurMoyenne(TEMPS_0 );
|
|
double epaisseur = epaisseur_0;double epaisseur_moyenne = 0.; // init
|
|
*/
|
|
|
|
// --- init pour l'init de la stabilisation de membrane-biel éventuelle
|
|
const Met_abstraite::Expli_t_tdt* 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
|
|
|
|
for (int ni =1; ni<= nbint; ni++) // boucle sur les pt d'integ
|
|
{ def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
|
|
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
|
|
TenseurHH & sigHH_t = *(ptIntegMeca.SigHH_t());
|
|
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
|
|
TenseurBB & epsBB = *(ptIntegMeca.EpsBB());
|
|
TenseurBB & DepsBB_ = *(ptIntegMeca.DepsBB());
|
|
EnergieMeca& energ = tab_energ(ni);
|
|
EnergieMeca& energ_t = tab_energ_t(ni);
|
|
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
|
|
if (loiTP != NULL) {sTP = tabSaveTP(ni);}; // 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(),ni);
|
|
DeuxDoubles delta_ener_bulk_vol; // init
|
|
|
|
Met_abstraite::Expli_t_tdt ex_inter; // une grandeur intermédiaire
|
|
// -- 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(ni), *def,tab_ddl,ptIntegMeca,d_epsBB,jacobien
|
|
,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
|
|
ex_inter= ex;
|
|
ex_final = &ex_inter; // stockage pour la stabilisation éventuelle de membraneBiel
|
|
|
|
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
|
|
if (bulk_viscosity)
|
|
delta_ener_bulk_vol=ModifContrainteBulk(TEMPS_tdt,*ex.gijHH_tdt,sigHH,DepsBB_,(*sig_bulkHH) );
|
|
}
|
|
else
|
|
{ const Met_abstraite::Expli& ex = loiComp->Cal_explicit_t
|
|
(tabSaveDon(ni), *def,tab_ddl,ptIntegMeca,d_epsBB,jacobien
|
|
,sTP,loiTP,dilatation,energ,energ_t,premier_calcul_meca_impli_expli);
|
|
ex_inter = ex.T_dans_tdt();
|
|
ex_final = &ex_inter; // stockage pour la stabilisation éventuelle de membraneBiel
|
|
// examen du cas où on désire utiliser la méthode d'atténuation des hautes fréquences avec le bulk viscosity
|
|
if (bulk_viscosity)
|
|
delta_ener_bulk_vol=ModifContrainteBulk(TEMPS_t,*ex.gijHH_t,sigHH,DepsBB_,(*sig_bulkHH) );
|
|
};
|
|
}
|
|
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= "<<ni << 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= "<<ni << endl;
|
|
throw (Err_inconnue_ElemMeca());
|
|
};
|
|
|
|
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_inter.jacobien_0))) {rap = Dabs(jacobien) / (*(ex_inter.jacobien_0));}
|
|
else {rap = (*(ex_inter.jacobien_0)) / Dabs(jacobien) ;}
|
|
};
|
|
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;
|
|
// là on ne peut plus rien faire donc on génère une exception
|
|
switch (pa.JabobienNegatif())
|
|
{ case 1:
|
|
{ cout << "\n on annulle sa contribution. Element nb: " << Num_elt() << " nb d'integ : " << ni;
|
|
break;
|
|
}
|
|
case 2:
|
|
{ // on génère une exception
|
|
throw ErrJacobienNegatif_ElemMeca();break;
|
|
}
|
|
// dans les autres cas on ne fait rien
|
|
};
|
|
};
|
|
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() << " nb d'integ : " << ni;
|
|
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() >= 6)
|
|
{// on affiche les infos des noeuds
|
|
for (int i=1;i<= tab_noeud.Taille();i++) tab_noeud(i)->Affiche();
|
|
// on affiche les vecteurs de bases à 0 et t
|
|
cout << "\n gi_0: " << *(ex_inter.giB_0)<< "\n gi_tdt: " << *(ex_inter.giB_tdt);
|
|
};
|
|
if (ParaGlob::NiveauImpression() >= 3)
|
|
{ cout << "\n *** attention la variation maximal du jacobien est atteinte !! *** "
|
|
<< "\n ele= "<< Num_elt() << " nbi= " << ni << " jacobien= " << jacobien << " jacobien_0= "
|
|
<< (*(ex_inter.jacobien_0)) << endl; };
|
|
// on génère une exception
|
|
throw ErrVarJacobienMini_ElemMeca();break;
|
|
}
|
|
else
|
|
{ double poid_jacobien= poids(ni) * jacobien;
|
|
// calcul des intégrales de volume et de volume + temps
|
|
// cas = 1 : pour un appel après un calcul implicit
|
|
// cas = 2 : pour un appel après un calcul explicit
|
|
Calcul_Integ_vol_et_temps(2,poid_jacobien,ni);
|
|
/* 19 nov: en fait il ne faut pas prendre en compte l'épaisseur méca: cf. doc théorique
|
|
//---cas éventuel de la prise en compte de la variation d'une épaisseur, ou section
|
|
switch (loiComp->Comportement_3D_CP_DP_1D())
|
|
{ case COMP_CONTRAINTES_PLANES :
|
|
// en contrainte plane on recalcule l'épaisseur
|
|
{epaisseur = epaisseur_0 * loiComp->HsurH0(tabSaveDon(ni));
|
|
epaisseur_moyenne += epaisseur;
|
|
poid_jacobien *= epaisseur;
|
|
break;
|
|
}
|
|
case COMP_DEFORMATIONS_PLANES :
|
|
// en deformation plane, l'épaisseur ne change pas
|
|
{epaisseur = epaisseur_0 ;
|
|
epaisseur_moyenne += epaisseur;
|
|
poid_jacobien *= epaisseur;
|
|
break;
|
|
}
|
|
default : // on n'a rien à faire
|
|
break;
|
|
};
|
|
*/
|
|
// 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,poid_jacobien);
|
|
// energie_totale += (energ-energ_t) * poid_jacobien;
|
|
// energie_totale += energ * poid_jacobien;
|
|
////debug
|
|
//cout << "\n debug: ElemMeca::Cal_explicit "
|
|
// << " energ visqueux = "<< energie_totale.DissipationVisqueuse() << " brut " << energ.DissipationVisqueuse()
|
|
// << " poid_jacobien = "<< poid_jacobien << endl;
|
|
//
|
|
////fin debug
|
|
|
|
volume += poid_jacobien;
|
|
ptIntegMeca.Volume_pti()=poid_jacobien;
|
|
// calcul du résidu si ce n'est pas une loi rien
|
|
if (!Loi_rien(loiComp->Id_comport()))
|
|
{E_elem_bulk_tdt += delta_ener_bulk_vol.un * poid_jacobien;
|
|
P_elem_bulk += delta_ener_bulk_vol.deux * poid_jacobien;
|
|
for (int j =1; j<= nbddl; j++)
|
|
{(*residu)(j) -= (sigHH && (*(d_epsBB(j)))) * poid_jacobien;
|
|
};
|
|
|
|
// --- cas éventuel d'une stabilisation membrane-biel ---
|
|
// ici il s'agit de la contribution précise à chaque pti
|
|
// pour une stabilisation via les normales aux pti
|
|
if (pt_StabMembBiel != NULL)
|
|
if (!(pt_StabMembBiel->Aa_calculer()))
|
|
Cal_explicit_StabMembBiel(ni,*ex_final,nbint,poids(ni),NULL);
|
|
|
|
////---debug
|
|
// #ifdef MISE_AU_POINT
|
|
// if (ParaGlob::NiveauImpression() > 2)
|
|
// if ((std::isinf((*residu)(j))||( std::isnan((*residu)(j)))))
|
|
// { cout << "\n ********** attention on a trouve une composante ("<<j<<") de residu infini ou nan !! = ("
|
|
// <<(*residu)(j)<<")********* " << endl;
|
|
// Sortie(2);
|
|
// };
|
|
// #endif
|
|
////--- fin debug
|
|
};
|
|
|
|
// --- cas éventuel d'une stabilisation membrane-biel ---
|
|
// modif: 10 janvier 2019 : je ne vois pas l'intérêt de faire le calcul à chaque pti
|
|
// car en fait le facteur de stabilisation est calculé à partir du résidu, or celui-ci
|
|
// n'est pas fini, et il commule à chaque nouveau pti, du coup le facteur est calculé à chaque pti
|
|
// via un résidu qui est transitoire et surtout qui ne contient pas que ce qui est calculé au pti !!
|
|
// donc ce n'est pas cohérent
|
|
// // ici il s'agit d'un calcul précis à chaque pti
|
|
// if (pt_StabMembBiel != NULL)
|
|
// if (pt_StabMembBiel->Aa_calculer())
|
|
// Cal_explicit_StabMembBiel(*ex_final);
|
|
};
|
|
if (bulk_viscosity) // dans le cas du bulk on retire la contrainte du bulk, qui est non physique
|
|
{ sigHH -= (*sig_bulkHH); } // 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 Second membre local: ? (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 Second membre local: elem= " << Num_elt_const()
|
|
<< ", maillage= " << Num_maillage();
|
|
cout << "\n second membre: " << (*residu);
|
|
};
|
|
};
|
|
#endif
|
|
/* 19 nov: en fait il ne faut pas prendre en compte l'épaisseur méca: cf. doc théorique
|
|
//---cas éventuel de la prise en compte de la variation d'une épaisseur, ou section
|
|
switch (loiComp->Comportement_3D_CP_DP_1D())
|
|
{ case COMP_CONTRAINTES_PLANES :
|
|
// en contrainte plane on recalcule l'épaisseur
|
|
{epaisseur_moyenne /= nbint;
|
|
Modifie_epaisseur_moyenne_tdt(epaisseur);
|
|
break;
|
|
}
|
|
default : // on n'a rien à faire
|
|
break;
|
|
};
|
|
*/
|
|
// --- intervention dans le cas d'une stabilisation d'hourglass
|
|
// stabilisation pour un calcul implicit
|
|
if (type_stabHourglass)
|
|
Cal_explicit_hourglass(atdt);
|
|
|
|
// --- 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)
|
|
{Cal_explicit_StabMembBiel(-1,*ex_final, nbint,volume,NULL);
|
|
};
|
|
#ifdef MISE_AU_POINT
|
|
if (ParaGlob::NiveauImpression() > 9)
|
|
{ if ((type_stabHourglass)||(pt_StabMembBiel != NULL))
|
|
{cout << "\n apres stabilisation: Second membre local: ? (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 Second membre local: elem= " << Num_elt_const()
|
|
<< ", maillage= " << Num_maillage();
|
|
cout << "\n second membre: " << (*residu);
|
|
};
|
|
};
|
|
};
|
|
#endif
|
|
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
};
|
|
|
|
// Calcul de la matrice géométrique et de la matrice initiale
|
|
// cette fonction est éventuellement appelée par les classes dérivées
|
|
// ddl represente les degres de liberte specifiques a l'element
|
|
// epsBB = deformation, sigHH = contrainte, d_epsbb = variation des def
|
|
// nbint = nb de pt d'integration , poids = poids d'integration
|
|
// cald_Dvirtuelle = indique si l'on doit calculer la dérivée de la vitesse de déformation virtuelle
|
|
void ElemMeca::Cal_matGeom_Init (Mat_pleine & matGeom, Mat_pleine & matInit,DdlElement & tab_ddl
|
|
,Tableau <TenseurBB *>& d_epsBB, Tableau < Tableau2 <TenseurBB *> > d2_epsBB_
|
|
,Tableau <TenseurHH *>& d_sigHH,int nbint,const Vecteur& poids
|
|
,const ParaAlgoControle & pa,bool cald_Dvirtuelle)
|
|
{ // le début du programme est semblable au calcul implicit
|
|
// seule la constitution des matrices finales est différente
|
|
int nbddl = tab_ddl.NbDdl();
|
|
double jacobien; // def du jacobien
|
|
Vecteur d_jacobien_tdt;
|
|
volume = 0. ; // init
|
|
|
|
for (int ni =1; ni<= nbint; ni++) // boucle sur les pt d'integ
|
|
{def->ChangeNumInteg(ni);def->Mise_a_jour_data_specif(tabSaveDefDon(ni));
|
|
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(ni);
|
|
TenseurHH & sigHH = *(ptIntegMeca.SigHH());
|
|
TenseurBB & DepsBB_tdt = *(ptIntegMeca.DepsBB());
|
|
EnergieMeca& energ = tab_energ(ni);
|
|
EnergieMeca& energ_t = tab_energ_t(ni);
|
|
CompThermoPhysiqueAbstraite::SaveResul* sTP=NULL; // les données spécifique thermo physiques
|
|
if (loiTP != NULL) {sTP = tabSaveTP(ni);}; // 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(),ni);
|
|
|
|
// -- 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
|
|
{loiComp->Cal_flamb_lin(tabSaveDon(ni), *def,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= "<<ni << 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= "<<ni << endl;
|
|
throw (Err_inconnue_ElemMeca());
|
|
};
|
|
|
|
// on calcul éventuellement la dérivée de la vitesse de déformation virtuelle
|
|
Tableau2 <TenseurBB *>& d2_epsBB = d2_epsBB_(ni); // pour simplifier
|
|
if (cald_Dvirtuelle)
|
|
def->Cal_var_def_virtuelle(false,d2_epsBB);
|
|
// calcul du volume
|
|
double poid_jacobien= poids(ni) * jacobien;
|
|
volume += poid_jacobien;ptIntegMeca.Volume_pti()=poid_jacobien;
|
|
for (int j =1; j<= nbddl; j++) // 1ere boucle sur les ddl
|
|
for (int i =1; i<= nbddl; i++) // 2ere boucle sur les ddl
|
|
{matInit(j,i) += ((*d_sigHH(i)) && (*(d_epsBB(j)))) * jacobien * poids(ni);
|
|
matGeom(j,i) += (sigHH && (*(d2_epsBB(j,i)))) * poid_jacobien;
|
|
}
|
|
}
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
};
|
|
|
|
// Calcul de la matrice masse selon différent choix donné par type_matrice_masse,
|
|
void ElemMeca::Cal_Mat_masse (DdlElement & tab_ddl,Enum_calcul_masse type_matrice_masse,
|
|
int nbint,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids)
|
|
{
|
|
int nbddl = tab_ddl.NbDdl();
|
|
double jacobien_0; // def du jacobien initial
|
|
// initialisation de la matrice masse
|
|
mat_masse->Initialise ();
|
|
|
|
// le calcul est fonction du type de matrice masse demandé
|
|
switch (type_matrice_masse)
|
|
{case MASSE_DIAG_COEF_EGAUX :
|
|
// même masse sur chaque noeud
|
|
{ // calcul du volume totale
|
|
double vol_totale =0.;
|
|
for (int ni =1; ni<= nbint; ni++) // boucle sur les pt d'integ
|
|
{ defMas->ChangeNumInteg(ni);
|
|
const Met_abstraite::Dynamiq& ex = defMas->Cal_matMass();
|
|
jacobien_0 = *ex.jacobien_0;
|
|
vol_totale += jacobien_0 * poids(ni);
|
|
}
|
|
// la masse élémentaire pour chaque noeud
|
|
double masse_elementaire = vol_totale * masse_volumique / tab_noeud.Taille();
|
|
// on initialise toutes les valeurs de la matrice masse à la valeur moyenne
|
|
mat_masse->Initialise (masse_elementaire );
|
|
break;
|
|
}
|
|
case MASSE_DIAG_COEF_VAR :
|
|
// masse diagonale répartie au prorata des fonctions d'interpolation
|
|
// (cf page 304, tome 2 Batoz)
|
|
{ // définition d'un vecteur qui contiendra l'intégral de chaque
|
|
// fonction d'interpolation (ici la masse est considérée constante)
|
|
Vecteur fonc2(nbne);
|
|
// calcul de fonc2 et du volume totale
|
|
double vol_totale =0.;
|
|
for (int ni =1; ni<= nbint; ni++) // boucle sur les pt d'integ
|
|
{ defMas->ChangeNumInteg(ni);
|
|
const Met_abstraite::Dynamiq& ex = defMas->Cal_matMass();
|
|
jacobien_0 = *ex.jacobien_0;
|
|
for (int ne =1;ne<=nbne;ne++)
|
|
fonc2(ne) += taphi(ni)(ne) * taphi(ni)(ne) * jacobien_0 * poids(ni);
|
|
vol_totale += jacobien_0 * poids(ni);
|
|
};
|
|
// calcul de la somme des composantes de fonc2
|
|
double somme_fonc2 = 0.;
|
|
for (int ne = 1;ne<=nbne;ne++) somme_fonc2 += fonc2(ne);
|
|
// calcul des termes diagonaux de la matrice de masse
|
|
int dimension = ParaGlob::Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'élément axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimension-1 coordonnées donc on décrémente
|
|
dimension--;
|
|
int ix=1;
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
{ double inter_toto = vol_totale * masse_volumique * fonc2(ne) / somme_fonc2;
|
|
for (int i=1;i<= dimension;i++,ix++)
|
|
(*mat_masse)(1,ix) = inter_toto;
|
|
};
|
|
break;
|
|
}
|
|
case MASSE_CONSISTANTE :
|
|
// matrice masse complète avec les mêmes fonctions d'interpolation
|
|
// que pour la raideur.
|
|
{
|
|
int dimension = ParaGlob::Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'élément axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimension-1 coordonnées donc on décrémente
|
|
dimension--;
|
|
for (int ni =1; ni<= nbint; ni++) // boucle sur les pt d'integ
|
|
{defMas->ChangeNumInteg(ni); // def du pt d'intégration
|
|
// récup des fonctions d'interpolation et du jacobien
|
|
const Met_abstraite::Dynamiq& ex = defMas->Cal_matMass();
|
|
double masse_i = (*ex.jacobien_0) * masse_volumique;
|
|
// calcul de la matrice
|
|
int ix=1;
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int a=1;a<= dimension;a++,ix++)
|
|
{ int jy=1;
|
|
for (int me =1; me<= nbne; me++)
|
|
for (int b=1;b<= dimension;b++,jy++)
|
|
if (a==b) // seule des termes de même indice sont non nulles
|
|
(*mat_masse)(ix,jy) += taphi(ni)(ne)* taphi(ni)(me) * masse_i * poids(ni);
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
default :
|
|
cout << "\nErreur : valeur incorrecte du type Enum_calcul_masse !\n";
|
|
cout << "ElemMeca::Cal_Mat_masse (... \n";
|
|
Sortie(1);
|
|
};
|
|
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
};
|
|
|
|
// ------------ calcul de second membre ou de l'action des efforts extérieures -----------
|
|
|
|
// initialisation éventuelle, nécessaire avant d'appliquer l'ensemble des charges
|
|
// par exemple des stockages intermédiaires
|
|
void ElemMeca::Initialisation_avant_chargement()
|
|
{ // init à 0 des efforts extérieurs éventuels
|
|
if (lesChargeExtSurEle != NULL)
|
|
lesChargeExtSurEle->Zero();
|
|
};
|
|
|
|
// calcul des seconds membres suivant les chargements
|
|
// cas d'un chargement surfacique de direction constante, sur les frontières des éléments
|
|
// force indique la force surfacique appliquée
|
|
// retourne le second membre résultant
|
|
// nSurf : le numéro de la surface externe
|
|
// calcul à tdt ou t suivant la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_surf_E (DdlElement & ,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & ,bool atdt)
|
|
{ // définition du vecteur de retour
|
|
Vecteur& SM = *((*res_extS)(nSurf));
|
|
Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesEffortsDirFixe() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesEffortsDirFixe_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesEffortsDirFixe_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_DirFixe= (*lesChargeExtSurEle->LesEffortsDirFixe())(nSurf); // pour simplifier
|
|
if (tab_DirFixe.Taille() != defS.Phi1_Taille())
|
|
tab_DirFixe.Change_taille(defS.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
int ni; // compteur globale de point d'integration
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
double jacobien;
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = NULL;
|
|
if (atdt)
|
|
{ const Met_abstraite::Expli_t_tdt& ex = defS.Cal_explicit_tdt(premier_calcul);
|
|
jacobien = (*ex.jacobien_tdt);
|
|
ex_expli_tdt = &ex;
|
|
}
|
|
else // normalement le cas à t n'est pas courant, donc on tolère une création
|
|
{ const Met_abstraite::Expli& ex = defS.Cal_explicit_t(premier_calcul);
|
|
jacobien = (*ex.jacobien_t);
|
|
ex_expli = &ex;
|
|
}
|
|
int ix=1; int dimf = force.Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
Coordonnee force_reelle(force); // init par défaut
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defS,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defS,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defS.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_reelle *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_reelle(3) *= tab_val(3);
|
|
case 2:force_reelle(2) *= tab_val(2);
|
|
case 1:force_reelle(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_surf_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " surfacique: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " surfacique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_DirFixe(ni)=force_reelle; // sauvegarde de l'effort
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* force_reelle(i) * (poids(ni) * jacobien);
|
|
}; // fin de la boucle sur les points d'intégraion
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
// calcul des seconds membres suivant les chargements
|
|
// cas d'un chargement surfacique de direction constante, sur les frontières des éléments
|
|
// force indique la force surfacique appliquée
|
|
// retourne le second membre résultant
|
|
// nSurf : le numéro de la surface externe
|
|
// -> implicite,
|
|
// pa : permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
// retourne le second membre et la matrice de raideur correspondant
|
|
Element::ResRaid ElemMeca::SMR_charge_surf_I (DdlElement & ddls,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & pa)
|
|
{ bool avec_raid = pa.Var_charge_externe(); // récup indicateur
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
int ni; // compteur globale de point d'integration
|
|
Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesEffortsDirFixe() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesEffortsDirFixe_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesEffortsDirFixe_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_DirFixe= (*lesChargeExtSurEle->LesEffortsDirFixe())(nSurf); // pour simplifier
|
|
if (tab_DirFixe.Taille() != defS.Phi1_Taille())
|
|
tab_DirFixe.Change_taille(defS.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
int nbddl = ddls.NbDdl();
|
|
Vecteur& SM = (*((*res_extS)(nSurf))); // " " "
|
|
Mat_pleine & KM = (*((*raid_extS)(nSurf))); // " " "
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en implicit
|
|
const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = &ex;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = NULL;
|
|
int ix=1; int dimf = force.Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
// on ne tient pas compte de la variation du point dans la fct nD
|
|
Coordonnee force_reelle(force); // init par défaut
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defS,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defS,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defS.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_reelle *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_reelle(3) *= tab_val(3);
|
|
case 2:force_reelle(2) *= tab_val(2);
|
|
case 1:force_reelle(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_surf_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " surfacique: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " surfacique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_DirFixe(ni)=force_reelle; // sauvegarde de l'effort
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) += taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
if (avec_raid)
|
|
for (int j =1; j<= nbddl; j++)
|
|
KM(ix,j) +=
|
|
taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j));
|
|
};
|
|
}
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = (*res_extS)(nSurf);
|
|
el.raid = (*raid_extS)(nSurf);
|
|
return el;
|
|
};
|
|
|
|
// calcul des seconds membres suivant les chargements
|
|
// cas d'un chargement pression, sur les frontières des éléments
|
|
// pression indique la pression appliquée
|
|
// retourne le second membre résultant
|
|
// nSurf : le numéro de la surface externe
|
|
// calcul à l'instant tdt ou t en fonction de la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_pres_E (DdlElement & ,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,double pression,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & ,bool atdt)
|
|
{
|
|
Deformation* definter=NULL; // variable de travail transitoire
|
|
Vecteur* SM_P = NULL;
|
|
if (!(ParaGlob::AxiSymetrie()))
|
|
{ definter = defSurf(nSurf); // le cas normal est non axisymétrique
|
|
SM_P = ((*res_extS)(nSurf));
|
|
}
|
|
else { definter = defArete(nSurf); // en axisymétrie, c'est une def d'arête
|
|
SM_P = ((*res_extA)(nSurf));
|
|
};
|
|
Deformation & defS = *definter; // pour simplifier l'écriture
|
|
Vecteur& SM = *SM_P; // " " "
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesPressionsExternes() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Pression_appliquee>& tab_Press= (*lesChargeExtSurEle->LesPressionsExternes())(nSurf); // pour simplifier
|
|
if (tab_Press.Taille() != defS.Phi1_Taille())
|
|
tab_Press.Change_taille(defS.Phi1_Taille());
|
|
|
|
// définition du vecteur de retour
|
|
int ni; // compteur globale de point d'integration
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
Met_abstraite::Expli ex;
|
|
if (atdt)
|
|
ex = ((defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t());
|
|
else
|
|
ex.operator=( defS.Cal_explicit_t(premier_calcul));
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = &ex;
|
|
// détermination de la normale à la surface
|
|
#ifdef MISE_AU_POINT
|
|
// test pour savoir si l'on a réellement affaire à une surface
|
|
if ((ex.giB_t->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3))
|
|
{ cout << "\n erreur, il doit y avoir deux vecteurs pour la surface de dimension 3"
|
|
<< "\n ici on a nbvec= " << ex.giB_t->NbVecteur() <<" vecteur(s) et on est en dimension "
|
|
<< " dim = " << ParaGlob::Dimension();
|
|
if (ParaGlob::NiveauImpression() > 2)
|
|
cout << "\n ElemMeca::SM_charge_pres (DdlElement & ...nbvec= ";
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// la normale vaut le produit vectoriel des 2 premiers vecteurs
|
|
Coordonnee normale = Util::ProdVec_coorBN( (*ex.giB_t)(1), (*ex.giB_t)(2));
|
|
|
|
// dans le cas où il y a une fonction nD, et qu'un de ces arguments
|
|
// est des coordonnées Xi on l'utilise
|
|
double pression_applique = pression;
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defS,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defS,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// if (pt_fonct->Depend_M())
|
|
// {// calcul de la pression au niveau de la position du point d'intégration
|
|
// // 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defS.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
|
|
// seule la première valeur est ok
|
|
pression_applique *= tab_val(1);
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " surfacique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
// calcul de la normale normée et pondérée de la pression
|
|
// la pression est positive qu'en elle appuis (diminution de volume)
|
|
normale.Normer();normale *= -pression_applique;
|
|
|
|
tab_Press(ni).P += normale;
|
|
int ix=1; int dimf = 3; // ne fonctionne qu'en dim 3
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_t));
|
|
// sauvegarde de la pression appliquée
|
|
tab_Press(ni).press += pression_applique;
|
|
|
|
//debug
|
|
//if (num_elt == 33839)
|
|
// cout << "\n debug: ElemMeca::SM_charge_pres_E "
|
|
// << " pression= " << pression << endl;
|
|
//if (tab_press_appliquee(ni).press_tdt > 120)
|
|
// cout << "\n debug: ElemMeca::SM_charge_pres_E "<< " on s'arrete ";
|
|
//fin debug
|
|
};
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
// calcul des seconds membres suivant les chargements
|
|
// cas d'un chargement pression, sur les frontières des éléments
|
|
// pression indique la pression appliquée
|
|
// retourne le second membre résultant
|
|
// nSurf : le numéro de la surface externe -> implicite,
|
|
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
Element::ResRaid ElemMeca::SMR_charge_pres_I(DdlElement & ddlS,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,double pression,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & pa)
|
|
{
|
|
int ni; // compteur globale de point d'integration
|
|
|
|
Deformation* definter=NULL; // variable de travail transitoire
|
|
Vecteur* SM_P = NULL; Mat_pleine* KM_P = NULL;
|
|
if (!(ParaGlob::AxiSymetrie()))
|
|
{ definter = defSurf(nSurf); // le cas normal est non axisymétrique
|
|
SM_P = ((*res_extS)(nSurf));
|
|
KM_P = ((*raid_extS)(nSurf));
|
|
}
|
|
else { definter = defArete(nSurf); // en axisymétrie, c'est une def d'arête
|
|
SM_P = ((*res_extA)(nSurf));
|
|
KM_P = ((*raid_extA)(nSurf));
|
|
};
|
|
Deformation & defS = *definter; // pour simplifier l'écriture
|
|
Vecteur& SM = *SM_P; // " " "
|
|
Mat_pleine & KM = *KM_P; // " " "
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesPressionsExternes() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Pression_appliquee>& tab_Press= (*lesChargeExtSurEle->LesPressionsExternes())(nSurf); // pour simplifier
|
|
if (tab_Press.Taille() != defS.Phi1_Taille())
|
|
tab_Press.Change_taille(defS.Phi1_Taille());
|
|
|
|
int nbddl = ddlS.NbDdl();
|
|
// controle
|
|
bool avec_raid = pa.Var_charge_externe();
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = &ex;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = NULL;
|
|
// détermination de la normale à la surface
|
|
#ifdef MISE_AU_POINT
|
|
// test pour savoir si l'on a réellement affaire à une surface
|
|
if ((ex.giB_tdt->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3))
|
|
{ cout << "\n erreur, il doit y avoir deux vecteurs pour la surface de dimension 3"
|
|
<< " ElemMeca::SM_charge_pres (DdlElement & ...nbvec= " << ex.giB_tdt->NbVecteur()
|
|
<< " dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// la normale vaut le produit vectoriel des 2 premiers vecteurs
|
|
CoordonneeB normale = Util::ProdVec_coorB( (*ex.giB_tdt)(1), (*ex.giB_tdt)(2));
|
|
// calcul de la variation de la normale
|
|
// 1) variation du produit vectoriel
|
|
Tableau <CoordonneeB > D_pasnormale =
|
|
Util::VarProdVect_coorB( (*ex.giB_tdt)(1),(*ex.giB_tdt)(2),(*ex.d_giB_tdt));
|
|
// 2) de la normale
|
|
Tableau <CoordonneeB> D_normale = Util::VarUnVect_coorB(normale,D_pasnormale,normale.Norme());
|
|
|
|
// dans le cas où il y a une fonction nD, et qu'un de ces arguments
|
|
// est des coordonnées Xi on l'utilise
|
|
double pression_applique = pression;
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defS,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defS,ex_impli,ex_expli_tdt,ex_expli);
|
|
// //---debug
|
|
// cout << "\n debug ElemMeca::SMR_charge_pres_I ";
|
|
// cout << li_quelc ;
|
|
// //-- fin debug
|
|
//
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// if (pt_fonct->Depend_M())
|
|
// {// calcul de la pression au niveau de la position du point d'intégration
|
|
// // 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defS.Position_tdt();
|
|
// // 2) on appelle la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// seule la première valeur est ok
|
|
pression_applique *= tab_val(1);
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " surfacique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
|
|
//*** important: pour l'instant on ne calcule pas la variation de la
|
|
// pression causée par la variation du point ==> ce qui est différent
|
|
// du calcul hydrostatique par exemple. Du coup, la convergence est
|
|
// sans doute un peu moins bonne pour des grands pas de temps ou des
|
|
// grandes variation de positions
|
|
|
|
// calcul de la normale normée et pondérée de la pression
|
|
// la pression est positive qu'en elle appuis (diminution de volume)
|
|
normale.Normer();normale *= -pression_applique;
|
|
tab_Press(ni).P += normale.Coor_const();
|
|
// également pondération de la variation de la
|
|
for (int ihi =1;ihi<= nbddl;ihi++)
|
|
D_normale(ihi) *= -pression_applique;
|
|
// sauvegarde de la pression appliquée
|
|
tab_Press(ni).press += pression_applique;
|
|
|
|
int ix=1; int dimf = 3; // ne fonctionne qu'en dim 3
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) += taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
if (avec_raid)
|
|
for (int j =1; j<= nbddl; j++)
|
|
KM(ix,j) -=
|
|
taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j))
|
|
+ taphi(ni)(ne)* D_normale(j)(i) * (poids(ni) * (*ex.jacobien_tdt)) ;
|
|
}
|
|
}
|
|
////-------- debug
|
|
//{cout << "\n debug: ElemMeca::SMR_charge_pres_I(..";
|
|
// cout << "\n effort externe de pression: Raideur et second membre : elem= " << Num_elt_const()
|
|
// << ", maillage= " << Num_maillage();
|
|
// cout << "\n raideur: ";
|
|
// KM.Affiche();
|
|
// cout << "\n second membre: " << SM;
|
|
//};
|
|
//
|
|
//
|
|
////------- fin debug
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = SM_P;
|
|
el.raid = KM_P;
|
|
return el;
|
|
};
|
|
|
|
// calcul des seconds membres suivant les chargements
|
|
// cas d'un chargement lineique, sur les aretes frontières des éléments
|
|
// force indique la force lineique appliquée
|
|
// retourne le second membre résultant
|
|
// nArete : le numéro de la surface externe
|
|
// calcul à l'instant tdt ou t en fonction de la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_line_E (DdlElement & ,int nArete
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & ,bool atdt)
|
|
{ // définition du vecteur de retour
|
|
Vecteur& SM = *((*res_extA)(nArete));
|
|
int ni; // compteur globale de point d'integration
|
|
Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesLineique() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesLineique_Change_taille(ElementGeometrique().Nbne());}
|
|
else
|
|
{lesChargeExtSurEle->LesLineique_Change_taille(ElementGeometrique().NbSe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_line= (*lesChargeExtSurEle->LesLineique())(nArete); // pour simplifier
|
|
if (tab_line.Taille() != defA.Phi1_Taille())
|
|
tab_line.Change_taille(defA.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
Met_abstraite::Expli ex;
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
if (atdt)
|
|
ex = (defA.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
|
|
else
|
|
ex = defA.Cal_explicit_t(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = &ex;
|
|
int ix=1; int dimf = force.Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
Coordonnee force_reelle(force); // init par défaut
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defA,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defA,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
try
|
|
{ Tableau <double> & tab_val =
|
|
pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_reelle *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_reelle(3) *= tab_val(3);
|
|
case 2:force_reelle(2) *= tab_val(2);
|
|
case 1:force_reelle(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_line_E(.. ";
|
|
ErrCalculFct_nD toto; throw (toto);
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " lineique: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
ErrCalculFct_nD toto; throw (toto);
|
|
Sortie(1);
|
|
};
|
|
}
|
|
catch (ErrCalculFct_nD excep)
|
|
{ if (ParaGlob::NiveauImpression() > 0)
|
|
cout << "\n *** erreur ElemMeca::SM_charge_line_E(.. "
|
|
<< " pb dans le calcul de la fonction nD associee ";
|
|
if (ParaGlob::param->ParaAlgoControleActifs().Cas_fctnD_charge() == 0)
|
|
{ErrCalculFct_nD toto; throw (toto);}
|
|
else // on neutralise la force
|
|
force_reelle *= 0.;
|
|
};
|
|
};
|
|
tab_line(ni)=force_reelle; // sauvegarde
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.jacobien_t));
|
|
}
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
// cas d'un chargement lineique, sur les arêtes frontières des éléments
|
|
// force indique la force lineique appliquée
|
|
// retourne le second membre résultant
|
|
// nArete : le numéro de l'arête externe -> implicite,
|
|
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
Element::ResRaid ElemMeca::SMR_charge_line_I(DdlElement & ddlA,int nArete
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & pa)
|
|
{ int ni; // compteur globale de point d'integration
|
|
Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
|
|
int nbddl = ddlA.NbDdl();
|
|
Vecteur& SM = (*((*res_extA)(nArete))); // " " "
|
|
Mat_pleine & KM = (*((*raid_extA)(nArete))); // " " "
|
|
// controle
|
|
bool avec_raid = pa.Var_charge_externe();
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesLineique() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesLineique_Change_taille(ElementGeometrique().Nbne());}
|
|
else
|
|
{lesChargeExtSurEle->LesLineique_Change_taille(ElementGeometrique().NbSe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_line= (*lesChargeExtSurEle->LesLineique())(nArete); // pour simplifier
|
|
if (tab_line.Taille() != defA.Phi1_Taille())
|
|
tab_line.Change_taille(defA.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en implicit
|
|
const Met_abstraite::Impli& ex = defA.Cal_implicit(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = &ex;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = NULL;
|
|
int ix=1; int dimf = force.Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
// on ne tient pas compte de la variation du point dans la fct nD
|
|
Coordonnee force_reelle(force); // init par défaut
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defA,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defA,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// {if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defA.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_reelle *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_reelle(3) *= tab_val(3);
|
|
case 2:force_reelle(2) *= tab_val(2);
|
|
case 1:force_reelle(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_line_I(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " lineique: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " lineique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_line(ni)=force_reelle; // sauvegarde
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) +=
|
|
taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
if (avec_raid)
|
|
for (int j =1; j<= nbddl; j++)
|
|
KM(ix,j) -=
|
|
taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j));
|
|
}
|
|
}
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = (*res_extA)(nArete);
|
|
el.raid = (*raid_extA)(nArete);
|
|
return el;
|
|
};
|
|
|
|
// cas d'un chargement lineique suiveur, sur les aretes frontières des éléments
|
|
// force indique la force lineique appliquée
|
|
// retourne le second membre résultant
|
|
// nArete : le numéro de la surface externe
|
|
// calcul à l'instant tdt ou t en fonction de la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_line_Suiv_E (DdlElement & ,int nArete
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & ,bool atdt)
|
|
{ // définition du vecteur de retour
|
|
Vecteur& SM = *((*res_extA)(nArete));
|
|
int ni; // compteur globale de point d'integration
|
|
Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesLineiqueSuiveuse() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesLineiqueSuiveuse_Change_taille(ElementGeometrique().Nbne());}
|
|
else
|
|
{lesChargeExtSurEle->LesLineiqueSuiveuse_Change_taille(ElementGeometrique().NbSe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_line= (*lesChargeExtSurEle->LesLineiqueSuiveuse())(nArete); // pour simplifier
|
|
if (tab_line.Taille() != defA.Phi1_Taille())
|
|
tab_line.Change_taille(defA.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
Met_abstraite::Expli ex;
|
|
if (atdt)
|
|
ex = (defA.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
|
|
else
|
|
ex = defA.Cal_explicit_t(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = &ex;
|
|
#ifdef MISE_AU_POINT
|
|
// test pour savoir si l'on est en dimension 2 ou en axy sinon pas possible
|
|
if ((ParaGlob::Dimension() != 2) & (!(ParaGlob::AxiSymetrie())))
|
|
{ cout << "\n erreur, pour l'instant seul la dimension 2 est permise avec les forces suiveuses"
|
|
<< "\n si pb contacter gerard Rio ! "
|
|
<< " ElemMeca::SM_charge_line_Suiv_E ( ..."
|
|
<< " \n dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// calcul du repère locale initiale orthonormée
|
|
// Coordonnee2 v1((*ex.giB_0)(1).Vect()); v1.Normer();
|
|
// on récupère que les deux premières coordonnées mais on travaille en 3D pour intégrer le cas axi
|
|
Coordonnee v1((*ex.giB_0).Coordo(1)); v1.Normer();
|
|
int dim = v1.Dimension();
|
|
Coordonnee v2(dim); v2(1)=-v1(2); v2(2)=v1(1);
|
|
// composante du vecteur force linéique dans le repère initiale
|
|
Coordonnee force_0(v1*force,v2*force);
|
|
//repère locale actuel orthonormée
|
|
Coordonnee giB_t = (*ex.giB_t)(1).Coor();
|
|
Coordonnee vf1(giB_t);
|
|
Coordonnee vf2(dim);vf2(1)=-vf1(2);vf2(2)=vf1(1);
|
|
// composantes actuelles de la force
|
|
Coordonnee force_t = force_0(1)*vf1 + force_0(2)*vf2;
|
|
// calcul du résidu
|
|
int ix=1; int dimf = 2 ; // ici uniquement en dim 2 sinon pb
|
|
// if(ParaGlob::AxiSymetrie())
|
|
// // cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// // dimf-1 coordonnées donc on décrémente
|
|
// dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defA,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defA,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// {if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defA.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_t *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_t(3) *= tab_val(3);
|
|
case 2:force_t(2) *= tab_val(2);
|
|
case 1:force_t(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_line_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " lineique suiveuse: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " lineique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_line(ni)=force_t; // sauvegarde
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_t));
|
|
}
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
// cas d'un chargement lineique suiveur, sur les arêtes frontières des éléments
|
|
// force indique la force lineique appliquée
|
|
// retourne le second membre résultant
|
|
// nArete : le numéro de l'arête externe -> implicite,
|
|
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
Element::ResRaid ElemMeca::SMR_charge_line_Suiv_I(DdlElement & ddlA,int nArete
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & pa)
|
|
{ int ni; // compteur globale de point d'integration
|
|
Deformation & defA = *defArete(nArete); // pour simplifier l'écriture
|
|
int nbddl = ddlA.NbDdl();
|
|
Vecteur& SM = (*((*res_extA)(nArete))); // " " "
|
|
Mat_pleine & KM = (*((*raid_extA)(nArete))); // " " "
|
|
// controle
|
|
bool avec_raid = pa.Var_charge_externe();
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesLineiqueSuiveuse() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesLineiqueSuiveuse_Change_taille(ElementGeometrique().Nbne());}
|
|
else
|
|
{lesChargeExtSurEle->LesLineiqueSuiveuse_Change_taille(ElementGeometrique().NbSe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_line= (*lesChargeExtSurEle->LesLineiqueSuiveuse())(nArete); // pour simplifier
|
|
if (tab_line.Taille() != defA.Phi1_Taille())
|
|
tab_line.Change_taille(defA.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defA.PremierPtInteg(), ni = 1;defA.DernierPtInteg();defA.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en implicit
|
|
const Met_abstraite::Impli& ex = defA.Cal_implicit(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = &ex;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = NULL;
|
|
#ifdef MISE_AU_POINT
|
|
// test pour savoir si l'on est en dimension 2 ou en axy sinon pas possible
|
|
if ((ParaGlob::Dimension() != 2) & (!(ParaGlob::AxiSymetrie())))
|
|
{ cout << "\n erreur, pour l'instant seul la dimension 2 est permise avec les forces suiveuses"
|
|
<< "\n si pb contacter gerard Rio ! "
|
|
<< " ElemMeca::SM_charge_line_Suiv_E ( ..."
|
|
<< " \n dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// calcul du repère locale initiale orthonormée
|
|
// Coordonnee2 v1((*ex.giB_0)(1).Vect()); v1.Normer();
|
|
// on récupère que les deux premières coordonnées mais on travaille en 3D pour intégrer le cas axi
|
|
Coordonnee v1((*ex.giB_0).Coordo(1)); v1.Normer();
|
|
int dim = v1.Dimension();
|
|
Coordonnee v2(dim); v2(1)=-v1(2); v2(2)=v1(1);
|
|
// composante du vecteur force linéique dans le repère initiale
|
|
Coordonnee force_0(v1*force,v2*force);
|
|
//repère locale actuel orthonormée
|
|
Coordonnee giB_tdt = (*ex.giB_tdt)(1).Coor();
|
|
Coordonnee vf1(giB_tdt);
|
|
double vf1Norme = vf1.Norme();vf1 /= vf1Norme;
|
|
Coordonnee vf2(dim);vf2(1)=-vf1(2);vf2(2)=vf1(1);
|
|
// composantes actuelles de la force
|
|
Coordonnee force_t = force_0(1)*vf1 + force_0(2)*vf2;
|
|
|
|
// calcul du résidu
|
|
int ix=1; int dimf = 2 ; // ici uniquement en dim 2 sinon pb
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
// on ne tient pas compte de la variation du point dans la fct nD
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defA,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defA,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// {if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defA.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_t *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_t(3) *= tab_val(3);
|
|
case 2:force_t(2) *= tab_val(2);
|
|
case 1:force_t(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_line_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " lineique suiveuse: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " lineique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_line(ni)=force_t; // sauvegarde
|
|
Coordonnee dvf1,dvf2,d_giB_tdt; // def de vecteurs intermediaires
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) += taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
if (avec_raid)
|
|
for (int j =1; j<= nbddl; j++)
|
|
{ // calcul des variations de la force
|
|
d_giB_tdt = ((((*ex.d_giB_tdt)(j))(1)).Coor());
|
|
dvf1 = Util::VarUnVect_coor // variation de vf1 par rapport au ddl j
|
|
(giB_tdt,d_giB_tdt,vf1Norme);
|
|
dvf2 = -(dvf1 * vf2) * vf1;
|
|
KM(ix,j) -= taphi(ni)(ne) * poids(ni)
|
|
* (force_t(i) *(*ex.d_jacobien_tdt)(j)
|
|
+ (*ex.jacobien_tdt) * ((force_0(1) * dvf1 + force_0(2) * dvf2))(i)
|
|
);
|
|
}
|
|
}
|
|
}
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = (*res_extA)(nArete);
|
|
el.raid = (*raid_extA)(nArete);
|
|
return el;
|
|
};
|
|
|
|
// cas d'un chargement surfacique suiveur, sur les surfaces de l'élément
|
|
// la direction varie selon le système suivant: on définit les coordonnées matérielles
|
|
// de la direction, ce qui sert ensuite à calculer les nouvelles directions. L'intensité
|
|
// elle est constante.
|
|
// force indique la force surfacique appliquée
|
|
// retourne le second membre résultant
|
|
// nSurf : le numéro de la surface externe
|
|
// calcul à l'instant tdt ou t en fonction de la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_surf_Suiv_E (DdlElement & ,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& forc,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & ,bool atdt)
|
|
{ // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
|
|
// *** sinon c'est intractable et long
|
|
// définition du vecteur de retour
|
|
Vecteur& SM = *((*res_extS)(nSurf));
|
|
double module_f = forc.Norme();
|
|
// dans le cas ou la force est de module nul, on retourne sans calcul
|
|
if (module_f < ConstMath::trespetit)
|
|
return SM;
|
|
#ifdef MISE_AU_POINT
|
|
// test : a priori non valide pour les éléments axi et pour un pb autre que 3D
|
|
if (( ParaGlob::Dimension() != 3) || ((ParaGlob::Dimension() == 3)&&(ParaGlob::AxiSymetrie())))
|
|
{ cout << "\n erreur, ce type de chargement n'est valide que pour la dimension 3 !! et non axi"
|
|
<< " ElemMeca::SM_charge_line_Suiv_E ( ..."
|
|
<< " \n dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
int ni; // compteur globale de point d'integration
|
|
Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesPressDir() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesPressDir_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesPressDir_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_Press= (*lesChargeExtSurEle->LesPressDir())(nSurf); // pour simplifier
|
|
if (tab_Press.Taille() != defS.Phi1_Taille())
|
|
tab_Press.Change_taille(defS.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
Met_abstraite::Expli ex;
|
|
if (atdt)
|
|
ex = (defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
|
|
else
|
|
ex = defS.Cal_explicit_t(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = NULL;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = &ex;
|
|
// détermination de la normale à la surface
|
|
// la normale vaut le produit vectoriel des 2 premiers vecteurs
|
|
Coordonnee normale_0 = Util::ProdVec_coorBN( (*ex.giB_0)(1), (*ex.giB_0)(2));
|
|
// on passe en BN pour dire que c'est du B en normal
|
|
const Coordonnee& giBN_t1 = (*ex.giB_t).Coordo(1); const Coordonnee& giBN_t2 = (*ex.giB_t).Coordo(2);
|
|
const Coordonnee& giB_01 = (*ex.giB_0).Coordo(1); const Coordonnee& giB_02 = (*ex.giB_0).Coordo(2);
|
|
Coordonnee normale_t = Util::ProdVec_coor( giBN_t1, giBN_t2);
|
|
|
|
// calcul de la normale normée
|
|
normale_0.Normer(); normale_t.Normer();
|
|
// calcul des composantes de direction de la force dans le repère local initial
|
|
Coordonnee dir_forceH_0(3);
|
|
// ici le fait d'utiliser la méthode .Vect() supprime la variance, mais justement
|
|
// on ne veut pas de vérification de variance qui est ici incohérente
|
|
dir_forceH_0(1) = giB_01 * forc; dir_forceH_0(2) = giB_02 * forc;
|
|
dir_forceH_0(3) = normale_0 * forc;
|
|
// calcul de la direction finale
|
|
Coordonnee dir_force_t = dir_forceH_0(1) * giBN_t1 + dir_forceH_0(2) * giBN_t2
|
|
+ dir_forceH_0(3) * normale_t;
|
|
dir_force_t.Normer();
|
|
// calcul de la force finale
|
|
Coordonnee force_t = (forc.Norme()) * dir_force_t;
|
|
// calcul du second membre
|
|
int ix=1; int dimf = 3; // ne fonctionne qu'en dim 3
|
|
// sauf pour les éléments axisymétriques pour lesquelles on ne considère que les deux premières composantes
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defS,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defS,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// {if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defS.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_t *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_t(3) *= tab_val(3);
|
|
case 2:force_t(2) *= tab_val(2);
|
|
case 1:force_t(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_surf_Suiv_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " surfacique suiveuse: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " surfacique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_Press(ni)=force_t; // sauvegarde de l'effort
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_t));
|
|
}
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
// idem SM_charge_surf_Suiv_E mais
|
|
// -> implicite,
|
|
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
// retourne le second membre et la matrice de raideur correspondant
|
|
Element::ResRaid ElemMeca::SMR_charge_surf_Suiv_I (DdlElement & ddls,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& forc,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & pa)
|
|
{ // *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
|
|
// *** sinon c'est intractable et long
|
|
double module_f = forc.Norme();
|
|
// dans le cas ou la force est de module nul, on retourne sans calcul
|
|
if (module_f < ConstMath::trespetit)
|
|
{ Element::ResRaid el;
|
|
el.res = (*res_extS)(nSurf);
|
|
el.raid = (*raid_extS)(nSurf);
|
|
return el;
|
|
}
|
|
// controle
|
|
bool avec_raid = pa.Var_charge_externe();
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
int ni; // compteur globale de point d'integration
|
|
Deformation & defS = *defSurf(nSurf); // pour simplifier l'écriture
|
|
int nbddl = ddls.NbDdl();
|
|
Vecteur& SM = (*((*res_extS)(nSurf))); // " " "
|
|
Mat_pleine & KM = (*((*raid_extS)(nSurf))); // " " "
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesPressDir() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesPressDir_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesPressDir_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Coordonnee>& tab_Press= (*lesChargeExtSurEle->LesPressDir())(nSurf); // pour simplifier
|
|
if (tab_Press.Taille() != defS.Phi1_Taille())
|
|
tab_Press.Change_taille(defS.Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en implicit
|
|
const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul);
|
|
const Met_abstraite::Impli* ex_impli = &ex;
|
|
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
|
|
const Met_abstraite::Expli* ex_expli = NULL;
|
|
#ifdef MISE_AU_POINT
|
|
// test : a priori non valide pour les éléments axi et pour un pb autre que 3D
|
|
if (( ParaGlob::Dimension() != 3) || ((ParaGlob::Dimension() == 3)&&(ParaGlob::AxiSymetrie())))
|
|
{ cout << "\n erreur, ce type de chargement n'est valide pour pour la dimension 3 !! et non axi"
|
|
<< " ElemMeca::SM_charge_line_Suiv_E ( ..."
|
|
<< " \n dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
};
|
|
// test pour savoir si l'on a réellement affaire à une surface
|
|
if (ex.giB_t->NbVecteur() != 2)
|
|
{ cout << "\n erreur, pb de nombre de vecteurs pour la surface !!"
|
|
<< " ElemMeca::SM_charge_surf_Suiv_I (.... " << ex.giB_t->NbVecteur()
|
|
<< " dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// détermination de la normale à la surface
|
|
// la normale vaut le produit vectoriel des 2 premiers vecteurs
|
|
Coordonnee normale_0 = Util::ProdVec_coorBN( (*ex.giB_0)(1), (*ex.giB_0)(2));
|
|
const Coordonnee& giB_tdt1 = (*ex.giB_tdt).Coordo(1); const Coordonnee& giB_tdt2 = (*ex.giB_tdt).Coordo(2);
|
|
const Coordonnee& giB_01 = (*ex.giB_0).Coordo(1); const Coordonnee& giB_02 = (*ex.giB_0).Coordo(2);
|
|
Coordonnee normale_t = Util::ProdVec_coor( giB_tdt1, giB_tdt2);
|
|
// calcul des normales normées
|
|
normale_0.Normer();
|
|
double norme_normarle_t = normale_t.Norme();normale_t /= norme_normarle_t;
|
|
// calcul des composantes de direction de la force dans le repère local initial
|
|
Coordonnee dir_forceH_0(3);
|
|
// on ne veut pas de vérification de variance
|
|
dir_forceH_0(1) = giB_01 * forc; dir_forceH_0(2) = giB_02 * forc;
|
|
dir_forceH_0(3) = normale_0 * forc;
|
|
// calcul de la direction finale
|
|
Coordonnee dir_f_t = dir_forceH_0(1) * giB_tdt1 + dir_forceH_0(2) * giB_tdt2
|
|
+ dir_forceH_0(3) * normale_t;
|
|
Coordonnee dir_f_t_normer = dir_f_t/(dir_f_t.Norme());
|
|
// calcul de la force finale
|
|
Coordonnee force_t = module_f * dir_f_t_normer;
|
|
// debug
|
|
//cout << "\n debug Element::ResRaid ElemMeca::SMR_charge_surf_Suiv_I "
|
|
// << "\n force= "<< force_t << endl;
|
|
//fin debug
|
|
// calcul du résidu et éventuellement de la raideur
|
|
int ix=1; int dimf = forc.Dimension(); // normalement ne fonctionne qu'en 3D
|
|
// sauf pour les éléments axisymétriques pour lesquelles on ne considère que les deux premières composantes
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
// on ne tient pas compte de la variation du point dans la fct nD
|
|
if (pt_fonct != NULL)
|
|
{ // ici on utilise les variables connues aux noeuds, ou calculées à partir de
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
|
|
// pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_enu_scal,defS,ex_impli,ex_expli_tdt,ex_expli)
|
|
);
|
|
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
|
|
// pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles_interpoler_ou_calculer
|
|
(absolue,TEMPS_tdt,li_quelc,defS,ex_impli,ex_expli_tdt,ex_expli);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
// {if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee & M = defS.Position_tdt();
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_t *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_t(3) *= tab_val(3);
|
|
case 2:force_t(2) *= tab_val(2);
|
|
case 1:force_t(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_surf_Suiv_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " surfacique suiveuse: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " surfacique: dependance autre que du point "<< endl;
|
|
// Sortie(1);
|
|
// };
|
|
};
|
|
tab_Press(ni)=force_t; // sauvegarde de l'effort
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* force_t(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
Coordonnee d_giB_tdt1,d_giB_tdt2,d_dir_f_t; // def de vecteurs intermediaires
|
|
Coordonnee d_dir_f_t_normer; // " "
|
|
Coordonnee D_normale_pasnormer,D_normale_t;
|
|
if (avec_raid)
|
|
{for (int ne =1; ne<= nbne; ne++)
|
|
for (int j =1; j<= nbddl; j++)
|
|
{ // calcul de la variation de la normale
|
|
// 1) variation du produit vectoriel pour la normale finale
|
|
D_normale_pasnormer = Util::ProdVec_coor(giB_tdt1,(*ex.d_giB_tdt)(j).Coordo(2))
|
|
+ Util::ProdVec_coor((*ex.d_giB_tdt)(j).Coordo(1),giB_tdt2);
|
|
// 2) de la normale finale
|
|
D_normale_t = Util::VarUnVect_coor(normale_t,D_normale_pasnormer,norme_normarle_t);
|
|
// calcul des variations de la force
|
|
d_giB_tdt1 = ((((*ex.d_giB_tdt)(j))(1)).Coor());
|
|
d_giB_tdt2 = ((((*ex.d_giB_tdt)(j))(2)).Coor());
|
|
d_dir_f_t = dir_forceH_0(1) * d_giB_tdt1 + dir_forceH_0(2) * d_giB_tdt2
|
|
+ dir_forceH_0(3) * D_normale_t;
|
|
d_dir_f_t_normer = Util::VarUnVect_coor( dir_f_t,d_dir_f_t,module_f);
|
|
for (int i=1;i<= dimf;i++)
|
|
{ ix = i+ (ne-1)*dimf;
|
|
KM(ix,j) -= taphi(ni)(ne) * poids(ni)
|
|
* (force_t(i) *(*ex.d_jacobien_tdt)(j)
|
|
+ (*ex.jacobien_tdt) * d_dir_f_t(i) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = (*res_extS)(nSurf);
|
|
el.raid = (*raid_extS)(nSurf);
|
|
return el;
|
|
};
|
|
|
|
// cas d'un chargement volumique, sur l'élément
|
|
// force indique la force volumique appliquée
|
|
// retourne le second membre résultant
|
|
// calcul à l'instant tdt ou t en fonction de la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_vol_E (DdlElement &
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & ,bool sur_volume_finale_,bool atdt )
|
|
{
|
|
// #ifdef MISE_AU_POINT
|
|
// axisymétrie: a priori ce n'est pas nécessaire, on met une erreur
|
|
if(ParaGlob::AxiSymetrie())
|
|
{ cout << "\n erreur, les charges volumique ne sont utilisable avec les elements axisymetriques"
|
|
<< "\n utilisez à la place les charges surfaciques !! "
|
|
<< " ElemMeca::SM_charge_vol_E ( ...";
|
|
Sortie(1);
|
|
};
|
|
// #endif
|
|
|
|
// définition du vecteur de retour
|
|
Vecteur& SM = * residu;
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->Force_volume() == NULL)
|
|
lesChargeExtSurEle->Force_volume_Change_taille(def->Phi1_Taille());
|
|
Tableau <Coordonnee>& tab_F_vol= (*lesChargeExtSurEle->Force_volume()); // pour simplifier
|
|
if (tab_F_vol.Taille() != def->Phi1_Taille())
|
|
tab_F_vol.Change_taille(def->Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
int ni; // compteur globale de point d'integration
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
Met_abstraite::Expli ex;
|
|
if (atdt)
|
|
ex = (def->Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
|
|
else
|
|
ex = def->Cal_explicit_t(premier_calcul);
|
|
// choix entre volume final ou initial
|
|
double jacobien = (*ex.jacobien_t);
|
|
if (!sur_volume_finale_) jacobien = (*ex.jacobien_0);
|
|
|
|
int ix=1; int dimf = force.Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
Coordonnee force_reelle(force); // init par défaut
|
|
if (pt_fonct != NULL)
|
|
{ // ici on peut utiliser toutes les variables connues de l'élément
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,ni,-1));
|
|
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,ni,-1);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_reelle *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_reelle(3) *= tab_val(3);
|
|
case 2:force_reelle(2) *= tab_val(2);
|
|
case 1:force_reelle(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_vol_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " volumique: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee * Mpt = NULL; // init
|
|
// // on regarde s'il s'agit du volume final ou initial
|
|
// if (sur_volume_finale_)
|
|
// { Mpt = &(def->Position_0()); }
|
|
// else
|
|
// { Mpt = &(def->Position_tdt());};
|
|
// const Coordonnee & M = *Mpt;
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// // deux cas: soit une fonction scalaire soit dimf composantes
|
|
// if (tab_val.Taille() == 1)
|
|
// {force_reelle *= tab_val(1);}
|
|
// else if (tab_val.Taille() == dimf)
|
|
// {switch (dimf)
|
|
// {case 3:force_reelle(3) *= tab_val(3);
|
|
// case 2:force_reelle(2) *= tab_val(2);
|
|
// case 1:force_reelle(1) *= tab_val(1);
|
|
// break;
|
|
// default: cout << "\n *** erreur ElemMeca::SM_charge_vol_E(.. ";
|
|
// Sortie(1);
|
|
// };
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " volumique: la dimension de la fct nD "<<tab_val.Taille()
|
|
// << " devrait etre soit 1, soit "<<dimf << endl;
|
|
// Sortie(1);
|
|
// };
|
|
// };
|
|
// };
|
|
tab_F_vol(ni) +=force_reelle; // on sauvegarde la force de volume
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* force_reelle(i) * (poids(ni) * jacobien);
|
|
}
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
// cas d'un chargement volumique, sur l'élément
|
|
// force indique la force volumique appliquée
|
|
// retourne le second membre résultant -> implicite,
|
|
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
void ElemMeca::SMR_charge_vol_I(DdlElement & tab_ddl
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids,const Coordonnee& force,Fonction_nD* pt_fonct
|
|
,const ParaAlgoControle & pa,bool sur_volume_finale_)
|
|
{
|
|
// #ifdef MISE_AU_POINT
|
|
// axisymétrie: a priori ce n'est pas nécessaire, on met une erreur
|
|
if(ParaGlob::AxiSymetrie())
|
|
{ cout << "\n erreur, les charges volumique ne sont pas utilisable avec les elements axisymetriques"
|
|
<< "\n utilisez à la place les charges surfaciques !! "
|
|
<< " ElemMeca::SMR_charge_vol_I ( ...";
|
|
Sortie(1);
|
|
};
|
|
// #endif
|
|
int ni; // compteur globale de point d'integration
|
|
int nbddl = tab_ddl.NbDdl();
|
|
Vecteur& SM = *residu;
|
|
Mat_pleine & KM = *raideur;
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->Force_volume() == NULL)
|
|
lesChargeExtSurEle->Force_volume_Change_taille(def->Phi1_Taille());
|
|
Tableau <Coordonnee>& tab_F_vol= (*lesChargeExtSurEle->Force_volume()); // pour simplifier
|
|
if (tab_F_vol.Taille() != def->Phi1_Taille())
|
|
tab_F_vol.Change_taille(def->Phi1_Taille(),Coordonnee(ParaGlob::Dimension()));
|
|
|
|
// controle
|
|
bool avec_raid = pa.Var_charge_externe();
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (def->PremierPtInteg(), ni = 1;def->DernierPtInteg();def->NevezPtInteg(),ni++)
|
|
{
|
|
// calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en implicit
|
|
const Met_abstraite::Impli& ex = def->Cal_implicit(premier_calcul);
|
|
int ix=1; int dimf = force.Dimension();
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
// prise en compte d'une dépendance à une fonction nD
|
|
// on ne tient pas compte de la variation du point dans la fct nD
|
|
Coordonnee force_reelle(force); // init par défaut
|
|
if (pt_fonct != NULL)
|
|
{ // ici on peut utiliser toutes les variables connues de l'élément
|
|
// on commence par récupérer les conteneurs des grandeurs à fournir
|
|
List_io <Ddl_enum_etendu>& li_enu_scal = pt_fonct->Li_enu_etendu_scalaire();
|
|
List_io <TypeQuelconque >& li_quelc = pt_fonct->Li_equi_Quel_evolue();
|
|
bool absolue = true; // on se place systématiquement en absolu
|
|
// on va utiliser la méhode Valeur_multi pour les grandeurs strictement scalaire
|
|
Tableau <double> val_ddl_enum(ElemMeca::Valeur_multi(absolue,TEMPS_tdt,li_enu_scal,ni,-1));
|
|
// on utilise la méthode des grandeurs évoluées pour les Coordonnees et Tenseur
|
|
Valeurs_Tensorielles(absolue, TEMPS_tdt,li_quelc,ni,-1);
|
|
// calcul de la valeur et retour dans tab_ret
|
|
Tableau <double> & tab_val = pt_fonct->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
|
|
|
|
// deux cas: soit une fonction scalaire soit dimf composantes
|
|
if (tab_val.Taille() == 1)
|
|
{force_reelle *= tab_val(1);}
|
|
else if (tab_val.Taille() == dimf)
|
|
{switch (dimf)
|
|
{case 3:force_reelle(3) *= tab_val(3);
|
|
case 2:force_reelle(2) *= tab_val(2);
|
|
case 1:force_reelle(1) *= tab_val(1);
|
|
break;
|
|
default: cout << "\n *** erreur ElemMeca::SM_charge_vol_E(.. ";
|
|
Sortie(1);
|
|
};
|
|
}
|
|
else
|
|
{cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
<< " volumique: la dimension de la fct nD "<<tab_val.Taille()
|
|
<< " devrait etre soit 1, soit "<<dimf << endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// {if (pt_fonct->Depend_M())
|
|
// {// 1) on récupère la position dans I_a du point d'intégration
|
|
// const Coordonnee * Mpt = NULL; // init
|
|
// // on regarde s'il s'agit du volume final ou initial
|
|
// if (sur_volume_finale_)
|
|
// { Mpt = &(def->Position_0()); }
|
|
// else
|
|
// { Mpt = &(def->Position_tdt());};
|
|
// const Coordonnee & M = *Mpt;
|
|
// // 2) on appel la fonction, qui ne doit dépendre que de M en var locales
|
|
// Tableau <Coordonnee> tab(1); tab(1)=M;
|
|
// Tableau <double> tab_val = pt_fonct->Val_FnD_Evoluee(NULL,&tab,NULL);
|
|
// // deux cas: soit une fonction scalaire soit dimf composantes
|
|
// if (tab_val.Taille() == 1)
|
|
// {force_reelle *= tab_val(1);}
|
|
// else if (tab_val.Taille() == dimf)
|
|
// {switch (dimf)
|
|
// {case 3:force_reelle(3) *= tab_val(3);
|
|
// case 2:force_reelle(2) *= tab_val(2);
|
|
// case 1:force_reelle(1) *= tab_val(1);
|
|
// break;
|
|
// default: cout << "\n *** erreur ElemMeca::SM_charge_vol_E(.. ";
|
|
// Sortie(1);
|
|
// };
|
|
// }
|
|
// else
|
|
// {cout << "\n *** erreur dans l'utilisation d'une fct nD pour un chargement"
|
|
// << " volumique: la dimension de la fct nD "<<tab_val.Taille()
|
|
// << " devrait etre soit 1, soit "<<dimf << endl;
|
|
// Sortie(1);
|
|
// };
|
|
// };
|
|
// };
|
|
tab_F_vol(ni) +=force_reelle; // on sauvegarde la force de volume
|
|
// on choisit entre volume final ou initial
|
|
if (sur_volume_finale_)
|
|
{for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) +=
|
|
taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
if (avec_raid)
|
|
for (int j =1; j<= nbddl; j++)
|
|
KM(ix,j) -=
|
|
taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j));
|
|
};
|
|
}
|
|
else // cas plus simple du volume initial
|
|
{for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) +=
|
|
taphi(ni)(ne)* force_reelle(i) * (poids(ni) * (*ex.jacobien_0));
|
|
// pas de contribution à la raideur car le jacobien ne bouge pas
|
|
};
|
|
};
|
|
}
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
};
|
|
|
|
// cas d'un chargement hydrostatique, sur les surfaces de l'élément
|
|
// la charge dépend de la hauteur à la surface libre du liquide déterminée par un point
|
|
// et une direction normale à la surface libre:
|
|
// nSurf : le numéro de la surface externe
|
|
// poidvol: indique le poids volumique du liquide
|
|
// M_liquide : un point de la surface libre
|
|
// dir_normal_liquide : direction normale à la surface libre
|
|
// retourne le second membre résultant
|
|
// calcul à l'instant tdt ou t en fonction de la variable atdt
|
|
Vecteur& ElemMeca::SM_charge_hydro_E (DdlElement & ,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids
|
|
,const Coordonnee& dir_normal_liquide,const double& poidvol
|
|
,const Coordonnee& M_liquide,bool sans_limitation
|
|
,const ParaAlgoControle & ,bool atdt)
|
|
{
|
|
// #ifdef MISE_AU_POINT
|
|
// // axisymétrie: pour l'instant non testé a priori avec les forces hydro
|
|
// if(ParaGlob::AxiSymetrie())
|
|
// { cout << "\n erreur, les charges hydro ne sont pas testees avec les elements axisymetriques"
|
|
// << " ElemMeca::SM_charge_hydro_E ( ...";
|
|
// Sortie(1);
|
|
// };
|
|
// #endif
|
|
|
|
// *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
|
|
// *** sinon c'est intractable et long
|
|
// définition du vecteur de retour
|
|
Deformation* definter=NULL; // variable de travail transitoire
|
|
Vecteur* SM_P = NULL;
|
|
if (!(ParaGlob::AxiSymetrie()))
|
|
{ definter = defSurf(nSurf); // le cas normal est non axisymétrique
|
|
SM_P = ((*res_extS)(nSurf));
|
|
}
|
|
else { definter = defArete(nSurf); // en axisymétrie, c'est une def d'arête
|
|
SM_P = ((*res_extA)(nSurf));
|
|
};
|
|
Deformation & defS = *definter; // pour simplifier l'écriture
|
|
Vecteur& SM = *SM_P; // " " "
|
|
|
|
// dans le cas ou le poid volumique est nul, on retourne sans calcul
|
|
if (Dabs(poidvol) < ConstMath::trespetit)
|
|
return SM;
|
|
int ni; // compteur globale de point d'integration
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesPressionsExternes() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Pression_appliquee>& tab_Press= (*lesChargeExtSurEle->LesPressionsExternes())(nSurf); // pour simplifier
|
|
if (tab_Press.Taille() != defS.Phi1_Taille())
|
|
tab_Press.Change_taille(defS.Phi1_Taille());
|
|
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul de la pression au niveau de la position du point d'intégration
|
|
// 1) on récupère la position dans I_a du point d'intégration
|
|
const Coordonnee & M = defS.Position_tdt();
|
|
// 2) calcul de la distance à la surface libre
|
|
Coordonnee MMp= M_liquide - M; double dis = dir_normal_liquide * MMp;
|
|
// 3) calcul de la pression effective uniquement si le point est dans l'eau
|
|
// en fait si le point est hors de l'eau il n'y a pas de pression hydrostatique
|
|
// sauf si sans_limitation est true
|
|
if ((dis > 0.) || sans_limitation)
|
|
{ double pression = -dis * poidvol;
|
|
// sauvegarde de la pression appliquée
|
|
tab_Press(ni).press += pression;
|
|
// -- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en explicit mais pour un calcul autre que mécanique
|
|
Met_abstraite::Expli ex;
|
|
if (atdt)
|
|
ex = (defS.Cal_explicit_tdt(premier_calcul)).Tdt_dans_t();
|
|
else
|
|
ex = defS.Cal_explicit_t(premier_calcul);
|
|
// détermination de la normale à la surface
|
|
#ifdef MISE_AU_POINT
|
|
// test pour savoir si l'on a réellement affaire à une surface
|
|
if ((ex.giB_t->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3))
|
|
{ cout << "\n erreur, il doit y avoir deux vecteurs pour la surface de dimension 3"
|
|
<< " ElemMeca::SM_charge_hydro_E (DdlElement & ...nbvec= " << ex.giB_t->NbVecteur()
|
|
<< " dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// la normale vaut le produit vectoriel des 2 premiers vecteurs
|
|
CoordonneeB normale = Util::ProdVec_coorB( (*ex.giB_t)(1), (*ex.giB_t)(2));
|
|
// calcul de la normale normée et pondérée de la pression
|
|
normale.Normer();normale *= -pression;
|
|
tab_Press(ni).P += normale.Coor_const();
|
|
int ix=1; int dimf = 3; // ne fonctionne qu'en dim 3
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
SM(ix) += taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_t));
|
|
}// fin test : point dans le liquide ou pas
|
|
} // fin boucle sur les points d'intégrations
|
|
// retour du second membre
|
|
return SM;
|
|
};
|
|
|
|
|
|
// idem SMR_charge_hydro_E mais
|
|
// -> implicite,
|
|
// pa: permet de déterminer si oui ou non on calcul la contribution à la raideur
|
|
// retourne le second membre et la matrice de raideur correspondant
|
|
Element::ResRaid ElemMeca::SMR_charge_hydro_I (DdlElement & ddls,int nSurf
|
|
,const Tableau <Vecteur>& taphi,int nbne
|
|
,const Vecteur& poids
|
|
,const Coordonnee& dir_normal_liquide,const double& poidvol
|
|
,const Coordonnee& M_liquide,bool sans_limitation
|
|
,const ParaAlgoControle & pa)
|
|
{
|
|
// #ifdef MISE_AU_POINT
|
|
// // axisymétrie: pour l'instant non testé a priori avec les forces hydro
|
|
// if(ParaGlob::AxiSymetrie())
|
|
// { cout << "\n erreur, les charges hydro ne sont pas testees avec les elements axisymetriques"
|
|
// << " ElemMeca::SM_charge_hydro_I ( ...";
|
|
// Sortie(1);
|
|
// };
|
|
// #endif
|
|
Deformation* definter=NULL; // variable de travail transitoire
|
|
Vecteur* SM_P = NULL; Mat_pleine* KM_P = NULL;
|
|
if (!(ParaGlob::AxiSymetrie()))
|
|
{ definter = defSurf(nSurf); // le cas normal est non axisymétrique
|
|
SM_P = ((*res_extS)(nSurf));
|
|
KM_P = ((*raid_extS)(nSurf));
|
|
}
|
|
else { definter = defArete(nSurf); // en axisymétrie, c'est une def d'arête
|
|
SM_P = ((*res_extA)(nSurf));
|
|
KM_P = ((*raid_extA)(nSurf));
|
|
};
|
|
Deformation & defS = *definter; // pour simplifier l'écriture
|
|
Vecteur& SM = *SM_P; // " " "
|
|
Mat_pleine & KM = *KM_P; // " " "
|
|
|
|
// *** dans l'ensemble du sp on utilise des vecteurs et on gére directement la variance
|
|
// *** sinon c'est intractable et long
|
|
// dans le cas ou le poid volumique est nul, on retourne sans calcul
|
|
if (Dabs(poidvol) < ConstMath::trespetit)
|
|
{ Element::ResRaid el;
|
|
el.res = SM_P;
|
|
el.raid = KM_P;
|
|
return el;
|
|
}
|
|
// controle
|
|
bool avec_raid = pa.Var_charge_externe();
|
|
// // dans le cas où la variation sur la raideur est demandé on s'assure que la variation
|
|
// // du jacobien est bien effective sinon on l'impose
|
|
// ParaAlgoControle pa_nevez(pa);
|
|
// if (!pa_nevez.Var_jacobien()) pa_nevez.Modif_Var_jacobien(true);
|
|
int ni; // compteur globale de point d'integration
|
|
|
|
// éventuellement dimensionnement de la sauvegarde
|
|
if (lesChargeExtSurEle == NULL)
|
|
lesChargeExtSurEle = new LesChargeExtSurElement;
|
|
if (lesChargeExtSurEle->LesPressionsExternes() == NULL)
|
|
{if(ParaGlob::AxiSymetrie())
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbSe());}
|
|
else
|
|
{lesChargeExtSurEle->LesPressionsExternes_Change_taille(ElementGeometrique().NbFe());}
|
|
};
|
|
Tableau <Pression_appliquee>& tab_Press= (*lesChargeExtSurEle->LesPressionsExternes())(nSurf); // pour simplifier
|
|
if (tab_Press.Taille() != defS.Phi1_Taille())
|
|
tab_Press.Change_taille(defS.Phi1_Taille());
|
|
|
|
int nbddl = ddls.NbDdl();
|
|
bool premier_calcul=true; // contrairement à la déformation, pas de sauvegarde
|
|
// donc il faut calculer tous les éléments de la métrique
|
|
for (defS.PremierPtInteg(), ni = 1;defS.DernierPtInteg();defS.NevezPtInteg(),ni++)
|
|
{
|
|
// calcul de la pression au niveau de la position du point d'intégration
|
|
// 1) on récupère la position dans I_a du point d'intégration
|
|
const Coordonnee & M = defS.Position_tdt();
|
|
// 2) calcul de la distance à la surface libre
|
|
Coordonnee MMp= M_liquide - M; double dis = dir_normal_liquide * MMp;
|
|
// 3) calcul de la pression effective uniquement si le point est dans l'eau
|
|
// en fait si le point est hors de l'eau il n'y a pas de pression hydrostatique
|
|
// sauf si sans_limitation est true
|
|
if ((dis > 0.) || sans_limitation)
|
|
{ double pression = -dis * poidvol;
|
|
// sauvegarde de la pression appliquée
|
|
tab_Press(ni).press += pression;
|
|
// dans le cas avec raideur on calcul la variation de la pression
|
|
const Tableau <Coordonnee> * d_M_pointe = NULL;
|
|
if (avec_raid)
|
|
d_M_pointe=&(defS.Der_Posi_tdt());
|
|
// --- calcul des éléments de la métrique, entre autre le jacobien, on utilise le même calcul
|
|
// que pour un calcul primaire en implicit
|
|
const Met_abstraite::Impli& ex = defS.Cal_implicit(premier_calcul);
|
|
#ifdef MISE_AU_POINT
|
|
// test pour savoir si l'on a réellement affaire à une surface
|
|
if ((ex.giB_t->NbVecteur() != 2) || ( ParaGlob::Dimension() != 3))
|
|
{ cout << "\n erreur, il doit y avoir deux vecteurs pour la surface de dimension 3"
|
|
<< " ElemMeca::SMR_charge_hydro_I (.... " << ex.giB_t->NbVecteur()
|
|
<< " dim = " << ParaGlob::Dimension();
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// la normale vaut le produit vectoriel des 2 premiers vecteurs
|
|
CoordonneeB normale = Util::ProdVec_coorB( (*ex.giB_tdt)(1), (*ex.giB_tdt)(2));
|
|
// calcul de la variation de la normale
|
|
// 1) variation du produit vectoriel
|
|
Tableau <CoordonneeB > D_pasnormale =
|
|
Util::VarProdVect_coorB( (*ex.giB_tdt)(1),(*ex.giB_tdt)(2),(*ex.d_giB_tdt));
|
|
// 2) de la normale
|
|
Tableau <CoordonneeB> D_normale = Util::VarUnVect_coorB(normale,D_pasnormale,normale.Norme());
|
|
// calcul de la normale normée et pondérée de la pression
|
|
normale.Normer();normale *= -pression;
|
|
tab_Press(ni).P += normale.Coor_const();
|
|
// également pondération de la variation de la
|
|
if (avec_raid)
|
|
{ for (int ihi =1;ihi<= nbddl;ihi++)
|
|
{ double d_pression = poidvol * (dir_normal_liquide * (*d_M_pointe)(ihi));
|
|
// le - devant dis disparait avec le - devant M
|
|
D_normale(ihi) = pression * D_normale(ihi) + d_pression * normale ;
|
|
}
|
|
}
|
|
int ix=1; int dimf = 3; // ne fonctionne qu'en dim 3
|
|
if(ParaGlob::AxiSymetrie())
|
|
// cas d'un chargement axisymétrique, dans ce cas on ne prend en compte que les
|
|
// dimf-1 coordonnées donc on décrémente
|
|
dimf--;
|
|
for (int ne =1; ne<= nbne; ne++)
|
|
for (int i=1;i<= dimf;i++,ix++)
|
|
{ SM(ix) += taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.jacobien_tdt));
|
|
// dans le cas avec_raideur on calcul la contribution à la raideur
|
|
if (avec_raid)
|
|
for (int j =1; j<= nbddl; j++)
|
|
KM(ix,j) -=
|
|
taphi(ni)(ne)* normale(i) * (poids(ni) * (*ex.d_jacobien_tdt)(j))
|
|
+ taphi(ni)(ne)* D_normale(j)(i) * (poids(ni) * (*ex.jacobien_tdt)) ;
|
|
}
|
|
}// fin test : point dans le liquide ou pas
|
|
} // fin boucle sur les points d'intégrations
|
|
|
|
// liberation des tenseurs intermediaires
|
|
LibereTenseur();
|
|
Element::ResRaid el;
|
|
el.res = SM_P;
|
|
el.raid = KM_P;
|
|
return el;
|
|
};
|
|
|
|
|