From 7aaccdc952f94271216de8286fdaa189728454aa Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 28 Aug 2016 14:32:37 +0000
Subject: [PATCH] faster removeDuplicateMeshVertices (single search per vertex
+ direct element modification)
---
Common/CommandLine.cpp | 2 +-
Geo/GModel.cpp | 90 +++++++++++++++++-------------------------
Geo/MHexahedron.h | 5 ++-
Geo/MLine.h | 2 +
Geo/MPrism.h | 5 ++-
Geo/MPyramid.h | 6 +--
Geo/MQuadrangle.h | 20 +++++-----
Geo/MTetrahedron.h | 2 +
Geo/MTriangle.h | 12 +++---
Geo/MVertexRTree.h | 15 +++----
doc/VERSIONS.txt | 3 +-
11 files changed, 79 insertions(+), 83 deletions(-)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 79e5e7c4ff..0083b29fec 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -660,7 +660,7 @@ void GetOptions(int argc, char *argv[])
#endif
// convert mesh to latest binary format
if(GModel::current()->getMeshStatus() > 0){
- CTX::instance()->mesh.mshFileVersion = 2.0;
+ CTX::instance()->mesh.mshFileVersion = 3.0;
CTX::instance()->mesh.binary = 1;
CreateOutputFile(fileName, FORMAT_MSH);
}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 8a55ab05be..04cd1e475c 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1721,79 +1721,59 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
std::vector<GEntity*> entities;
getEntities(entities);
-
- std::vector<MVertex*> vertices;
+ std::map<int, MVertex*> vertices;
+ MVertexRTree pos(eps);
+ std::map<MVertex*,MVertex*> newVertex;
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];
- vertices.push_back(new MVertex(v->x(), v->y(), v->z()));
+ GEntity* ge = entities[i];
+ for(unsigned int j = 0; j < ge->mesh_vertices.size(); j++){
+ MVertex *v = ge->mesh_vertices[j];
+ MVertex *v2 = pos.insert(v);
+ if(v2){
+ newVertex[v] = v2; // v should be removed
+ vertices[v2->getNum()] = v2;
+ }
+ else{
+ vertices[v->getNum()] = v;
+ }
}
}
- MVertexRTree pos(eps);
- int num = pos.insert(vertices);
+ int num = (int)newVertex.size();
Msg::Info("Found %d duplicate vertices ", num);
if(!num){
- for(unsigned int i = 0; i < vertices.size(); i++)
- delete vertices[i];
Msg::Info("No duplicate vertices found");
return 0;
}
- std::map<MVertex*,MVertex*> newVertex;
- std::map<int, std::vector<MElement*> > elements[11];
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::vector<MVertex*> verts;
+ GEntity* ge = entities[i];
+ // clear list of vertices owned by entity
+ ge->mesh_vertices.clear();
+ // replace vertices in element
+ for(unsigned int j = 0; j < ge->getNumMeshElements(); j++){
+ MElement *e = ge->getMeshElement(j);
for(int k = 0; k < e->getNumVertices(); k++){
- MVertex *v = e->getVertex(k);
- MVertex *v2 = pos.find(v->x(), v->y(), v->z());
- if(v2) {
- verts.push_back(v2);
- newVertex[v] = v2;
- }
+ std::map<MVertex*, MVertex*>::iterator it = newVertex.find(e->getVertex(k));
+ if(it != newVertex.end())
+ e->setVertex(k, it->second);
}
- if((int)verts.size() == e->getNumVertices()){
- MElementFactory factory;
- MElement *e2 = factory.create(e->getTypeForMSH(), verts, e->getNum(),
- e->getPartition());
- switch(e2->getType()){
- case TYPE_PNT: elements[0][entities[i]->tag()].push_back(e2); break;
- case TYPE_LIN: elements[1][entities[i]->tag()].push_back(e2); break;
- case TYPE_TRI: elements[2][entities[i]->tag()].push_back(e2); break;
- case TYPE_QUA: elements[3][entities[i]->tag()].push_back(e2); break;
- case TYPE_TET: elements[4][entities[i]->tag()].push_back(e2); break;
- case TYPE_HEX: elements[5][entities[i]->tag()].push_back(e2); break;
- case TYPE_PRI: elements[6][entities[i]->tag()].push_back(e2); break;
- case TYPE_PYR: elements[7][entities[i]->tag()].push_back(e2); break;
- case TYPE_TRIH: elements[8][entities[i]->tag()].push_back(e2); break;
- case TYPE_POLYG: elements[9][entities[i]->tag()].push_back(e2); break;
- case TYPE_POLYH: elements[10][entities[i]->tag()].push_back(e2); break;
- }
- }
- else
- Msg::Error("Could not recreate element %d", e->getNum());
}
- }
-
- // replace vertices in the periodic copies
- for(unsigned int i = 0; i < entities.size(); i++){
- GEntity* ge = entities[i];
+ // replace vertices in periodic copies
std::map<MVertex*,MVertex*>& corrVtcs = ge->correspondingVertices;
std::map<MVertex*,MVertex*>::iterator cIter;
- for (cIter=newVertex.begin();cIter!=newVertex.end();++cIter) {
+ for (cIter = newVertex.begin(); cIter != newVertex.end(); ++cIter) {
MVertex* oldTgt = cIter->first;
MVertex* newTgt = cIter->second;
- std::map<MVertex*,MVertex*>::iterator cvIter = corrVtcs.find(oldTgt);
+ std::map<MVertex*, MVertex*>::iterator cvIter = corrVtcs.find(oldTgt);
if (cvIter != corrVtcs.end()) {
MVertex* src = cvIter->second;
corrVtcs.erase(cvIter);
corrVtcs[newTgt] = src;
}
}
- for (cIter=corrVtcs.begin();cIter!=corrVtcs.end();++cIter) {
+ for (cIter = corrVtcs.begin(); cIter != corrVtcs.end(); ++cIter) {
MVertex* oldSrc = cIter->second;
std::map<MVertex*,MVertex*>::iterator nIter = newVertex.find(oldSrc);
if (nIter != newVertex.end()) {
@@ -1804,14 +1784,18 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
}
}
- for(unsigned int i = 0; i < entities.size(); i++)
- entities[i]->deleteMesh();
-
- for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
- _storeElementsInEntities(elements[i]);
+ destroyMeshCaches();
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertices);
+ // delete unused vertices
+ std::vector<MVertex*> to_delete;
+ for(std::map<MVertex *, MVertex*>::iterator it = newVertex.begin();
+ it != newVertex.end(); it++)
+ to_delete.push_back(it->first);
+ for(unsigned int i = 0; i < to_delete.size(); i++)
+ delete to_delete[i];
+
if(num)
Msg::Info("Removed %d duplicate mesh %s", num, num > 1 ? "vertices" : "vertex");
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index c410956441..7caea46813 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -243,6 +243,7 @@ class MHexahedron20 : public MHexahedron {
virtual int getNumVertices() const { return 20; }
virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 8) _v[num] = v; else _vs[num - 8] = v; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[20] = {0, 8, 1, 11, 2, 13, 3, 9, 10, 12,
@@ -371,7 +372,7 @@ class MHexahedron27 : public MHexahedron {
virtual int getNumVertices() const { return 27; }
virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }
-
+ virtual void setVertex(int num, MVertex *v){ if(num < 8) _v[num] = v; else _vs[num - 8] = v; }
virtual MVertex *getVertexINP(int num)
{
static const int map[27] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 11, 13, 9, 16, 18, 19,
@@ -502,7 +503,7 @@ class MHexahedronN : public MHexahedron {
virtual int getNumVertices() const { return 8 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
virtual const MVertex *getVertex(int num) const { return num < 8 ? _v[num] : _vs[num - 8]; }
-
+ virtual void setVertex(int num, MVertex *v){ if(num < 8) _v[num] = v; else _vs[num - 8] = v; }
virtual int getNumEdgeVertices() const { return 12 * (_order - 1); }
virtual int getNumFaceVertices() const
{
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 0e26cab51e..ad7a9a5722 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -126,6 +126,7 @@ class MLine3 : public MLine {
virtual int getNumVertices() const { return 3; }
virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
virtual const MVertex *getVertex(int num) const{ return num < 2 ? _v[num] : _vs[num - 2]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 2) _v[num] = v; else _vs[num - 2] = v; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[3] = {0, 2, 1};
@@ -183,6 +184,7 @@ class MLineN : public MLine {
virtual int getNumVertices() const { return _vs.size() + 2; }
virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
virtual const MVertex *getVertex(int num) const{ return num < 2 ? _v[num] : _vs[num - 2]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 2) _v[num] = v; else _vs[num - 2] = v; }
virtual int getNumEdgeVertices() const { return _vs.size(); }
virtual int getNumEdgesRep(bool curved);
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n);
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 1ecc03c95f..1d4bbf29f1 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -242,6 +242,7 @@ class MPrism15 : public MPrism {
virtual int getNumVertices() const { return 15; }
virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
virtual const MVertex *getVertex(int num) const { return num < 6 ? _v[num] : _vs[num - 6]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 6) _v[num] = v; else _vs[num - 6] = v; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[15] = {0, 6, 1, 9, 2, 7, 8, 10, 11, 3, 12, 4, 14, 5, 13};
@@ -346,6 +347,7 @@ class MPrism18 : public MPrism {
virtual int getNumVertices() const { return 18; }
virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num - 6]; }
virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num - 6]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 6) _v[num] = v; else _vs[num - 6] = v; }
virtual int getNumEdgeVertices() const { return 9; }
virtual int getNumFaceVertices() const { return 3; }
virtual int getNumEdgesRep(bool curved);
@@ -422,7 +424,8 @@ class MPrismN : public MPrism {
virtual int getPolynomialOrder() const { return _order; }
virtual int getNumVertices() const { return 6+_vs.size(); }
virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num-6]; }
- virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num-6]; }
+ virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num - 6]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 6) _v[num] = v; else _vs[num - 6] = v; }
virtual int getNumEdgeVertices() const { return 9*(_order-1); }
virtual int getNumFaceVertices() const
{
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 05903346f3..57508dbf8d 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -243,10 +243,8 @@ class MPyramidN : public MPyramid {
virtual int getPolynomialOrder() const { return _order; }
virtual int getNumVertices() const { return 5 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
- virtual const MVertex *getVertex(int num) const
- {
- return num < 5 ? _v[num] : _vs[num - 5];
- }
+ virtual const MVertex *getVertex(int num) const { return num < 5 ? _v[num] : _vs[num - 5]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 5) _v[num] = v; else _vs[num - 5] = v; }
virtual int getNumEdgeVertices() const { return 8 * (_order - 1); }
virtual int getNumFaceVertices() const
{
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index f81b261147..4a522f28cd 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -139,10 +139,10 @@ class MQuadrangle : public MElement {
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
}
// reorient the quadrangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
-
+
virtual bool isInside(double u, double v, double w) const
{
double tol = _isInsideTolerance;
@@ -205,6 +205,7 @@ class MQuadrangle8 : public MQuadrangle {
virtual int getNumVertices() const { return 8; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual const MVertex *getVertex(int num) const { return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 4) _v[num] = v; else _vs[num - 4] = v; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[8] = {0, 4, 1, 5, 2, 6, 3, 7};
@@ -250,10 +251,10 @@ class MQuadrangle8 : public MQuadrangle {
}
// reorient the quadrangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
-
+
virtual void getNode(int num, double &u, double &v, double &w) const
{
num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
@@ -298,6 +299,7 @@ class MQuadrangle9 : public MQuadrangle {
virtual int getNumVertices() const { return 9; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual const MVertex *getVertex(int num) const { return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 4) _v[num] = v; else _vs[num - 4] = v; }
virtual MVertex *getVertexDIFF(int num)
{
static const int map[9] = {0, 2, 8, 6, 1, 5, 7, 3, 4};
@@ -339,10 +341,10 @@ class MQuadrangle9 : public MQuadrangle {
}
// reorient the quadrangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
-
+
virtual void getNode(int num, double &u, double &v, double &w) const
{
num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
@@ -461,12 +463,12 @@ class MQuadrangleN : public MQuadrangle {
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
-
+
// reorient the quadrangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
-
+
virtual void getNode(int num, double &u, double &v, double &w) const
{
num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index ec41b1a8cf..66b56bd0e6 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -229,6 +229,7 @@ class MTetrahedron10 : public MTetrahedron {
virtual int getNumVertices() const { return 10; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual const MVertex *getVertex(int num) const { return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 4) _v[num] = v; else _vs[num - 4] = v; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[10] = {0, 4, 1, 5, 2, 6, 7, 9, 8, 3};
@@ -336,6 +337,7 @@ class MTetrahedronN : public MTetrahedron {
virtual int getNumVertices() const { return 4 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual const MVertex *getVertex(int num) const{ return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 4) _v[num] = v; else _vs[num - 4] = v; }
virtual int getNumEdgeVertices() const { return 6 * (_order - 1); }
virtual int getNumFaceVertices() const
{
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index d0318b44c7..f053a68d22 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -125,12 +125,12 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
-
+
// reorient the triangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
-
+
virtual void getNode(int num, double &u, double &v, double &w) const
{
w = 0.;
@@ -200,6 +200,7 @@ class MTriangle6 : public MTriangle {
virtual int getNumVertices() const { return 6; }
virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
virtual const MVertex *getVertex(int num) const { return num < 3 ? _v[num] : _vs[num - 3]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 3) _v[num] = v; else _vs[num - 3] = v; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[6] = {0, 3, 1, 4, 2, 5};
@@ -244,7 +245,7 @@ class MTriangle6 : public MTriangle {
num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
}
// reorient the triangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
};
@@ -287,6 +288,7 @@ class MTriangleN : public MTriangle {
virtual int getNumVertices() const { return 3 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
virtual const MVertex *getVertex(int num) const { return num < 3 ? _v[num] : _vs[num - 3]; }
+ virtual void setVertex(int num, MVertex *v){ if(num < 3) _v[num] = v; else _vs[num - 3] = v; }
virtual int getNumFaceVertices() const
{
if (getIsAssimilatedSerendipity())
@@ -355,7 +357,7 @@ class MTriangleN : public MTriangle {
num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
}
// reorient the triangle to conform with other face
- // orientation computed with MFace based on this face with respect to other
+ // orientation computed with MFace based on this face with respect to other
// in computeCorrespondence
virtual void reorient(int rotation, bool swap);
};
diff --git a/Geo/MVertexRTree.h b/Geo/MVertexRTree.h
index 166eef90cb..effa3ca7d6 100644
--- a/Geo/MVertexRTree.h
+++ b/Geo/MVertexRTree.h
@@ -34,7 +34,7 @@ class MVertexRTree{
_rtree->RemoveAll();
delete _rtree;
}
- int insert(MVertex *v, bool warnIfExists=false)
+ MVertex *insert(MVertex *v, bool warnIfExists=false)
{
MVertex *out;
double _min[3] = {v->x() - _tol, v->y() - _tol, v->z() - _tol};
@@ -43,18 +43,19 @@ class MVertexRTree{
_rtree->Insert(_min, _max, v);
return 0;
}
- else if(warnIfExists){
- Msg::Warning("Vertex %d (%.16g, %.16g, %.16g) already exists in the "
- "mesh with tolerance %g", v->getNum(),
- v->x(), v->y(), v->z(), _tol);
+ else{
+ if(warnIfExists)
+ Msg::Warning("Vertex %d (%.16g, %.16g, %.16g) already exists in the "
+ "mesh with tolerance %g", v->getNum(),
+ v->x(), v->y(), v->z(), _tol);
+ return out;
}
- return 1; // one vertex not inserted
}
int insert(std::vector<MVertex*> &v, bool warnIfExists=false)
{
int num = 0;
for(unsigned int i = 0; i < v.size(); i++)
- num += insert(v[i], warnIfExists);
+ num += (insert(v[i], warnIfExists) ? 1 : 0);
return num; // number of vertices not inserted
}
MVertex *find(double x, double y, double z)
diff --git a/doc/VERSIONS.txt b/doc/VERSIONS.txt
index c9c52ca789..dfad83e3dc 100644
--- a/doc/VERSIONS.txt
+++ b/doc/VERSIONS.txt
@@ -1,5 +1,6 @@
2.13.3: new Tochnog file format export; added ability to remove last command in
-scripts generated interactively.
+scripts generated interactively; small ONELAB usability improvements; faster
+"Coherence Mesh".
2.13.2 (August 18, 2016)): small improvements (scale labels, periodic and
high-order meshes) and bug fixes.
--
GitLab