diff --git a/Fltk/Makefile b/Fltk/Makefile
index 98fe1f830f1e9e72cb130210cdb34f2754ef914f..ba29d76f8e3d0a174875cb5261e9e063c6b8f4e6 100644
--- a/Fltk/Makefile
+++ b/Fltk/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.71 2005-07-03 08:02:24 geuzaine Exp $
+# $Id: Makefile,v 1.72 2005-08-19 14:07:31 remacle Exp $
 #
 # Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 #
@@ -23,7 +23,7 @@ include ../variables
 
 LIB     = ../lib/libGmshFltk.a
 INCLUDE = -I../Common -I../DataStr -I../Graphics -I../Geo -I../Mesh\
-          -I../Numeric -I../Parser -I../Fltk -I../Plugin -I../utils/solvers
+          -I../Numeric -I../Parser -I../Fltk -I../Plugin -I../utils/solvers -I../ANN/include
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE}
 
 SRC = Main.cpp \
diff --git a/Graphics/Makefile b/Graphics/Makefile
index 06644cc247f6f7d1281777c724aabd2992da6281..57db350a9afa3c76b16032862cfe486cef419a53 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.71 2005-06-03 17:32:29 geuzaine Exp $
+# $Id: Makefile,v 1.72 2005-08-19 14:07:32 remacle Exp $
 #
 # Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 #
@@ -23,7 +23,7 @@ include ../variables
 
 LIB     = ../lib/libGmshGraphics.a
 INCLUDE = -I../Common -I../DataStr -I../Geo -I../Graphics\
-          -I../Fltk -I../Mesh -I../Numeric -I../Parser -I../Plugin
+          -I../Fltk -I../Mesh -I../Numeric -I../Parser -I../Plugin -I../ANN/include
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE}
 
 SRC = Draw.cpp \
diff --git a/Mesh/2D_Mesh.cpp b/Mesh/2D_Mesh.cpp
index cb520c0f371b9a1fdee6297a32a49c07ba98a83e..2fcc67e0e7c79ae048f3ca6175bab769b06d2af7 100644
--- a/Mesh/2D_Mesh.cpp
+++ b/Mesh/2D_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_Mesh.cpp,v 1.78 2005-06-10 20:59:15 geuzaine Exp $
+// $Id: 2D_Mesh.cpp,v 1.79 2005-08-19 14:07:32 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -850,12 +850,18 @@ void Maillage_Surface(void *data, void *dum)
 
   Msg(STATUS3, "Meshing surface %d", s->Num);
 
