cleanup previous commit + move "carve hole" routines in separate file

Geo/MVertex.cpp

-// $Id: MVertex.cpp,v 1.13 2007-04-21 19:40:00 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.14 2007-07-31 22:09:11 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -25,6 +25,16 @@
 int MVertex::_globalNum = 0;
 double MVertexLessThanLexicographic::tolerance = 1.e-6;
 
+bool MVertexLessThanLexicographic::operator()(const MVertex *v1, const MVertex *v2) const
+{
+  if(v1->x() - v2->x() >  tolerance) return true;
+  if(v1->x() - v2->x() < -tolerance) return false;
+  if(v1->y() - v2->y() >  tolerance) return true;
+  if(v1->y() - v2->y() < -tolerance) return false;
+  if(v1->z() - v2->z() >  tolerance) return true;
+  return false;
+}
+
 void MVertex::writeMSH(FILE *fp, bool binary, double scalingFactor)
 {
   if(_num < 0) return; // negative vertices are never saved
@@ -133,3 +143,14 @@ void MVertex::writeBDF(FILE *fp, int format, double scalingFactor)
     fprintf(fp, "*N%-6d%-16.9G\n", _num, z1);
   }
 }
+
+std::set<MVertex*, MVertexLessThanLexicographic>::iterator 
+MVertex::linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+{
+  double tol = MVertexLessThanLexicographic::tolerance;
+  for(std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.begin();
+      it != pos.end(); ++it){
+    if(distance(*it) < tol) return it;
+  }
+  return pos.end();
+}

Geo/MVertex.h

@@ -21,10 +21,17 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 #include <stdio.h>
-#include <algorithm>
+#include <set>
 #include "SPoint3.h"
 
 class GEntity;
