diff --git a/Common/Context.h b/Common/Context.h
index b7a21401b76def9be073398e0c87b4d358159ec8..97e0d304ac09df8a3a770700d574568a70e9609f 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -136,7 +136,8 @@ public :
 
   int forced_bbox; // dynamic variable tracking if the bbox is currently imposed
 
-  int enable_mouse_selection; // enable selection using the mouse
+  int mouse_selection; // enable selection using the mouse
+  int mouse_hover_meshes; // enable mouse hover on meshes
   int pick_elements; // allow individual element picking
 
   int expert_mode; // to disable some warnings for beginners
@@ -178,7 +179,8 @@ public :
     int dual;
     int light, light_two_side, light_lines;
     int format, nb_smoothing, algo2d, algo3d, algo_recombine;
-    int order, second_order_linear, second_order_incomplete, smooth_internal_edges, c1_continuity;
+    int order, second_order_linear, second_order_incomplete;
+    int smooth_internal_edges, c1_continuity;
     int min_circ_points;
     int constrained_bgmesh, lc_from_curvature;
     double normals, tangents, explode;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 843f7400df55cbdbcf69f8c714f4d4a2f44aa399..20ec21abbc071d022b08b46e62ba903defecf802 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -671,8 +671,10 @@ StringXNumber GeneralOptions_Number[] = {
     "Minimum model coordinate along the Y-axis (read-only)" }, 
   { F,   "MinZ" , opt_general_zmin , 0. , 
     "Minimum model coordinate along the Z-axis (read-only)" }, 
