diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index f252cb603ec4c16c3b107e793c17a1b591dbc270..04465bd228b91b4fe6f6fa0300ba878eb82b7971 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -28,7 +28,7 @@
#include "Context.h"
#include "OS.h"
#include "GEdgeLoop.h"
-
+#include "MVertexPositionSet.h"
#include "OpenFile.h"
#include "CreateFile.h"
@@ -956,12 +956,15 @@ void GModel::_associateEntityWithMeshVertices()
void GModel::_storeVerticesInEntities(std::map<int, MVertex*> &vertices)
{
- std::map<int, MVertex*>::const_iterator it = vertices.begin();
+ std::map<int, MVertex*>::iterator it = vertices.begin();
for(; it != vertices.end(); ++it){
MVertex *v = it->second;
GEntity *ge = v->onWhat();
if(ge) ge->mesh_vertices.push_back(v);
- else delete v; // we delete all unused vertices
+ else{
+ delete v; // we delete all unused vertices
+ it->second = 0;
+ }
}
}
@@ -972,7 +975,10 @@ void GModel::_storeVerticesInEntities(std::vector<MVertex*> &vertices)
if(v){ // the vector is allowed to have null entries
GEntity *ge = v->onWhat();
if(ge) ge->mesh_vertices.push_back(v);
- else delete v; // we delete all unused vertices
+ else{
+ delete v; // we delete all unused vertices
+ vertices[i] = 0;
+ }
}
}
}
@@ -986,65 +992,51 @@ void GModel::checkMeshCoherence(double tolerance)
SBoundingBox3d bbox = bounds();
double lc = bbox.empty() ? 1. : norm(SVector3(bbox.max(), bbox.min()));
+ double eps = lc * tolerance;
std::vector<GEntity*> entities;
getEntities(entities);
// check for duplicate mesh vertices
{
- double old_tol = MVertexLessThanLexicographic::tolerance;
- MVertexLessThanLexicographic::tolerance = tolerance * lc;
- std::set<MVertex*, MVertexLessThanLexicographic> pos;
+ std::vector<MVertex*> vertices;
+ for(unsigned int i = 0; i < entities.size(); i++)
+ vertices.insert(vertices.end(), entities[i]->mesh_vertices.begin(),
+ entities[i]->mesh_vertices.end());
+ MVertexPositionSet pos(vertices);
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ pos.find(vertices[i]->x(), vertices[i]->y(), vertices[i]->z(), eps);
int num = 0;
- for(unsigned int i = 0; i < entities.size(); i++){
- for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
- MVertex *v = entities[i]->mesh_vertices[j];
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.find(v);
- if(it == pos.end()){
- pos.insert(v);
- }
- else{
- Msg::Info("Vertices %d and %d have identical position (%g, %g, %g)",
- (*it)->getNum(), v->getNum(), v->x(), v->y(), v->z());
- num++;
- }
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ if(!vertices[i]->getIndex()){
+ Msg::Info("Duplicate vertex at (%.16g,%.16g,%.16g)",
+ vertices[i]->x(), vertices[i]->y(), vertices[i]->z());
+ num++;
}
- }
- if(num) Msg::Warning("%d duplicate vertices", num);
- MVertexLessThanLexicographic::tolerance = old_tol;
+ if(num) Msg::Warning("%d duplicate vert%s", num, num > 1 ? "ices" : "ex");
}
// check for duplicate elements
{
- double old_tol = MElementLessThanLexicographic::tolerance;
- MElementLessThanLexicographic::tolerance = tolerance * lc;
- std::set<MElement*, MElementLessThanLexicographic> pos;
- int num = 0;
- for(unsigned int i = 0; i < entities.size(); i++){
+ std::vector<MVertex*> vertices;
+ for(unsigned int i = 0; i < entities.size(); i++)
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- std::set<MElement*, MElementLessThanLexicographic>::iterator it = pos.find(e);
- if(it == pos.end()){
- pos.insert(e);
- }
- else{
- std::ostringstream sstream;
- sstream << "Element " << e->getNum() << " [ ";
- for (int k = 0; k < e->getNumVertices(); k++)
- sstream << e->getVertex(k)->getNum() << " ";
- sstream << "] on entity " << entities[i]->tag()
- << " has same barycenter as element " << (*it)->getNum()
- << " [ ";
- for (int k = 0; k < (*it)->getNumVertices(); k++)
- sstream << (*it)->getVertex(k)->getNum() << " ";
- sstream << "]";
- Msg::Info("%s", sstream.str().c_str());
- num++;
- }
+ SPoint3 p = entities[i]->getMeshElement(j)->barycenter();
+ vertices.push_back(new MVertex(p.x(), p.y(), p.z()));
}
+ MVertexPositionSet pos(vertices);
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ pos.find(vertices[i]->x(), vertices[i]->y(), vertices[i]->z(), eps);
+ int num = 0;
+ for(unsigned int i = 0; i < vertices.size(); i++){
+ if(!vertices[i]->getIndex()){
+ Msg::Info("Duplicate element with barycenter (%.16g,%.16g,%.16g)",
+ vertices[i]->x(), vertices[i]->y(), vertices[i]->z());
+ num++;
+ }
+ delete vertices[i];
}
- if(num) Msg::Warning("%d duplicate elements", num);
- MElementLessThanLexicographic::tolerance = old_tol;
+ if(num) Msg::Warning("%d duplicate element%s", num, num > 1 ? "s" : "");
}
}
@@ -1052,29 +1044,27 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
{
SBoundingBox3d bbox = bounds();
double lc = bbox.empty() ? 1. : norm(SVector3(bbox.max(), bbox.min()));
+ double eps = lc * tolerance;
std::vector<GEntity*> entities;
getEntities(entities);
- double old_tol = MVertexLessThanLexicographic::tolerance;
- MVertexLessThanLexicographic::tolerance = tolerance * lc;
- std::set<MVertex*, MVertexLessThanLexicographic> pos;
-
- for(unsigned int i = 0; i < entities.size(); i++){
+ std::vector<MVertex*> vertices;
+ for(unsigned int i = 0; i < entities.size(); i++)
for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
MVertex *v = entities[i]->mesh_vertices[j];
- MVertex w(v->x(), v->y(), v->z());
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.find(&w);
- if(it == pos.end())
- pos.insert(new MVertex(v->x(), v->y(), v->z()));
+ vertices.push_back(new MVertex(v->x(), v->y(), v->z()));
}
- }
+ MVertexPositionSet pos(vertices);
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ pos.find(vertices[i]->x(), vertices[i]->y(), vertices[i]->z(), eps);
+ int num = 0;
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ if(!vertices[i]->getIndex()) num++;
- int diff = getNumMeshVertices() - pos.size();
- if(!diff){
- for(std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.begin();
- it != pos.end(); it++)
- delete *it;
+ if(!num){
+ for(unsigned int i = 0; i < vertices.size(); i++)
+ delete vertices[i];
Msg::Info("No duplicate vertices found");
return 0;
}
@@ -1086,17 +1076,8 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
std::vector<MVertex*> verts;
for(int k = 0; k < e->getNumVertices(); k++){
MVertex *v = e->getVertex(k);
- std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.find(v);
- if(it == pos.end()){
- Msg::Info("Linear search for (%.16g, %.16g, %.16g)", v->x(), v->y(), v->z());
- it = v->linearSearch(pos);
- if(it == pos.end()){
- Msg::Error("Could not find unique vertex (%.16g, %.16g, %.16g)",
- v->x(), v->y(), v->z());
- break;
- }
- }
- verts.push_back(*it);
+ MVertex *v2 = pos.find(v->x(), v->y(), v->z(), eps);
+ if(v2) verts.push_back(v2);
}
if((int)verts.size() == e->getNumVertices()){
MElementFactory factory;
@@ -1115,27 +1096,22 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
case TYPE_POLYH: elements[9][entities[i]->tag()].push_back(e2); break;
}
}
+ else
+ Msg::Error("Could not recreate element %d", e->getNum());
}
}
for(unsigned int i = 0; i < entities.size(); i++)
entities[i]->deleteMesh();
- std::vector<MVertex*> vertices;
- for(std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.begin();
- it != pos.end(); it++)
- vertices.push_back(*it);
-
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
_storeElementsInEntities(elements[i]);
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertices);
- MVertexLessThanLexicographic::tolerance = old_tol;
-
- Msg::Info("Removed %d duplicate mesh %s", diff, diff > 1 ? "vertices" : "vertex");
+ Msg::Info("Removed %d duplicate mesh %s", num, num > 1 ? "vertices" : "vertex");
- return diff;
+ return num;
}
void GModel::createTopologyFromMesh()
@@ -1652,13 +1628,13 @@ GEntity *GModel::extrude(GEntity *e, std::vector<double> p1, std::vector<double>
return _factory->extrude(this, e, p1, p2);
return 0;
}
+
GEntity *GModel::addPipe(GEntity *e, std::vector<GEdge *> edges){
if(_factory)
return _factory->addPipe(this,e,edges);
return 0;
}
-
GEntity *GModel::addSphere(double cx, double cy, double cz, double radius)
{
if(_factory) return _factory->addSphere(this, cx, cy, cz, radius);
@@ -1910,10 +1886,9 @@ static void glueFacesInRegions(GModel *model,
std::list<GFace*> old = r->faces(), fnew;
for (std::list<GFace*>::iterator fit = old.begin(); fit != old.end(); fit++){
std::map<GFace*, GFace*>::iterator itR = Duplicates2Unique.find(*fit);
- if (itR == Duplicates2Unique.end())
- {
- Msg::Fatal("error in the gluing process");
- }
+ if (itR == Duplicates2Unique.end()){
+ Msg::Fatal("Error in the gluing process");
+ }
GFace *temp = itR->second;;
fnew.push_back(temp);
if (temp != *fit){
@@ -1948,6 +1923,7 @@ void GModel::glue(double eps)
glueFacesInRegions(this, Unique2Duplicates, Duplicates2Unique);
}
}
+
void GModel::insertRegion(GRegion *r)
{
regions.insert(r);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index d5591b259e61433ea546e8db22a3b8d0ce07d0da..138dd33a0200da3978c192792b5b4d85a5339e45 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1043,7 +1043,7 @@ int GModel::readSTL(const std::string &name, double tolerance)
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]));
+ points[i][j].z()));
MVertexPositionSet pos(vertices);
for(unsigned int i = 0; i < points.size(); i ++){
@@ -1054,14 +1054,13 @@ int GModel::readSTL(const std::string &name, double tolerance)
double y = points[i][j + k].y();
double z = points[i][j + k].z();
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]));
}
}
- for(unsigned int i = 0; i < vertices.size(); i++)
- if(!vertices[i]->getIndex()) delete vertices[i];
+ _associateEntityWithMeshVertices();
+ _storeVerticesInEntities(vertices); // will delete unused vertices
fclose(fp);
return 1;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 3076ad52c8ecc496694879a5dd2fac1d284a4b26..6bea2b0079bfc37d6b8be3893a9905b58ec79ce1 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -27,7 +27,6 @@
int MElement::_globalNum = 0;
double MElement::_isInsideTolerance = 1.e-6;
-double MElementLessThanLexicographic::tolerance = 1.e-6;
MElement::MElement(int num, int part) : _visible(1)
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3131ee236a41c6f04408e57a1c855ba8107a0c64..7f0992c0af1877fe9c8420f19cc4d8e92b6e3360 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -309,22 +309,6 @@ class MElement
std::map<MElement*, MElement*> &newDomains);
};
-class MElementLessThanLexicographic{
- public:
- static double tolerance;
- bool operator()(MElement *e1, MElement *e2) const
- {
- SPoint3 b1 = e1->barycenter();
- SPoint3 b2 = e2->barycenter();
- if(b1.x() - b2.x() > tolerance) return true;
- if(b1.x() - b2.x() < -tolerance) return false;
- if(b1.y() - b2.y() > tolerance) return true;
- if(b1.y() - b2.y() < -tolerance) return false;
- if(b1.z() - b2.z() > tolerance) return true;
- return false;
- }
-};
-
class MElementFactory{
public:
MElement *create(int type, std::vector<MVertex*> &v, int num=0, int part=0,
diff --git a/Geo/MVertexPositionSet.h b/Geo/MVertexPositionSet.h
index c3ee5e076f5c1d20da396274615d7ae867ccf8b4..91770e913fad2de7ae7ef63e0d1766ab60e5167e 100644
--- a/Geo/MVertexPositionSet.h
+++ b/Geo/MVertexPositionSet.h
@@ -14,10 +14,10 @@
#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.
+// to a prescribed tolerance). 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;
@@ -54,14 +54,16 @@ class MVertexPositionSet{
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)
+ if(_index[i] >= 0 && sqrt(_dist[i]) < tolerance &&
+ _vertices[_index[i]]->getIndex() < 0)
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);
+ Msg::Error("Could not find vertex (%g,%g,%g) (tol %g)",
+ x, y, z, tolerance);
return 0;
}
};
diff --git a/benchmarks/2d/wing-splines.geo b/benchmarks/2d/wing-splines.geo
index a6eeab357b6fdb9c543f084781a70a7617c231cc..0aea8841a1e4b89961e10d51b4c16029446fc027 100644
--- a/benchmarks/2d/wing-splines.geo
+++ b/benchmarks/2d/wing-splines.geo
@@ -1,7 +1,7 @@
scale = .01 ;
-lc_wing = 0.0001 * scale ;
+lc_wing = 0.001 * scale ;
lc_box = 10 * scale ;
Point(3895) = {1.177410e-02*scale,-2.768003e-03*scale,0,lc_wing};
diff --git a/benchmarks/stl/car_body_with_antenna.geo b/benchmarks/stl/car_body_with_antenna.geo
index 211b30e07ce8ebc8d831b9b5fffa132c09a23e14..7e41369609856bc69bb92ef6fc4df776b9d3cdaa 100644
--- a/benchmarks/stl/car_body_with_antenna.geo
+++ b/benchmarks/stl/car_body_with_antenna.geo
@@ -10,6 +10,8 @@ Point(p2) = { -0.000210105 , 0.000769998 , 1.30307 + ll , lc1 };
l1 = newl;
Line(l1) = { p1 , p2 };
+
Mesh 2;
Coherence Mesh;
Save "car_body_with_antenna.msh";
+