-  if(MeshDiscreteSurface(s)){
+  int tag = MeshDiscreteSurface(s);
+
+  if(tag == 1){
     Tree_Action(THEM->Points, PutVertex_OnSurf);
     Tree_Action(s->Vertices, PutVertex_OnSurf);
     Tree_Action(s->Vertices, Add_In_Mesh);
     return;
   }
+  else if (tag == 2)
+    {
+      return;
+    }
 
   THESUPPORT = s->Support;
   THESURFACE = s;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 8929b7a6de2613e659ccadb2275e7fac0adfa92a..4e2bd936d05ee325061b748b040c6be38b3f8f5f 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -12,6 +12,7 @@
   
 
 */
+
 double BDS_Quadric:: normalCurv ( double x, double y, double z ) const
 {
     // K = div n = div ( Grad(f) / ||Grad(f)||)
@@ -89,6 +90,67 @@ void BDS_Quadric::projection ( double xa, double ya, double za,
    }
 }
 
+void BDS_GeomEntity::getClosestTriangles ( double x, double y, double z, 
+					   std::list<BDS_Triangle*> &l )
+  
+{
+#ifdef HAVE_ANN_
+  l.clear();
+
+  if (kdTree == 0)
+    {
+      l = t;
+      return;
+    }
+
+  queryPt[0] = x;
+  queryPt[1] = y;
+  queryPt[2] = z;
+  double eps;
+  kdTree->annkSearch(                                             // search
+		     queryPt,                                     // query point
+		     1,                                           // number of near neighbors
+		     nnIdx,                                       // nearest neighbors (returned)
+		     dists,                                       // distance (returned)
+		     eps);                                        // error bound
+  if (nnIdx[0] < sP.size())
+    {
+      sP[nnIdx[0]]->getTriangles (l);
+    }
+  else if (nnIdx[0] <  sP.size() + sE.size() )
+    {
+
+      BDS_Edge *e = sE[nnIdx[0]- sP.size()];
+      std::list<BDS_Triangle *> l1;
+      //e->p1->getTriangles (l);
+      //e->p2->getTriangles (l1);
+      //l.insert(l.begin(),l1.begin(),l1.end());
+            for (int i=0;i<e->numfaces();i++)
+      	{
+      	  l.push_back(e->faces(i));
+      	}
+    }
+  else
+    {
+      std::list<BDS_Triangle *> l1,l2;
+      BDS_Point *pts[3];
+      sT[nnIdx[0] - sP.size() - sE.size()]->getNodes (pts); 
+      pts[0]->getTriangles (l);
+      pts[0]->getTriangles (l1);
+      pts[0]->getTriangles (l2);
+      l.insert(l.begin(),l1.begin(),l1.end());
+      l.insert(l.begin(),l2.begin(),l2.end());
+      
+      //      l.push_back( sT[nnIdx[0] - sP.size() - sE.size()]);
+    } 
+#else
+  {
+    l = t;
+  }
+#endif
+}
+
+
 BDS_Vector::BDS_Vector (const BDS_Point &p2,const BDS_Point &p1)
     : x(p2.X-p1.X),y(p2.Y-p1.Y),z(p2.Z-p1.Z)
 {    
@@ -134,6 +196,50 @@ double dist_droites_gauches(BDS_Point *p1, BDS_Point *p2,
     
 }
 
+bool proj_point_triangle ( double xa, double ya, double za,
+			   BDS_Triangle *t,
+			   double &x, double &y, double &z)
+{
+  const double eps_prec = 1.e-10;
+  BDS_Vector n = t->N();
+  double mat[3][3];
+  double b[3];
+  double res[2];
+  double det;
+  BDS_Point *pts[3];
+  t->getNodes (pts); 
+
+  mat[0][0] = pts[1]->X - pts[0]->X;
+  mat[0][1] = pts[2]->X - pts[0]->X;
+  mat[0][2] = -n.x;
+
+  mat[1][0] = pts[1]->Y - pts[0]->Y;
+  mat[1][1] = pts[2]->Y - pts[0]->Y;
+  mat[1][2] = -n.y;
+
+  mat[2][0] = pts[1]->Z - pts[0]->Z;
+  mat[2][1] = pts[2]->Z - pts[0]->Z;
+  mat[2][2] = -n.z;
+
+  b[0] = xa - pts[0]->X;
+  b[1] = ya - pts[0]->Y;
+  b[2] = za - pts[0]->Z; 
+  if(!sys3x3(mat, b, res, &det)) return false; 
+
+  /* coordonnees dans la face */
+  if(res[0] >= 1.0 + eps_prec || res[0] <= -eps_prec)
+    return false;
+  if(res[1] >= 1.0 + eps_prec || res[1] <= -eps_prec)
+    return false;
+  if(res[1] >= 1. + eps_prec - res[0])
+    return false;
+  
+  x = xa + res[2] * n.x;
+  y = ya + res[2] * n.y;
+  z = za + res[2] * n.z;
+  return true;
+}
+
 void proj_point_plane ( double xa, double ya, double za,
 			BDS_Point *p1, BDS_Point *p2, BDS_Point *p3 ,
 			double &x, double &y, double &z)
@@ -562,6 +668,90 @@ void BDS_Mesh :: reverseEngineerCAD ( )
 	}
     }
 }
