From 164b5be98ca29482a215f78fa2cfe01adbe02615 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 21 Nov 2013 15:09:34 +0000
Subject: [PATCH] improve reproductibility of 2D mesh algorithm (compare
vertices and triangles by number, not by pointer)
---
Geo/MEdge.h | 19 +++++++-------
Mesh/meshGFace.cpp | 42 +++++++++++++++---------------
Mesh/meshGFaceOptimize.h | 55 +++++++++++++++++++---------------------
3 files changed, 56 insertions(+), 60 deletions(-)
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index 94e31eed59..86c9857c27 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -25,7 +25,7 @@ class MEdge {
MEdge(MVertex *v0, MVertex *v1)
{
_v[0] = v0; _v[1] = v1;
- if(_v[1] < _v[0]) {
+ if(_v[1]->getNum() < _v[0]->getNum()) {
_si[0] = 1;
_si[1] = 0;
}
@@ -41,13 +41,13 @@ class MEdge {
inline MVertex *getMaxVertex() const { return _v[int(_si[1])]; }
SVector3 scaledTangent() const
{
- return SVector3(_v[1]->x() - _v[0]->x(),
+ return SVector3(_v[1]->x() - _v[0]->x(),
_v[1]->y() - _v[0]->y(),
_v[1]->z() - _v[0]->z());
}
SVector3 tangent() const
{
- SVector3 t(_v[1]->x() - _v[0]->x(),
+ SVector3 t(_v[1]->x() - _v[0]->x(),
_v[1]->y() - _v[0]->y(),
_v[1]->z() - _v[0]->z());
t.normalize();
@@ -55,12 +55,12 @@ class MEdge {
}
double length() const
{
- SVector3 t(_v[1]->x() - _v[0]->x(),
+ SVector3 t(_v[1]->x() - _v[0]->x(),
_v[1]->y() - _v[0]->y(),
_v[1]->z() - _v[0]->z());
return t.norm();
}
- SVector3 normal() const
+ SVector3 normal() const
{
// this computes one of the normals to the edge
SVector3 t = tangent(), ex(0., 0., 0.);
@@ -98,7 +98,7 @@ inline bool operator!=(const MEdge &e1, const MEdge &e2)
return (e1.getMinVertex() != e2.getMinVertex() ||
e1.getMaxVertex() != e2.getMaxVertex());
}
-
+
struct Equal_Edge : public std::binary_function<MEdge, MEdge, bool> {
bool operator()(const MEdge &e1, const MEdge &e2) const
{
@@ -109,12 +109,11 @@ struct Equal_Edge : public std::binary_function<MEdge, MEdge, bool> {
struct Less_Edge : public std::binary_function<MEdge, MEdge, bool> {
bool operator()(const MEdge &e1, const MEdge &e2) const
{
- if(e1.getMinVertex() < e2.getMinVertex()) return true;
- if(e1.getMinVertex() > e2.getMinVertex()) return false;
- if(e1.getMaxVertex() < e2.getMaxVertex()) return true;
+ if(e1.getMinVertex()->getNum() < e2.getMinVertex()->getNum()) return true;
+ if(e1.getMinVertex()->getNum() > e2.getMinVertex()->getNum()) return false;
+ if(e1.getMaxVertex()->getNum() < e2.getMaxVertex()->getNum()) return true;
return false;
}
-
};
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 1590d0d43d..a7e6b12dbf 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -2509,9 +2509,9 @@ void directions_storage(GFace* gf){
MElementOctree* octree;
std::set<MVertex*> vertices;
std::set<MVertex*>::iterator it;
-
- vertices.clear();
-
+
+ vertices.clear();
+
for(i=0;i<gf->getNumMeshElements();i++){
element = gf->getMeshElement(i);
for(j=0;j<element->getNumVertices();j++){
@@ -2519,18 +2519,18 @@ void directions_storage(GFace* gf){
vertices.insert(vertex);
}
}
-
+
backgroundMesh::set(gf);
octree = backgroundMesh::current()->get_octree();
-
+
gf->storage1.clear();
gf->storage2.clear();
gf->storage3.clear();
gf->storage4.clear();
-
+
for(it=vertices.begin();it!=vertices.end();it++){
ok = 0;
-
+
if(!gf->getCompound()){
if(gf->geomType()==GEntity::CompoundSurface){
ok = translate(gf,octree,*it,SPoint2(0.0,0.0),v1,v2);
@@ -2539,16 +2539,16 @@ void directions_storage(GFace* gf){
ok = improved_translate(gf,*it,v1,v2);
}
}
-
+
if(ok){
gf->storage1.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
gf->storage2.push_back(v1);
gf->storage3.push_back(v2);
reparamMeshVertexOnFace(*it,gf,point);
gf->storage4.push_back(backgroundMesh::current()->operator()(point.x(),point.y(),0.0));
- }
+ }
}
-
+
backgroundMesh::unset();
}
@@ -2565,7 +2565,7 @@ bool translate(GFace* gf,MElementOctree* octree,MVertex* vertex,SPoint2 corr,SVe
SPoint2 point;
GPoint gp1;
GPoint gp2;
-
+
ok = true;
k = 0.0001;
reparamMeshVertexOnFace(vertex,gf,point);
@@ -2576,12 +2576,12 @@ bool translate(GFace* gf,MElementOctree* octree,MVertex* vertex,SPoint2 corr,SVe
delta_x = k*size*cos(angle);
delta_y = k*size*sin(angle);
-
+
x1 = x + delta_x;
y1 = y + delta_y;
x2 = x + delta_y;
y2 = y - delta_x;
-
+
if(!inside_domain(octree,x1,y1)){
x1 = x - delta_x;
y1 = y - delta_y;
@@ -2592,9 +2592,9 @@ bool translate(GFace* gf,MElementOctree* octree,MVertex* vertex,SPoint2 corr,SVe
y2 = y + delta_x;
if(!inside_domain(octree,x2,y2)) ok = false;
}
-
+
ok = true; //?
-
+
if(ok){
gp1 = gf->point(x1,y1);
gp2 = gf->point(x2,y2);
@@ -2616,27 +2616,27 @@ bool improved_translate(GFace* gf,MVertex* vertex,SVector3& v1,SVector3& v2){
SVector3 normal;
SVector3 basis_u,basis_v;
Pair<SVector3,SVector3> derivatives;
-
+
reparamMeshVertexOnFace(vertex,gf,point);
x = point.x();
y = point.y();
-
+
angle = backgroundMesh::current()->getAngle(x,y,0.0);
derivatives = gf->firstDer(point);
-
+
s1 = derivatives.first();
s2 = derivatives.second();
normal = crossprod(s1,s2);
-
+
basis_u = s1;
basis_u.normalize();
basis_v = crossprod(normal,basis_u);
basis_v.normalize();
-
+
v1 = basis_u*cos(angle) + basis_v*sin(angle);
v2 = crossprod(v1,normal);
v2.normalize();
-
+
return 1;
}
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index a866812ec0..8b7644cfe7 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -13,13 +13,10 @@
#include "meshGFaceDelaunayInsertion.h"
#include "STensor3.h"
-
-
class GFace;
class GVertex;
class MVertex;
-
class edge_angle {
public :
MVertex *v1, *v2;
@@ -28,10 +25,10 @@ class edge_angle {
bool operator < (const edge_angle &other) const
{
return other.angle < angle;
- }
+ }
};
-typedef std::map<MVertex*, std::vector<MElement*> > v2t_cont;
+typedef std::map<MVertex*, std::vector<MElement*>, MVertexLessThanNum > v2t_cont;
typedef std::map<MEdge, std::pair<MElement*, MElement*>, Less_Edge> e2t_cont;
template <class T> void buildVertexToElement(std::vector<T*> &eles, v2t_cont &adj){
@@ -51,6 +48,7 @@ template <class T> void buildVertexToElement(std::vector<T*> &eles, v2t_cont &ad
}
}
}
+
template <class T> void buildEdgeToElement(std::vector<T*> &eles, e2t_cont &adj);
void buildVertexToTriangle(std::vector<MTriangle*> &, v2t_cont &adj);
@@ -71,18 +69,18 @@ void edgeSwappingLawson(GFace *gf);
enum swapCriterion {SWCR_DEL, SWCR_QUAL, SWCR_NORM, SWCR_CLOSE};
enum splitCriterion {SPCR_CLOSE, SPCR_QUAL, SPCR_ALLWAYS};
-int edgeSwapPass(GFace *gf,
- std::set<MTri3*, compareTri3Ptr> &allTris,
+int edgeSwapPass(GFace *gf,
+ std::set<MTri3*, compareTri3Ptr> &allTris,
const swapCriterion &cr, bidimMeshData &DATA);
void removeFourTrianglesNodes(GFace *gf, bool replace_by_quads);
void removeThreeTrianglesNodes(GFace *gf);
-void buildMeshGenerationDataStructures(GFace *gf,
+void buildMeshGenerationDataStructures(GFace *gf,
std::set<MTri3*, compareTri3Ptr> &AllTris,
bidimMeshData & data);
void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,bidimMeshData &DATA);
void computeEquivalences(GFace *gf,bidimMeshData &DATA);
-void recombineIntoQuads(GFace *gf,
- bool topologicalOpti = true,
+void recombineIntoQuads(GFace *gf,
+ bool topologicalOpti = true,
bool nodeRepositioning = true);
//used for meshGFaceRecombine development
@@ -120,7 +118,7 @@ struct swapquad{
std::sort(v, v + 4);
}
};
-
+
struct RecombineTriangle
{
MElement *t1, *t2;
@@ -132,14 +130,14 @@ struct RecombineTriangle
{
n1 = me.getVertex(0);
n2 = me.getVertex(1);
-
+
if(t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
else if(t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
else if(t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
if(t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
else if(t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
else if(t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
-
+
MQuadrangle q (n1,n3,n2,n4);
angle = q.etaShapeMeasure();
@@ -162,21 +160,21 @@ struct RecombineTriangle
/******************class Temporary******************/
/***************************************************/
-class Temporary{
- private :
- static double w1,w2,w3;
- static std::vector<SVector3> gradients;
- void read_data(std::string);
- static SVector3 compute_normal(MElement*);
- static SVector3 compute_other_vector(MElement*);
- static SVector3 compute_gradient(MElement*);
- public :
- Temporary();
- ~Temporary();
- void quadrilaterize(std::string,double,double,double);
- static double compute_total_cost(double,double);
- static void select_weights(double,double,double);
- static double compute_alignment(const MEdge&,MElement*,MElement*);
+class Temporary{
+ private :
+ static double w1,w2,w3;
+ static std::vector<SVector3> gradients;
+ void read_data(std::string);
+ static SVector3 compute_normal(MElement*);
+ static SVector3 compute_other_vector(MElement*);
+ static SVector3 compute_gradient(MElement*);
+ public :
+ Temporary();
+ ~Temporary();
+ void quadrilaterize(std::string,double,double,double);
+ static double compute_total_cost(double,double);
+ static void select_weights(double,double,double);
+ static double compute_alignment(const MEdge&,MElement*,MElement*);
};
/***************************************************/
@@ -184,4 +182,3 @@ class Temporary{
/***************************************************/
#endif
-
--
GitLab