From 399374aab2e71d00ed7e95a583ba02654e16f1c6 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 23 Aug 2010 16:22:02 +0000
Subject: [PATCH] new MVertexPositionSet using kd-tree
---
Geo/CMakeLists.txt | 2 +-
Geo/GModelIO_Mesh.cpp | 32 +++++++------
Geo/Geo.cpp | 9 ++--
Geo/MVertex.cpp | 17 ++-----
Geo/MVertexOctree.cpp | 76 -------------------------------
Geo/MVertexOctree.h | 26 -----------
Geo/MVertexPositionSet.h | 84 +++++++++++++++++++++++++++++++++++
Mesh/Field.cpp | 27 +++++------
Mesh/meshGRegionCarveHole.cpp | 2 +-
Plugin/NearestNeighbor.cpp | 2 +-
10 files changed, 122 insertions(+), 155 deletions(-)
delete mode 100644 Geo/MVertexOctree.cpp
delete mode 100644 Geo/MVertexOctree.h
create mode 100644 Geo/MVertexPositionSet.h
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index f9be3537e7..b35eb7beb1 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -30,7 +30,7 @@ set(SRC
findLinks.cpp
SOrientedBoundingBox.cpp
GeomMeshMatcher.cpp
- MVertex.cpp MVertexOctree.cpp
+ MVertex.cpp
MFace.cpp
MElement.cpp MElementOctree.cpp
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 207e81cb7c..d5591b259e 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -28,6 +28,7 @@
#include "discreteEdge.h"
#include "discreteFace.h"
#include "discreteRegion.h"
+#include "MVertexPositionSet.h"
void GModel::_storePhysicalTagsInEntities(int dim,
std::map<int, std::map<int, std::string> > &map)
@@ -1036,11 +1037,15 @@ int GModel::readSTL(const std::string &name, double tolerance)
add(face);
}
- // create (unique) vertices and triangles
- double lc = norm(SVector3(bbox.max(), bbox.min()));
- double old_tol = MVertexLessThanLexicographic::tolerance;
- MVertexLessThanLexicographic::tolerance = lc * tolerance;
- std::set<MVertex*, MVertexLessThanLexicographic> vertices;
+ // create triangles using unique vertices
+ double eps = norm(SVector3(bbox.max(), bbox.min())) * tolerance;
+ std::vector<MVertex*> vertices;
+ for(unsigned int i = 0; i < points.size(); i++)
+ for(unsigned int j = 0; j < points[i].size(); j++)
+ vertices.push_back(new MVertex(points[i][j].x(), points[i][j].y(),
+ points[i][j].z(), faces[i]));
+ MVertexPositionSet pos(vertices);
+
for(unsigned int i = 0; i < points.size(); i ++){
for(unsigned int j = 0; j < points[i].size(); j += 3){
MVertex *v[3];
@@ -1048,22 +1053,15 @@ int GModel::readSTL(const std::string &name, double tolerance)
double x = points[i][j + k].x();
double y = points[i][j + k].y();
double z = points[i][j + k].z();
- MVertex w(x, y, z);
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = vertices.find(&w);
- if(it != vertices.end()) {
- v[k] = *it;
- }
- else {
- v[k] = new MVertex(x, y, z, faces[i]);
- vertices.insert(v[k]);
- faces[i]->mesh_vertices.push_back(v[k]);
- }
+ v[k] = pos.find(x, y, z, eps);
+ faces[i]->mesh_vertices.push_back(v[k]);
}
faces[i]->triangles.push_back(new MTriangle(v[0], v[1], v[2]));
}
}
-
- MVertexLessThanLexicographic::tolerance = old_tol;
+
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ if(!vertices[i]->getIndex()) delete vertices[i];
fclose(fp);
return 1;
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index f1e8918901..565aded654 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -2707,10 +2707,11 @@ static void ReplaceDuplicatePoints()
{
// FIXME: This routine is in fact logically wrong (the
// compareTwoPoints function used in the avl tree is not a
- // appropriate comparison function). The fix is simple (use a
- // multi-d tree, e.g., MVertexOctree), but fixing the routine would
- // break backward compatibility with old .geo files. This will be
- // fixed in the new abstract GModel CAD creation routines.
+ // appropriate comparison function). The fix is simple (use a multi
+ // dimensional tree, e.g., MVertexPositionSet), but fixing the
+ // routine would break backward compatibility with old .geo
+ // files. This will be fixed in the new abstract GModel CAD creation
+ // routines.
Vertex *v, **pv, **pv2;
Curve *c;
Surface *s;
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index f67b7374cf..36f18b6f3d 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -221,21 +221,10 @@ void MVertex::writeDIFF(FILE *fp, bool binary, double scalingFactor)
std::set<MVertex*, MVertexLessThanLexicographic>::iterator
MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos)
{
- double tol = MVertexLessThanLexicographic::tolerance;
- double dmin = 1e22;
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator itmin = pos.end();
for(std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.begin();
- it != pos.end(); ++it){
- double d = distance(*it);
- if(d < dmin){
- dmin = d;
- itmin = it;
- }
- if(distance(*it) < tol) return it;
- }
- Msg::Warning("Could not find point: returning closest (dist = %g)", dmin);
- //return pos.end();
- return itmin;
+ it != pos.end(); ++it)
+ if(distance(*it) < MVertexLessThanLexicographic::tolerance) return it;
+ return pos.end();
}
static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> ¶ms)
diff --git a/Geo/MVertexOctree.cpp b/Geo/MVertexOctree.cpp
deleted file mode 100644
index 9dcc936a81..0000000000
--- a/Geo/MVertexOctree.cpp
+++ /dev/null
@@ -1,76 +0,0 @@
-// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#include "Octree.h"
-#include "GModel.h"
-#include "MVertex.h"
-#include "MVertexOctree.h"
-#include "SBoundingBox3d.h"
-
-double MVertexOctree::tolerance = 1.e-6;
-
-static void MVertexBB(void *a, double *min, double *max)
-{
- MVertex *v = (MVertex*)a;
- double tol = MVertexOctree::tolerance;
- min[0] = v->x() - tol;
- max[0] = v->x() + tol;
- min[1] = v->y() - tol;
- max[1] = v->y() + tol;
- min[2] = v->z() - tol;
- max[2] = v->z() + tol;
-}
-
-static void MVertexCentroid(void *a, double *x)
-{
- MVertex *v = (MVertex*)a;
- x[0] = v->x();
- x[1] = v->y();
- x[2] = v->z();
-}
-
-static int MVertexInEle(void *a, double *x)
-{
- MVertex *v = (MVertex*)a;
- double d = sqrt((x[0] - v->x()) * (x[0] - v->x()) +
- (x[1] - v->y()) * (x[1] - v->y()) +
- (x[2] - v->z()) * (x[2] - v->z()));
- return (d < MVertexOctree::tolerance);
-}
-
-MVertexOctree::MVertexOctree(GModel *m, double tol)
-{
- tolerance = tol;
- SBoundingBox3d bb = m->bounds();
- bb *= 1.2;
- double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
- double size[3] = {bb.max().x() - bb.min().x(),
- bb.max().y() - bb.min().y(),
- bb.max().z() - bb.min().z()};
- const int maxElePerBucket = 100; // memory vs. speed trade-off
- _octree = Octree_Create(maxElePerBucket, min, size,
- MVertexBB, MVertexCentroid, MVertexInEle);
-}
-
-MVertexOctree::~MVertexOctree()
-{
- Octree_Delete(_octree);
-}
-
-void MVertexOctree::insert(MVertex *v)
-{
- Octree_Insert(v, _octree);
-}
-
-void MVertexOctree::finalize()
-{
- Octree_Arrange(_octree);
-}
-
-MVertex *MVertexOctree::find(double x, double y, double z)
-{
- double P[3] = {x, y, z};
- return (MVertex*)Octree_Search(P, _octree);
-}
diff --git a/Geo/MVertexOctree.h b/Geo/MVertexOctree.h
deleted file mode 100644
index 79a7806eb6..0000000000
--- a/Geo/MVertexOctree.h
+++ /dev/null
@@ -1,26 +0,0 @@
-// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
-//
-// See the LICENSE.txt file for license information. Please report all
-// bugs and problems to <gmsh@geuz.org>.
-
-#ifndef _MVERTEX_OCTREE_
-#define _MVERTEX_OCTREE_
-
-class Octree;
-class GModel;
-class MVertex;
-
-class MVertexOctree{
- private:
- Octree *_octree;
- public:
- static double tolerance;
- public:
- MVertexOctree(GModel *, double);
- ~MVertexOctree();
- void insert(MVertex *);
- void finalize();
- MVertex *find(double x, double y, double z);
-};
-
-#endif
diff --git a/Geo/MVertexPositionSet.h b/Geo/MVertexPositionSet.h
new file mode 100644
index 0000000000..c3ee5e076f
--- /dev/null
+++ b/Geo/MVertexPositionSet.h
@@ -0,0 +1,84 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _MVERTEX_POSITION_SET_
+#define _MVERTEX_POSITION_SET_
+
+#include <vector>
+#include "GmshMessage.h"
+#include "MVertex.h"
+
+#if defined(HAVE_ANN)
+#include "ANN/ANN.h"
+
+// Stores MVertices in a kd-tree so we can query unique vertices (up
+// to a tolerance d). The constructor tags all the vertices with 0;
+// find() tags the returned vertex with -1; if no negatively-tagged
+// vertex exists, find() returns the closest vertex up to the
+// prescribed tolerance.
+class MVertexPositionSet{
+ private:
+ ANNkd_tree *_kdtree;
+ ANNpointArray _zeronodes;
+ ANNidxArray _index;
+ ANNdistArray _dist;
+ int _maxDuplicates;
+ std::vector<MVertex*> &_vertices;
+ public:
+ MVertexPositionSet(std::vector<MVertex*> &vertices, int maxDuplicates=10)
+ : _vertices(vertices), _maxDuplicates(maxDuplicates)
+ {
+ int totpoints = vertices.size();
+ _zeronodes = annAllocPts(totpoints, 3);
+ for(int i = 0; i < totpoints; i++){
+ vertices[i]->setIndex(0);
+ _zeronodes[i][0] = vertices[i]->x();
+ _zeronodes[i][1] = vertices[i]->y();
+ _zeronodes[i][2] = vertices[i]->z();
+ }
+ _kdtree = new ANNkd_tree(_zeronodes, totpoints, 3);
+ _index = new ANNidx[_maxDuplicates];
+ _dist = new ANNdist[_maxDuplicates];
+ }
+ ~MVertexPositionSet()
+ {
+ delete _kdtree;
+ annDeallocPts(_zeronodes);
+ delete [] _index;
+ delete [] _dist;
+ }
+ MVertex *find(double x, double y, double z, double tolerance)
+ {
+ double xyz[3] = {x, y, z};
+ _kdtree->annkSearch(xyz, _maxDuplicates, _index, _dist);
+ for(int i = 0; i < _maxDuplicates; i++){
+ if(_index[i] >= 0 && _vertices[_index[i]]->getIndex() < 0 && sqrt(_dist[i]) < tolerance)
+ return _vertices[_index[i]];
+ }
+ if(_index[0] >= 0 && sqrt(_dist[0]) < tolerance){
+ _vertices[_index[0]]->setIndex(-1);
+ return _vertices[_index[0]];
+ }
+ Msg::Error("Could not find point (%g,%g,%g)", x, y, z);
+ return 0;
+ }
+};
+
+#else
+
+class MVertexPositionSet{
+ public:
+ MVertexPositionSet(std::vector<MVertex*> &vertices, int maxDuplicates=10){}
+ ~MVertexPositionSet(){}
+ MVertex *find(double x, double y, double z, double tolerance)
+ {
+ Msg::Error("Gmsh must be compiled with ANN to use MVertexPositionSet");
+ return 0;
+ }
+};
+
+#endif
+
+#endif
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 5a9a5d9cfe..468f166012 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -665,7 +665,6 @@ class GradientField : public Field
}
};
-
class CurvatureField : public Field
{
int iField;
@@ -937,7 +936,8 @@ class MathEvalExpressionAniso
variables[1] = "y";
variables[2] = "z";
i = 3;
- for(std::set<int>::iterator it = _fields[iFunction].begin(); it != _fields[iFunction].end(); it++){
+ for(std::set<int>::iterator it = _fields[iFunction].begin();
+ it != _fields[iFunction].end(); it++){
std::ostringstream sstream;
sstream << "F" << *it;
variables[i++] = sstream.str();
@@ -962,7 +962,8 @@ class MathEvalExpressionAniso
values[1] = y;
values[2] = z;
int i = 3;
- for(std::set<int>::iterator it = _fields[iFunction].begin(); it != _fields[iFunction].end(); it++){
+ for(std::set<int>::iterator it = _fields[iFunction].begin();
+ it != _fields[iFunction].end(); it++){
Field *field = GModel::current()->getFields()->get(*it);
values[i++] = field ? (*field)(x, y, z) : MAX_LC;
}
@@ -1008,7 +1009,6 @@ class MathEvalField : public Field
}
};
-
class MathEvalFieldAniso : public Field
{
MathEvalExpressionAniso expr;
@@ -1209,7 +1209,6 @@ class MinAnisoField : public Field
}
metr = v;
}
-
double operator() (double x, double y, double z, GEntity *ge=0)
{
double v = MAX_LC;
@@ -1225,7 +1224,6 @@ class MinAnisoField : public Field
}
};
-
class MinField : public Field
{
std::list<int> idlist;
@@ -1375,10 +1373,12 @@ class AttractorField : public Field
"is replaced by NNodesByEdge equidistant nodes and the distance from those "
"nodes is computed.";
}
- std::pair<AttractorInfo,SPoint3> getAttractorInfo() const {
- return std::make_pair(_infos[index[0]], SPoint3 (zeronodes[index[0]][0],zeronodes[index[0]][1],zeronodes[index[0]][2] ) );
+ std::pair<AttractorInfo,SPoint3> getAttractorInfo() const
+ {
+ return std::make_pair(_infos[index[0]], SPoint3(zeronodes[index[0]][0],
+ zeronodes[index[0]][1],
+ zeronodes[index[0]][2]));
}
-
virtual double operator() (double X, double Y, double Z, GEntity *ge=0)
{
if(update_needed) {
@@ -1389,7 +1389,7 @@ class AttractorField : public Field
int totpoints = nodes_id.size() + n_nodes_by_edge * edges_id.size() +
n_nodes_by_edge * n_nodes_by_edge * faces_id.size();
if(totpoints){
- zeronodes = annAllocPts(totpoints, 4);
+ zeronodes = annAllocPts(totpoints, 3);
_infos.resize(totpoints);
}
int k = 0;
@@ -1519,8 +1519,8 @@ public:
options["hfar"] = new FieldOptionDouble
(hfar, "Element size far from the wall");
}
-
- virtual double operator() (double x, double y, double z, GEntity *ge=0){
+ virtual double operator() (double x, double y, double z, GEntity *ge=0)
+ {
Field *field = GModel::current()->getFields()->get(iField);
if(!field || iField == id) {
return MAX_LC;
@@ -1530,7 +1530,6 @@ public:
// double lc = hwall * pow (ratio, dist / hwall);
return std::min (hfar,lc);
}
-
virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
{
Field *field = GModel::current()->getFields()->get(iField);
@@ -1637,8 +1636,6 @@ public:
};
#endif
-
-
template<class F> class FieldFactoryT : public FieldFactory {
public:
Field * operator()() { return new F; }
diff --git a/Mesh/meshGRegionCarveHole.cpp b/Mesh/meshGRegionCarveHole.cpp
index 837ed08ffe..ffe26fa699 100644
--- a/Mesh/meshGRegionCarveHole.cpp
+++ b/Mesh/meshGRegionCarveHole.cpp
@@ -84,7 +84,7 @@ void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces
numnodes += gf->mesh_vertices.size();
}
- ANNpointArray kdnodes = annAllocPts(numnodes, 4);
+ ANNpointArray kdnodes = annAllocPts(numnodes, 3);
int k = 0;
for(unsigned int i = 0; i < surfaces.size(); i++){
GFace *gf = m->getFaceByTag(surfaces[i]);
diff --git a/Plugin/NearestNeighbor.cpp b/Plugin/NearestNeighbor.cpp
index db789e314d..1b598d9a61 100644
--- a/Plugin/NearestNeighbor.cpp
+++ b/Plugin/NearestNeighbor.cpp
@@ -55,7 +55,7 @@ PView *GMSH_NearestNeighborPlugin::execute(PView *v)
}
#if defined(HAVE_ANN)
- ANNpointArray zeronodes = annAllocPts(totpoints, 4);
+ ANNpointArray zeronodes = annAllocPts(totpoints, 3);
int k = 0, step = 0;
for(int ent = 0; ent < data1->getNumEntities(step); ent++){
for(int ele = 0; ele < data1->getNumElements(step, ent); ele++){
--
GitLab