diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 3b659a7ef9a2660769aad6d1157fc22ec8c2d5e2..59e42ee8d4106ca9218b5703b5069a07a47d79f6 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -988,7 +988,7 @@ StringXNumber MeshOptions_Number[] = {
     "Global scaling factor applied to the saved mesh" },
   { F|O, "SecondOrderIncomplete" , opt_mesh_second_order_incomplete , 1. ,
     "Create incomplete second order elements? (8-node quads, 20-node hexas, etc.)" },
-  { F|O, "SecondOrderLinear" , opt_mesh_second_order_linear , 1. ,
+  { F|O, "SecondOrderLinear" , opt_mesh_second_order_linear , 0. ,
     "Should second order vertices simply be created by linear interpolation?" },
   { F|O, "Smoothing" , opt_mesh_nb_smoothing , 1. ,
     "Number of smoothing steps applied to the final mesh" },
@@ -1486,7 +1486,7 @@ StringXColor GeometryOptions_Color[] = {
 
 StringXColor MeshOptions_Color[] = {
   { F|O, "Points" , opt_mesh_color_points , 
-    {0, 0, 128, 255}, {0, 0, 128, 255}, {0, 0, 0, 255},
+    {0, 0, 255, 255}, {0, 0, 255, 255}, {0, 0, 0, 255},
     "Mesh node color" },
   { F|O, "Lines" , opt_mesh_color_lines , 
     {0, 0, 0, 255}, {0, 0, 0, 255}, {0, 0, 0, 255},
diff --git a/Common/Options.cpp b/Common/Options.cpp
index f38860f7114a1295e7b42c81363c5abaacca5965..457f3b4921a00e533335304b4f24dc608c271069 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.310 2006-09-08 02:39:42 geuzaine Exp $
+// $Id: Options.cpp,v 1.311 2006-09-22 19:28:49 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -4873,7 +4873,10 @@ double opt_mesh_use_cut_plane(OPT_ARGS_NUM)
     CTX.mesh.use_cut_plane = (int)val;
 #if defined(HAVE_FLTK)
     if(WID){
-      double val1 = CTX.lc;
+      double val1 = 0;
+      for(int i = 0; i < 3; i++)
+	val1 = std::max(val1, std::max(fabs(CTX.min[i]), fabs(CTX.max[i])));
+      val1 *= 1.5;
       WID->mesh_value[17]->step(val1/200.);
       WID->mesh_value[17]->minimum(-val1);
       WID->mesh_value[17]->maximum(val1);
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index f38ba9323fd887b1c627f52ce7dce2546d5902cf..debda568f6480e9e43096fdc9274fa49b2a4e9f3 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.548 2006-09-08 02:39:42 geuzaine Exp $
+// $Id: GUI.cpp,v 1.549 2006-09-22 19:28:50 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -4180,7 +4180,10 @@ void GUI::reset_clip_browser()
     clip_value[i]->minimum(-1.0);
     clip_value[i]->maximum(1.0);
   }
-  double val1 = CTX.lc;
+  double val1 = 0;
+  for(int i = 0; i < 3; i++)
+    val1 = std::max(val1, std::max(fabs(CTX.min[i]), fabs(CTX.max[i])));
+  val1 *= 1.5;
   for(int i = 3; i < 10; i++){
     clip_value[i]->step(val1/200.);
     clip_value[i]->minimum(-val1);
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index ce8da2b242d98d11dd0b68dd58e1fd8f25c5f933..f9b6ce720fbb1dd2517e3b0158d9dcaea1a92432 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -77,6 +77,15 @@ class MVertex{
   virtual bool getParameter(int i, double &par){ return false; }
   virtual bool setParameter(int i, double par){ return false; }
 
+  // measure distance to another vertex
+  double distance(MVertex *v)
+  {
+    double dx = _x - v->x();
+    double dy = _y - v->y();
+    double dz = _z - v->z();
+    return sqrt(dx * dx + dy * dy + dz * dz);
+  }
+
   // Get the data associated with this vertex
   virtual void *getData(){ return 0; }
 
diff --git a/Geo/fourierModel.cpp b/Geo/fourierModel.cpp
index 274ca40edb0e192dc482228a998bf7ef0660c617..18967f8e149b6b62c7524112c6208aee6aa57efd 100644
--- a/Geo/fourierModel.cpp
+++ b/Geo/fourierModel.cpp
@@ -3,23 +3,6 @@
 #include "Context.h"
 #include "Views.h"
 
-/*
-  Investigate the following approach:
-
-  compute and store the pou of each overlapping patch in the nodes of
-  all the patches
-
-  for each pair of overlapping patches, find the line pou1=pou2 by
-  interpolation on the overlapping grids
-
-  compute the intersections of these lines
-
-  this should define a non-overlapping partitioning of the grid, which
-  could be used as the boundary constrain for the unstructured algo
-
-  Would that work??
-*/
-
 extern Context_T CTX;
 
 #if defined(HAVE_FOURIER_MODEL)
@@ -353,6 +336,183 @@ void meshGrout(GFace *gf, std::vector<MVertex*> &loop, std::vector<MVertex*> &ho
   delete grout;
 }
 
+void meshGroutWithHole(GFace *gf,
+		       std::vector<std::vector<std::vector<MVertex*> > > &inside,
+		       std::vector<std::vector<std::vector<MVertex*> > > &outside)
+{		       
+  std::vector<MVertex*> hole, loop;
+
+  // create hole
+  SPoint2 ic(0., 0.);
+  for(unsigned int i = 0; i < inside[0][0].size(); i++){
+    hole.push_back(inside[0][0][i]);
+    double u, v;
+    inside[0][0][i]->getParameter(0, u);
+    inside[0][0][i]->getParameter(1, v);
+    ic += SPoint2(u, v);
+  }
+  ic *= 1. / (double)inside[0][0].size();
+  hole.push_back(hole[0]);
+
+  // order exterior parts and create exterior loop
+  std::vector<std::pair<double, int> > angle;
+  std::map<int, std::vector<MVertex*> > outside_flat;
+  int part = 0;
+  for(unsigned int i = 0; i < outside.size(); i++){
+    for(unsigned int j = 0; j < outside[i].size(); j++){
+      for(unsigned int k = 0; k < outside[i][j].size(); k++){
+	outside_flat[part].push_back(outside[i][j][k]);
+      }
+      part++;
+    }
+  }
+  std::map<int, std::vector<MVertex*> >::iterator it = outside_flat.begin();
+  for(; it != outside_flat.end(); it++){
+    SPoint2 oc(0., 0.);
+    for(unsigned int i = 0; i < it->second.size(); i++){
+      double u, v;
+      it->second[i]->getParameter(0, u);
+      it->second[i]->getParameter(1, v);
+      oc += SPoint2(u, v);
+    }
+    oc *= 1. / (double)it->second.size();
+    double a = atan2(oc[1] - ic[1], oc[0] - ic[0]);
+    angle.push_back(std::make_pair(a, it->first));
+  }
+  std::sort(angle.begin(), angle.end());
+  for(unsigned int i = 0; i < angle.size(); i++){
+    for(unsigned int j = 0; j < outside_flat[angle[i].second].size(); j++)
+      loop.push_back(outside_flat[angle[i].second][j]);
+  }
+  loop.push_back(hole[0]);
+
+  // mesh the grout
+  meshGrout(gf, loop, hole);
+}
+
+bool onHardEdge(GFace *gf, MVertex *vertex)
+{
+  std::vector<int> edges;
+  FM->GetEdge(gf->tag(), edges);
+  if(edges.empty()) return false;
+  double u, v;
+  vertex->getParameter(0, u);
+  vertex->getParameter(1, v);
+  for(unsigned int i = 0; i < edges.size(); i++){
+    switch(edges[i]){
+    case 0: if(u == 0.) return true;
+    case 1: if(u == 1.) return true;
+    case 2: if(v == 0.) return true;
+    case 3: if(v == 1.) return true;
+    }
+  }
+  return false;
+}
+
+bool removeHardEdges(GFace *gf, std::vector<MVertex*> &loop,
+		     std::vector<std::vector<MVertex*> > &subloops)
+{
+  std::vector<MVertex*> tmp;
+  tmp.push_back(loop[0]);
+  for(unsigned int i = 1; i < loop.size() - 1; i++){
+    if(onHardEdge(gf, loop[i - 1]) && 
+       onHardEdge(gf, loop[i]) &&
+       onHardEdge(gf, loop[i + 1])){
+      // skip + create new path
+    }
+    else{
+      tmp.push_back(loop[i]);
+    }
+  }
+  tmp.push_back(loop[0]);
+  if(tmp.size() != loop.size()){
+    Msg(INFO, "Generating subloops");
+    
+
+    return true;
+  }
+  return false;
+}
+
+/*
+  std::vector<std::vector<MVertex*> > self, other;
+  
+  self.resize(1);
+  for(unsigned int i = 0; i < input.size(); i++){
+    if(input[i]->onWhat() == gf){
+      if(i && input[i]->onWhat() != input[i - 1]->onWhat())
+	self.resize(self.size() + 1);
+      self[self.size() - 1].push_back(input[i]);
+    }
+  }
+
+  other.resize(1);
+  for(unsigned int i = 0; i < input.size(); i++){
+    if(input[i]->onWhat() != gf){
+      if(i && input[i]->onWhat() != input[i - 1]->onWhat())
+	other.resize(other.size() + 1);
+      other[other.size() - 1].push_back(input[i]);
+    }
+  }
+
+  if(self.size() == 1 || other.size() == 0){
+    // nothing special to do
+    output.resize(1);
+    for(unsigned int i = 0; i < input.size(); i++){
+      output[0].push_back(input[i]);
+    }
+  }
+  else if(self.size() == other.size()){
+
+    for(unsigned int i = 0; i < self.size(); i++)
+      debugVertices(self[i], "self", false, i);
+    for(unsigned int i = 0; i < other.size(); i++)
+      debugVertices(other[i], "other", false, i);
+
+
+    output.resize(self.size());
+    for(unsigned int i = 0; i < self.size(); i++){
+      // add self pnts
+      for(unsigned int j = 0; j < self[i].size(); j++)
+	output[i].push_back(self[i][j]);
+      // check which other is closest
+      int which = 0;
+      double dist = 1e20;
+      for(unsigned int j = 0; j < other.size(); j++){
+	if(output[i][output[i].size() - 1]->distance(other[j][0]) < dist)
+	  which = j;
+      }
+      for(unsigned int j = 0; j < other[which].size(); j++)
+	output[i].push_back(other[which][j]);
+      output[i].push_back(output[i][0]);
+    }
+  }
+  else{
+    Msg(GERROR, "General subloop creation not implemented yet");
+  }
+}
+*/
+
+void meshGroutWithoutHole(GFace *gf,
+			  std::vector<std::vector<MVertex*> > &inside)
+{
+  std::vector<MVertex*> hole, tmp;
+  std::vector<std::vector<MVertex*> > loops;
+
+  for(unsigned int i = 0; i < inside.size(); i++)
+    for(unsigned int j = 0; j < inside[i].size(); j++)
+      tmp.push_back(inside[i][j]);
+  tmp.push_back(tmp[0]);
+
+  if(removeHardEdges(gf, tmp, loops)){
+    for(unsigned int i = 0; i < loops.size(); i++)
+      meshGrout(gf, loops[i], hole);
+  }
+  else{
+    meshGrout(gf, tmp, hole);
+  }
+}
+
 class createGrout{
 public:
   void operator() (GFace *gf)
@@ -392,79 +552,17 @@ public:
     for(unsigned int i = 0; i < outside.size(); i++)
       Msg(INFO, "   outside loop %d intersect has %d part(s)", i, (int)outside[i].size());
 
-    std::vector<MVertex*> hole, loop;
-    
     if(inside.size() == 1 && inside[0].size() == 1){
-      Msg(INFO, "CASE 1: MESH AROUND ONE-PART HOLE");
-      // create hole
-      SPoint2 ic(0., 0.);
-      {
-	for(unsigned int i = 0; i < inside[0][0].size(); i++){
-	  hole.push_back(inside[0][0][i]);
-	  double u, v;
-	  inside[0][0][i]->getParameter(0, u);
-	  inside[0][0][i]->getParameter(1, v);
-	  ic += SPoint2(u, v);
-	}
-	ic *= 1. / (double)inside[0][0].size();
-	hole.push_back(hole[0]);
-      }
-      // order exterior parts and create exterior loop
-      {
-	std::vector<std::pair<double, int> > angle;
-	std::map<int, std::vector<MVertex*> > outside_flat;
-	int part = 0;
-	for(unsigned int i = 0; i < outside.size(); i++){
-	  for(unsigned int j = 0; j < outside[i].size(); j++){
-	    for(unsigned int k = 0; k < outside[i][j].size(); k++){
-	      outside_flat[part].push_back(outside[i][j][k]);
-	    }
-	    part++;
-	  }
-	}
-	std::map<int, std::vector<MVertex*> >::iterator it = outside_flat.begin();
-	for(; it != outside_flat.end(); it++){
-	  SPoint2 oc(0., 0.);
-	  for(unsigned int i = 0; i < it->second.size(); i++){
-	    double u, v;
-	    it->second[i]->getParameter(0, u);
-	    it->second[i]->getParameter(1, v);
-	    oc += SPoint2(u, v);
-	  }
-	  oc *= 1. / (double)it->second.size();
-	  double a = atan2(oc[1] - ic[1], oc[0] - ic[0]);
-	  angle.push_back(std::make_pair(a, it->first));
-	}
-	std::sort(angle.begin(), angle.end());
-	for(unsigned int i = 0; i < angle.size(); i++){
-	  for(unsigned int j = 0; j < outside_flat[angle[i].second].size(); j++)
-	    loop.push_back(outside_flat[angle[i].second][j]);
-	}
-	loop.push_back(hole[0]);
-      }
-      // mesh the grout
-      meshGrout(gf, loop, hole);
+      Msg(INFO, "*** CASE 1: grout has a single hole in one part ***");
+      meshGroutWithHole(gf, inside, outside);
     }
     else if(!outside.size()){
-      Msg(INFO, "CASE 2: MESH SIMPLE HOLES");
-      hole.clear();
-      for(unsigned int i = 0; i < inside.size(); i++){
-	loop.clear();
-	for(unsigned int j = 0; j < inside[i].size(); j++){
-	  for(unsigned int k = 0; k < inside[i][j].size(); k++){
-	    loop.push_back(inside[i][j][k]);
-	  }
-	}
-	loop.push_back(loop[0]);
-	meshGrout(gf, loop, hole);
-      }
+      Msg(INFO, "*** CASE 2: grout has no outside contributions ***");
+      for(unsigned int i = 0; i < inside.size(); i++)
+	meshGroutWithoutHole(gf, inside[i]);
     }
     else{
-      // num individual meshes = num parts with onWhat() == gf!
-      // for each one
-      //   - find other parts that are "close" (using POUs?)
-      //   - sort parts w.r.t barycenter of each group      
-      Msg(GERROR, "This case is not implemented yet");
+      Msg(WARNING, "*** TODO: UNKNOWN CASE ***");
     }
   }
 };
@@ -496,6 +594,15 @@ fourierModel::fourierModel(const std::string &name)
 
   // remove any duplicate vertices on hard edges
 
+  // Here's an alternative approach that might be worth investigating:
+  // - compute and store the pou of each overlapping patch in the nodes of
+  //   all the patches
+  // - for each pair of overlapping patches, find the line pou1=pou2 by
+  //   interpolation on the overlapping grids
+  // - compute the intersections of these lines
+  // This should define a non-overlapping partitioning of the grid, which
+  // could be used as the boundary constrain for the unstructured algo
+
   CTX.mesh.changed = ENT_ALL;
 }
 
diff --git a/doc/TODO b/doc/TODO
index 97a260c683e43ef062aa63c73a5350739a5685d6..48ddfd5c272b727b43c7f08911f4e52b28f7ebcd 100644
--- a/doc/TODO
+++ b/doc/TODO
@@ -1,4 +1,9 @@
-$Id: TODO,v 1.19 2006-09-14 15:48:35 geuzaine Exp $
+$Id: TODO,v 1.20 2006-09-22 19:28:51 geuzaine Exp $
+
+********************************************************************
+
+- option to fl_beep() when mesh is done
+- warning+error counter (display counters on Exit())
 
 ********************************************************************