Herezh_dev/Elements/Mecanique/ElemMeca.cc
Gérard Rio ece68c42db V 7.012 : résolution du pb de l'utilisation de l'accélération dans une fonction nD
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
2023-06-01 08:47:54 +02:00

3582 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;
};