+void BDS_Mesh :: createSearchStructures ( ) 
+{
+#ifdef HAVE_ANN_
+
+  printf("creating the ANN search structure\n");
+
+    for (std::set<BDS_GeomEntity*,GeomLessThan>::iterator it = geom.begin();
+	 it != geom.end();
+	 ++it)
+    {
+	if ((*it)->classif_degree == 2)
+	{
+	    std::set<BDS_Point*> pts;
+	    std::set<BDS_Edge*> eds;
+	    
+	    if ((*it)->t.size() > 10)
+	      {
+		std::list<BDS_Triangle*>::iterator tit  = (*it)->t.begin();
+		std::list<BDS_Triangle*>::iterator tite = (*it)->t.end();
+		while (tit!=tite)
+		  {
+		    BDS_Point *nod[3];
+		    (*tit)->getNodes (nod);
+		    pts.insert (nod[0]);
+		    pts.insert (nod[1]);
+		    pts.insert (nod[2]);
+		    eds.insert((*tit)->e1);
+		    eds.insert((*tit)->e2);
+		    eds.insert((*tit)->e3);
+		    tit++;
+		  }
+		
+		printf("copying %d nodes %d edges and %d triangles\n",pts.size(),eds.size(),(*it)->t.size());
+		(*it)->sT.resize((*it)->t.size());
+		(*it)->sE.resize(eds.size());
+		(*it)->sP.resize(pts.size());
+		std::copy (  (*it)->t.begin(), (*it)->t.end(), (*it)->sT.begin() );
+		std::copy ( pts.begin(), pts.end(), (*it)->sP.begin()  );
+		std::copy ( eds.begin(), eds.end(), (*it)->sE.begin()  );
+		
+		int maxPts = (*it)->sT.size() + (*it)->sE.size() + (*it)->sP.size();
+		
+		printf("%d pts found\n",maxPts);
+		
+		const int dim = 3;
+		(*it)->queryPt = annAllocPt(dim);                              // allocate query point
+		(*it)->dataPts = annAllocPts(maxPts, dim);                     // allocate data points
+		(*it)->nnIdx = new ANNidx[1];                                  // allocate near neigh indices
+		(*it)->dists = new ANNdist[1];                                 // allocate near neighbor dists  
+		
+		int I = 0;
+		
+		for (int i=0;i<(*it)->sP.size();i++)
+		  {
+		    (*it)->dataPts[I][0] = (*it)->sP[i]->X;
+		    (*it)->dataPts[I][1] = (*it)->sP[i]->Y;
+		    (*it)->dataPts[I][2] = (*it)->sP[i]->Z;
+		    I++;
+		  }
+		for (int i=0;i<(*it)->sE.size();i++)
+		  {
+		    (*it)->dataPts[I][0] = 0.5 * ((*it)->sE[i]->p1->X+(*it)->sE[i]->p2->X);
+		    (*it)->dataPts[I][1] = 0.5 * ((*it)->sE[i]->p1->Y+(*it)->sE[i]->p2->Y);
+		    (*it)->dataPts[I][2] = 0.5 * ((*it)->sE[i]->p1->Z+(*it)->sE[i]->p2->Z);
+		    I++;
+		  }
+		for (int i=0;i<(*it)->sT.size();i++)
+		  {
+		    const BDS_Vector cog = (*it)->sT[i]->cog();
+		    (*it)->dataPts[I][0] = cog.x;
+		    (*it)->dataPts[I][1] = cog.y;
+		    (*it)->dataPts[I][2] = cog.z;
+		    I++;
+		  } 
+		printf("building kd Tree for surface %d -- %d points\n",(*it)->classif_tag,maxPts);
+		(*it)->kdTree = new ANNkd_tree(                                        // build search structure
+					       (*it)->dataPts,                                // the data points
+					       maxPts,                                 // number of points
+					       dim);
+	      }
+	}
+    }
+#endif
+}
 
 void BDS_Point :: compute_curvature ( )
 {
@@ -717,7 +907,7 @@ void BDS_Mesh :: classify ( double angle, int NB_T )
 
     // color plane surfaces
 
-    if (NB_T > 0)color_plane_surf (1.e-4, NB_T);
+    //    if (NB_T > 0)color_plane_surf (1.e-4, NB_T);
 
 
     {
@@ -948,8 +1138,9 @@ void BDS_Mesh :: classify ( double angle, int NB_T )
     }
     printf(" reverse engineering surfaces\n");
     reverseEngineerCAD ( ) ;
-
-  printf("  end classifying %d edgetags %d vertextags\n",edgetag-1,vertextag-1);
+    printf(" creating search  structures\n");
+    createSearchStructures ( ) ;
+    printf("  end classifying %d edgetags %d vertextags\n",edgetag-1,vertextag-1);
 }
 double PointLessThanLexicographic::t = 0;
 double BDS_Vector::t = 0;
