diff --git a/Common/Context.h b/Common/Context.h
index ecfbd3c2aab55badd3afc8cb11c85b7e48551aae..114928da9346a5b49cf08cad70d6959b6f7e752e 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -197,7 +197,7 @@ public :
     int smooth_normals;
     double angle_smooth_normals;
     double stl_distance_tol, dihedral_angle_tol;
-    int edge_prolongation_threshold;
+    int edge_prolongation_threshold, do_not_coarsen;
     double  nb_elem_per_rc ,   min_elem_size_fact , target_elem_size_fact, beta_smooth_metric;
   } mesh;
 
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 269143727232b2963fd214566374c60fe5755418..7c431b15c0a3a5af02d1fe616cef4d94bc9e13a2 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.465 2005-11-02 19:31:49 geuzaine Exp $
+// $Id: GUI.cpp,v 1.466 2005-11-03 13:44:02 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -2276,6 +2276,12 @@ void GUI::create_option_window()
       mesh_value[25]->step(.01);
       mesh_value[25]->align(FL_ALIGN_RIGHT); 
 
+      mesh_value[26] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 8 * BH, IW, BH, "Allow coarsening of the initial mesh");
+      mesh_value[26]->minimum(0);
+      mesh_value[26]->maximum(1);
+      mesh_value[26]->step(1);
+      mesh_value[26]->align(FL_ALIGN_RIGHT); 
+
       o->end();
     }
 
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 9518f3f4af4dffc97d358fbb1d868815701c4fa0..fbf8c1bdf5b69de9f50451ad92650ff800732fb5 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.42 2005-11-03 03:55:14 geuzaine Exp $
+// $Id: BDS.cpp,v 1.43 2005-11-03 13:44:03 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -777,7 +777,7 @@ void BDS_Mesh::createSearchStructures()
 
   Msg(INFO, "Creating the ANN search structure");
 
-  const double LC_SEARCH = LC * 1e-2;
+  const double LC_SEARCH = LC * 1e-3;
 
   for(std::set < BDS_GeomEntity *, GeomLessThan >::iterator it = geom.begin();
       it != geom.end(); ++it) {
@@ -1742,37 +1742,42 @@ template < class IT > void DESTROOOY(IT beg, IT end)
 void BDS_Mesh::cleanup()
 {
   {
-    for(std::list < BDS_Triangle * >::iterator it = triangles.begin();
-        it != triangles.end(); it++) {
-      while(it != triangles.end() && (*it)->deleted) {
-        delete *it;
-        it = triangles.erase(it);
-      }
-    }
-  }
-  {
-    for(std::list < BDS_Edge * >::iterator it = edges.begin();
-        it != edges.end(); it++) {
-      while(it != edges.end() && (*it)->deleted) {
-        delete *it;
-        it = edges.erase(it);
-      }
-    }
-  }
-  // DESTROOOY(edges_to_delete.begin(),edges_to_delete.end());
-  // edges_to_delete.clear();
+    for (std::list<BDS_Triangle*> :: iterator it = triangles.begin();
+	 it != triangles.end();
+	 it++)
+      {
+	while (it != triangles.end() && (*it)->deleted)
+	  {
+	    delete *it;
+	    it = triangles.erase (it);
+	  }
+      }	   
+  }
+  { 
+    for (std::list<BDS_Edge*> :: iterator it = edges.begin();
+	 it != edges.end();
+	 it++)
+      {
+	while (it != edges.end() && (*it)->deleted)
+	  {
+	    delete *it;
+	    it = edges.erase (it);
+	  }	
+      }	   
+  } 
 }
 
-BDS_Mesh::~BDS_Mesh()
+BDS_Mesh ::~ BDS_Mesh ()
 {
-  // DESTROOOY(geom.begin(),geom.end());
-  // DESTROOOY(points.begin(),points.end());
-  // cleanup();
-  // DESTROOOY(edges.begin(),edges.end());
-  // DESTROOOY(triangles.begin(),triangles.end());
-  // DESTROOOY(tets.begin(),tets.end());
+   DESTROOOY ( geom.begin(),geom.end());
+   DESTROOOY ( points.begin(),points.end());
+   cleanup();
+   DESTROOOY ( edges.begin(),edges.end());
+   DESTROOOY ( triangles.begin(),triangles.end());
+   DESTROOOY ( tets.begin(),tets.end());
 }
 