+  { F|O, "MouseHoverMeshes" , opt_general_mouse_hover_meshes , 0. ,
+    "Enable mouse hover on meshes" },
   { F|O, "MouseSelection" , opt_general_mouse_selection , 1. ,
-    "Mouse hover and selection mode (0=none, 1=hover and select geometry but only select mesh, 2=hover and select geometry and mesh)" },
+    "Enable mouse selection" },
 
   { F|O, "NoPopup" , opt_general_nopopup , 0. , 
     "Disable interactive dialog windows in scripts (and use default values instead)" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index fc0d535f5a9f4848c8fef9202bb167b6854e3648..3601f2247c024b9225ac6c031a2d97cbfa320cd8 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.350 2007-08-28 08:38:47 geuzaine Exp $
+// $Id: Options.cpp,v 1.351 2007-09-03 20:09:13 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -2801,28 +2801,32 @@ double opt_general_orthographic(OPT_ARGS_NUM)
 double opt_general_mouse_selection(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
-    CTX.enable_mouse_selection = (int)val;
+    CTX.mouse_selection = (int)val;
 #if defined(HAVE_FLTK)
   if(WID && (action & GMSH_GUI)) {
-    switch(CTX.enable_mouse_selection){
-    case 1:
+    if(CTX.mouse_selection){
       if(!CTX.batch) Msg(STATUS2N, "Mouse selection ON");
       WID->g_status_butt[9]->color(FL_BACKGROUND_COLOR);
-      break;
-    case 2:
-      if(!CTX.batch) Msg(STATUS2N, "Mouse selection ON (with mesh hover)");
-      WID->g_status_butt[9]->color(FL_GREEN);
-      break;
-    case 0:
-    default:
+    }
+    else{
       if(!CTX.batch) Msg(STATUS2N, "Mouse selection OFF");
       WID->g_status_butt[9]->color(FL_RED);
-      break;
     }
     WID->g_status_butt[9]->redraw();
   }
 #endif
-  return CTX.enable_mouse_selection;
+  return CTX.mouse_selection;
+}
+
+double opt_general_mouse_hover_meshes(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.mouse_hover_meshes = (int)val;
+#if defined(HAVE_FLTK)
+  if(WID && (action & GMSH_GUI))
+    WID->gen_butt[11]->value(CTX.mouse_hover_meshes);
+#endif
+  return CTX.mouse_hover_meshes;
 }
 
 double opt_general_fast_redraw(OPT_ARGS_NUM)
diff --git a/Common/Options.h b/Common/Options.h
index 1b6992b9ab6e63e44a1fb4d0f3f65836e2ecc411..729ad253bbfc27046e9bf0abc67a88ed4c6a91b3 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -293,6 +293,7 @@ double opt_general_tooltips(OPT_ARGS_NUM);
 double opt_general_confirm_overwrite(OPT_ARGS_NUM);
 double opt_general_orthographic(OPT_ARGS_NUM);
 double opt_general_mouse_selection(OPT_ARGS_NUM);
+double opt_general_mouse_hover_meshes(OPT_ARGS_NUM);
 double opt_general_draw_bounding_box(OPT_ARGS_NUM);
 double opt_general_fast_redraw(OPT_ARGS_NUM);
 double opt_general_xmin(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index e30273017a22f218d156f4438d6d407ca03a5deb..09c8d33a0909ce1169db98ba077d379c568ecdfc 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.539 2007-08-22 21:02:11 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.540 2007-09-03 20:09:13 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -480,20 +480,12 @@ void status_xyz1p_cb(CALLBACK_ARGS)
     WID->create_message_window();
   }
   else if(!strcmp(str, "S")){ // mouse selection
-    if(Fl::event_state(FL_SHIFT)){
-      // mouse hover and select for geometry and mesh
-      opt_general_mouse_selection(0, GMSH_SET | GMSH_GUI, 2);
-    }
-    else if(CTX.enable_mouse_selection){
-      // mouse does nothing
+    if(CTX.mouse_selection){
       opt_general_mouse_selection(0, GMSH_SET | GMSH_GUI, 0);
       WID->g_opengl_window->cursor(FL_CURSOR_DEFAULT, FL_BLACK, FL_WHITE);
     }
-    else{
-      // mouse hover and select for geometry, but mouse select only
-      // for mesh (default, for performance reasons)
+    else
       opt_general_mouse_selection(0, GMSH_SET | GMSH_GUI, 1);
-    }
   }
   WID->update_manip_window();
 }
@@ -942,6 +934,7 @@ void general_options_ok_cb(CALLBACK_ARGS)
   opt_general_axes_auto_position(0, GMSH_SET, WID->gen_butt[0]->value());
   opt_general_small_axes(0, GMSH_SET, WID->gen_butt[1]->value());
   opt_general_fast_redraw(0, GMSH_SET, WID->gen_butt[2]->value());
+  opt_general_mouse_hover_meshes(0, GMSH_SET, WID->gen_butt[11]->value());
   if(opt_general_double_buffer(0, GMSH_GET, 0) != WID->gen_butt[3]->value())
     opt_general_double_buffer(0, GMSH_SET, WID->gen_butt[3]->value());
   if(opt_general_antialiasing(0, GMSH_GET, 0) != WID->gen_butt[12]->value())
@@ -3904,7 +3897,7 @@ void mesh_optimize_cb(CALLBACK_ARGS)
     return;
   }
   CTX.threads_lock = 1;
-  OptimizeMesh();
+  OptimizeMesh(GModel::current());
   CTX.threads_lock = 0;
   CTX.mesh.changed = ENT_LINE | ENT_SURFACE | ENT_VOLUME;
   Draw();
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index ba4e9dabb44596f8a3e08352773835c29e8bfd1c..f94966fe45bb0a279eded5bfa8646f5c9ca1d247 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.633 2007-08-28 08:38:47 geuzaine Exp $
+// $Id: GUI.cpp,v 1.634 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1843,30 +1843,34 @@ void GUI::create_option_window()
       gen_butt[2]->type(FL_TOGGLE_BUTTON);
       gen_butt[2]->callback(general_options_ok_cb, (void*)"fast_redraw");
 
-      gen_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 5 * BH, BW, BH, "Enable double buffering");
+      gen_butt[11] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 5 * BH, BW, BH, "Enable mouse hover over meshes");
+      gen_butt[11]->type(FL_TOGGLE_BUTTON);
+      gen_butt[11]->callback(general_options_ok_cb);
+
+      gen_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 6 * BH, BW, BH, "Enable double buffering");
       gen_butt[3]->type(FL_TOGGLE_BUTTON);
       gen_butt[3]->callback(general_options_ok_cb);
 
-      gen_butt[12] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 6 * BH, BW, BH, "Enable antialiasing");
+      gen_butt[12] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 7 * BH, BW, BH, "Enable antialiasing");
       gen_butt[12]->type(FL_TOGGLE_BUTTON);
       gen_butt[12]->callback(general_options_ok_cb);
 
-      gen_butt[5] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 7 * BH, BW, BH, "Use trackball rotation mode instead of Euler angles");
+      gen_butt[5] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 8 * BH, BW, BH, "Use trackball rotation mode instead of Euler angles");
       gen_butt[5]->type(FL_TOGGLE_BUTTON);
       gen_butt[5]->callback(general_options_ok_cb);
 
-      gen_butt[15] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 8 * BH, BW, BH, "Rotate around pseudo center of mass");
+      gen_butt[15] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 9 * BH, BW, BH, "Rotate around pseudo center of mass");
       gen_butt[15]->type(FL_TOGGLE_BUTTON);
       gen_butt[15]->callback(general_options_ok_cb, (void*)"rotation_center");
 