@@ -1476,7 +1667,7 @@ BDS_Mesh ::~ BDS_Mesh ()
     DESTROOOY ( triangles.begin(),triangles.end());
 }
 
-bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
+bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord, BDS_Mesh *geom )
 {
 
 
@@ -1615,6 +1806,9 @@ bool BDS_Mesh ::split_edge ( BDS_Edge *e, double coord)
     triangles.push_back (t3); 
     triangles.push_back (t4); 
 
+    snap_point (mid, geom);
+
+
     return true;
 }
 
@@ -1841,54 +2035,88 @@ bool BDS_Mesh ::collapse_edge ( BDS_Edge *e, BDS_Point *p, const double eps)
 
 */
 
-#ifdef _HAVE_ANN
-#include <ANN/ANN.h>
-#endif
-
-void project_point_on_a_list_of_triangles ( BDS_Point *p , 
+bool project_point_on_a_list_of_triangles ( BDS_Point *p , 
 					    const std::list<BDS_Triangle*> &t,
 					    double &X, double &Y, double &Z)
 {
-    int p_classif = p->g->classif_tag;
-    BDS_Vector N = p->N();
-    std::list<BDS_Triangle *>::const_iterator it = t.begin();
-    std::list<BDS_Triangle *>::const_iterator ite = t.end() ;
-    double dist = 1.e22;
-    double XX = p->X;
-    double YY = p->Y;
-    double ZZ = p->Z;
-
-//    printf("looking at %d triangles\n",t.size());
-
-    while (it != ite)
+  bool global_ok = false;
+  
+  int p_classif = p->g->classif_tag;
+  std::list<BDS_Triangle *>::const_iterator it = t.begin();
+  std::list<BDS_Triangle *>::const_iterator ite = t.end() ;
+  double dist2 = 1.e22;
+  double XX = p->X;
+  double YY = p->Y;
+  double ZZ = p->Z;
+  
+  while (it != ite)
     { 
-	if ((*it)->g->classif_tag == p_classif && fabs(N.angle((*it)->N())) < (M_PI/6))
-	{
-	    BDS_Point *pts[3];
-	    (*it)->getNodes (pts); 
+      if ((*it)->g->classif_tag == p_classif)
+	{	  
+	  {
 	    double xp,yp,zp;
-	    proj_point_plane ( p->X,p->Y,p->Z,pts[0],pts[1],pts[2],xp,yp,zp);
-	    BDS_Point ppp (0,xp,yp,zp);
-	    double a1 = surface_triangle ( &ppp, pts[0],pts[1] );
-	    double a2 = surface_triangle ( &ppp, pts[1],pts[2] );
-	    double a3 = surface_triangle ( &ppp, pts[2],pts[0] );
-	    double a  = surface_triangle ( pts[0], pts[1],pts[2] );
-	    if ( fabs (a-a1-a2-a3) < 1.e-2 * a)
-	    {
-		double d = sqrt ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
-		if (d < dist )
-		{
-//		    printf("one found %g %g %g %g\n",a,a1,a2,a3);
+	    bool ok = proj_point_triangle ( p->X,p->Y,p->Z,*it,xp,yp,zp);
+	    if (ok)
+	      {
+		global_ok = true;
+		double d2 = ((xp-X)*(xp-X)+(yp-Y)*(yp-Y)+(zp-Z)*(zp-Z));
+		if (d2 < dist2 )
+		  {
+		    //		    printf("one found among %d\n",t.size());
 		    XX = xp; YY = yp; ZZ = zp;
-		    dist = d;
-		}
-	    }
+		    dist2 = d2;
+		  }
+	      }
+	  }
 	}
-	++it;
+      ++it;
     }
+
     X = XX;
     Y = YY;
     Z = ZZ;
+    return global_ok;
+}
+
+void BDS_Mesh :: snap_point ( BDS_Point *p , BDS_Mesh *geom_mesh )
+{
+    if (p->g->surf)
+    {
+	p->g->surf->projection ( p->X,p->Y,p->Z,p->X,p->Y,p->Z);
+    }
+    else if (p->g && p->g->classif_degree == 2 && geom_mesh)
+	{
+	  std::list<BDS_Triangle*> l;
+	  BDS_GeomEntity *gg = geom_mesh->get_geom(p->g->classif_tag,p->g->classif_degree);
+	  gg->getClosestTriangles (p->X,p->Y,p->Z,l);
+	  bool ok = project_point_on_a_list_of_triangles ( p , l, p->X,p->Y,p->Z);
+	  
+	  if (!ok)
+	    {
+	      SNAP_FAILURE++;
+	    }
+	  else
+	    {
+	      SNAP_SUCCESS++;
+	    }
+	  
+	}
+    else
+      {
+	return;
+      }
+    
+    {
+      std::list<BDS_Triangle *> t;
+      p->getTriangles (t); 	
+      std::list<BDS_Triangle *>::iterator tit  = t.begin(); 	
+      std::list<BDS_Triangle *>::iterator tite =  t.end();
+      while (tit!=tite)
+	{
+	  (*tit)->_update();
+	  ++tit;
+	}      
+    }
 }
 
 bool BDS_Mesh ::smooth_point ( BDS_Point *p , BDS_Mesh *geom_mesh )
@@ -1915,38 +2143,12 @@ bool BDS_Mesh ::smooth_point ( BDS_Point *p , BDS_Mesh *geom_mesh )
     Y /= p->edges.size();
     Z /= p->edges.size();
 
-    if (p->g->surf)
-    {
-	p->g->surf->projection ( X,Y,Z,p->X,p->Y,p->Z);
-//	printf("coucou\n");
-    }
-    else
-    {
-      //	return false;
-	p->X = X;
-	p->Y = Y;
-	p->Z = Z;
-	//	std::list<BDS_Triangle *>t;
-	//	p->getTriangles (t);
-	//	project_point_on_a_list_of_triangles ( p , t, p->X,p->Y,p->Z);	
-	if (p->g && geom_mesh)
-	  {
-	    BDS_GeomEntity *gg = geom_mesh->get_geom(p->g->classif_tag,p->g->classif_degree);
-	    project_point_on_a_list_of_triangles ( p , gg->t, p->X,p->Y,p->Z);
-	  }
-    }
-    
-    {
-      std::list<BDS_Triangle *> t;
-      p->getTriangles (t); 	
-      std::list<BDS_Triangle *>::iterator tit  = t.begin(); 	
-      std::list<BDS_Triangle *>::iterator tite =  t.end();
-      while (tit!=tite)
-	{
-	  (*tit)->_update();
-	  ++tit;
-	}      
-    }
+    p->X = X;
+    p->Y = Y;
+    p->Z = Z;
+
+    snap_point ( p, geom_mesh );
+
     return true;
 }
 
@@ -2045,6 +2247,8 @@ void BDS_Mesh :: compute_metric_edge_lengths (const BDS_Metric & metric)
 int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh) 
 {
     int nb_modif = 0;
+    SNAP_SUCCESS = 0;
+    SNAP_FAILURE = 0;
 
     BDS_Metric metric ( l , LC/100 , LC );
     //    printf("METRIC %g %g %g\n",LC,metric._min,metric._max);
@@ -2073,7 +2277,7 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
 	    {
 	      double length = (*it)->length();
 	      if (!(*it)->deleted && length > (*it)->target_length / 0.7 ){
-		split_edge (*it, 0.5 );
+		split_edge (*it, 0.5,geom_mesh );
 		nb_modif++;
 	      }
 	    }
@@ -2180,6 +2384,7 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
     cleanup();  
     if (smooth)
     {    
+      printf("smoothing %d points\n",points.size());
 	std::set<BDS_Point*, PointLessThan>::iterator it   = points.begin();
 	std::set<BDS_Point*, PointLessThan>::iterator ite  = points.end();
 	while (it != ite)
@@ -2188,6 +2393,9 @@ int BDS_Mesh :: adapt_mesh ( double l, bool smooth, BDS_Mesh *geom_mesh)
 	    ++it;
 	}
     }
