rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/Elements/Mecanique/SFE/SfeMembT3.cc
Gérard Rio 9692dbd130 intégration du répertoire Mecanique: 
- contient les éléments finis, métriques associées, déformations ...
intégration du réperoire Géométrie:
- contient les géométries 1D 2D et 3D, les frontières des éléments géométriques
2021-09-27 12:42:13 +02:00

803 lines
40 KiB
C++
Raw Blame History

// FICHIER : SfeMembT.cc
// CLASSE : SfeMembT
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2021 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
# include <iostream>
using namespace std; //introduces namespace std
#include <stdlib.h>
#include "Sortie.h"
#include "MathUtil.h"
#include "SfeMembT.h"
#include "TypeConsTens.h"
#include "FrontPointF.h"
#include "Loi_Umat.h"
#include "ExceptionsElemMeca.h"
#include "TypeQuelconqueParticulier.h"
#include "Coordonnee2.h"
#include "Coordonnee3.h"
#include "Util.h"
// calcul de la nouvelle épaisseur (sans raideur) avec métrique en explicite
// cas où l'épaisseur est stocké dans l'élément
// mise à jour du volume au pti
void SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(const double& epaisseur0, const ParaAlgoControle &
,const double & epaisseur_t, PtIntegMecaInterne & ptIntegMeca
,double & epaisseur_tdt, const Met_abstraite::Expli_t_tdt& ex)
{ // on commence par calculer la trace du tenseur des contraintes
double traceSig = (ptIntegMeca.SigHH_const() * (* ex.gijBB_tdt)).Trace();
double traceSig_t = (ptIntegMeca.SigHH_t_const() * (* ex.gijBB_t)).Trace();
const double& troisK = 3. * ptIntegMeca.ModuleCompressibilite_const(); // récup de la compressibilité
// calcul de la nouvelle épaisseur
//--- vérification et gestion des cas particuliers ---
bool cas_particulier=false;
{// on suppose un module de compressibilité positif sinon ce n'est pas normal
if (troisK < 0.)
{if (ParaGlob::NiveauImpression()>2)
cout << "\n erreur*** on a trouve une compressibilite negative = " << troisK/3.
<< "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(...";
// on ne change pas l'épaisseur
cas_particulier = true;
};
// classiquement on suppose que troisK doit-être bien supérieur à traceSig:
// on considère un facteur 10 arbitraire
// si ce n'est pas le cas, on risque d'avoir un calcul d'épaisseur erroné
if (Dabs(traceSig) > 10. * Dabs(troisK/3.))
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " le coef de compressibilite = "<<troisK/3
<< " est par rapport a delta la trace de sigma = "
<< traceSig << " 10 fois inferieur !! ce n'est pas normal a priori "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(..."
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
double delta_traceSig = traceSig-traceSig_t;
// maintenant on regarde si le module peut conduire à une épaisseur négative
if (delta_traceSig > (troisK/3. - ConstMath::petit))
// dans le cas où delta traceSig > troisK, cela va donner une épaisseur négative !!
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " coef de compressibilite = "<<troisK/3
<< " est inferieur a delta trace de sigma = "
<< delta_traceSig
<< " ce qui va donner une epaisseur finale negative !! "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(..."
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
// maintenant on regarde si le module n'est pas trop petit
// si oui il y a aussi un pb
if ((Dabs(troisK) <= ConstMath::unpeupetit)
&& (Dabs(delta_traceSig) > 10. * Dabs(troisK/3.))
)
// on ne va pas pouvoir peut-être calculer une nouvelle épaisseur
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " coef de compressibilite = "<<troisK/3
<< " et delta trace de sigma = "
<< delta_traceSig
<< " sont de valeurs trop faibles (et un peu etrange) !! "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurAtdtExp_et_vol_pti(..."
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
};
if (cas_particulier)
{ epaisseur_tdt = epaisseur_t;
ptIntegMeca.Volume_pti() *= epaisseur_t;
}
else // sinon cas sans pb
{
// epaisseur_tdt = (epaisseur0 * (*(ex.jacobien_0))) / (*(ex.jacobien_tdt))
// * troisK / (troisK - traceSig);
epaisseur_tdt = (epaisseur_t * (*(ex.jacobien_t))) / (*(ex.jacobien_tdt))
* troisK / (troisK - traceSig+traceSig_t);
// calcul du volume au pti
ptIntegMeca.Volume_pti() *= epaisseur_tdt ;
};
};
// calcul de la nouvelle épaisseur (sans raideur) avec métrique en implicite
// cas où l'épaisseur est stocké dans l'élément
// mise à jour du volume au pti
void SfeMembT::CalEpaisseurAtdtImp_et_vol_pti(const double& epaisseur0, const ParaAlgoControle &
,const double & epaisseur_t, PtIntegMecaInterne & ptIntegMeca
,double & epaisseur_tdt, const Met_abstraite::Impli& ex)
{ // on commence par calculer la trace du tenseur des contraintes
// double traceSig = (ptIntegMeca.SigHH_const() * (* ex.gijBB_tdt)).Trace();
// double traceSig_t = (ptIntegMeca.SigHH_t_const() * (* ex.gijBB_t)).Trace();
double traceSig = (ptIntegMeca.SigHH_const() && (* ex.gijBB_tdt));
double traceSig_t = (ptIntegMeca.SigHH_t_const() && (* ex.gijBB_t));
const double& troisK = 3. * ptIntegMeca.ModuleCompressibilite_const(); // récup de la compressibilité
// calcul de la nouvelle épaisseur
//--- vérification et gestion des cas particuliers ---
bool cas_particulier=false;
{// on suppose un module de compressibilité positif sinon ce n'est pas normal
if (troisK < 0.)
{if (ParaGlob::NiveauImpression()>2)
cout << "\n erreur*** on a trouve une compressibilite negative = " << troisK
<< "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) ";
// on ne change pas l'épaisseur
cas_particulier = true;
};
// classiquement on suppose que troisK doit-être bien supérieur à traceSig:
// on considère un facteur 10 arbitraire
// si ce n'est pas le cas, on risque d'avoir un calcul d'épaisseur erroné
if (Dabs(traceSig) > 10. * Dabs(troisK/3.))
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " le coef de compressibilite = "<<troisK/3
<< " est par rapport a delta la trace de sigma = "
<< traceSig << " 10 fois inferieur !! ce n'est pas normal a priori "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) "
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
double delta_traceSig = traceSig-traceSig_t;
// maintenant on regarde si le module peut conduire à une épaisseur négative
if (delta_traceSig > (troisK/3. - ConstMath::petit))
// dans le cas où delta traceSig > troisK, cela va donner une épaisseur négative !!
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " coef de compressibilite = "<<troisK/3
<< " est inferieur a delta trace de sigma = "
<< delta_traceSig
<< " ce qui va donner une epaisseur finale negative !! "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) "
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
// maintenant on regarde si le module n'est pas trop petit
// si oui il y a aussi un pb
if ((Dabs(troisK) <= ConstMath::unpeupetit)
&& (Dabs(delta_traceSig) > 10. * Dabs(troisK/3.))
)
// on ne va pas pouvoir peut-être calculer une nouvelle épaisseur
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " coef de compressibilite = "<<troisK/3
<< " et delta trace de sigma = "
<< delta_traceSig
<< " sont de valeurs trop faibles (et un peu etrange) !! "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurAtdtImp_et_vol_pti...) "
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
};
if (cas_particulier)
{ epaisseur_tdt = epaisseur_t;
ptIntegMeca.Volume_pti() *= epaisseur_t;
}
else // sinon cas sans pb
{
epaisseur_tdt = (epaisseur_t * (*(ex.jacobien_t))) / (*(ex.jacobien_tdt))
* troisK / (troisK - traceSig+traceSig_t);
// calcul du volume au pti
ptIntegMeca.Volume_pti() *= epaisseur_tdt ;
};
};
// calcul de la nouvelle épaisseur moyenne finale (sans raideur)
// mise à jour des volumes aux pti
// ramène l'épaisseur moyenne calculée à atdt
const double& SfeMembT::CalEpaisseurMoyenne_et_vol_pti(bool atdt)
{ // .. en bouclant sur les pt d'integ enregistré ..
// -- on récupère et on calcule les jacobiens moyens à 0 et final
double jacobien_moy_fin = 0.; // init
double jacobien_moy_ini = 0.;
// -- de même on récupère et on calcul la trace moyenne de la contrainte
double traceSig_moy = 0.;double traceSig_moy_ini = 0.;
// -- de même on récupère et on calcul le module de compressibilité moyen
double troisK_moy = 0.;
if (donnee_specif.epais != NULL)
{Epai& epais = *(donnee_specif.epais); // pour simplifier
int indice = 1;
SfeMembT::DonnComSfe* CoSfe = unefois->doCoMemb; // pour simplifier l'écriture
int nbtotsurf =CoSfe->eleCentre->Nbi();
int nbtotepais =(CoSfe->segment).Nbi();
int nbi = nbtotsurf * nbtotepais;
for (int nisur =1; nisur<= nbtotsurf; nisur++) // boucle sur les pt d'integ de surface
for (int niepais =1; niepais<= nbtotepais; niepais++, indice++) // pt d'integ d'epaisseur
// for (int i=1;i<= nombre->nbi;i++)
{ // cas de la compressibilité
const double& troisK = 3. * (*lesPtIntegMecaInterne)(indice).ModuleCompressibilite_const();
troisK_moy += troisK;
// cas des jacobiens
const double& jacobien_0 = *(tabSaveDefDon(indice)->Meti_00().jacobien_);
if (atdt)
{ const double& jacobien_ini = *(tabSaveDefDon(indice)->Meti_t().jacobien_);
jacobien_moy_ini += jacobien_ini;
const double& jacobien_fin = *(tabSaveDefDon(indice)->Meti_tdt().jacobien_);
jacobien_moy_fin += jacobien_fin;
// cas de la trace de sigma
const double traceSig = (*lesPtIntegMecaInterne)(indice).SigHH_const() && (*(tabSaveDefDon(indice)->Meti_tdt().gijBB_));
traceSig_moy += traceSig;
const double traceSig_ini = (*lesPtIntegMecaInterne)(indice).SigHH_t_const() && (*(tabSaveDefDon(indice)->Meti_t().gijBB_));
traceSig_moy_ini += traceSig_ini;
(*lesPtIntegMecaInterne)(indice).Volume_pti() *= (epais.epaisseur_t * jacobien_ini) / jacobien_fin
* troisK / (troisK - traceSig+traceSig_ini); // la pression = - traceSig/3 !!!
}
else
{ const double& jacobien_ini = *(tabSaveDefDon(indice)->Meti_00().jacobien_);
jacobien_moy_ini += jacobien_ini;
const double& jacobien_fin = *(tabSaveDefDon(indice)->Meti_t().jacobien_);
jacobien_moy_fin += jacobien_fin;
// cas de la trace de sigma
double traceSig = (*lesPtIntegMecaInterne)(indice).SigHH_const() && (*(tabSaveDefDon(indice)->Meti_t().gijBB_));
traceSig_moy += traceSig;
(*lesPtIntegMecaInterne)(indice).Volume_pti() *= (epais.epaisseur0 * jacobien_ini) / jacobien_fin
* troisK / (troisK - traceSig);
};
};
jacobien_moy_ini /= nbi;
jacobien_moy_fin /= nbi;
traceSig_moy_ini /= nbi;
traceSig_moy /= nbi;
troisK_moy /= nbi;
// d'où le calcul de la nouvelle épaisseur en utilisant la relation:
// (V-V_0)/V = trace(sigma)/3 /K_moy
//--- vérification et gestion des cas particuliers ---
bool cas_particulier=false;
{// on suppose un module de compressibilité positif sinon ce n'est pas normal
if (troisK_moy < 0.)
{if (ParaGlob::NiveauImpression()>2)
cout << "\n erreur*** on a trouve une compressibilite moyenne negative = "
<< troisK_moy/3.
<< "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(...";
// on ne change pas l'épaisseur
cas_particulier = true;
};
// classiquement on suppose que troisK_moy doit-être bien supérieur à traceSig_moy:
// on considère un facteur 10 arbitraire
// si ce n'est pas le cas, on risque d'avoir un calcul d'épaisseur erroné
double delta_traceSig = traceSig_moy-traceSig_moy_ini;
if (Dabs(delta_traceSig) > 10. * Dabs(troisK_moy))
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " le coef de compressibilite moyen = "<<troisK_moy/3.
<< " est par rapport a delta la trace de sigma = "
<< delta_traceSig << " 10 fois inferieur !! ce n'est pas normal a priori "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(... "
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
// maintenant on regarde si le module peut conduire à une épaisseur négative
if (delta_traceSig > (troisK_moy/3. - ConstMath::petit))
// dans le cas où delta_traceSig > troisK_moy, cela va donner une épaisseur négative !!
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " coef de compressibilite moyen = "<<troisK_moy/3.
<< " est inferieur a delta la trace de sigma = "
<< delta_traceSig
<< " ce qui va donner une epaisseur finale negative !! "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(... "
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
// maintenant on regarde si le module n'est pas trop petit
// si oui il y a aussi un pb
if ((Dabs(troisK_moy) <= ConstMath::unpeupetit)
&& (Dabs(traceSig_moy) > 10. * Dabs(troisK_moy/3.))
)
// on ne va pas pouvoir peut-être calculer une nouvelle épaisseur
{if (ParaGlob::NiveauImpression()>2)
cout << "\n *** pb dans le calcul de la nouvelle epaisseur: "
<< " coef de compressibilite moyen = "<<troisK_moy/3
<< " et delta trace de sigma = "
<< delta_traceSig
<< " sont de valeurs trop faibles (et un peu etrange) !! "
<< " on continue sans changer l'epaisseur ... "
<< "\n SfeMembT::CalEpaisseurMoyenne_et_vol_pti(... "
<< endl;
// on ne change pas l'épaisseur
cas_particulier = true;
};
};
// si cas particulier: on ne peut pas calculer, on reste sur les valeurs précédentes
if (cas_particulier)
{if (atdt)
{ // on met à jour l'épaisseur
epais.epaisseur_tdt = epais.epaisseur_t;
// puis retour
return epais.epaisseur_tdt;
}
else
{// retour
return epais.epaisseur_t;
};
}
else // sinon c'est ok
{if (atdt)
{epais.epaisseur_tdt = (epais.epaisseur_t * jacobien_moy_ini) / jacobien_moy_fin
* troisK_moy / (troisK_moy - traceSig_moy+traceSig_moy_ini);
//epais.epaisseur_tdt = (epais.epaisseur0 * jacobien_moy_0) / jacobien_moy_fin
// * troisK_moy / (troisK_moy + traceSig_moy);
//--debug
//if (num_elt==1)
// { Noeud* noe = tab_noeud(1);
// double R_0 = noe->Coord0().Norme(); double R = noe->Coord2().Norme();
//cout << "\n e0= " << epais.epaisseur_tdt<< " troisK_moy=" << troisK_moy << " traceSig_moy=" << traceSig_moy
// << " J0= " << jacobien_moy_0 << " J= " << jacobien_moy_fin << " R_0 " << R_0 << " R= " << R;
// };
//-- fin debug
return epais.epaisseur_tdt;
}
else
{epais.epaisseur_t = (epais.epaisseur0 * jacobien_moy_ini) / jacobien_moy_fin
* troisK_moy / (troisK_moy - traceSig_moy);
return epais.epaisseur_t;
};
};
};
// sinon c'est une erreur
cout << "\n *** erreur dans SfeMembT::CalEpaisseurMoyenne_et_vol_pti( ";
Sortie(1);
};
//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
bool SfeMembT::Delta_Jacobien_anormal(const Deformation::SaveDefResul* don
,int nisur,int niepais)
{// on récupère le jacobien de la surface médiane sans courbure
const DeformationSfe1::SaveDefResulSfe1* donsfe = ((DeformationSfe1::SaveDefResulSfe1*) don);
#ifdef MISE_AU_POINT
if (donsfe->Meti_a_tdt().jacobien_ == NULL)
{cout << "\n*** erreur, le jacobien de la facette n'est pas attribue ! "
<< " pti de surface "<<nisur
<< "\n element SFE: num= " << Num_elt_const()
<< ", num maillage= " << Num_maillage()
<< "\n SfeMembT::Delta_Jacobien_anormal(..."<<flush;
Sortie(1);
};
if (donsfe->Meti_tdt().giB_ == NULL)
{cout << "\n*** erreur, la bas giB_tdt n'est pas attribue ! "
<< " pti de surface "<<nisur << " et au pti d'epaisseur " << niepais
<< "\n element SFE: num= " << Num_elt_const()
<< ", num maillage= " << Num_maillage()
<< "\n SfeMembT::Delta_Jacobien_anormal(..."<<flush;
Sortie(1);
};
#endif
double jaco_init = *(donsfe->Meti_a_tdt().jacobien_);
// on récupère la base au pti courant
const BaseB * giB_ = donsfe->Meti_tdt().giB_;
// on calcule le produit mixte pour les gi finaux
double jacobien_final = Util::ProduitMixte(*giB_);
double ratio = jacobien_final/(jaco_init+ConstMath::pasmalpetit);
double limite = ParaGlob::param->ParaAlgoControleActifs().Ratio_maxi_jacoMembrane_jacoPti();
if ((ratio > limite) || (ratio < 1./(limite+ConstMath::pasmalpetit)))
cout << "\n *** attention *** le rapport entre le jacobien initial ("<< jaco_init
<< ") et le jacobien ("<< jacobien_final
<< ") calcule au pti de surface "<<nisur << " et au pti d'epaisseur " << niepais
<< " est en dehors de l'intervalle ["<<1./(limite+ConstMath::pasmalpetit) <<","<<limite<<"] "
<< "( peut-etre que les courbure-epaisseur sont trop importantes localement ? )"
<< "\n element SFE: num= " << Num_elt_const()
<< ", num maillage= " << Num_maillage()
<< flush;
};
// retourne la liste abond้e de tous les donn้es particuli่res interne actuellement utilis้s
// par l'้l้ment (actif ou non), sont exclu de cette liste les donn้es particuli่res des noeuds
// reli้s เ l'้l้ment
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
List_io <TypeQuelconque> SfeMembT::Les_types_particuliers_internes(bool absolue) const
{ // on commence par récupérer la liste général provenant d'ElemMeca
List_io <TypeQuelconque> ret = ElemMeca::Les_types_particuliers_internes(absolue);
// ensuite on va ajouter les données particulières aux sfe
Grandeur_scalaire_double grand_courant; // def d'une grandeur courante
if (donnee_specif.epais != NULL)
{ // $$$ cas de l'้paisseur initiale
TypeQuelconque typQ1(EPAISSEUR_MOY_INITIALE,SIG11,grand_courant);
ret.push_back(typQ1);
// $$$ cas de l'้paisseur finale
TypeQuelconque typQ2(EPAISSEUR_MOY_FINALE,SIG11,grand_courant);
ret.push_back(typQ2);
};
// $$$ cas des efforts normaux classiques, linéique
{TypeQuelconque typQ1(NN_11,SIG11,grand_courant); ret.push_back(typQ1);};
{TypeQuelconque typQ1(NN_22,SIG11,grand_courant); ret.push_back(typQ1);};
{TypeQuelconque typQ1(NN_12,SIG11,grand_courant); ret.push_back(typQ1);};
// $$$ cas des moments linéiques
{TypeQuelconque typQ1(MM_11,SIG11,grand_courant); ret.push_back(typQ1);};
{TypeQuelconque typQ1(MM_22,SIG11,grand_courant); ret.push_back(typQ1);};
{TypeQuelconque typQ1(MM_12,SIG11,grand_courant); ret.push_back(typQ1);};
// $$$ cas du tenseur de courbure
// on ne l'exprime jamais (pour l'instant) dans le repère absolu)
// mais on sort toutes les infos: repère local et courbures principales
Tenseur2BB tens; // un tenseur typique
Grandeur_TenseurBB grand2(tens); // def d'une grandeur courante TenseurBB
TypeQuelconque typQ3(TENSEUR_COURBURE,SIG11,grand2);
ret.push_back(typQ3);
// $$$ cas du vecteur contenant les courbures principales
Tableau <double> t_d(2); // un tableau qui sert d'indicateur de type
Tab_Grandeur_scalaire_double grand3(t_d); // def d'une grandeur courante tab de scalaire
TypeQuelconque typQ4(COURBURES_PRINCIPALES,SIG11,grand3);
ret.push_back(typQ4);
// $$$ cas des directions principales de courbure (dans le rep่re local
Coordonnee3 v_c_princ; // un vecteur qui sert d'indicateur de type
Tab_Grandeur_Coordonnee grand4(v_c_princ,2); // def d'une grandeur courante coordonn้e
TypeQuelconque typQ5(DIRECTIONS_PRINC_COURBURE,SIG11,grand4);
ret.push_back(typQ5);
// $$$ cas du repère local orthonormé d'expression
Coordonnee3 v_rep;
Tab_Grandeur_Coordonnee grand5(v_rep,3); // def d'une grandeur courante coordonn้e
TypeQuelconque typQ6(REPERE_LOCAL_ORTHO,SIG11,grand5);
ret.push_back(typQ6);
// liberation des tenseurs intermediaires
LibereTenseur();
return ret;
};
// récupération de grandeurs particulières au numéro d'ordre = iteg
// celles-ci peuvent être quelconques
// en retour liTQ est modifié et contiend les infos sur les grandeurs particulières
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
void SfeMembT::Grandeur_particuliere (bool absolue,List_io<TypeQuelconque>& liTQ,int iteg)
{ // *** pour l'instant on ne sort rien de particulier dans le repère absolue !! cf. les grandeurs
// on balaie la liste transmise pour les grandeurs propres
List_io<TypeQuelconque>::iterator il,ilfin = liTQ.end();
// on fait un premier passage pour voir s'il y a une demande d'infos relatives à la courbure
bool besoin_tenseur_courbure=false;TenseurBB * curbBB_tdt = NULL;
bool besoin_courbures_principales=false;Coordonnee2 val_princ;
bool besoin_dir_princ_courb=false; Tableau <Coordonnee3 > t_vec_princ;
bool besoin_rep_local_ortho=false; Tableau <Coordonnee3 > t_ortho_local;
for (il=liTQ.begin();il!=ilfin;il++)
{TypeQuelconque& tipParticu = (*il); // pour simplifier
if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
switch (tipParticu.EnuTypeQuelconque().EnumTQ())
{ case EPAISSEUR_MOY_INITIALE: (*il).Inactive(); break; // on inactive la grandeur quelconque
case EPAISSEUR_MOY_FINALE: (*il).Inactive(); break; // on inactive la grandeur quelconque
// --- cas des éléments liés à la courbure on inactive et on crée des booléen adoc
case TENSEUR_COURBURE:
{ (*il).Inactive();besoin_tenseur_courbure=true;break;}
case COURBURES_PRINCIPALES:
{ (*il).Inactive();besoin_courbures_principales=true;break;}
case DIRECTIONS_PRINC_COURBURE:
{ (*il).Inactive();besoin_dir_princ_courb=true;break;}
case REPERE_LOCAL_ORTHO:
{ (*il).Inactive();besoin_rep_local_ortho=true;break;}
default: break; // on ne fait rien
};
};
// dans le cas où on demande des infos relatives à la courbure on les récupères
if ((besoin_tenseur_courbure)||(besoin_courbures_principales)||(besoin_dir_princ_courb)
||(besoin_rep_local_ortho))
{// ----- def de grandeurs de travail
// def de la dimension des tenseurs
int dim = lesPtIntegMecaInterne->DimTens();
PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(iteg);
// on recupère le tableau pour la lecture des coordonnées des tenseurs
int nbcompo = ParaGlob::NbCompTens();
const Tableau2 <int> & OrdreCont = OrdreContrainteR(nbcompo);
// définition des grandeurs qui sont indépendante de la boucle sur les ddl_enum_etendue
def->ChangeNumInteg(iteg); // on change le num้ro de point d'int้gration courant
// def du cas de sortie
int cas=0;
if ((unefois->dualSortSfe == 0) && (unefois->CalimpPrem != 0))
{ cas=1;unefois->dualSortSfe = 1;
}
else if ((unefois->dualSortSfe == 1) && (unefois->CalimpPrem != 0))
{ cas = 11;}
else if ((unefois->dualSortSfe == 0) && (unefois->CalResPrem_tdt != 0))
{ cas=2;unefois->dualSortSfe = 1;
}
else if ((unefois->dualSortSfe == 1) && (unefois->CalResPrem_tdt != 0))
{ cas = 12;}
// sinon problème
else
{ cout << "\n warning: les grandeurs ne sont pas calculees : il faudrait au moins un pas de calcul"
<< " pour inialiser les conteneurs des tenseurs resultats ";
if (ParaGlob::NiveauImpression() >= 4)
cout << "\n cas non pr้vu, unefois->dualSortSfe= " << unefois->dualSortSfe
<< " unefois->CalimpPrem= " << unefois->CalimpPrem
<< "\n SfeMembT::Grandeur_particuliere(List_io<TypeQuelconque>&...";
Sortie(1);
};
// choix d'init en fonction de la valeur de "cas"
if (((cas == 1) || (cas == 2)))
{ // cas d'une premiere initialisation
Tableau<Enum_variable_metrique> tab(5);
tab(1) = iM0; tab(2) = iMt; tab(3) = iMtdt;
tab(4) = igiH_0; tab(5) = igijHH_0;
met->PlusInitVariables(tab) ;
};
//él้ments de m้trique et matrices de passage
TenseurHH* gijHH;TenseurBB* gijBB;
Mat_pleine jB0(dim,dim),jBfin(dim,dim);
if ((cas==1) || (cas==11))
// calcul d'une ortho interessance de visualisation des tenseurs
// cas de tenseur 3D -> Ia, cas 1D on prend un vecteur norme collineaire
// avec g1, dans le cas 2D
// la nouvelle base jB est calculee dans def par projection de "Ia" sur Galpha
// le resultat est une matrice de passage utilisable pour le changement de base
// jB0 a t=0, B pour les tenseurs BB, jH0 idem pour les tenseurs HH
// resultat a t+dt
{const Met_abstraite::InfoImp& ex = def->RemontImp(absolue,jB0,jBfin);
gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt;
}
else if ((cas==2) || (cas==12))
// resultat a t
{const Met_abstraite::InfoExp_tdt& ex= def->RemontExp_tdt(absolue,jB0,jBfin);
gijHH = ex.gijHH_tdt;gijBB = ex.gijBB_tdt;
};
// maintenant on récupère le tenseur de courbure
const Met_Sfe1::Courbure_t_tdt& curb = ((DeformationSfe1*)def)->RecupCourbure0_t_tdt();
// passage en absolu si nécessaire
if (besoin_tenseur_courbure)
{Tenseur_ns2BB& cur_ns_bBB_tdt = *((Tenseur_ns2BB*) curb.curbBB_tdt) ;
// on le symétrise
curbBB_tdt = NevezTenseurBB(dim);
(*curbBB_tdt).Coor(1,1) = cur_ns_bBB_tdt(1,1);(*curbBB_tdt).Coor(2,2) = cur_ns_bBB_tdt(2,2);
(*curbBB_tdt).Coor(1,2) = 0.5*(cur_ns_bBB_tdt(1,2) + cur_ns_bBB_tdt(2,1));
curbBB_tdt->ChBase(jBfin); // le tenseur de courbure dans le rep่re local
};
// --- bases ortho locales (du plan de la facette ---
if (besoin_dir_princ_courb || besoin_rep_local_ortho)
{ t_ortho_local.Change_taille(3);
t_ortho_local(1) = jBfin(1,1) * curb.aiB_tdt->Coordo(1) + jBfin(2,1) * curb.aiB_tdt->Coordo(2);
t_ortho_local(2) = jBfin(1,2) * curb.aiB_tdt->Coordo(1) + jBfin(2,2) * curb.aiB_tdt->Coordo(2);
t_ortho_local(3) = Util::ProdVec_coor(t_ortho_local(1),t_ortho_local(2));
};
// --- les directions principales et/ou les valeurs principales
int ret=0;
if (besoin_dir_princ_courb)
{ Tenseur_ns2BB& cur_ns_bBB_tdt = *((Tenseur_ns2BB*) curb.curbBB_tdt) ;
// on le sym้trise
Tenseur2BB curb_symBB_tdt;
curb_symBB_tdt.Coor(1,1) = cur_ns_bBB_tdt(1,1);curb_symBB_tdt.Coor(2,2) = cur_ns_bBB_tdt(2,2);
curb_symBB_tdt.Coor(1,2) = 0.5*(cur_ns_bBB_tdt(1,2) + cur_ns_bBB_tdt(2,1));
// on cr้e le BH correspondant
Tenseur2BH curbBH = curb_symBB_tdt * (*gijHH);
Mat_pleine mat(2,2);
val_princ = curbBH.ValPropre(ret,mat);
// // calcul des vecteurs principaux dans la base locale
// Coordonnee3 V_1=mat(1,1) * curb.aiB_tdt->Coordo(1) + mat(2,1) * curb.aiB_tdt->Coordo(2);
// Coordonnee3 V_2=mat(1,2) * curb.aiB_tdt->Coordo(1) + mat(2,2) * curb.aiB_tdt->Coordo(2);
// t_vec_princ.Change_taille(2);
// t_vec_princ(1)(1)= V_1 * t_ortho_local(1); t_vec_princ(1)(2)= V_1 * t_ortho_local(2);
// t_vec_princ(2)(1)= V_2 * t_ortho_local(1); t_vec_princ(2)(2)= V_2 * t_ortho_local(2);
// non finalement on les exprimes dans la base globale sinon on ne peut pas visualiser sans manip
t_vec_princ.Change_taille(2);
// s'il y a eu un pb dans le calcul des valeurs propres ou vecteurs propres on évite
if (ret != -1)
{t_vec_princ(1)=mat(1,1) * curb.aiH_tdt->Coordo(1) + mat(2,1) * curb.aiH_tdt->Coordo(2);
t_vec_princ(2)=mat(1,2) * curb.aiB_tdt->Coordo(1) + mat(2,2) * curb.aiB_tdt->Coordo(2);
t_vec_princ(1).Normer();t_vec_princ(2).Normer();
};
}
else if (besoin_courbures_principales)
// si on arrive ici c'est que l'on n'a pas besoin des directions principales mais
// seulement des valeurs propres de courbure
{ Tenseur_ns2BB& cur_ns_bBB_tdt = *((Tenseur_ns2BB*) curb.curbBB_tdt) ;
// on le sym้trise
Tenseur2BB curb_symBB_tdt;
curb_symBB_tdt.Coor(1,1) = cur_ns_bBB_tdt(1,1);curb_symBB_tdt.Coor(2,2) = cur_ns_bBB_tdt(2,2);
curb_symBB_tdt.Coor(1,2) = 0.5*(cur_ns_bBB_tdt(1,2) + cur_ns_bBB_tdt(2,1));
// on cr้e le BH correspondant
Tenseur2BH curbBH = curb_symBB_tdt * (*gijHH);
val_princ = curbBH.ValPropre(ret);
};
if (ret == -1)
{ cout << "\n pb dans le calcul de valeurs ou vecteurs propres du tenseur de courbure !";
if (ParaGlob::NiveauImpression() > 6 )
cout << "\n SfeMembT::Grandeur_particuliere (...";
};
};
// ------ maintenant on fait la sortie réelle -----
// on commence par appeler la fonction de la classe mère
// il n'y aura pas de calcul des grandeurs inactivées
// absolue: indique si oui ou non on sort les tenseurs dans la base absolue ou une base particulière
ElemMeca::Grandeur_particuliere (absolue,liTQ,iteg);
// puis les grandeurs sp้cifiques
for (il=liTQ.begin();il!=ilfin;il++)
{TypeQuelconque& tipParticu = (*il); // pour simplifier
if (tipParticu.EnuTypeQuelconque().Nom_vide()) // veut dire que c'est un enum pur
switch (tipParticu.EnuTypeQuelconque().EnumTQ())
{
// 1) -----cas de l'้paisseur moyenne initiale, ici elle ne d้pend pas du point d'int้gration
case EPAISSEUR_MOY_INITIALE:
{if (donnee_specif.epais != NULL)
{ *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=donnee_specif.epais->epaisseur0;
(*il).Active();
}; break;
}
// 2) -----cas de l'้paisseur moyenne finale, ici elle ne d้pend pas du point d'int้gration
case EPAISSEUR_MOY_FINALE: (*il).Inactive(); break; // on inactive la grandeur quelconque
{ if (donnee_specif.epais != NULL)
{ *((Grandeur_scalaire_double*) tipParticu.Grandeur_pointee())=donnee_specif.epais->epaisseur_tdt;
(*il).Active();
}; break;
}
// 3) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas du tenseur de courbure ฐฐฐฐฐฐฐฐฐฐ
case TENSEUR_COURBURE:
{ ((Grandeur_TenseurBB*) tipParticu.Grandeur_pointee())->RefConteneurTenseur() = *curbBB_tdt;
(*il).Active();break;
}
// 4) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas du vecteur contenant les courbures principales ฐฐฐฐฐฐฐฐฐฐ
case COURBURES_PRINCIPALES:
{ (*((Tab_Grandeur_scalaire_double*) tipParticu.Grandeur_pointee()))(1) = val_princ(1);
(*((Tab_Grandeur_scalaire_double*) tipParticu.Grandeur_pointee()))(2) = val_princ(2);
(*il).Active();break;
}
// 5) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas des directions principales de courbure (dans le rep่re local ฐฐฐฐฐฐฐฐฐฐ
case DIRECTIONS_PRINC_COURBURE:
{ (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(1) = t_vec_princ(1);
(*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(2) = t_vec_princ(2);
(*il).Active();break;
}
// 6) ฐฐฐฐฐฐฐฐฐฐฐฐฐ cas du rep่re local orthonorm้ d'expression ฐฐฐฐฐฐฐฐฐฐ
case REPERE_LOCAL_ORTHO:
{ (*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(1) = t_ortho_local(1);
(*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(2) = t_ortho_local(2);
(*((Tab_Grandeur_Coordonnee*) tipParticu.Grandeur_pointee()))(3) = t_ortho_local(3);
(*il).Active();break;
};
default: break; // on ne fait rien
};
};
// liberation des tenseurs intermediaires
LibereTenseur();
if (curbBB_tdt != NULL) delete curbBB_tdt;
};
// calcul éventuel de la normale à un noeud
// ce calcul existe pour les éléments 2D, 1D axi, et aussi pour les éléments 1D
// qui possède un repère d'orientation
// en retour coor = la normale si coor.Dimension() est = à la dimension de l'espace
// si le calcul n'existe pas --> coor.Dimension() = 0
// ramène un entier :
// == 1 : calcul normal
// == 0 : problème de calcul -> coor.Dimension() = 0
// == 2 : indique que le calcul n'est pas licite pour le noeud passé en paramètre
// c'est le cas par exemple des noeuds exterieurs pour les éléments SFE
// mais il n'y a pas d'erreur, c'est seulement que l'élément n'est pas ad hoc pour
// calculer la normale à ce noeud là
// temps: indique à quel moment on veut le calcul
// pour des éléments particulier (ex: SFE) la méthode est surchargée
int SfeMembT::CalculNormale_noeud(Enum_dure temps,const Noeud& noe,Coordonnee& coor)
{ // === par rapport à la méthode d'ElemMeca, ici on ne s'occupe des noeuds de l'élément central
int retour = 1; // on part avec l'idée que tout va être ok
Enum_type_geom enutygeom = Type_geom_generique(this->Id_geometrie());
// if ((enutygeom == LIGNE) || (enutygeom == SURFACE))
if (enutygeom != SURFACE) // dans une première étape on ne s'occupe que des surfaces
{coor.Libere(); // pas de calcul possible
retour = 0;
}
else // sinon le calcul est possible
{ // on commence par repérer le noeud dans la numérotation locale
int nuoe=0;
// int borne_nb_noeud=nombre->nbnce+1; // ici on ne prend en compte que les noeuds centraux
int borne_nb_noeud=tab_noeud.Taille()+1;
for (int i=1;i< borne_nb_noeud;i++)
{Noeud& noeloc = *tab_noeud(i);
if ( (noe.Num_noeud() == noeloc.Num_noeud())
&& (noe.Num_Mail() == noeloc.Num_Mail())
)
{nuoe = i; break;
};
};
// on ne continue que si on a trouvé le noeud
if (nuoe != 0)
{ // il y a deux cas: soit le noeud fait partie des noeuds centraux ou non
if (nuoe <= nombre->nbnce) // cas où le noeud fait partie du centre
{ElemGeomC0& elemgeom = ElementGeometrique(); // récup de la géométrie
// récup des coordonnées locales du noeuds
const Coordonnee& theta_noeud = elemgeom.PtelemRef()(nuoe);
// récup des phi et dphi au noeud
const Vecteur & phi = elemgeom.Phi(theta_noeud);
const Mat_pleine& dphi = elemgeom.Dphi(theta_noeud);
switch (temps)
{case TEMPS_0 :
{const BaseB& baseB = met->BaseNat_0(tab_noeud,dphi,phi);
coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
coor.Normer();
break;
}
case TEMPS_t :
{const BaseB& baseB = met->BaseNat_t(tab_noeud,dphi,phi);
coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
coor.Normer();
break;
}
case TEMPS_tdt :
{const BaseB& baseB = met->BaseNat_tdt(tab_noeud,dphi,phi);
coor = Util::ProdVec_coor( baseB.Coordo(1), baseB.Coordo(2));
coor.Normer();
break;
}
default :
cout << "\nErreur : valeur incorrecte du temps demande = "
<< Nom_dure(temps) << " !\n";
cout << "\n SfeMembT::CalculNormale_noeud(Enum_dure temps... \n";
retour = 0;
Sortie(1);
};
}
else // sinon le noeud fait partie des noeuds externes
// on considère que la normale n'est pas à calculer
{coor.Libere(); // pas de calcul possible
retour = 2;
};
}
else
{cout << "\n *** erreur le noeud demande num= "<<noe.Num_noeud()
<< " du maillage "<< noe.Num_Mail()
<< " ne fait pas parti de l'element "
<< num_elt << " du maillage "<< noe.Num_Mail()
<< " on ne peut pas calculer la normale au noeud !!"
<< "\n SfeMembT::CalculNormale_noeud(...";
retour = 0;
Sortie(1);
}
};
// retour
return retour;
};
Reference in a new issue View git blame Copy permalink