diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 448766489261a5f6641d747dfe608746050431d2..e63b430dda9f0f10d4401ef575d6331690b2dc0e 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $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();
+}
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 157a84195a0d11a21727277893b8a59ab23cebd0..895e1d12f2e8a73e6e3010eccb8ce514ad039c8c 100644
--- a/Geo/MVertex.h
+++ b/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
diff --git a/Geo/Makefile b/Geo/Makefile
index 32b798047dff8480fa4e911f9e9ad8e8a2d05939..627a7815f44731177e0026948984b34026118766 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $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 \
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 836d395680b669c2f98ad07712e29480b968410a..ae455eb57f295b0a2c29552caba097ec0fdf5a03 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $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 \
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index da58bf00ef66fca4e6b9d0ed72537bd20c17c062..df7171e49c0baa15f48cf8eee24c89ac014a1296 100644
--- a/Mesh/meshGFace.h
+++ b/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
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 759f50b8b7d14c4ec3163adbc4b81eeebd651c4c..b29c958b154e3cd9d69633b04209ed1caec75614 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -1,4 +1,4 @@
-// $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());
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index b59fbb4320cc31fd836bf2c1f1466a0c58dcc2c5..cc9081fb479777035708fbbf115a2bd3eb332d92 100644
--- a/Mesh/meshGRegion.h
+++ b/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
diff --git a/Mesh/meshGRegionCarveHole.cpp b/Mesh/meshGRegionCarveHole.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d66438cc5117dce6004e9837413616c074e63cc
--- /dev/null
+++ b/Mesh/meshGRegionCarveHole.cpp
@@ -0,0 +1,152 @@
+// $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
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 27ee65851d49fdb3c2c1615f6f2a424afad8f0f5..b2e0933c9030819b68f25ea69dad1ad793f94fca 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -1,4 +1,4 @@
-// $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