+class MVertex;
+
+class MVertexLessThanLexicographic{
+ public:
+  static double tolerance;
+  bool operator()(const MVertex *v1, const MVertex *v2) const;
+};
 
 // A mesh vertex.
 class MVertex{
@@ -90,6 +97,10 @@ class MVertex{
     return sqrt(dx * dx + dy * dy + dz * dz);
   }
 
+  // linear coordinate search for the vertex in a set
+  std::set<MVertex*, MVertexLessThanLexicographic>::iterator 
+  linearSearch(std::set<MVertex*, MVertexLessThanLexicographic> &pos);
+
   // Get the data associated with this vertex
   virtual void *getData(){ return 0; }
 
@@ -143,18 +154,4 @@ class MDataFaceVertex : public MFaceVertex{
   virtual void *getData(){ return (void*)&_data; }
 };
 
-class MVertexLessThanLexicographic{
- public:
-  static double tolerance;
-  bool operator()(const MVertex *v1, const MVertex *v2) const
-  {
-    if(v1->x() - v2->x() >  tolerance) return true;
-    if(v1->x() - v2->x() < -tolerance) return false;
-    if(v1->y() - v2->y() >  tolerance) return true;
-    if(v1->y() - v2->y() < -tolerance) return false;
-    if(v1->z() - v2->z() >  tolerance) return true;
-    return false;
-  }
-};
-
 #endif

Geo/Makefile

-# $Id: Makefile,v 1.153 2007-07-31 20:07:38 geuzaine Exp $
+# $Id: Makefile,v 1.154 2007-07-31 22:09:11 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -223,8 +223,7 @@ GModelIO_F.o: GModelIO_F.cpp GModel.h GVertex.h GEntity.h Range.h \
   ExtrudeParams.h ../Common/SmoothData.h GFace.h GEdgeLoop.h Pair.h \
   GRegion.h ../Common/Message.h ../Post/Views.h ../Post/ColorTable.h \
   ../Common/VertexArray.h ../Post/AdaptiveViews.h ../Common/GmshMatrix.h \
-  FFace.h FEdge.h FVertex.h ../Mesh/meshGFace.h ../Geo/MVertex.h \
-  GModelIO_F.h
+  FFace.h FEdge.h FVertex.h ../Mesh/meshGFace.h GModelIO_F.h
 GModelIO_CGNS.o: GModelIO_CGNS.cpp GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h \
   SPoint2.h GEdge.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \

Mesh/Makefile

-# $Id: Makefile,v 1.177 2007-07-31 20:07:38 geuzaine Exp $
+# $Id: Makefile,v 1.178 2007-07-31 22:09:11 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -42,6 +42,7 @@ SRC = Generator.cpp \
         meshGRegionDelaunayInsertion.cpp \
         meshGRegionTransfinite.cpp \
         meshGRegionExtruded.cpp \
+        meshGRegionCarveHole.cpp \
         DivideAndConquer.cpp \
         BackgroundMesh.cpp \
         BoundaryLayer.cpp \
@@ -150,9 +151,9 @@ meshGEdgeExtruded.o: meshGEdgeExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
   ../Common/Message.h
-meshGFace.o: meshGFace.cpp meshGFace.h ../Geo/MVertex.h ../Geo/SPoint3.h \
-  meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
-  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+meshGFace.o: meshGFace.cpp meshGFace.h meshGFaceDelaunayInsertion.h \
+  ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
+  ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Geo/SPoint3.h ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h \
   ../Geo/SVector3.h ../Numeric/Numeric.h ../Common/Context.h \
   ../DataStr/List.h DivideAndConquer.h BackgroundMesh.h ../Geo/GVertex.h \
@@ -170,19 +171,18 @@ meshGFace.o: meshGFace.cpp meshGFace.h ../Geo/MVertex.h ../Geo/SPoint3.h \
   ../Common/OS.h BDS.h ../Post/Views.h ../Post/ColorTable.h \
   ../Post/AdaptiveViews.h ../Common/GmshMatrix.h
 meshGFaceTransfinite.o: meshGFaceTransfinite.cpp meshGFace.h \
-  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/GVertex.h ../Geo/GEntity.h \
-  ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
-  ../Geo/SPoint3.h ../Common/GmshDefines.h ../Geo/MVertex.h \
-  ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
-  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h \
-  ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/MVertex.h ../Geo/MEdge.h \
-  ../Geo/MVertex.h ../Geo/SVector3.h ../Common/Hash.h ../Geo/MFace.h \
-  ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/Numeric.h \
-  ../Common/Context.h ../DataStr/List.h ../Geo/ExtrudeParams.h \
-  ../Common/SmoothData.h ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h \
-  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h \
-  ../Common/Message.h
+  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
+  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h \
+  ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h \
+  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+  ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Geo/GFace.h \
+  ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
+  ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Geo/ExtrudeParams.h ../Common/Message.h
 meshGFaceExtruded.o: meshGFaceExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Common/SmoothData.h ../Numeric/Numeric.h ../Geo/GModel.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
@@ -254,19 +254,19 @@ meshGRegionDelaunayInsertion.o: meshGRegionDelaunayInsertion.cpp \
   ../Geo/GEntity.h ../Geo/MElement.h ../Geo/ExtrudeParams.h \
   ../Geo/SBoundingBox3d.h ../Common/Message.h
 meshGRegionTransfinite.o: meshGRegionTransfinite.cpp meshGFace.h \
-  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/GFace.h ../Geo/GPoint.h \
-  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
-  ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
-  ../Geo/GEntity.h ../Geo/MVertex.h ../Geo/GPoint.h ../Geo/SPoint2.h \
-  ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/MElement.h ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h \
-  ../Geo/SVector3.h ../Numeric/Numeric.h ../Common/Context.h \
-  ../DataStr/List.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
-  ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
-  ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
-  ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Common/Message.h
+  ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+  ../Common/GmshDefines.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/MVertex.h \
+  ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h \
+  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+  ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Geo/MElement.h \
+  ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \
+  ../Geo/ExtrudeParams.h ../Common/Message.h
 meshGRegionExtruded.o: meshGRegionExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Common/SmoothData.h ../Numeric/Numeric.h ../Geo/GModel.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
@@ -282,6 +282,21 @@ meshGRegionExtruded.o: meshGRegionExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
   meshGFace.h meshGRegion.h ../Common/Message.h
+meshGRegionCarveHole.o: meshGRegionCarveHole.cpp ../Geo/GModel.h \
+  ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
+  ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h \
+  ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h \
+  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+  ../Geo/ExtrudeParams.h ../Common/SmoothData.h ../Geo/GFace.h \
+  ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
+  ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
+  ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
+  ../Common/Message.h
 DivideAndConquer.o: DivideAndConquer.cpp ../Common/Gmsh.h \
   ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \

Mesh/meshGFace.h

@@ -22,9 +22,9 @@
 
 #include <vector>
 #include <set>
-#include "MVertex.h"
 
 class GFace;
+class MVertex;
 
 // Create the mesh of the face
 class meshGFace {
@@ -59,7 +59,4 @@ int MeshTransfiniteSurface(GFace *gf);
 int MeshExtrudedSurface(GFace *gf, 
 			std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges=0);
 
-std::set<MVertex*, MVertexLessThanLexicographic>::iterator 
-linearFind(std::set<MVertex*, MVertexLessThanLexicographic> &pos, MVertex *p);
-
 #endif

Mesh/meshGFaceExtruded.cpp

-// $Id: meshGFaceExtruded.cpp,v 1.20 2007-07-31 20:07:38 geuzaine Exp $
+// $Id: meshGFaceExtruded.cpp,v 1.21 2007-07-31 22:09:11 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -60,26 +60,6 @@ void createQuaTri(std::vector<MVertex*> &v, GFace *to,
   }
 }
 