+    
+    printf("%d snaps have succeeded , %d have failed\n",SNAP_SUCCESS,SNAP_FAILURE);
+
     return nb_modif;
 }
 
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 19231803763d1ddf980b55b0566ff69256694e51..98036b0a134860a3b03eee8cf30d8ecd1627c5f9 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -10,6 +10,9 @@
 #include <list>
 #include <math.h>
 #include <algorithm>
+#ifdef HAVE_ANN_
+#include "ANN/ANN.h"
+#endif
 
 class BDS_Edge;
 class BDS_Triangle;
@@ -62,11 +65,22 @@ public:
     int classif_tag;
     int classif_degree;
 
+#ifdef HAVE_ANN_
+    ANNpointArray           dataPts;                                // data points
+    ANNpoint                queryPt;                                // query point
+    ANNidxArray             nnIdx;                                  // near neighbor indices
+    ANNdistArray            dists;                                  // near neighbor distances
+    ANNkd_tree*             kdTree;                                 // search structure
+    std::vector<BDS_Triangle *> sT;
+    std::vector<BDS_Edge *> sE;
+    std::vector<BDS_Point *> sP;
+#endif
     std::list<BDS_Triangle *> t;
     std::list<BDS_Edge *>     e;
     BDS_Point   *p;
     BDS_Surface *surf;