-      gen_push_butt[0] = new Fl_Button(L + 2 * IW - 2 * WB, 2 * WB + 9 * BH, BB, BH, "Select");
+      gen_push_butt[0] = new Fl_Button(L + 2 * IW - 2 * WB, 2 * WB + 10 * BH, BB, BH, "Select");
       gen_push_butt[0]->callback(general_options_rotation_center_select_cb);
 
-      gen_value[8] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 9 * BH, IW / 3, BH);
+      gen_value[8] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 10 * BH, IW / 3, BH);
       gen_value[8]->callback(general_options_ok_cb, (void*)"rotation_center_coord");
-      gen_value[9] = new Fl_Value_Input(L + 2 * WB + IW / 3, 2 * WB + 9 * BH, IW / 3, BH);
+      gen_value[9] = new Fl_Value_Input(L + 2 * WB + IW / 3, 2 * WB + 10 * BH, IW / 3, BH);
       gen_value[9]->callback(general_options_ok_cb, (void*)"rotation_center_coord");
-      gen_value[10] = new Fl_Value_Input(L + 2 * WB + 2 * IW / 3, 2 * WB + 9 * BH, IW / 3, BH, "Rotation center");
+      gen_value[10] = new Fl_Value_Input(L + 2 * WB + 2 * IW / 3, 2 * WB + 10 * BH, IW / 3, BH, "Rotation center");
       gen_value[10]->align(FL_ALIGN_RIGHT);
       gen_value[10]->callback(general_options_ok_cb, (void*)"rotation_center_coord");
 
diff --git a/Fltk/Opengl_Window.cpp b/Fltk/Opengl_Window.cpp
index 4bcd1a027febfd37f34e9616c36f8a0f5ab14ce1..f258819b8a2506294f5e1b6327cab6463bd2602a 100644
--- a/Fltk/Opengl_Window.cpp
+++ b/Fltk/Opengl_Window.cpp
@@ -1,4 +1,4 @@
-// $Id: Opengl_Window.cpp,v 1.78 2006-12-17 22:23:16 geuzaine Exp $
+// $Id: Opengl_Window.cpp,v 1.79 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -120,7 +120,7 @@ void Opengl_Window::draw()
     glLoadIdentity();
     glColor3d(1., 1., 1.);
     glDisable(GL_DEPTH_TEST);
-    if(SelectionMode && CTX.enable_mouse_selection){
+    if(SelectionMode && CTX.mouse_selection){
       glEnable(GL_LINE_STIPPLE);
       glLineStipple(1, 0x0F0F);
     }
@@ -139,7 +139,7 @@ void Opengl_Window::draw()
       if(!i) lasso.set();
     }
     glDisable(GL_BLEND);
-    if(SelectionMode && CTX.enable_mouse_selection)
+    if(SelectionMode && CTX.mouse_selection)
       glDisable(GL_LINE_STIPPLE);
     glEnable(GL_DEPTH_TEST);
   }
