From 9c5da6b2161cd6df8fab9da077012d5c7d0ff339 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 3 Aug 2010 07:34:22 +0000
Subject: [PATCH] pp+MVertexOctree
---
Geo/CMakeLists.txt | 5 +--
Geo/MElementOctree.h | 5 +++
Geo/MVertex.cpp | 11 ++++-
Geo/MVertexOctree.cpp | 72 +++++++++++++++++++++++++++++++
Geo/MVertexOctree.h | 19 ++++++++
Mesh/meshGFaceExtruded.cpp | 33 ++++++++------
Mesh/meshGRegionExtruded.cpp | 55 +++++++++++++----------
benchmarks/extrude/degenerate.geo | 2 +-
8 files changed, 161 insertions(+), 41 deletions(-)
create mode 100644 Geo/MVertexOctree.cpp
create mode 100644 Geo/MVertexOctree.h
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 70310e4ae6..f9be3537e7 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -30,10 +30,9 @@ set(SRC
findLinks.cpp
SOrientedBoundingBox.cpp
GeomMeshMatcher.cpp
- MVertex.cpp
+ MVertex.cpp MVertexOctree.cpp
MFace.cpp
- MElement.cpp
- MElementOctree.cpp
+ MElement.cpp MElementOctree.cpp
MLine.cpp MTriangle.cpp MQuadrangle.cpp MTetrahedron.cpp
MHexahedron.cpp MPrism.cpp MPyramid.cpp MElementCut.cpp
MZone.cpp MZoneBoundary.cpp
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index 8e883f996f..35606fcd53 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -1,3 +1,8 @@
+// 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 _MELEMENT_OCTREE_
#define _MELEMENT_OCTREE_
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 314c69861a..f67b7374cf 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -222,11 +222,20 @@ 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;
}
- return pos.end();
+ Msg::Warning("Could not find point: returning closest (dist = %g)", dmin);
+ //return pos.end();
+ return itmin;
}
static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> ¶ms)
diff --git a/Geo/MVertexOctree.cpp b/Geo/MVertexOctree.cpp
new file mode 100644
index 0000000000..0dedfcd98f
--- /dev/null
+++ b/Geo/MVertexOctree.cpp
@@ -0,0 +1,72 @@
+// 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"
+
+static void MVertexBB(void *a, double *min, double *max)
+{
+ MVertex *v = (MVertex*)a;
+ double tol = MVertexLessThanLexicographic::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 < MVertexLessThanLexicographic::tolerance);
+}
+
+MVertexOctree::MVertexOctree(GModel *m)
+{
+ SBoundingBox3d bb = m->bounds();
+ 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
new file mode 100644
index 0000000000..e94e5b9104
--- /dev/null
+++ b/Geo/MVertexOctree.h
@@ -0,0 +1,19 @@
+#ifndef _MVERTEX_OCTREE_
+#define _MVERTEX_OCTREE_
+
+class Octree;
+class GModel;
+class MVertex;
+
+class MVertexOctree{
+ private:
+ Octree *_octree;
+ public:
+ MVertexOctree(GModel *);
+ ~MVertexOctree();
+ void insert(MVertex *);
+ void finalize();
+ MVertex *find(double x, double y, double z);
+};
+
+#endif
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 8e7e254c0e..00c8761086 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -12,45 +12,50 @@
#include "Context.h"
#include "GmshMessage.h"
-static void addTriangle(MVertex* v1, MVertex* v2, MVertex* v3, GFace *to, MElement* source) {
+static void addTriangle(MVertex* v1, MVertex* v2, MVertex* v3,
+ GFace *to, MElement* source)
+{
MTriangle* newTri = new MTriangle(v1, v2, v3);
to->triangles.push_back(newTri);
- to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newTri);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem(source, (MElement*)newTri);
}
-static void addQuadrangle(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, GFace *to, MElement* source) {
+static void addQuadrangle(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
+ GFace *to, MElement* source)
+{
MQuadrangle* newQuad = new MQuadrangle(v1, v2, v3, v4);
to->quadrangles.push_back(newQuad);
- to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newQuad);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem(source, (MElement*)newQuad);
}
static void createQuaTri(std::vector<MVertex*> &v, GFace *to,
- std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges,MLine* source)
+ std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges,
+ MLine* source)
{
ExtrudeParams *ep = to->meshAttributes.extrude;
if(v[0] == v[1] || v[1] == v[3])
- addTriangle(v[0], v[3], v[2],to,source);
+ addTriangle(v[0], v[3], v[2], to, source);
else if(v[0] == v[2] || v[2] == v[3])
- addTriangle(v[0], v[1], v[3],to,source);
+ addTriangle(v[0], v[1], v[3], to, source);
else if(v[0] == v[3] || v[1] == v[2])
Msg::Error("Uncoherent extruded quadrangle in surface %d", to->tag());
else{
if(ep->mesh.Recombine){
- addQuadrangle(v[0], v[1], v[3], v[2],to,source);
+ addQuadrangle(v[0], v[1], v[3], v[2], to, source);
}
else if(!constrainedEdges){
- addTriangle(v[0], v[1], v[3],to,source);
- addTriangle(v[0], v[3], v[2],to,source);
+ addTriangle(v[0], v[1], v[3], to, source);
+ addTriangle(v[0], v[3], v[2], to, source);
}
else{
std::pair<MVertex*, MVertex*> p(std::min(v[1], v[2]), std::max(v[1], v[2]));
if(constrainedEdges->count(p)){
- addTriangle(v[2], v[1], v[0],to,source);
- addTriangle(v[2], v[3], v[1],to,source);
+ addTriangle(v[2], v[1], v[0], to, source);
+ addTriangle(v[2], v[3], v[1], to, source);
}
else{
- addTriangle(v[2], v[3], v[0],to,source);
- addTriangle(v[0], v[3], v[1],to,source);
+ addTriangle(v[2], v[3], v[0], to, source);
+ addTriangle(v[0], v[3], v[1], to, source);
}
}
}
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 72dcff796e..f4d658099d 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -17,27 +17,37 @@
#include "Context.h"
#include "GmshMessage.h"
-static void addTetrahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, GRegion *to, MElement* source) {
+static void addTetrahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
+ GRegion *to, MElement* source)
+{
MTetrahedron* newElem = new MTetrahedron(v1, v2, v3, v4);
to->tetrahedra.push_back(newElem);
- to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem(source, (MElement*)newElem);
}
-static void addPyramid(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, MVertex* v5, GRegion *to, MElement* source) {
+static void addPyramid(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
+ MVertex* v5, GRegion *to, MElement* source)
+{
MPyramid* newElem = new MPyramid(v1, v2, v3, v4, v5);
to->pyramids.push_back(newElem);
- to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem(source, (MElement*)newElem);
}
-static void addPrism(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, MVertex* v5, MVertex* v6, GRegion *to, MElement* source) {
+static void addPrism(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
+ MVertex* v5, MVertex* v6, GRegion *to, MElement* source)
+{
MPrism* newElem = new MPrism(v1, v2, v3, v4, v5, v6);
to->prisms.push_back(newElem);
- to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem(source, (MElement*)newElem);
}
-static void addHexahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4, MVertex* v5, MVertex* v6, MVertex* v7, MVertex* v8, GRegion *to, MElement* source) {
+
+static void addHexahedron(MVertex* v1, MVertex* v2, MVertex* v3, MVertex* v4,
+ MVertex* v5, MVertex* v6, MVertex* v7, MVertex* v8,
+ GRegion *to, MElement* source)
+{
MHexahedron* newElem = new MHexahedron(v1, v2, v3, v4, v5, v6, v7, v8);
to->hexahedra.push_back(newElem);
- to->meshAttributes.extrude->elementMap.addExtrudedElem((MElement*)source,(MElement*)newElem);
+ to->meshAttributes.extrude->elementMap.addExtrudedElem(source, (MElement*)newElem);
}
static void createPriPyrTet(std::vector<MVertex*> &v, GRegion *to, MElement* source)
@@ -50,22 +60,22 @@ static void createPriPyrTet(std::vector<MVertex*> &v, GRegion *to, MElement* sou
if(j == 2) {
if(dup[0] == 0 && dup[1] == 1)
- addTetrahedron(v[0], v[1], v[2], v[5],to,source);
+ addTetrahedron(v[0], v[1], v[2], v[5], to, source);
else if(dup[0] == 1 && dup[1] == 2)
- addTetrahedron(v[0], v[1], v[2], v[3],to,source);
+ addTetrahedron(v[0], v[1], v[2], v[3], to, source);
else
- addTetrahedron(v[0], v[1], v[2], v[4],to,source);
+ addTetrahedron(v[0], v[1], v[2], v[4], to, source);
}
else if(j == 1) {
if(dup[0] == 0)
- addPyramid(v[1], v[4], v[5], v[2], v[0],to,source);
+ addPyramid(v[1], v[4], v[5], v[2], v[0], to, source);
else if(dup[0] == 1)
- addPyramid(v[0], v[2], v[5], v[3], v[1],to,source);
+ addPyramid(v[0], v[2], v[5], v[3], v[1], to, source);
else
- addPyramid(v[0], v[1], v[4], v[3], v[2],to,source);
+ addPyramid(v[0], v[1], v[4], v[3], v[2], to, source);
}
else {
- addPrism(v[0], v[1], v[2], v[3], v[4], v[5],to,source);
+ addPrism(v[0], v[1], v[2], v[3], v[4], v[5], to, source);
if(j) Msg::Error("Degenerated prism in extrusion of volume %d", to->tag());
}
}
@@ -80,26 +90,27 @@ static void createHexPri(std::vector<MVertex*> &v, GRegion *to, MElement* source
if(j == 2) {
if(dup[0] == 0 && dup[1] == 1)
- addPrism(v[0], v[3], v[7], v[1], v[2], v[6],to,source);
+ addPrism(v[0], v[3], v[7], v[1], v[2], v[6], to, source);
else if(dup[0] == 1 && dup[1] == 2)
- addPrism(v[0], v[1], v[4], v[3], v[2], v[7],to,source);
+ addPrism(v[0], v[1], v[4], v[3], v[2], v[7], to, source);
else if(dup[0] == 2 && dup[1] == 3)
- addPrism(v[0], v[3], v[4], v[1], v[2], v[5],to,source);
+ addPrism(v[0], v[3], v[4], v[1], v[2], v[5], to, source);
else if(dup[0] == 0 && dup[1] == 3)
- addPrism(v[0], v[1], v[5], v[3], v[2], v[6],to,source);
+ addPrism(v[0], v[1], v[5], v[3], v[2], v[6], to, source);
else
Msg::Error("Uncoherent hexahedron in extrusion of volume %d", to->tag());
}
else {
- addHexahedron(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7],to,source);
+ addHexahedron(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], to, source);
if(j) Msg::Error("Degenerated hexahedron in extrusion of volume %d", to->tag());
}
}
-static void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, GRegion *to, MElement* source)
+static void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4, GRegion *to,
+ MElement* source)
{
if(v1 != v2 && v1 != v3 && v1 != v4 && v2 != v3 && v2 != v4 && v3 != v4)
- addTetrahedron(v1, v2, v3, v4,to,source);
+ addTetrahedron(v1, v2, v3, v4, to, source);
}
static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
diff --git a/benchmarks/extrude/degenerate.geo b/benchmarks/extrude/degenerate.geo
index 8688e37f8d..9e0ceb013a 100644
--- a/benchmarks/extrude/degenerate.geo
+++ b/benchmarks/extrude/degenerate.geo
@@ -63,4 +63,4 @@ Extrude Surface { 21, {0,1,0} , {0,0,0} , Pi/2 } {
};
-Symmetry { 0,1,0,0 } { Duplicata{ Surface{158}; } }
+//Symmetry { 0,1,0,0 } { Duplicata{ Surface{158}; } }
--
GitLab