-// FIXME: this is a temporary hack to try to circumvent problems
-// experienced with pos.find(). We need to find a better solution.
-std::set<MVertex*, MVertexLessThanLexicographic>::iterator 
-linearFind(std::set<MVertex*, MVertexLessThanLexicographic> &pos, MVertex *p)
-{
-  double eps = MVertexLessThanLexicographic::tolerance;
-  Msg(INFO, "Trying linear find for point %.16g %.16g %.16g (eps = %.16g)",
-      p->x(), p->y(), p->z(), eps);
-  for(std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.begin();
-      it != pos.end(); ++it){
-    MVertex *v = *it;
-    double dx = v->x() - p->x();
-    double dy = v->y() - p->y();
-    double dz = v->z() - p->z();
-    if(sqrt(dx * dx + dy * dy + dz * dz) < eps) return it;
-  }
-  return pos.end();
-}
-	   
-
 void extrudeMesh(GEdge *from, GFace *to,
 		 std::set<MVertex*, MVertexLessThanLexicographic> &pos,
 		 std::set<std::pair<MVertex*, MVertex*> > *constrainedEdges)
@@ -124,7 +104,10 @@ void extrudeMesh(GEdge *from, GFace *to,
 	for(int p = 0; p < 4; p++){
 	  MVertex tmp(x[p], y[p], z[p], 0, -1);
 	  itp = pos.find(&tmp);
-	  if(itp == pos.end()) itp = linearFind(pos, &tmp);
+	  if(itp == pos.end()){ // FIXME: workaround
+	    Msg(INFO, "Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
+	    itp = tmp.linearSearch(pos);
+	  }
 	  if(itp == pos.end()){
 	    Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
 		tmp.x(), tmp.y(), tmp.z(), to->tag());
@@ -164,7 +147,10 @@ void copyMesh(GFace *from, GFace *to,
       ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
 		  tmp.x(), tmp.y(), tmp.z());
       itp = pos.find(&tmp);
-      if(itp == pos.end()) itp = linearFind(pos, &tmp);
+      if(itp == pos.end()){ // FIXME: workaround
+	Msg(INFO, "Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
+	itp = tmp.linearSearch(pos);
+      }
       if(itp == pos.end()) {
 	Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d",
 	    tmp.x(), tmp.y(), tmp.z(), to->tag());
@@ -182,7 +168,10 @@ void copyMesh(GFace *from, GFace *to,
       ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
 		  tmp.x(), tmp.y(), tmp.z());
       itp = pos.find(&tmp);
-      if(itp == pos.end()) itp = linearFind(pos, &tmp);
+      if(itp == pos.end()){ // FIXME: workaround
+	Msg(INFO, "Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
+	itp = tmp.linearSearch(pos);
+      }
       if(itp == pos.end()) {
 	Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g) in surface %d", 
 	    tmp.x(), tmp.y(), tmp.z(), to->tag());

Mesh/meshGRegion.h

@@ -53,5 +53,6 @@ class deMeshGRegion {
 void MeshDelaunayVolume(std::vector<GRegion*> &delaunay);
 int MeshTransfiniteVolume(GRegion *gr);
 int SubdivideExtrudedMesh(GModel *m);
+void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces);
 
 #endif

Mesh/meshGRegionCarveHole.cpp

+// $Id: meshGRegionCarveHole.cpp,v 1.1 2007-07-31 22:09:11 geuzaine Exp $
+//
+// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+// 
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <set>
+#include "GModel.h"
+#include "Message.h"
+
+#if !defined(HAVE_ANN_)
+
+void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces)
+{
+  Msg(GERROR, "Gmsh must be compiled with ANN support to carve holes in meshes");
+}
+
+#else
+
+#include "ANN/ANN.h"
+
+template <class T>
+void carveHole(std::vector<T*> &elements, double distance, ANNkd_tree *kdtree)
+{
+  // delete all elements that have at least one vertex closer than
+  // 'distance' from the carving surface vertices
+  ANNidxArray index = new ANNidx[1];
+  ANNdistArray dist = new ANNdist[1];
+  std::vector<T*> temp;
+  for(unsigned int i = 0; i < elements.size(); i++){
+    for(int j = 0; j < elements[i]->getNumVertices(); j++){
+      MVertex *v = elements[i]->getVertex(j);
+      double xyz[3] = {v->x(), v->y(), v->z()};
+      kdtree->annkSearch(xyz, 1, index, dist);
+      double d = sqrt(dist[0]);
+      if(d < distance){
+	delete elements[i];
+	break;
+      }
+      else if(j == elements[i]->getNumVertices() - 1){
+        temp.push_back(elements[i]);
+      }
+    }
+  }
+  elements = temp;
+  delete [] index;
+  delete [] dist;
+}
+
+template <class T>
+void addFaces(std::vector<T*> &elements, std::set<MFace, Less_Face> &faces)
+{
+  for(unsigned int i = 0; i < elements.size(); i++){
+    for(int j = 0; j < elements[i]->getNumFaces(); j++){
+      MFace f = elements[i]->getFace(j);
+      std::set<MFace, Less_Face>::iterator it = faces.find(f);
+      if(it == faces.end())
+	faces.insert(f);
+      else
+	faces.erase(it);
+    }
+  }
+}
+
+void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces)
+{
+  Msg(INFO, "Carving hole %d from surface %d at distance %g", num, surfaces[0], distance);
+  GModel *m = gr->model();
+
+  // add all points from carving surfaces into kdtree
+  int numnodes = 0;
+  for(unsigned int i = 0; i < surfaces.size(); i++){
+    GFace *gf = m->faceByTag(surfaces[i]);
+    if(!gf){
+      Msg(GERROR, "Unknown carving surface %d", surfaces[i]);
+      return;
+    }
+    numnodes += gf->mesh_vertices.size();
+  }
+
+  ANNpointArray kdnodes = annAllocPts(numnodes, 4);
+  int k = 0;
+  for(unsigned int i = 0; i < surfaces.size(); i++){
+    GFace *gf = m->faceByTag(surfaces[i]);
+    for(unsigned int j = 0; j < gf->mesh_vertices.size(); j++){
+      kdnodes[k][0] = gf->mesh_vertices[j]->x();
+      kdnodes[k][1] = gf->mesh_vertices[j]->y();
+      kdnodes[k][2] = gf->mesh_vertices[j]->z();
+      k++;
+    }
+  }
+  ANNkd_tree *kdtree = new ANNkd_tree(kdnodes, numnodes, 3);
+
+  // remove the volume elements that are within 'distance' of the
+  // carved surface
+  carveHole(gr->tetrahedra, distance, kdtree);
+  carveHole(gr->hexahedra, distance, kdtree);
+  carveHole(gr->prisms, distance, kdtree);
+  carveHole(gr->pyramids, distance, kdtree);
+
+  delete kdtree;
+  annDeallocPts(kdnodes);
+
+  // TODO: remove any interior elements left inside the carved surface
+  // (could shoot a line from each element's barycenter and count
+  // intersections o see who's inside)
+
+  // generate discrete boundary mesh of the carved hole
+  GFace *gf = m->faceByTag(num);
+  if(!gf) return;
+  std::set<MFace, Less_Face> faces;
+  std::list<GFace*> f = gr->faces();
+  for(std::list<GFace*>::iterator it = f.begin(); it != f.end(); it++){
+    addFaces((*it)->triangles, faces);
+    addFaces((*it)->quadrangles, faces);
+  }
+  addFaces(gr->tetrahedra, faces);
+  addFaces(gr->hexahedra, faces);
+  addFaces(gr->prisms, faces);
+  addFaces(gr->pyramids, faces);
+
+  std::set<MVertex*> verts;
+  for(std::set<MFace, Less_Face>::iterator it = faces.begin(); it != faces.end(); it++){
+    for(int i = 0; i < it->getNumVertices(); i++){
+      it->getVertex(i)->setEntity(gf);
+      verts.insert(it->getVertex(i));
+    }
+    if(it->getNumVertices() == 3)
+      gf->triangles.push_back(new MTriangle(it->getVertex(0), it->getVertex(1),
+					    it->getVertex(2)));
+    else if(it->getNumVertices() == 4)
+      gf->quadrangles.push_back(new MQuadrangle(it->getVertex(0), it->getVertex(1),
+						it->getVertex(2), it->getVertex(3)));
+  }
+}
+
+#endif

Mesh/meshGRegionExtruded.cpp

-// $Id: meshGRegionExtruded.cpp,v 1.18 2007-07-31 20:07:38 geuzaine Exp $
+// $Id: meshGRegionExtruded.cpp,v 1.19 2007-07-31 22:09:11 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -111,7 +111,10 @@ int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
   for(int p = 0; p < 2 * n; p++){
     MVertex tmp(x[p], y[p], z[p], 0, -1);
     itp = pos.find(&tmp);
-    if(itp == pos.end()) itp = linearFind(pos, &tmp);
+    if(itp == pos.end()){ // FIXME: workaround
+      Msg(INFO, "Linear search for (%.16g, %.16g, %.16g)", tmp.x(), tmp.y(), tmp.z());
+      itp = tmp.linearSearch(pos);
+    }
     if(itp == pos.end())
       Msg(GERROR, "Could not find extruded vertex (%.16g, %.16g, %.16g)",
 	  tmp.x(), tmp.y(), tmp.z());
@@ -187,131 +190,6 @@ void insertAllVertices(GRegion *gr,
   }
 }
 
-#if defined(HAVE_ANN_)
-
-#include "ANN/ANN.h"
-
-template <class T>
-void carveHole(std::vector<T*> &elements, double distance, ANNkd_tree *kdtree)
-{
-  // delete all elements that have at least one vertex closer than
-  // 'distance' from the carving surface vertices
-  ANNidxArray index = new ANNidx[1];
-  ANNdistArray dist = new ANNdist[1];
-  std::vector<T*> temp;
-  for(unsigned int i = 0; i < elements.size(); i++){
-    for(int j = 0; j < elements[i]->getNumVertices(); j++){
-      MVertex *v = elements[i]->getVertex(j);
-      double xyz[3] = {v->x(), v->y(), v->z()};
-      kdtree->annkSearch(xyz, 1, index, dist);
-      double d = sqrt(dist[0]);
-      if(d < distance){
-	delete elements[i];
-	break;
-      }
-      else if(j == elements[i]->getNumVertices() - 1){
-        temp.push_back(elements[i]);
-      }
-    }
-  }
-  elements = temp;
-  delete [] index;
-  delete [] dist;
-}
-
-#endif
-
-template <class T>
-void addFaces(std::vector<T*> &elements, std::set<MFace, Less_Face> &faces)
-{
-  for(unsigned int i = 0; i < elements.size(); i++){
-    for(int j = 0; j < elements[i]->getNumFaces(); j++){
-      MFace f = elements[i]->getFace(j);
-      std::set<MFace, Less_Face>::iterator it = faces.find(f);
-      if(it == faces.end())
-	faces.insert(f);
-      else
-	faces.erase(it);
-    }
-  }
-}
-
-void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces)
-{
-#if !defined(HAVE_ANN_)
-  Msg(GERROR, "Gmsh must be compiled with ANN support to carve holes in extruded meshes");
-#else
-  Msg(INFO, "Carving hole %d from surface %d at distance %g", num, surfaces[0], distance);
-  GModel *m = gr->model();
-
-  // add all points from carving surfaces into kdtree
-  int numnodes = 0;
-  for(unsigned int i = 0; i < surfaces.size(); i++){
-    GFace *gf = m->faceByTag(surfaces[i]);
-    if(!gf){
-      Msg(GERROR, "Unknown carving surface %d", surfaces[i]);
-      return;
-    }
-    numnodes += gf->mesh_vertices.size();
-  }
-
-  ANNpointArray kdnodes = annAllocPts(numnodes, 4);
-  int k = 0;
-  for(unsigned int i = 0; i < surfaces.size(); i++){
-    GFace *gf = m->faceByTag(surfaces[i]);
-    for(unsigned int j = 0; j < gf->mesh_vertices.size(); j++){
-      kdnodes[k][0] = gf->mesh_vertices[j]->x();
-      kdnodes[k][1] = gf->mesh_vertices[j]->y();
-      kdnodes[k][2] = gf->mesh_vertices[j]->z();
-      k++;
-    }
-  }
-  ANNkd_tree *kdtree = new ANNkd_tree(kdnodes, numnodes, 3);
-
-  // remove the volume elements that are within 'distance' of the
-  // carved surface
-  carveHole(gr->tetrahedra, distance, kdtree);
-  carveHole(gr->hexahedra, distance, kdtree);
-  carveHole(gr->prisms, distance, kdtree);
-  carveHole(gr->pyramids, distance, kdtree);
-
-  delete kdtree;
-  annDeallocPts(kdnodes);
-
-  // TODO: remove any interior elements left inside the carved surface
-  // (could shoot a line from each element's barycenter and count
-  // intersections o see who's inside)
-
-  // generate discrete boundary mesh of the carved hole
-  GFace *gf = m->faceByTag(num);
-  if(!gf) return;
-  std::set<MFace, Less_Face> faces;
-  std::list<GFace*> f = gr->faces();
-  for(std::list<GFace*>::iterator it = f.begin(); it != f.end(); it++){
-    addFaces((*it)->triangles, faces);
-    addFaces((*it)->quadrangles, faces);
-  }
-  addFaces(gr->tetrahedra, faces);
-  addFaces(gr->hexahedra, faces);
-  addFaces(gr->prisms, faces);
-  addFaces(gr->pyramids, faces);
-
-  std::set<MVertex*> verts;
-  for(std::set<MFace, Less_Face>::iterator it = faces.begin(); it != faces.end(); it++){
-    for(int i = 0; i < it->getNumVertices(); i++){
-      it->getVertex(i)->setEntity(gf);
-      verts.insert(it->getVertex(i));
-    }
-    if(it->getNumVertices() == 3)
-      gf->triangles.push_back(new MTriangle(it->getVertex(0), it->getVertex(1),
-					    it->getVertex(2)));
-    else if(it->getNumVertices() == 4)
-      gf->quadrangles.push_back(new MQuadrangle(it->getVertex(0), it->getVertex(1),
-						it->getVertex(2), it->getVertex(3)));
-  }
-#endif
-}
-
 void meshGRegionExtruded::operator() (GRegion *gr) 
 {  
   if(gr->geomType() == GEntity::DiscreteVolume) return;
@@ -320,7 +198,6 @@ void meshGRegionExtruded::operator() (GRegion *gr)
 
   if(!ep || !ep->mesh.ExtrudeMesh || ep->geo.Mode != EXTRUDED_ENTITY) return;
 
-  // Send a messsage to the GMSH environment
   Msg(STATUS2, "Meshing volume %d (extruded)", gr->tag());
 
   // destroy the mesh if it exists