@@ -202,7 +202,7 @@ int Opengl_Window::handle(int event)
       }
       else if(LassoMode) {
         LassoMode = false;
-	if(SelectionMode && CTX.enable_mouse_selection){
+	if(SelectionMode && CTX.mouse_selection){
 	  WID->try_selection = 2; // will try to select multiple entities
 	  WID->try_selection_xywh[0] = (int)(click.win[0] + curr.win[0])/2;
 	  WID->try_selection_xywh[1] = (int)(click.win[1] + curr.win[1])/2;
@@ -213,7 +213,7 @@ int Opengl_Window::handle(int event)
 	  lasso_zoom(click, curr);
 	}
       }
-      else if(CTX.enable_mouse_selection){
+      else if(CTX.mouse_selection){
         WID->try_selection = 1; // will try to select clicked entity
 	WID->try_selection_xywh[0] = (int)curr.win[0];
 	WID->try_selection_xywh[1] = (int)curr.win[1];
@@ -231,7 +231,7 @@ int Opengl_Window::handle(int event)
       }
       else if(LassoMode) {
 	LassoMode = false;
-	if(SelectionMode && CTX.enable_mouse_selection){
+	if(SelectionMode && CTX.mouse_selection){
 	  WID->try_selection = -2; // will try to unselect multiple entities
 	  WID->try_selection_xywh[0] = (int)(click.win[0] + curr.win[0])/2;
 	  WID->try_selection_xywh[1] = (int)(click.win[1] + curr.win[1])/2;
@@ -242,7 +242,7 @@ int Opengl_Window::handle(int event)
 	  lasso_zoom(click, curr);
 	}
       }
-      else if(CTX.enable_mouse_selection){
+      else if(CTX.mouse_selection){
         WID->try_selection = -1; // will try to unselect clicked entity
 	WID->try_selection_xywh[0] = (int)curr.win[0];
 	WID->try_selection_xywh[1] = (int)curr.win[1];
@@ -388,12 +388,12 @@ int Opengl_Window::handle(int event)
 	std::vector<GFace*> faces;
 	std::vector<GRegion*> regions;
 	std::vector<MElement*> elements;
-	bool something = ProcessSelectionBuffer(WID->selection, false,
-						CTX.enable_mouse_selection > 1,
-						(int)curr.win[0], (int)curr.win[1], 5, 5, 
-						vertices, edges, faces, regions,
-						elements);
-	if((WID->selection == ENT_ALL && something) ||
+	bool res = ProcessSelectionBuffer(WID->selection, false, 
+					  CTX.mouse_hover_meshes, 
+					  (int)curr.win[0], (int)curr.win[1], 5, 5, 
+					  vertices, edges, faces, regions,
+					  elements);
+	if((WID->selection == ENT_ALL && res) ||
 	   (WID->selection == ENT_POINT && vertices.size()) ||
 	   (WID->selection == ENT_LINE && edges.size()) || 
 	   (WID->selection == ENT_SURFACE && faces.size()) ||
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 300dd998cec817fb5fc7cae0c5ac3c34d148142d..25d1f22f99b0012476eef32a171491f0931b1df2 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.28 2007-08-29 07:57:39 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.29 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -376,6 +376,14 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
     else
       return InterpolateCubicSpline(v, t, c->mat, 0, t1, t2);
 
+  case MSH_SEGM_BND_LAYER:
+    Msg(GERROR, "Cannot interpolate boundary layer curve");
+    return V;
+
+  case MSH_SEGM_DISCRETE:
+    Msg(GERROR, "Cannot interpolate discrete curve");
+    return V;
+
   default:
     Msg(GERROR, "Unknown curve type in interpolation");
     return V;
@@ -600,6 +608,18 @@ Vertex InterpolateSurface(Surface * s, double u, double v, int derivee, int u_v)
 	T.Pos.Z = 0.;
       return T;
     }
+  case MSH_SURF_BND_LAYER:
+    {
+      Msg(GERROR, "Cannot interpolate boundary layer surface");
+      Vertex T(0., 0., 0.);
+      return T;
+    }    
+  case MSH_SURF_DISCRETE:
+    {
+      Msg(GERROR, "Cannot interpolate discrete surface");
+      Vertex T(0., 0., 0.);
+      return T;
+    }    
   default:
     {
       Msg(GERROR, "Unknown surface type in interpolation");
diff --git a/Geo/Makefile b/Geo/Makefile
index 7384705c6754e6a07c46e591f7b99c5e84514a84..0c9401378ebbd0ee3ba95c41025a0bc528c8b0c6 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.159 2007-09-01 14:27:55 geuzaine Exp $
+# $Id: Makefile,v 1.160 2007-09-03 20:09:14 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -81,7 +81,8 @@ GVertex.o: GVertex.cpp GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h ../Common/GmshDefines.h MVertex.h GPoint.h SPoint2.h \
   GFace.h GEdgeLoop.h GEdge.h SVector3.h MElement.h MEdge.h \
   ../Common/Hash.h MFace.h ../Numeric/Numeric.h ../Common/Context.h \
-  ../DataStr/List.h ExtrudeParams.h ../Common/SmoothData.h Pair.h
+  ../DataStr/List.h ExtrudeParams.h ../Common/SmoothData.h Pair.h \
+  ../Common/Message.h
 GEdge.o: GEdge.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 MFace.h \
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 22eb6c40827acf7484848aaf181da5a9d9b01374..d0df3bd6aa515b473305c58e0351682537c9e339 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.138 2007-08-21 19:05:39 geuzaine Exp $
+// $Id: Geom.cpp,v 1.139 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -308,13 +308,24 @@ class drawGFace {
       std::list<GEdge*> edges = f->edges();
       SBoundingBox3d bb;
       for(std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); it++){
-	Range<double> t_bounds = (*it)->parBounds(0);
-	GPoint p[3] = {(*it)->point(t_bounds.low()),
-		       (*it)->point(0.5 * (t_bounds.low() + t_bounds.high())),
-		       (*it)->point(t_bounds.high())};
-	for(int i = 0; i < 3; i++){
-	  SPoint2 uv = f->parFromPoint(SPoint3(p[i].x(), p[i].y(), p[i].z()));
-	  bb += SPoint3(uv.x(), uv.y(), 0.);
+	GEdge *ge = *it;
+	if(ge->geomType() == GEntity::DiscreteCurve || 
+	   ge->geomType() == GEntity::BoundaryLayerCurve){
+	  // don't try again
+	  f->cross.clear();
+	  f->cross.push_back(SPoint3(0., 0., 0.));
+	  CTX.threads_lock = 0;
+	  return;
+	}
+	else{
+	  Range<double> t_bounds = ge->parBounds(0);
+	  GPoint p[3] = {ge->point(t_bounds.low()),
+			 ge->point(0.5 * (t_bounds.low() + t_bounds.high())),
+			 ge->point(t_bounds.high())};
+	  for(int i = 0; i < 3; i++){
+	    SPoint2 uv = f->parFromPoint(SPoint3(p[i].x(), p[i].y(), p[i].z()));
+	    bb += SPoint3(uv.x(), uv.y(), 0.);
+	  }
 	}
       }
       bb *= 1.1;
diff --git a/Mesh/BoundaryLayer.cpp b/Mesh/BoundaryLayer.cpp
index 969cbe761af96647aef3f2be810a7a5e98ae9292..f50b50eea56c3e01a331e292c9a6119ed99cf233 100644
--- a/Mesh/BoundaryLayer.cpp
+++ b/Mesh/BoundaryLayer.cpp
@@ -1,4 +1,4 @@
-// $Id: BoundaryLayer.cpp,v 1.3 2007-09-03 12:00:28 geuzaine Exp $
+// $Id: BoundaryLayer.cpp,v 1.4 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -46,74 +46,73 @@ static void addExtrudeNormals(std::vector<T*> &elements)
   }
 }
 
-int MeshBoundaryLayerFaces(GModel *m)
+int Mesh2DWithBoundaryLayers(GModel *m)
 {
-  bool haveBoundaryLayers = false;
+  ExtrudeParams *ep = 0;
+
+  std::set<GFace*> sourceFaces, otherFaces;
+  std::set<GEdge*> sourceEdges, otherEdges;
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
     GFace *gf = *it;
     if(gf->geomType() == GEntity::BoundaryLayerSurface){
-      haveBoundaryLayers = true;
-      break;
+      ep = gf->meshAttributes.extrude;
+      if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
+	GFace *from = m->faceByTag(std::abs(ep->geo.Source));
+	if(!from){
+	  Msg(GERROR, "Unknown source face %d for boundary layer", ep->geo.Source);
+	  return 0;
+	}
+	sourceFaces.insert(from);
+	std::list<GEdge*> e = from->edges();
+	sourceEdges.insert(e.begin(), e.end());
+      }
     }
   }
-  if(!haveBoundaryLayers) return 0;
+  if(sourceFaces.empty()) return 0;
 
-  // make sure the surface mesh is oriented correctly (normally we do
-  // this only after the 3D mesh is done; but here it's critical since
-  // we use the normals for the extrusion)
-  std::for_each(m->firstFace(), m->lastFace(), orientMeshGFace());
+  // compute set of non-source edges and faces
+  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++)
+    if(sourceEdges.find(*it) == sourceEdges.end()) otherEdges.insert(*it);
+  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++)
+    if(sourceFaces.find(*it) == sourceFaces.end()) otherFaces.insert(*it);
 
-  // compute a normal field for the extrusion
+  // mesh source surfaces
+  std::for_each(sourceFaces.begin(), sourceFaces.end(), meshGFace());  
+  
+  // make sure the source surfaces for the boundary layers are
+  // oriented correctly (normally we do this only after the 3D mesh is
+  // done; but here it's critical since we use the normals for the
+  // extrusion)
+  std::for_each(sourceFaces.begin(), sourceFaces.end(), orientMeshGFace());
+
+  // compute a normal field for all the source faces
   if(ExtrudeParams::normals) delete ExtrudeParams::normals;
   ExtrudeParams::normals = new smooth_data();
-  ExtrudeParams *myep = 0;
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
+  for(std::set<GFace*>::iterator it = sourceFaces.begin(); 
+      it != sourceFaces.end(); it++){
     GFace *gf = *it;
-    if(gf->geomType() == GEntity::BoundaryLayerSurface){
-      ExtrudeParams *ep = myep = gf->meshAttributes.extrude;
-      if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
-	if(ep->mesh.ViewIndex >= 0 && ep->mesh.ViewIndex < List_Nbr(CTX.post.list)){ 
-	  // use external vector point post-pro view to get normals
-	  // FIXME: should use an octree on a general view instead
-	  xyzv::eps = 1.e-4;
-	  Post_View *v = *(Post_View**)List_Pointer(CTX.post.list, ep->mesh.ViewIndex);
-	  if(v->NbVP){
-	    int nb = List_Nbr(v->VP) / v->NbVP;
-	    for(int i = 0; i < List_Nbr(v->VP); i += nb){
-	      double *data = (double*)List_Pointer_Fast(v->VP, i);
-	      ExtrudeParams::normals->add(data[0], data[1], data[2], 3, &data[3]);
-	    }
-	  }
-	}
-	else{ 
-	  // compute smooth normal field from surfaces
-	  GFace *from = gf->model()->faceByTag(std::abs(ep->geo.Source));
-	  if(!from){
-	    Msg(GERROR, "Unknown source face %d for boundary layer", ep->geo.Source);
-	    continue;
-	  }
-	  addExtrudeNormals(from->triangles);
-	  addExtrudeNormals(from->quadrangles);
-	}
-      }
-    }
+    addExtrudeNormals(gf->triangles);
+    addExtrudeNormals(gf->quadrangles);
   }
   ExtrudeParams::normals->normalize();
-  //ExtrudeParams::normals->exportview("normals.pos");
-  if(!myep) return 0;
 
-  // set the position of bounding points (FIXME: should check
-  // coherence of all extrude parameters)
+  // set the position of boundary layer points using the smooth normal
+  // field
   for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); it++){
     GVertex *gv = *it;
     if(gv->geomType() == GEntity::BoundaryLayerPoint){
       GPoint p = gv->point();
-      myep->Extrude(myep->mesh.NbLayer - 1, myep->mesh.NbElmLayer[myep->mesh.NbLayer - 1],
-		    p.x(), p.y(), p.z());
+      ep->Extrude(ep->mesh.NbLayer - 1, ep->mesh.NbElmLayer[ep->mesh.NbLayer - 1],
+		  p.x(), p.y(), p.z());
       gv->setPosition(p);
     }
   }
