Skip to content

Upgrade beta law #173

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lelandaisb merged 3 commits into from
Nov 26, 2025
Merged

Upgrade beta law #173

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
634b3b2
Add init with first edge for beta law
mpoudot Nov 14, 2025
40af2cc
Add test for beta law
mpoudot Nov 20, 2025
04403d8
Merge branch 'main' into upgrade-beta-law
mpoudot Nov 20, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
95 changes: 45 additions & 50 deletions src/Core/Topo/EdgeMeshingPropertyBeta.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,30 +17,27 @@ namespace Mgx3D {
/*----------------------------------------------------------------------------*/
namespace Topo {
/*----------------------------------------------------------------------------*/
#define _DEBUG_beta
/*----------------------------------------------------------------------------*/EdgeMeshingPropertyBeta::
EdgeMeshingPropertyBeta(int nb, double beta, bool isDirect, bool initWithFirstEdge, double meshingEdgeLength)
: CoEdgeMeshingProperty(nb, beta_resserrement, isDirect)
//#define _DEBUG_beta
/*----------------------------------------------------------------------------*/
EdgeMeshingPropertyBeta::
EdgeMeshingPropertyBeta(int nbBras, double beta, bool isDirect,
bool initWithFirstEdge, double meshingEdgeLength)
: CoEdgeMeshingProperty(nbBras, beta_resserrement, isDirect)
, m_beta(beta)
, m_arm1(meshingEdgeLength)
, m_initWithArm1(initWithFirstEdge)
{
if (!m_initWithArm1 && m_beta<=1.0){
if (!m_initWithArm1 && (m_beta<1.00001 || m_beta>1.01)){
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
messErr << "EdgeMeshingPropertyBeta, la valeur de beta doit être légèrement supérieur à 1, actuellement : "<<m_beta;
messErr << "EdgeMeshingPropertyBeta, la valeur de beta doit être comprise entre 1 et 1.01, actuellement : "<<m_beta;
throw TkUtil::Exception(messErr);
}
if (!m_initWithArm1 && m_beta>2.0){
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
messErr << "EdgeMeshingPropertyBeta, la longueur de beta n'a pas été testé jusque là "<<m_beta;
throw TkUtil::Exception(messErr);
}
if (m_initWithArm1 && m_arm1 < 0){
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
messErr << "EdgeMeshingPropertyGeometric, la longueur du premier bras ne doit pas être négative : "<<m_arm1;
throw TkUtil::Exception(messErr);
}
if (nb<1)
if (nbBras < 1)
throw TkUtil::Exception (TkUtil::UTF8String ("EdgeMeshingPropertyBeta, le nombre de bras doit être au moins de 1", TkUtil::Charset::UTF_8));
}
/*----------------------------------------------------------------------------*/
Expand Down Expand Up @@ -82,63 +79,61 @@ bool EdgeMeshingPropertyBeta::operator == (const CoEdgeMeshingProperty& cedp) co
return CoEdgeMeshingProperty::operator == (cedp);
}
/*----------------------------------------------------------------------------*/

void EdgeMeshingPropertyBeta::setNbEdges(const int nb)
{
CoEdgeMeshingProperty::setNbEdges(nb);
}
/*----------------------------------------------------------------------------*/
void EdgeMeshingPropertyBeta::initCoeff(double length)
{
//if (m_initWithArm1){
// c'est l'occasion de calculer la raison
computeBeta(0.0334298);
//}
if (m_initWithArm1){
// on calcule le coefficient de resserement
m_beta = computeBeta(m_arm1);
}
initCoeff();
}

double EdgeMeshingPropertyBeta::
nextCoeff()
void EdgeMeshingPropertyBeta::initCoeff()
{
m_dernierIndice+=1;

//#ifdef _DEBUG
if (m_dernierIndice>m_nb_edges){
#ifdef _DEBUG_beta
std::cout<<"EdgeMeshingPropertyBeta::initCoeff() "<<std::endl;
#endif
if (m_beta==0.0){
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
messErr<<"EdgeMeshingPropertyBeta::nextCoeff est en dehors des limites: dernierIndice : "
<<(long)m_dernierIndice<<", nb de bras : "<<(long)m_nb_edges;
messErr << "Erreur interne: EdgeMeshingPropertyBeta, le coefficient de resserrement n'a pas été initialisé";
throw TkUtil::Exception(messErr);
}
//#endif

m_dernierIndice = 0;
}

double EdgeMeshingPropertyBeta::
nextCoeff()
{
m_dernierIndice+=1;

double feta = 0.0;
double eta = 0.0;
if (m_sens){
eta = ((double)m_dernierIndice)/((double)m_nb_edges);
std::cout << "eta = " << eta << std::endl;
feta=resserre(eta, m_beta);
std::cout << "feta = " << feta << std::endl;
}
else {
eta = ((double)m_nb_edges-m_dernierIndice)/((double)m_nb_edges);
feta = 1.0-resserre(eta, m_beta);
}

//#ifdef _DEBUG2
std::cout<<"EdgeMeshingPropertyBeta::nextCoeff (m_beta="<<m_beta<<") return "<<feta<<" pour eta = "<<eta<<std::endl;
//#endif
#ifdef _DEBUG_beta
std::cout<<"EdgeMeshingPropertyBeta::nextCoeff (m_beta="<<m_beta<<") return "<<feta<<std::endl;
#endif

return feta;
}
/*----------------------------------------------------------------------------*/
double EdgeMeshingPropertyBeta::resserre(double eta, double beta)
{
std::cout << "beta = " << beta << std::endl;
double ratio = (beta + 1.0) / (beta - 1.0);
std::cout << "ratio = " << ratio << std::endl;
double zlog = std::log(ratio);
std::cout << "zlog = " << zlog << std::endl;
double puiss = zlog * (1.0 - eta);
puiss = std::exp(puiss);
std::cout << "puiss = " << puiss << std::endl;
double rapp = (1.0 - puiss) / (1.0 + puiss);
std::cout << "rapp = " << rapp << std::endl;
double feta = 1.0 + beta * rapp;
std::cout << "feta = " << feta << std::endl;
return feta;
double p = std::log((beta + 1.0)/(beta - 1.0)) * (1.0 - eta);
return 1.0 + beta * (1.0 - std::exp(p)) / (1.0 + std::exp(p));
}
/*----------------------------------------------------------------------------*/
TkUtil::UTF8String EdgeMeshingPropertyBeta::
Expand Down Expand Up @@ -187,18 +182,18 @@ computeBeta(const double lg)
std::cout << "eta = " << eta << std::endl;
}


// on limite la recherche
double beta1 = 1.0000001;
double beta2 = 2.0;
double beta1 = 1.00001;
double beta2 = 1.01;
double lg1 = resserre(eta, beta1);
double lg2 = resserre(eta, beta2);

std::cout << "lg1 = " << lg1 << ", lg2 = " << lg2 << ", lg = " << lg << std::endl;
if ((lg1<lg && lg2<lg) || (lg1>lg && lg2>lg)){
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
messErr << "EdgeMeshingPropertyBeta, on ne peut pas trouver de beta pour lg "
<<lg<<" et nombre de bras de "<<(short)m_nb_edges;
throw TkUtil::Exception(messErr);
throw TkUtil::Exception(messErr);
}
uint iter = 0;
bool sens = (lg1<lg2);
Expand All @@ -208,8 +203,8 @@ computeBeta(const double lg)

double lg_iter = resserre(eta, beta);
#ifdef _DEBUG_beta
std::cout<<"iter "<<iter<<std::endl;
std::cout<<" beta "<<beta<<", lg_iter "<<lg_iter<<std::endl;
//std::cout<<"iter "<<iter<<std::endl;
//std::cout<<" beta "<<beta<<", lg_iter "<<lg_iter<<std::endl;
#endif

if (sens){
Expand Down
19 changes: 6 additions & 13 deletions src/Core/Topo/EdgeMeshingPropertyGeometric.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,8 @@ namespace Topo {
//#define _DEBUG_raison
/*----------------------------------------------------------------------------*/
EdgeMeshingPropertyGeometric::
EdgeMeshingPropertyGeometric(int nbBras, double raison, bool isDirect, bool initWithFirstEdge, double meshingEdgeLength)
EdgeMeshingPropertyGeometric(int nbBras, double raison, bool isDirect,
bool initWithFirstEdge, double meshingEdgeLength)
: CoEdgeMeshingProperty(nbBras, geometrique, isDirect)
, m_raison(raison)
, m_arm1(meshingEdgeLength)
Expand All @@ -39,7 +40,7 @@ EdgeMeshingPropertyGeometric(int nbBras, double raison, bool isDirect, bool init
messErr << "EdgeMeshingPropertyGeometric, la longueur du premier bras ne doit pas être négative : "<<m_arm1;
throw TkUtil::Exception(messErr);
}
if (nbBras<1)
if (nbBras < 1)
throw TkUtil::Exception (TkUtil::UTF8String ("EdgeMeshingPropertyGeometric, le nombre de bras doit être au moins de 1", TkUtil::Charset::UTF_8));
// calcul de la somme à l'aide de la raison
if (!m_initWithArm1)
Expand Down Expand Up @@ -138,7 +139,7 @@ void EdgeMeshingPropertyGeometric::initSomme()
void EdgeMeshingPropertyGeometric::initCoeff(double length)
{
if (m_initWithArm1){
// c'est l'occasion de calculer la raison
// on calcule la raison
m_raison = computeRaison(length/m_arm1);
initSomme();
}
Expand Down Expand Up @@ -181,14 +182,6 @@ nextCoeff()
{
m_dernierIndice+=1;

#ifdef _DEBUG
if (m_dernierIndice>m_nb_edges){
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
messErr<<"EdgeMeshingPropertyGeometric::nextCoeff est en dehors des limites: dernierIndice : "
<<(long)m_dernierIndice<<", nb de bras : "<<(long)m_nb_edges;
throw TkUtil::Exception(messErr);
}
#endif
if (m_sens){
m_dernierCoeff *= m_raison;
m_dernierSommeCoeff += m_dernierCoeff;
Expand Down Expand Up @@ -280,8 +273,8 @@ computeRaison(const double lg)

double lg_iter = computePolynome(raison);
#ifdef _DEBUG_raison
std::cout<<"iter "<<iter<<std::endl;
std::cout<<" raison "<<raison<<", lg_iter "<<lg_iter<<std::endl;
//std::cout<<"iter "<<iter<<std::endl;
//std::cout<<" raison "<<raison<<", lg_iter "<<lg_iter<<std::endl;
#endif

if (sens){
Expand Down
31 changes: 22 additions & 9 deletions src/Core/protected/Topo/EdgeMeshingPropertyBeta.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,15 @@ class EdgeMeshingPropertyBeta : public CoEdgeMeshingProperty {
public:

/*------------------------------------------------------------------------*/
/** Constructeur avec un nombre de bras spécifié et la valeur de beta
*
/** Constructeur avec un nombre de bras spécifié,
* la valeur de beta
* et la longueur du premier bras
* \param beta coeff de resserrement , beta=1+epsilon ; epsilon> 0;
* en pratique beta peut varier de 1.01 (peu resséré) à 1.00001 (très resséré)
* Si la longueur du premier bras est non nulle, alors beta est calculé
*/
EdgeMeshingPropertyBeta(int nb, double beta, bool isDirect=true, bool initWithFirstEdge=false, double meshingEdgeLength=0.0);

/// Pour les pythoneries, cast par opérateur = :
EdgeMeshingPropertyBeta (const CoEdgeMeshingProperty&);

Expand All @@ -43,19 +46,31 @@ class EdgeMeshingPropertyBeta : public CoEdgeMeshingProperty {
/// Accesseur sur le nom de la méthode de maillage
virtual std::string getMeshLawName() const {return "Beta";}

/// change le nombre de bras demandé pour un côté
virtual void setNbEdges(const int nb);

/// le ratio entre deux bras successifs
double getBeta() const {return m_beta;}

#ifndef SWIG
bool initWithFirstEdge() const {return m_initWithArm1;}

/// longueur de la première arête
double getFirstEdgeLength() const {return m_arm1;}
void setFirstEdgeLength(const double & arm1) {m_arm1 = arm1;}

#ifndef SWIG
/** Création d'un clone, on copie toutes les informations */
virtual CoEdgeMeshingProperty* clone() {return new EdgeMeshingPropertyBeta(*this);}

/// retourne le coefficient suivant pour les noeuds internes
virtual double nextCoeff();
// initialisation avant appel à nextCoeff
virtual void initCoeff();

/// initialisation avant appel à nextCoeff pour le cas où la longueur est utile
virtual void initCoeff(double length);

/// retourne le coefficient suivant pour les noeuds internes
virtual double nextCoeff();

/*------------------------------------------------------------------------*/
/// Script pour la commande de création Python
virtual TkUtil::UTF8String getScriptCommand() const;
Expand All @@ -72,9 +87,7 @@ class EdgeMeshingPropertyBeta : public CoEdgeMeshingProperty {
/** f(eta) : est en sortie de resserre, ça peut être vu comme la
position du point après resserrement (feta=0 --> paroi,
feta=1--> frontière amont)

transformation de (0,1) sur (0,1) permettant de resserrer
le maillage au voisinage du corps
transformation de (0,1) sur (0,1) permettant de resserrer le maillage au voisinage du corps
*/
double resserre(double eta, double beta);

Expand All @@ -83,8 +96,8 @@ class EdgeMeshingPropertyBeta : public CoEdgeMeshingProperty {

/// valeur de beta pour le resserrement (de la forme 1+epsilon)
double m_beta;

double m_arm1;

/// vrai si c'est la longueur du premier bras qui permet de déterminer beta
bool m_initWithArm1;
};
Expand Down
9 changes: 5 additions & 4 deletions src/Core/protected/Topo/EdgeMeshingPropertyGeometric.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,11 +31,14 @@ class EdgeMeshingPropertyGeometric : public CoEdgeMeshingProperty {
* Si la longueur du premier bras est non nulle, alors le ratio est calculé
* à l'aide de la longueur de l'arête (projetée s'il y a lieu)
*/
EdgeMeshingPropertyGeometric(int nb, double raison, bool isDirect=true, bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
EdgeMeshingPropertyGeometric(int nb, double raison, bool isDirect=true,
bool initWithFirstEdge=false, double meshingEdgeLength=0.0);

/// idem avec un centre pour découpage polaire
EdgeMeshingPropertyGeometric(int nb, double raison,
Utils::Math::Point polar_center, bool isDirect=true, bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
Utils::Math::Point polar_center, bool isDirect=true,
bool initWithFirstEdge=false, double meshingEdgeLength=0.0);

/// Pour les pythoneries, cast par opérateur = :
EdgeMeshingPropertyGeometric (const CoEdgeMeshingProperty&);
/*------------------------------------------------------------------------*/
Expand All @@ -51,14 +54,12 @@ class EdgeMeshingPropertyGeometric : public CoEdgeMeshingProperty {

/// le ratio entre deux bras successifs
double getRatio() const {return m_raison;}

void setRatio(const double & ratio) {m_raison = ratio;}

bool initWithFirstEdge() const {return m_initWithArm1;}

/// longueur de la première arête
double getFirstEdgeLength() const {return m_arm1;}

void setFirstEdgeLength(const double & arm1) {m_arm1 = arm1;}

#ifndef SWIG
Expand Down
Loading