Skip to content

Commit b9ac706

Browse files
mpoudotlelandaisb
mpoudot
authored and
lelandaisb
committed
Add init with first edge for beta law
1 parent 3f2de38 commit b9ac706

File tree

6 files changed

+165
-93
lines changed

6 files changed

+165
-93
lines changed

‎src/Core/Topo/EdgeMeshingPropertyBeta.cpp‎

Lines changed: 45 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -17,30 +17,27 @@ namespace Mgx3D {
1717
/*----------------------------------------------------------------------------*/
1818
namespace Topo {
1919
/*----------------------------------------------------------------------------*/
20-
#define _DEBUG_beta
21-
/*----------------------------------------------------------------------------*/EdgeMeshingPropertyBeta::
22-
EdgeMeshingPropertyBeta(int nb, double beta, bool isDirect, bool initWithFirstEdge, double meshingEdgeLength)
23-
: CoEdgeMeshingProperty(nb, beta_resserrement, isDirect)
20+
//#define _DEBUG_beta
21+
/*----------------------------------------------------------------------------*/
22+
EdgeMeshingPropertyBeta::
23+
EdgeMeshingPropertyBeta(int nbBras, double beta, bool isDirect,
24+
bool initWithFirstEdge, double meshingEdgeLength)
25+
: CoEdgeMeshingProperty(nbBras, beta_resserrement, isDirect)
2426
, m_beta(beta)
2527
, m_arm1(meshingEdgeLength)
2628
, m_initWithArm1(initWithFirstEdge)
2729
{
28-
if (!m_initWithArm1 && m_beta<=1.0){
30+
if (!m_initWithArm1 && (m_beta<1.00001 || m_beta>1.01)){
2931
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
30-
messErr << "EdgeMeshingPropertyBeta, la valeur de beta doit être légèrement supérieur à 1, actuellement : "<<m_beta;
32+
messErr << "EdgeMeshingPropertyBeta, la valeur de beta doit être comprise entre 1 et 1.01, actuellement : "<<m_beta;
3133
throw TkUtil::Exception(messErr);
3234
}
33-
if (!m_initWithArm1 && m_beta>2.0){
34-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
35-
messErr << "EdgeMeshingPropertyBeta, la longueur de beta n'a pas été testé jusque là "<<m_beta;
36-
throw TkUtil::Exception(messErr);
37-
}
3835
if (m_initWithArm1 && m_arm1 < 0){
3936
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
4037
messErr << "EdgeMeshingPropertyGeometric, la longueur du premier bras ne doit pas être négative : "<<m_arm1;
4138
throw TkUtil::Exception(messErr);
4239
}
43-
if (nb<1)
40+
if (nbBras < 1)
4441
throw TkUtil::Exception (TkUtil::UTF8String ("EdgeMeshingPropertyBeta, le nombre de bras doit être au moins de 1", TkUtil::Charset::UTF_8));
4542
}
4643
/*----------------------------------------------------------------------------*/
@@ -82,63 +79,61 @@ bool EdgeMeshingPropertyBeta::operator == (const CoEdgeMeshingProperty& cedp) co
8279
return CoEdgeMeshingProperty::operator == (cedp);
8380
}
8481
/*----------------------------------------------------------------------------*/
85-
82+
void EdgeMeshingPropertyBeta::setNbEdges(const int nb)
83+
{
84+
CoEdgeMeshingProperty::setNbEdges(nb);
85+
}
86+
/*----------------------------------------------------------------------------*/
8687
void EdgeMeshingPropertyBeta::initCoeff(double length)
8788
{
88-
//if (m_initWithArm1){
89-
// c'est l'occasion de calculer la raison
90-
computeBeta(0.0334298);
91-
//}
89+
if (m_initWithArm1){
90+
// on calcule le coefficient de resserement
91+
m_beta = computeBeta(m_arm1);
92+
}
93+
initCoeff();
9294
}
9395

94-
double EdgeMeshingPropertyBeta::
95-
nextCoeff()
96+
void EdgeMeshingPropertyBeta::initCoeff()
9697
{
97-
m_dernierIndice+=1;
98-
99-
//#ifdef _DEBUG
100-
if (m_dernierIndice>m_nb_edges){
98+
#ifdef _DEBUG_beta
99+
std::cout<<"EdgeMeshingPropertyBeta::initCoeff() "<<std::endl;
100+
#endif
101+
if (m_beta==0.0){
101102
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
102-
messErr<<"EdgeMeshingPropertyBeta::nextCoeff est en dehors des limites: dernierIndice : "
103-
<<(long)m_dernierIndice<<", nb de bras : "<<(long)m_nb_edges;
103+
messErr << "Erreur interne: EdgeMeshingPropertyBeta, le coefficient de resserrement n'a pas été initialisé";
104104
throw TkUtil::Exception(messErr);
105105
}
106-
//#endif
106+
107+
m_dernierIndice = 0;
108+
}
109+
110+
double EdgeMeshingPropertyBeta::
111+
nextCoeff()
112+
{
113+
m_dernierIndice+=1;
107114

108115
double feta = 0.0;
109116
double eta = 0.0;
110117
if (m_sens){
111118
eta = ((double)m_dernierIndice)/((double)m_nb_edges);
112-
std::cout << "eta = " << eta << std::endl;
113119
feta=resserre(eta, m_beta);
114-
std::cout << "feta = " << feta << std::endl;
115120
}
116121
else {
117122
eta = ((double)m_nb_edges-m_dernierIndice)/((double)m_nb_edges);
118123
feta = 1.0-resserre(eta, m_beta);
119124
}
120125

121-
//#ifdef _DEBUG2
122-
std::cout<<"EdgeMeshingPropertyBeta::nextCoeff (m_beta="<<m_beta<<") return "<<feta<<" pour eta = "<<eta<<std::endl;
123-
//#endif
126+
#ifdef _DEBUG_beta
127+
std::cout<<"EdgeMeshingPropertyBeta::nextCoeff (m_beta="<<m_beta<<") return "<<feta<<std::endl;
128+
#endif
129+
124130
return feta;
125131
}
126132
/*----------------------------------------------------------------------------*/
127133
double EdgeMeshingPropertyBeta::resserre(double eta, double beta)
128134
{
129-
std::cout << "beta = " << beta << std::endl;
130-
double ratio = (beta + 1.0) / (beta - 1.0);
131-
std::cout << "ratio = " << ratio << std::endl;
132-
double zlog = std::log(ratio);
133-
std::cout << "zlog = " << zlog << std::endl;
134-
double puiss = zlog * (1.0 - eta);
135-
puiss = std::exp(puiss);
136-
std::cout << "puiss = " << puiss << std::endl;
137-
double rapp = (1.0 - puiss) / (1.0 + puiss);
138-
std::cout << "rapp = " << rapp << std::endl;
139-
double feta = 1.0 + beta * rapp;
140-
std::cout << "feta = " << feta << std::endl;
141-
return feta;
135+
double p = std::log((beta + 1.0)/(beta - 1.0)) * (1.0 - eta);
136+
return 1.0 + beta * (1.0 - std::exp(p)) / (1.0 + std::exp(p));
142137
}
143138
/*----------------------------------------------------------------------------*/
144139
TkUtil::UTF8String EdgeMeshingPropertyBeta::
@@ -187,18 +182,18 @@ computeBeta(const double lg)
187182
std::cout << "eta = " << eta << std::endl;
188183
}
189184

190-
191185
// on limite la recherche
192-
double beta1 = 1.0000001;
193-
double beta2 = 2.0;
186+
double beta1 = 1.00001;
187+
double beta2 = 1.01;
194188
double lg1 = resserre(eta, beta1);
195189
double lg2 = resserre(eta, beta2);
196190

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

209204
double lg_iter = resserre(eta, beta);
210205
#ifdef _DEBUG_beta
211-
std::cout<<"iter "<<iter<<std::endl;
212-
std::cout<<" beta "<<beta<<", lg_iter "<<lg_iter<<std::endl;
206+
//std::cout<<"iter "<<iter<<std::endl;
207+
//std::cout<<" beta "<<beta<<", lg_iter "<<lg_iter<<std::endl;
213208
#endif
214209

215210
if (sens){

‎src/Core/Topo/EdgeMeshingPropertyGeometric.cpp‎

Lines changed: 6 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@ namespace Topo {
2020
//#define _DEBUG_raison
2121
/*----------------------------------------------------------------------------*/
2222
EdgeMeshingPropertyGeometric::
23-
EdgeMeshingPropertyGeometric(int nbBras, double raison, bool isDirect, bool initWithFirstEdge, double meshingEdgeLength)
23+
EdgeMeshingPropertyGeometric(int nbBras, double raison, bool isDirect,
24+
bool initWithFirstEdge, double meshingEdgeLength)
2425
: CoEdgeMeshingProperty(nbBras, geometrique, isDirect)
2526
, m_raison(raison)
2627
, m_arm1(meshingEdgeLength)
@@ -39,7 +40,7 @@ EdgeMeshingPropertyGeometric(int nbBras, double raison, bool isDirect, bool init
3940
messErr << "EdgeMeshingPropertyGeometric, la longueur du premier bras ne doit pas être négative : "<<m_arm1;
4041
throw TkUtil::Exception(messErr);
4142
}
42-
if (nbBras<1)
43+
if (nbBras < 1)
4344
throw TkUtil::Exception (TkUtil::UTF8String ("EdgeMeshingPropertyGeometric, le nombre de bras doit être au moins de 1", TkUtil::Charset::UTF_8));
4445
// calcul de la somme à l'aide de la raison
4546
if (!m_initWithArm1)
@@ -138,7 +139,7 @@ void EdgeMeshingPropertyGeometric::initSomme()
138139
void EdgeMeshingPropertyGeometric::initCoeff(double length)
139140
{
140141
if (m_initWithArm1){
141-
// c'est l'occasion de calculer la raison
142+
// on calcule la raison
142143
m_raison = computeRaison(length/m_arm1);
143144
initSomme();
144145
}
@@ -181,14 +182,6 @@ nextCoeff()
181182
{
182183
m_dernierIndice+=1;
183184

184-
#ifdef _DEBUG
185-
if (m_dernierIndice>m_nb_edges){
186-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
187-
messErr<<"EdgeMeshingPropertyGeometric::nextCoeff est en dehors des limites: dernierIndice : "
188-
<<(long)m_dernierIndice<<", nb de bras : "<<(long)m_nb_edges;
189-
throw TkUtil::Exception(messErr);
190-
}
191-
#endif
192185
if (m_sens){
193186
m_dernierCoeff *= m_raison;
194187
m_dernierSommeCoeff += m_dernierCoeff;
@@ -280,8 +273,8 @@ computeRaison(const double lg)
280273

281274
double lg_iter = computePolynome(raison);
282275
#ifdef _DEBUG_raison
283-
std::cout<<"iter "<<iter<<std::endl;
284-
std::cout<<" raison "<<raison<<", lg_iter "<<lg_iter<<std::endl;
276+
//std::cout<<"iter "<<iter<<std::endl;
277+
//std::cout<<" raison "<<raison<<", lg_iter "<<lg_iter<<std::endl;
285278
#endif
286279

287280
if (sens){

‎src/Core/protected/Topo/EdgeMeshingPropertyBeta.h‎

Lines changed: 22 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -26,12 +26,15 @@ class EdgeMeshingPropertyBeta : public CoEdgeMeshingProperty {
2626
public:
2727

2828
/*------------------------------------------------------------------------*/
29-
/** Constructeur avec un nombre de bras spécifié et la valeur de beta
30-
*
29+
/** Constructeur avec un nombre de bras spécifié,
30+
* la valeur de beta
31+
* et la longueur du premier bras
3132
* \param beta coeff de resserrement , beta=1+epsilon ; epsilon> 0;
3233
* en pratique beta peut varier de 1.01 (peu resséré) à 1.00001 (très resséré)
34+
* Si la longueur du premier bras est non nulle, alors beta est calculé
3335
*/
3436
EdgeMeshingPropertyBeta(int nb, double beta, bool isDirect=true, bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
37+
3538
/// Pour les pythoneries, cast par opérateur = :
3639
EdgeMeshingPropertyBeta (const CoEdgeMeshingProperty&);
3740

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

49+
/// change le nombre de bras demandé pour un côté
50+
virtual void setNbEdges(const int nb);
51+
4652
/// le ratio entre deux bras successifs
4753
double getBeta() const {return m_beta;}
4854

49-
#ifndef SWIG
55+
bool initWithFirstEdge() const {return m_initWithArm1;}
56+
57+
/// longueur de la première arête
58+
double getFirstEdgeLength() const {return m_arm1;}
59+
void setFirstEdgeLength(const double & arm1) {m_arm1 = arm1;}
60+
61+
#ifndef SWIG
5062
/** Création d'un clone, on copie toutes les informations */
5163
virtual CoEdgeMeshingProperty* clone() {return new EdgeMeshingPropertyBeta(*this);}
5264

53-
/// retourne le coefficient suivant pour les noeuds internes
54-
virtual double nextCoeff();
65+
// initialisation avant appel à nextCoeff
66+
virtual void initCoeff();
5567

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

71+
/// retourne le coefficient suivant pour les noeuds internes
72+
virtual double nextCoeff();
73+
5974
/*------------------------------------------------------------------------*/
6075
/// Script pour la commande de création Python
6176
virtual TkUtil::UTF8String getScriptCommand() const;
@@ -72,9 +87,7 @@ class EdgeMeshingPropertyBeta : public CoEdgeMeshingProperty {
7287
/** f(eta) : est en sortie de resserre, ça peut être vu comme la
7388
position du point après resserrement (feta=0 --> paroi,
7489
feta=1--> frontière amont)
75-
76-
transformation de (0,1) sur (0,1) permettant de resserrer
77-
le maillage au voisinage du corps
90+
transformation de (0,1) sur (0,1) permettant de resserrer le maillage au voisinage du corps
7891
*/
7992
double resserre(double eta, double beta);
8093

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

8497
/// valeur de beta pour le resserrement (de la forme 1+epsilon)
8598
double m_beta;
86-
8799
double m_arm1;
100+
88101
/// vrai si c'est la longueur du premier bras qui permet de déterminer beta
89102
bool m_initWithArm1;
90103
};

‎src/Core/protected/Topo/EdgeMeshingPropertyGeometric.h‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -31,11 +31,14 @@ class EdgeMeshingPropertyGeometric : public CoEdgeMeshingProperty {
3131
* Si la longueur du premier bras est non nulle, alors le ratio est calculé
3232
* à l'aide de la longueur de l'arête (projetée s'il y a lieu)
3333
*/
34-
EdgeMeshingPropertyGeometric(int nb, double raison, bool isDirect=true, bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
34+
EdgeMeshingPropertyGeometric(int nb, double raison, bool isDirect=true,
35+
bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
3536

3637
/// idem avec un centre pour découpage polaire
3738
EdgeMeshingPropertyGeometric(int nb, double raison,
38-
Utils::Math::Point polar_center, bool isDirect=true, bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
39+
Utils::Math::Point polar_center, bool isDirect=true,
40+
bool initWithFirstEdge=false, double meshingEdgeLength=0.0);
41+
3942
/// Pour les pythoneries, cast par opérateur = :
4043
EdgeMeshingPropertyGeometric (const CoEdgeMeshingProperty&);
4144
/*------------------------------------------------------------------------*/
@@ -51,14 +54,12 @@ class EdgeMeshingPropertyGeometric : public CoEdgeMeshingProperty {
5154

5255
/// le ratio entre deux bras successifs
5356
double getRatio() const {return m_raison;}
54-
5557
void setRatio(const double & ratio) {m_raison = ratio;}
5658

5759
bool initWithFirstEdge() const {return m_initWithArm1;}
5860

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

6465
#ifndef SWIG

0 commit comments

Comments
 (0)