-  
+
+  // remesh non-source edges (since they might have been modified by
+  // the change in boundary layer points)
+  std::for_each(otherFaces.begin(), otherFaces.end(), deMeshGFace());
+  std::for_each(otherEdges.begin(), otherEdges.end(), meshGEdge());
+
   // mesh the curves bounding the boundary layers by extrusion using
   // the smooth normal field
   for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
@@ -126,6 +125,9 @@ int MeshBoundaryLayerFaces(GModel *m)
     }
   }
 
+  // mesh non-source surfaces
+  std::for_each(otherFaces.begin(), otherFaces.end(), meshGFace());  
+
   // mesh the surfaces bounding the boundary layers by extrusion using
   // the smooth normal field
   for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
diff --git a/Mesh/BoundaryLayer.h b/Mesh/BoundaryLayer.h
index e7cce6b78b3563c0b37fbce9b83e9c1b288a9cf6..2691ac51f50cac0a5b81966fe3d6128aea7bbe8c 100644
--- a/Mesh/BoundaryLayer.h
+++ b/Mesh/BoundaryLayer.h
@@ -22,6 +22,6 @@
 
 #include "GModel.h"
 
-int MeshBoundaryLayerFaces(GModel *m);
+int Mesh2DWithBoundaryLayers(GModel *m);
 
 #endif
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 3dcf3e08c1d1de0c770b6f01abf82d2d4c62ceec..e22f5a98994f6d3d93b6d799faf4d834a662cc14 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -1,4 +1,4 @@
-// $Id: Generator.cpp,v 1.121 2007-08-21 19:05:40 geuzaine Exp $
+// $Id: Generator.cpp,v 1.122 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -180,9 +180,8 @@ void GetStatistics(double stat[50], double quality[3][100])
 
 }
 