-
+    void getClosestTriangles ( double x, double y, double z, 
+			  std::list<BDS_Triangle*> &l );
     inline bool operator <  ( const BDS_GeomEntity & other ) const
 	{
 	    if (classif_degree < other.classif_degree)return true;
@@ -74,10 +88,25 @@ public:
 	    if (classif_tag < other.classif_tag)return true;
 	    return false;
 	}
-    BDS_GeomEntity (int a, int b)
-	: classif_tag (a),classif_degree(b),p(0),surf(0)
-	{
-	}
+    BDS_GeomEntity (int a, int b)  
+      : classif_tag (a),classif_degree(b),p(0),surf(0)
+      {
+#ifdef HAVE_ANN_
+	kdTree = 0;
+#endif
+      }
+    ~BDS_GeomEntity ()  
+      {
+#ifdef HAVE_ANN_
+	if (kdTree)
+	  {
+	    delete [] nnIdx;                                                    // clean things up
+	    delete [] dists;
+	    delete kdTree;
+	  }
+#endif
+      }
+
 };
 
 
@@ -446,7 +475,8 @@ class EdgeLessThan
 
 class BDS_Mesh 
 {    
-    int MAXPOINTNUMBER;
+    int MAXPOINTNUMBER,SNAP_SUCCESS,SNAP_FAILURE;
+    
  public:
     double Min[3],Max[3],LC;
 
@@ -472,12 +502,14 @@ class BDS_Mesh
     BDS_GeomEntity *get_geom  (int p1, int p2);
     bool swap_edge ( BDS_Edge *);
     bool collapse_edge ( BDS_Edge *, BDS_Point*, const double eps);
+    void snap_point   ( BDS_Point* , BDS_Mesh *geom = 0);
     bool smooth_point   ( BDS_Point* , BDS_Mesh *geom = 0);
     bool smooth_point_b ( BDS_Point* );
-    bool split_edge ( BDS_Edge *, double coord);
+    bool split_edge ( BDS_Edge *, double coord, BDS_Mesh *geom = 0);
     void classify ( double angle, int nb = -1); 
     void color_plane_surf ( double eps , int nb);
     void reverseEngineerCAD ( ) ;
+    void createSearchStructures ( ) ;
     int adapt_mesh(double,bool smooth = false,BDS_Mesh *geom = 0); 
     void compute_metric_edge_lengths (const BDS_Metric & metric);
     void cleanup();
@@ -489,5 +521,5 @@ class BDS_Mesh
     bool read_vrml ( const char *filename);
     void save_gmsh_format (const char *filename);
 };
