2021-09-27 12:42:13 +02:00
|
|
|
// 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.
|
|
|
|
|
//
|
2023-05-03 17:23:49 +02:00
|
|
|
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
|
2021-09-27 12:42:13 +02:00
|
|
|
// 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
|
2023-05-03 17:23:49 +02:00
|
|
|
// ??? const double& jacobien_0 = *(tabSaveDefDon(indice)->Meti_00().jacobien_);
|
2021-09-27 12:42:13 +02:00
|
|
|
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);
|
2023-05-03 17:23:49 +02:00
|
|
|
return ConstMath::trespetit; // taire le compilo
|
2021-09-27 12:42:13 +02:00
|
|
|
};
|
|
|
|
|
|
|
|
|
|
//test si le jacobien due aux gi finaux est très différent du jacobien de la facette
|
2023-05-03 17:23:49 +02:00
|
|
|
void SfeMembT::Delta_Jacobien_anormal(const Deformation::SaveDefResul* don
|
2021-09-27 12:42:13 +02:00
|
|
|
,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();
|
2023-05-03 17:23:49 +02:00
|
|
|
// PtIntegMecaInterne & ptIntegMeca = (*lesPtIntegMecaInterne)(iteg);
|
2021-09-27 12:42:13 +02:00
|
|
|
// on recupère le tableau pour la lecture des coordonnées des tenseurs
|
2023-05-03 17:23:49 +02:00
|
|
|
// int nbcompo = ParaGlob::NbCompTens();
|
|
|
|
|
// const Tableau2 <int> & OrdreCont = OrdreContrainteR(nbcompo);
|
2021-09-27 12:42:13 +02:00
|
|
|
// 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
|
2023-05-03 17:23:49 +02:00
|
|
|
const Vecteur & phi = elemgeom.Phi_point(theta_noeud);
|
|
|
|
|
const Mat_pleine& dphi = elemgeom.Dphi_point(theta_noeud);
|
2021-09-27 12:42:13 +02:00
|
|
|
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;
|
|
|
|
|
|
|
|
|
|
};
|