+
 bool BDS_Mesh::split_edge(BDS_Edge * e, double coord, BDS_Mesh * geom)
 {
   /*
@@ -2054,9 +2059,12 @@ bool BDS_Mesh::collapse_edge(BDS_Edge * e, BDS_Point * p, const double eps)
         // this does not concern the triangles with the small edge that
         // are collapsed too
         double angle = n1.angle(n2);
-        if(fabs(angle) > M_PI / 2)
+        if(fabs(angle) > M_PI / 12)
           return false;
-        if(s_after < 1.e-2 * s_before) {
+        if(s_after+s_before < 3.e-1 * fabs(s_before-s_after)) {
+          return false;
+        }
+        if(s_after < 1.e-1 * s_before) {
           return false;
         }
         gs[nt] = t->g;
@@ -2374,6 +2382,21 @@ void BDS_Mesh::compute_metric_edge_lengths(const BDS_Metric & metric)
 }
 
 
+void BDS_Mesh::applyOptimizationPatterns()
+{
+  std::set<BDS_Point*, PointLessThan>::iterator it   = points.begin();
+  std::set<BDS_Point*, PointLessThan>::iterator ite  = points.end();
+  while (it != ite)
+    {
+      if ((*it)->edges.size() == 3 && (*it)->g->classif_degree == 2)
+	{
+	  BDS_Edge *e = *((*it)->edges.begin());
+	  if (!e->deleted)collapse_edge ( e , *it , .1);
+	}
+      ++it;
+    }
+}
+
 int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
                          BDS_Mesh * geom_mesh)
 {
@@ -2410,93 +2433,127 @@ int BDS_Mesh::adapt_mesh(const BDS_Metric & metric, bool smooth,
   cleanup();
   small_to_long = edges;
 
-  // collapse 
+  // split edges
   {
-    std::list < BDS_Edge * >::iterator it = small_to_long.begin();
-    std::list < BDS_Edge * >::iterator ite = small_to_long.end();
-    //  compute_metric_edge_lengths (metric);
-
-    while(it != ite) {
-      double length = (*it)->length();
-      if(!(*it)->deleted && length < (*it)->target_length * 0.7) {
-        if(nb_modif % 2) {
-          if(collapse_edge(*it, (*it)->p1, 0.1))
-            nb_modif++;
-        }
-        else {
-          if(collapse_edge(*it, (*it)->p2, 0.1))
-            nb_modif++;
-        }
+    std::list<BDS_Edge*>::iterator it = small_to_long.begin();
+    std::list<BDS_Edge*>::iterator ite  = small_to_long.end();
+    compute_metric_edge_lengths (metric);
+    
+    while (it != ite)
+      {
+	if ((*it)->numfaces()==2)
+	  {
+	    double length = (*it)->length();
+	    if (!(*it)->deleted && length > (*it)->target_length / 0.7 ){
+	      split_edge (*it, 0.5,geom_mesh );
+	      //split_edge (*it, 0.5, 0  );
+	      nb_modif++;
+	    }
+	  }
+	++it;
       }
-      ++it;
     }
+  
+  // re-create small_to_long
+  cleanup();    
+  small_to_long = edges;
+  
+  // collapse 
+  {    	
+    std::list<BDS_Edge*>::iterator it = small_to_long.begin();
+    std::list<BDS_Edge*>::iterator ite  = small_to_long.end();
+    //	compute_metric_edge_lengths (metric);
+    
+    while (it != ite)
+      {
+	double length = (*it)->length();
+	if (!(*it)->deleted && length < (*it)->target_length * 0.7 ){
+	  if (nb_modif %2 )
+	    {
+	      if (collapse_edge (*it, (*it)->p1, 0.1 ))
+		nb_modif++;
+	    }
+	  else
+	    {
+	      if (collapse_edge (*it, (*it)->p2, 0.1 ))
+		nb_modif++;
+	    }
+	}
+	++it;
+      }
   }
-  cleanup();
+  cleanup();  
   small_to_long = edges;
-
-  {
-    std::list < BDS_Edge * >::iterator it = small_to_long.begin();
-    std::list < BDS_Edge * >::iterator ite = small_to_long.end();
-    while(it != ite) {
-      if(!(*it)->deleted && (*it)->numfaces() == 2) {
-        BDS_Point *op[2];
-        (*it)->oppositeof(op);
-
-        double a1 = surface_triangle((*it)->p1, (*it)->p2, op[0]);
-        double a2 = surface_triangle((*it)->p1, (*it)->p2, op[1]);
-        double b1 = surface_triangle((*it)->p1, op[0], op[1]);
-        double b2 = surface_triangle((*it)->p2, op[0], op[1]);
-
-        double cb1[3], cb2[3];
-        normal_triangle((*it)->p1, op[0], op[1], cb1);
-        normal_triangle((*it)->p2, op[0], op[1], cb2);
-
-        double prosc = cb1[0] * cb2[0] + cb1[1] * cb2[1] + cb1[2] * cb2[2];
-
-        if(fabs(a1 + a2 - b1 - b2) < 0.1 * (a1 + a2) && prosc < 0)      //&& fabs (prosc) > 0.7)
-        {
-          double qa1 = quality_triangle((*it)->p1, (*it)->p2, op[0]);
-          double qa2 = quality_triangle((*it)->p1, (*it)->p2, op[1]);
-          double qb1 = quality_triangle((*it)->p1, op[0], op[1]);
-          double qb2 = quality_triangle((*it)->p2, op[0], op[1]);
-
-          double d = dist_droites_gauches((*it)->p1, (*it)->p2,
-                                          op[0], op[1]);
-          BDS_Vector v1(*((*it)->p1), *((*it)->p2));
-          BDS_Vector v2(*(op[0]), *(op[1]));
-
-          double dd = sqrt(v1 * v1 + v2 * v2);
-
-          double qa = (qa1 < qa2) ? qa1 : qa2;
-          double qb = (qb1 < qb2) ? qb1 : qb2;
-          // printf("qa %g qb %g ..\n",qa,qb);
-          if(qb > qa && d < 0.01 * dd) {
-            // printf("swop ..\n");
-            nb_modif++;
-            swap_edge(*it);
-          }
-        }
+  
+  {    
+    std::list<BDS_Edge*>::iterator it = small_to_long.begin();
+    std::list<BDS_Edge*>::iterator ite  = small_to_long.end();
+    while (it != ite)
+      {
+	if (!(*it)->deleted && (*it)->numfaces()==2)
+	  {
+	    BDS_Point *op[2];
+	    (*it)->oppositeof (op);
+	    
+	    double a1 = surface_triangle ( (*it)->p1 , (*it)->p2 , op[0] );
+	    double a2 = surface_triangle ( (*it)->p1 , (*it)->p2 , op[1] );
+	    double b1 = surface_triangle ( (*it)->p1 , op[0] , op[1] );
+	    double b2 = surface_triangle ( (*it)->p2 , op[0] , op[1] );
+	    
+	    double cb1[3],cb2[3];
+	    normal_triangle ( (*it)->p1 , op[0] , op[1],cb1 );
+	    normal_triangle ( (*it)->p2 , op[0] , op[1],cb2 );
+	    
+	    double prosc = cb1[0]*cb2[0]+cb1[1]*cb2[1]+cb1[2]*cb2[2];
+	    
+	    if (fabs(a1+a2-b1-b2) < 0.1 * (a1+a2))
+	      {
+		double qa1 = quality_triangle ( (*it)->p1 , (*it)->p2 , op[0] );
+		double qa2 = quality_triangle ( (*it)->p1 , (*it)->p2 , op[1] );
+		double qb1 = quality_triangle ( (*it)->p1 , op[0] , op[1] );
+		double qb2 = quality_triangle ( (*it)->p2 , op[0] , op[1] );
+		
+		double d = dist_droites_gauches((*it)->p1, (*it)->p2, 
+						op[0],op[1]);
+		BDS_Vector v1 (*((*it)->p1), *((*it)->p2));
+		BDS_Vector v2 (*(op[0]),*(op[1]));
+		
+		double dd = sqrt (v1*v1 + v2*v2);
+
+		double qa = (qa1<qa2)?qa1:qa2; 
+		double qb = (qb1<qb2)?qb1:qb2; 
+		//		  printf("qa %g qb %g ..\n",qa,qb);
+		if (qb > qa && d < 0.01 * dd)
+		  {
+		    nb_modif++;
+		    swap_edge ( *it );
+		  }
+	      }
+	  }
+	++it;
       }
-      ++it;
-    }
   }
-  cleanup();
-  if(smooth) {
-    Msg(INFO, "smoothing %d points", points.size());
-    std::set < BDS_Point *, PointLessThan >::iterator it = points.begin();
-    std::set < BDS_Point *, PointLessThan >::iterator ite = points.end();
-    while(it != ite) {
-      smooth_point(*it, geom_mesh);
-      ++it;
-    }
+  cleanup();  
+  if (smooth){
+    Msg(INFO,"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)
+      {
+	smooth_point(*it,geom_mesh);
+	++it;
+      }
   }
-
-  Msg(INFO, "%d snaps have succeeded , %d have failed", SNAP_SUCCESS,
-      SNAP_FAILURE);
-  // outputScalarField (triangles,"b.pos");
+  
+  Msg(INFO,"%d snaps have succeeded , %d have failed\n",SNAP_SUCCESS,SNAP_FAILURE);
+    //    outputScalarField (triangles,"b.pos");
+  applyOptimizationPatterns();
+  cleanup();  
   return nb_modif;
 }
 
+
+
 BDS_Mesh::BDS_Mesh(const BDS_Mesh & other)
 {
   MAXPOINTNUMBER = other.MAXPOINTNUMBER;
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 77c7f2e2dea75b9b651263326b2847a69a8d4e7a..70a8606b9919faafa64b59e46b368c4696db1d4d 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -119,11 +119,13 @@ public:
   ~BDS_GeomEntity()  
   {
 #ifdef HAVE_ANN_
-    if(kdTree){
-      delete [] nnIdx;                                                    // clean things up
-      delete [] dists;
-      delete kdTree;
-    }
+	if (kdTree)
+	  {
+	    annDeallocPts(dataPts);
+	    delete [] nnIdx;                                                    // clean things up
+	    delete [] dists;	    
+	    delete kdTree;
+	  }
 #endif
   }
 };
@@ -638,6 +640,7 @@ public:
   bool read_mesh(const char *filename);
   bool read_vrml(const char *filename);
   void save_gmsh_format(const char *filename);
+  void applyOptimizationPatterns();
 };
 
 bool project_point_on_a_list_of_triangles(BDS_Point *p , const std::list<BDS_Triangle*> &t,
diff --git a/Mesh/DiscreteSurface.cpp b/Mesh/DiscreteSurface.cpp
index c058d3a37c52f4c6f0b08c9e1cb5e8f53b702b61..ef52ef61048db86187dc0eb81e6df614600adbb5 100644
--- a/Mesh/DiscreteSurface.cpp
+++ b/Mesh/DiscreteSurface.cpp
@@ -1,4 +1,4 @@
-// $Id: DiscreteSurface.cpp,v 1.32 2005-11-03 03:55:15 geuzaine Exp $
+// $Id: DiscreteSurface.cpp,v 1.33 2005-11-03 13:44:03 remacle Exp $
 //
 // Copyright (C) 1997-2005 C. Geuzaine, J.-F. Remacle
 //
@@ -214,6 +214,7 @@ void BDS_To_Mesh_2(Mesh * m)
   Msg(STATUS2, "Moving surface mesh into old structure");
 
   Tree_Action(m->Vertices, Free_Vertex);
+
   Tree_Delete(m->Vertices);
   m->Vertices = Tree_Create(sizeof(Vertex *), compareVertex);
 
@@ -229,7 +230,6 @@ void BDS_To_Mesh_2(Mesh * m)
       ++it;
     }
   }
-
   {
     std::list < BDS_Edge * >::iterator it = m->bds_mesh->edges.begin();
     std::list < BDS_Edge * >::iterator ite = m->bds_mesh->edges.end();
@@ -247,7 +247,6 @@ void BDS_To_Mesh_2(Mesh * m)
       ++it;
     }
   }
-
   {
     std::list < BDS_Triangle * >::iterator it =
       m->bds_mesh->triangles.begin();
@@ -274,7 +273,6 @@ void BDS_To_Mesh_2(Mesh * m)
       ++it;
     }
   }
-
   {
     std::list < BDS_Tet * >::iterator it = m->bds_mesh->tets.begin();
     std::list < BDS_Tet * >::iterator ite = m->bds_mesh->tets.end();
@@ -391,7 +389,7 @@ int ReMesh(Mesh * M)
 
 int MeshDiscreteSurface(Surface * s)
 {
-  const int NITER = 10;
+  static  int NITER = 10;
   if(THEM->bds) {
     BDS_Metric metric(THEM->bds->LC / CTX.mesh.target_elem_size_fact,
                       THEM->bds->LC / CTX.mesh.min_elem_size_fact,
@@ -414,6 +412,7 @@ int MeshDiscreteSurface(Surface * s)
 
       double t2 = Cpu();
       Msg(STATUS2, "Remesh 2D complete (%g s)", t2 - t1);
+      //      NITER++;
       return 1;
     }
     return 2;