-void project_point_on_a_list_of_triangles ( BDS_Point *p , const std::list<BDS_Triangle*> &t,
+bool project_point_on_a_list_of_triangles ( BDS_Point *p , const std::list<BDS_Triangle*> &t,
 					    double &X, double &Y, double &Z);	   
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index 05d68834405b812ec685e893716152960e7d6e27..35d7ade82d7a0e86ad026a3e81cdce13243e1179 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.20 2005-07-12 15:01:17 remacle Exp $
+// $Id: DiscreteSurface.cpp,v 1.21 2005-08-19 14:07:34 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -178,30 +178,33 @@ void BDS_To_Mesh(Mesh *m)
 
 int MeshDiscreteSurface(Surface *s)
 { 
-    if(s->bds){
+
+  const int NITER = 20;
+  if(s->bds){
     Msg(STATUS2, "Discrete Surface Mesh Generator...");
-	// s->bds is the discrete surface that defines the geometry
-	if(!THEM->bds_mesh){
-	    THEM->bds_mesh = new BDS_Mesh (*(THEM->bds));
-	    int iter = 0;
-	    while(iter < 20 && THEM->bds_mesh->adapt_mesh(CTX.mesh.lc_factor * THEM->bds->LC, 
+    // s->bds is the discrete surface that defines the geometry
+    if(!THEM->bds_mesh){
+      THEM->bds_mesh = new BDS_Mesh (*(THEM->bds));
+      int iter = 0;
+      while(iter < NITER && THEM->bds_mesh->adapt_mesh(CTX.mesh.lc_factor * THEM->bds->LC, 
 						    true, THEM->bds)){
-		Msg(STATUS2, "Iteration %2d/20 done (%d triangles)\n",iter, THEM->bds_mesh->triangles.size());
-		iter ++;
-	    }
-	    BDS_To_Mesh_2(THEM);
-	    Msg(STATUS2, "Mesh has %d vertices (%d)\n",Tree_Nbr(THEM->Vertices),THEM->bds->points.size());
-//	    THEM->bds_mesh->save_gmsh_format ( "3.msh" );
-	}
-	return 1;
+	Msg(STATUS2, "Iteration %2d/%d done (%d triangles)\n",iter, NITER,THEM->bds_mesh->triangles.size());
+	iter ++;
+      }
+      BDS_To_Mesh_2(THEM);
+      Msg(STATUS2, "Mesh has %d vertices (%d)\n",Tree_Nbr(THEM->Vertices),THEM->bds->points.size());
+      //	    THEM->bds_mesh->save_gmsh_format ( "3.msh" ); 
+      return 1;
     }
-    else if(s->Typ == MSH_SURF_DISCRETE){
-	// nothing to do: we suppose that the surface is represented by
-	// a mesh that will not be modified
-	return 1;
-    }
-    else
-	return 0;
+    return 2;
+  }
+  else if(s->Typ == MSH_SURF_DISCRETE){
+    // nothing to do: we suppose that the surface is represented by
+    // a mesh that will not be modified
+    return 2;
+  }
+  else
+    return 0;
 }
 
 int MeshDiscreteCurve(Curve *c)
diff --git a/Mesh/Makefile b/Mesh/Makefile
index a35696474f5f57e6a2428a0c6bb1db0d2be875da..f955b0b1f24a07f6bd65f14aad3e22aa5154a753 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.92 2005-07-15 10:31:06 geuzaine Exp $
+# $Id: Makefile,v 1.93 2005-08-19 14:07:34 remacle Exp $
 #
 # Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 #
@@ -24,7 +24,7 @@ include ../variables
 LIB     = ../lib/libGmshMesh.a
 INCLUDE = -I../Numeric -I../NR -I../Common -I../DataStr -I../Geo -I../Mesh\
           -I../Graphics -I../Parser -I../Fltk -I../Triangle -I../Tetgen\
-          -I../Netgen -I../Netgen/libsrc/include  -I../Netgen/libsrc/interface
+          -I../Netgen -I../Netgen/libsrc/include  -I../Netgen/libsrc/interface -I../ANN/include
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE}
 
 SRC = 1D_Mesh.cpp \
diff --git a/Parser/Makefile b/Parser/Makefile
index 48546d61a81a9cec613c5815432164809a645ded..875bef0aac261945126f28fa95cfa8573db6b5f4 100644
--- a/Parser/Makefile
+++ b/Parser/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.73 2005-06-03 17:32:30 geuzaine Exp $
+# $Id: Makefile,v 1.74 2005-08-19 14:07:35 remacle Exp $
 #
 # Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 #
@@ -23,7 +23,7 @@ include ../variables
 
 LIB     = ../lib/libGmshParser.a
 INCLUDE = -I../Common -I../DataStr -I../Geo -I../Graphics\
-          -I../Mesh -I../Numeric -I../Fltk -I../Plugin -I../Parallel
+          -I../Mesh -I../Numeric -I../Fltk -I../Plugin -I../Parallel -I../ANN/include
 CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE}
 
 SRC = Gmsh.tab.cpp\
diff --git a/configure b/configure
index 3ba70a49f13f44c30902bdd2f83ba586f130c16a..33da8eb0fc85864fdc4533ac34a42d9376d28be3 100755
--- a/configure
+++ b/configure
@@ -851,6 +851,7 @@ Optional Features:
   --enable-gsl            use GSL as numerical toolkit (default=yes)
   --enable-gui            build the graphical user interface (default=yes)
   --enable-parallel       enable parallel version (default=no)
+  --enable-ann            compile ANN if available (default=yes)
   --enable-triangle       compile Triangle if available (default=yes)
   --enable-netgen         compile Netgen if available (default=yes)
   --enable-tetgen         compile Tetgen if available (default=yes)
@@ -1366,6 +1367,11 @@ fi;
 if test "${enable_parallel+set}" = set; then
   enableval="$enable_parallel"
 
+fi;
+# Check whether --enable-ann or --disable-ann was given.
+if test "${enable_ann+set}" = set; then
+  enableval="$enable_ann"
+
 fi;
 # Check whether --enable-triangle or --disable-triangle was given.
 if test "${enable_triangle+set}" = set; then
@@ -3778,6 +3784,56 @@ else
   fi
 fi
 
+
+echo "$as_me:$LINENO: checking for ./ANN/include/ANN/ANN.h" >&5
+echo $ECHO_N "checking for ./ANN/include/ANN/ANN.h... $ECHO_C" >&6
+if test "${ac_cv_file___ANN_include_ANN_ANN_h+set}" = set; then
+  echo $ECHO_N "(cached) $ECHO_C" >&6
+else
+  test "$cross_compiling" = yes &&
+  { { echo "$as_me:$LINENO: error: cannot check for file existence when cross compiling" >&5
+echo "$as_me: error: cannot check for file existence when cross compiling" >&2;}
+   { (exit 1); exit 1; }; }
+if test -r "./ANN/include/ANN/ANN.h"; then
+  ac_cv_file___ANN_include_ANN_ANN_h=yes
+else
+  ac_cv_file___ANN_include_ANN_ANN_h=no
+fi
+fi
+echo "$as_me:$LINENO: result: $ac_cv_file___ANN_include_ANN_ANN_h" >&5
+echo "${ECHO_T}$ac_cv_file___ANN_include_ANN_ANN_h" >&6
+if test $ac_cv_file___ANN_include_ANN_ANN_h = yes; then
+  ANN="yes"
+else
+  ANN="no"
+fi
+
+if test "x${ANN}" = "xyes"; then
+  if test "x$enable_ann" != "xno"; then
+     GMSH_DIRS="${GMSH_DIRS} ANN"
+     GMSH_LIBS="${GMSH_LIBS} -lGmshANN"
+     FLAGS="-DHAVE_ANN_ ${FLAGS}"
+     echo "********************************************************************"
+     echo "You are building a version of Gmsh that contains ANN, the Approximate"
+     echo "Nearest Neighbor library."
+     echo "Please note that by doing so, you agree with ANN's licensing"
+     echo "requirements stated in ./ANN/Copyright.txt."
+     echo "To disable ANN, run configure again with the --disable-ann"
+     echo "option."
+     echo "********************************************************************"
+  fi
+else
+  if test "x$enable_ann" != "xno"; then
+     echo "********************************************************************"
+     echo "If you want to use ANN for doing fast geometrical searchs in the"
+     echo "STL mesher, please download ANN from the"
+     echo "author's web site at http://www.cs.umd.edu/~mount/ANN/,"
+     echo "unpack the archive and copy both src and include directories in the ./ANN subdirectory."
+     echo "Then run ./configure again."
+     echo "********************************************************************"
+  fi
+fi
+
 echo "$as_me:$LINENO: checking for ./Netgen/libsrc/meshing/meshclass.cpp" >&5
 echo $ECHO_N "checking for ./Netgen/libsrc/meshing/meshclass.cpp... $ECHO_C" >&6
 if test "${ac_cv_file___Netgen_libsrc_meshing_meshclass_cpp+set}" = set; then