pp+MVertexOctree

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

Geo/MElementOctree.h

+// 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_
 

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> &params)

Geo/MVertexOctree.cpp

+// 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);
+}

Geo/MVertexOctree.h

+#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

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);
       }
     }
   }

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, 

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}; } }