-bool TooManyElements(int dim)
+bool TooManyElements(GModel *m, int dim)
 {
-  GModel *m = GModel::current();
   if(CTX.expert_mode || !m->numVertex()) return false;
 
   // try to detect obvious mistakes in characteristic lenghts (one of
@@ -201,13 +200,12 @@ bool TooManyElements(int dim)
   return false;
 }
 
-void Mesh1D()
+void Mesh1D(GModel *m)
 {
-  if(TooManyElements(1)) return;
+  if(TooManyElements(m, 1)) return;
   Msg(STATUS1, "Meshing 1D...");
   double t1 = Cpu();
 
-  GModel *m = GModel::current();
   std::for_each(m->firstEdge(), m->lastEdge(), meshGEdge());
 
   double t2 = Cpu();
@@ -216,9 +214,9 @@ void Mesh1D()
   Msg(STATUS1, "Mesh");
 }
 
-void Mesh2D()
+void Mesh2D(GModel *m)
 {
-  if(TooManyElements(2)) return;
+  if(TooManyElements(m, 2)) return;
 
   if(CTX.mesh.algo2d == ALGO_2D_DELAUNAY && !CTX.expert_mode){
     if(!GetBinaryAnswer("The 2D Delaunay algorithm is still highly experimental\n"
@@ -232,13 +230,11 @@ void Mesh2D()
   Msg(STATUS1, "Meshing 2D...");
   double t1 = Cpu();
 
-  GModel *m = GModel::current();
-  std::for_each(m->firstFace(), m->lastFace(), meshGFace());  
-
-  // boundary layers are special: their definition (including vertices
-  // and curve meshes) relies on the surface mesh--and it is global
-  // since we use a smooth normal field
-  MeshBoundaryLayerFaces(m);
+  // boundary layers are special: their generation (including vertices
+  // and curve meshes) is global as it depends on a smooth normal
+  // field generated from the surface mesh of the source surfaces
+  if(!Mesh2DWithBoundaryLayers(m))
+    std::for_each(m->firstFace(), m->lastFace(), meshGFace());  
 
   double t2 = Cpu();
   CTX.mesh_timer[1] = t2 - t1;
@@ -254,14 +250,12 @@ void FindConnectedRegions(std::vector<GRegion*> &delaunay,
   connected.push_back(delaunay);
 }
 
-void Mesh3D()
+void Mesh3D(GModel *m)
 {
-  if(TooManyElements(3)) return;
+  if(TooManyElements(m, 3)) return;
   Msg(STATUS1, "Meshing 3D...");
   double t1 = Cpu();
 
-  GModel *m = GModel::current();
-
   // mesh the extruded volumes first
   std::for_each(m->firstRegion(), m->lastRegion(), meshGRegionExtruded());
 
@@ -287,12 +281,11 @@ void Mesh3D()
   Msg(STATUS1, "Mesh");
 }
 
-void OptimizeMesh()
+void OptimizeMesh(GModel *m)
 {
   Msg(STATUS1, "Optimizing 3D...");
   double t1 = Cpu();
 
-  GModel *m = GModel::current();
   std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegion());
 
   double t2 = Cpu();
@@ -319,18 +312,18 @@ void GenerateMesh(int ask)
   if(ask == 1 || (ask > 1 && old < 1)) {
     std::for_each(m->firstRegion(), m->lastRegion(), deMeshGRegion());
     std::for_each(m->firstFace(), m->lastFace(), deMeshGFace());
-    Mesh1D();
+    Mesh1D(m);
   }
 
   // 2D mesh
   if(ask == 2 || (ask > 2 && old < 2)) {
     std::for_each(m->firstRegion(), m->lastRegion(), deMeshGRegion());
-    Mesh2D();
+    Mesh2D(m);
   }
 
   // 3D mesh
   if(ask == 3) {
-    Mesh3D();
+    Mesh3D(m);
   }
 
   // Orient the surface mesh so that it matches the geometry
@@ -339,7 +332,7 @@ void GenerateMesh(int ask)
   
   // Optimize quality
   if(m->getMeshStatus() == 3 && CTX.mesh.optimize)
-    OptimizeMesh();
+    OptimizeMesh(m);
   
   // Create high order elements
   if(m->getMeshStatus() && CTX.mesh.order > 1) 
diff --git a/Mesh/Generator.h b/Mesh/Generator.h
index 0faf6a651adf3324d22a91954721eb4656ae809b..ad54f320e1c0b87690f6dfd8616162f0bcc942de 100644
--- a/Mesh/Generator.h
+++ b/Mesh/Generator.h
@@ -20,8 +20,10 @@
 // 
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
+class GModel;
+
 void GetStatistics(double stat[50], double quality[3][100]=0);
 void GenerateMesh(int dimension);
-void OptimizeMesh();
+void OptimizeMesh(GModel *m);
 
 #endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 3324147d221ff5bd6e49ea1588a00ecf30a3278d..63586b8452b1e9cc9ee09b4a735360993f4b25e1 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.13 2007-04-26 09:47:38 remacle Exp $
+// $Id: HighOrder.cpp,v 1.14 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -217,7 +217,9 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
       MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);            
       double u0 = 0., u1 = 0.;
       bool reparamOK = true;
-      if(c1 || (!linear && ge->geomType() != GEntity::DiscreteCurve)){
+      if(c1 || (!linear && 
+		ge->geomType() != GEntity::DiscreteCurve &&
+		ge->geomType() != GEntity::BoundaryLayerCurve)){
 	reparamOK &= reparamOnEdge(v0, ge, u0);
 	reparamOK &= reparamOnEdge(v1, ge, u1);
       }
@@ -279,7 +281,9 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
       MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
       SPoint2 p0, p1;
       bool reparamOK = true;
-      if(!linear && gf->geomType() != GEntity::DiscreteSurface){
+      if(!linear && 
+	 gf->geomType() != GEntity::DiscreteSurface &&
+	 gf->geomType() != GEntity::BoundaryLayerSurface){
 	reparamOK &= reparamOnFace(v0, gf, p0);
 	reparamOK &= reparamOnFace(v1, gf, p1);
       }
@@ -352,7 +356,9 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
     else{
       SPoint2 p0, p1, p2, p3;
       bool reparamOK = true;
-      if(!linear && gf->geomType() != GEntity::DiscreteSurface){
+      if(!linear && 
+	 gf->geomType() != GEntity::DiscreteSurface &&
+	 gf->geomType() != GEntity::BoundaryLayerSurface){
 	reparamOK &= reparamOnFace(p[0], gf, p0);
 	reparamOK &= reparamOnFace(p[1], gf, p1);
 	reparamOK &= reparamOnFace(p[2], gf, p2);
diff --git a/Mesh/meshGEdgeExtruded.cpp b/Mesh/meshGEdgeExtruded.cpp
index 58f6560f1bc7040366e7562719af875ca666c11b..fb396836ae766d4fcce63c2ed12318ac5bb0e3eb 100644
--- a/Mesh/meshGEdgeExtruded.cpp
+++ b/Mesh/meshGEdgeExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdgeExtruded.cpp,v 1.6 2007-01-16 11:31:41 geuzaine Exp $
+// $Id: meshGEdgeExtruded.cpp,v 1.7 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2d61c1af1a9e78ed323ad14339f71d72e4d9df13..d0e94805a0f6bf5385fd0f04487c1067180669f7 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.82 2007-07-26 16:28:27 anand Exp $
+// $Id: meshGFace.cpp,v 1.83 2007-09-03 20:09:14 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -612,8 +612,6 @@ bool recover_medge ( BDS_Mesh *m, GEdge *ge)
 bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
 {
 
-  //if (gf->tag() != 16)return true;
-
   typedef std::set<MVertex*> v_container ;
   v_container all_vertices;
   std::map<int, MVertex*>numbered_vertices;
@@ -648,24 +646,28 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
 
   v_container::iterator itv = all_vertices.begin();
 
-  //  FILE *fdeb = fopen("debug.dat","w");
-  //  fprintf(fdeb,"surface %d\n" ,gf->tag());
+  //char tmp[256]; sprintf(tmp, "surf%d.pos", gf->tag());
+  //FILE *fdeb = fopen(tmp,"w");
+  //fprintf(fdeb,"View \"surf%d\"{\n" ,gf->tag());
+
   int count = 0;
   SBoundingBox3d bbox;
   while(itv != all_vertices.end()){
     MVertex *here = *itv;
     SPoint2 param;
-    if(here->onWhat()->dim() == 0){
+    if(here->onWhat()->geomType() == GEntity::DiscreteCurve ||
+       here->onWhat()->geomType() == GEntity::BoundaryLayerCurve){
+      param = gf->parFromPoint(SPoint3(here->x(), here->y(), here->z()));
+    }
+    else if(here->onWhat()->dim() == 0){
       GVertex *gv = (GVertex*)here->onWhat();
-      if (gv->edges().size() == 1)
-	{
-	  GEdge *ge = *(gv->edges().begin());
-	  Range<double> bb = ge->parBounds(0);
-	  param = ge->reparamOnFace(gf, bb.low(), 1);	  
-	}
+      if(gv->edges().size() == 1){
+	GEdge *ge = *(gv->edges().begin());
+	Range<double> bb = ge->parBounds(0);
+	param = ge->reparamOnFace(gf, bb.low(), 1);	  
+      }
       else
 	param = gv->reparamOnFace(gf,1);
-      
     }
     else if(here->onWhat()->dim() == 1){
       GEdge *ge = (GEdge*)here->onWhat();
@@ -680,7 +682,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
       else
 	param = gf->parFromPoint(SPoint3(here->x(), here->y(), here->z()));
     }
-    // fprintf(fdeb,"%d %g %g %g\n" ,here->getNum(),here->x(),here->y(),here->z());
+    //fprintf(fdeb,"SP(%g,%g,%g){%d};\n" ,here->x(),here->y(),here->z(),here->getNum());
+    //fprintf(fdeb,"SP(%g,%g,0){%d};\n" ,param.x(),param.y(),here->getNum());
     U_[count] = param.x();
     V_[count] = param.y();
     (*itv)->setNum(count);
@@ -690,7 +693,8 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
     ++itv;
   }
   
-  //  fclose (fdeb);
+  //fprintf(fdeb,"};\n");
+  //fclose(fdeb);
 
   // compute the bounding box in parametric space
   // I do not have SBoundingBox, so I use a 3D one...