diff --git a/Common/Context.h b/Common/Context.h
index 627a195d4c1cf41d4c4f53dd3b3c2fe1f045cc5e..3277b14eba14ca751ea30127024b15d53ca0c109 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -140,6 +140,8 @@ public :
 
   int expert_mode; // to disable some warnings for beginners
   int printing; // dynamic: equal to 1 while gmsh is printing
+  
+  int hide_unselected; // hide all unselected entities
 
   // geometry options 
   struct{
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 3741abb7d9f7333e736913c2dd91bd6edf0c77e1..0c8fd1a9a66227beff784709f90cfd310947311f 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.346 2007-06-22 08:07:45 geuzaine Exp $
+// $Id: Options.cpp,v 1.347 2007-07-26 13:10:47 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -152,6 +152,7 @@ void Init_Options(int num)
   CTX.gl_font_enum = -1;
 #endif
   CTX.forced_bbox = 0;
+  CTX.hide_unselected = 0;
 }
 
 void ReInit_Options(int num)
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 06e911c21ba0d0b2d36b85cf72eb5d30fc8a2db7..eb6c2d909dcfa5a8f1fecbdc610a79383e4cb611 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.534 2007-07-22 15:48:28 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.535 2007-07-26 13:10:47 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -3693,11 +3693,6 @@ void mesh_3d_cb(CALLBACK_ARGS)
   Msg(STATUS2N, " ");
 }
 
-void mesh_edit_cb(CALLBACK_ARGS)
-{
-  WID->set_context(menu_mesh_edit, 0);
-}
-
 void mesh_delete_cb(CALLBACK_ARGS)
 {
   WID->set_context(menu_mesh_delete, 0);
diff --git a/Fltk/Callbacks.h b/Fltk/Callbacks.h
index 343b53923cdf4e2f11ccef721e4fb7cc1ddd2684..c68d6600afefb0716c2f43c87f61c66cc097aedb 100644
--- a/Fltk/Callbacks.h
+++ b/Fltk/Callbacks.h
@@ -274,7 +274,6 @@ void mesh_define_cb(CALLBACK_ARGS);
 void mesh_1d_cb(CALLBACK_ARGS);
 void mesh_2d_cb(CALLBACK_ARGS); 
 void mesh_3d_cb(CALLBACK_ARGS); 
-void mesh_edit_cb(CALLBACK_ARGS);
 void mesh_delete_cb(CALLBACK_ARGS);
 void mesh_delete_parts_cb(CALLBACK_ARGS);
 void mesh_inspect_cb(CALLBACK_ARGS);
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index f7425c13a86ac3ca7669a7e946aea55c2c654ced..d2c79d56bdb2713726a4d20c8948635ec0cbac5b 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.625 2007-07-25 15:48:32 geuzaine Exp $
+// $Id: GUI.cpp,v 1.626 2007-07-26 13:10:47 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -281,7 +281,7 @@ Context_Item menu_mesh[] = {
   {"1Mesh", NULL} ,
   {"Define",       (Fl_Callback *)mesh_define_cb} ,
   {"Inspect",      (Fl_Callback *)mesh_inspect_cb} , 
-  {"Edit",         (Fl_Callback *)mesh_edit_cb} , 
+  {"Delete",       (Fl_Callback *)mesh_delete_cb} , 
   {"1D",           (Fl_Callback *)mesh_1d_cb} ,
   {"2D",           (Fl_Callback *)mesh_2d_cb} , 
   {"3D",           (Fl_Callback *)mesh_3d_cb} , 
@@ -290,26 +290,10 @@ Context_Item menu_mesh[] = {
 #if defined(HAVE_NETGEN)
   {"Optimize quality", (Fl_Callback *)mesh_optimize_cb} , 
 #endif
+  {"Reparameterize",   (Fl_Callback *)mesh_parameterize_cb} , 
   {"Save",         (Fl_Callback *)mesh_save_cb} ,
   {0} 
 };  
-    Context_Item menu_mesh_edit[] = {
-      {"1Mesh>Edit", NULL} ,
-      {"Delete", (Fl_Callback *)mesh_delete_cb} ,
-      //{"Update edges",   (Fl_Callback *)mesh_update_edges_cb} ,
-      {"Reparameterize", (Fl_Callback *)mesh_parameterize_cb} ,
-      //{"Remesh 2D",      (Fl_Callback *)mesh_remesh_cb} , 
-      {0} 
-    };  
-        Context_Item menu_mesh_delete[] = {
-          {"1Mesh>Edit>Delete", NULL} ,
-          {"Elements", (Fl_Callback *)mesh_delete_parts_cb, (void*)"elements"} ,
-	  {"Point",   (Fl_Callback *)geometry_elementary_delete_point_cb} ,
-          {"Lines",    (Fl_Callback *)mesh_delete_parts_cb, (void*)"lines"} ,
-          {"Surfaces", (Fl_Callback *)mesh_delete_parts_cb, (void*)"surfaces"} ,
-          {"Volumes",  (Fl_Callback *)mesh_delete_parts_cb, (void*)"volumes"} ,
-          {0} 
-        };  
     Context_Item menu_mesh_define[] = {
       {"1Mesh>Define", NULL} ,
       {"Characteristic length", (Fl_Callback *)mesh_define_length_cb  } ,
@@ -324,6 +308,14 @@ Context_Item menu_mesh[] = {
 	  {"Volume",  (Fl_Callback *)mesh_define_transfinite_volume_cb} , 
 	  {0} 
 	};  
+    Context_Item menu_mesh_delete[] = {
+      {"1Mesh>Edit>Delete", NULL} ,
+      {"Elements", (Fl_Callback *)mesh_delete_parts_cb, (void*)"elements"} ,
+      {"Lines",    (Fl_Callback *)mesh_delete_parts_cb, (void*)"lines"} ,
+      {"Surfaces", (Fl_Callback *)mesh_delete_parts_cb, (void*)"surfaces"} ,
+      {"Volumes",  (Fl_Callback *)mesh_delete_parts_cb, (void*)"volumes"} ,
+      {0} 
+    };  
 
 // FIXME: should create MAXSOLVERS items...
 Context_Item menu_solver[] = {
diff --git a/Fltk/GUI_Projection.cpp b/Fltk/GUI_Projection.cpp
index 9706436b5de81f1a24426ab671e42f513b3854f3..8bf0566ba0d797d87d89058b12268235e294916d 100644
--- a/Fltk/GUI_Projection.cpp
+++ b/Fltk/GUI_Projection.cpp
@@ -5,132 +5,293 @@
 #include "Draw.h"
 #include "Options.h"
 #include "Context.h"
+#include "SelectBuffer.h"
 #include "GUI.h"
 #include "Shortcut_Window.h"
 
+#include <FL/Fl_Toggle_Button.H>
+#include <FL/Fl_Round_Button.H>
+
 extern GModel *GMODEL;
 extern Context_T CTX;
 
-int projection_editor(char *title, projectionFace *p)
+void select_cb(Fl_Widget *w, void *data);
+void browse_cb(Fl_Widget *w, void *data);
+void update_cb(Fl_Widget *w, void *data);
+void close_cb(Fl_Widget *w, void *data);
+void action_cb(Fl_Widget *w, void *data);
+void compute_cb(Fl_Widget *w, void *data);
+
+class uvPlot : public Fl_Window{
+ private:
+  std::vector<double> _u, _v;
+  void draw()
+  {
+    // draw background
+    fl_color(FL_WHITE);
+    fl_rectf(0, 0, w(), h());
+    // draw points
+    fl_color(FL_BLACK);
+    if(_u.size() != _v.size()) return;
+    for(unsigned int i = 0; i < _u.size(); i++)
+      fl_rect(_u[i] * w(), _v[i] * h(), 3, 3);
+  }
+ public:
+  uvPlot(int x, int y, int w, int h, const char *l=0) : Fl_Window(x, y, w, h, l){}
+  void fill(std::vector<double> &u, std::vector<double> &v){ _u = u; _v = v; redraw(); }
+};
+
+class projectionEditor{
+ private:
+  std::vector<projectionFace*> _faces;
+  Fl_Window *_window;
+  Fl_Value_Input *_input[20];
+  Fl_Hold_Browser *_browser;
+  Fl_Round_Button *_select[3];
+  uvPlot *_uvPlot;
+ public:
+  projectionEditor(std::vector<projectionFace*> &faces) : _faces(faces)
+  {
+    const int BH = 2 * GetFontSize() + 1, BB = 7 * GetFontSize(), WB = 7;
+    const int width = 3.5 * BB, height = 18 * BH;
+    _window = new Dialog_Window(width, height, "Reparameterize");
+    new Fl_Box(WB, WB + BH, 0.5 * BB, BH, "Select:");
+    Fl_Group *o = new Fl_Group(WB, WB, 2 * BB, 3 * BH);
+    _select[0] = new Fl_Round_Button(2 * WB + 0.5 * BB, WB, BB, BH, "Points");
+    _select[0]->callback(select_cb, (void*)"Points");
+    _select[0]->value(1);
+    _select[1] = new Fl_Round_Button(2 * WB + 0.5 * BB, WB + BH, BB, BH, "Elements");
+    _select[1]->callback(select_cb, (void*)"Elements");
+    _select[2] = new Fl_Round_Button(2 * WB + 0.5 * BB, WB + 2 * BH, BB, BH, "Surfaces");
+    _select[2]->callback(select_cb, (void*)"Surfaces");
+    for(int i = 0; i < 3; i++) 
+      _select[i]->type(FL_RADIO_BUTTON);
+    o->end();
+    Fl_Toggle_Button *b1 = new Fl_Toggle_Button(width - WB - 1.5 * BB, WB, 1.5 * BB, BH, 
+						"Hide unselected");
+    b1->callback(action_cb, (void*)"Hide");
+    Fl_Button *b2 = new Fl_Button(width - WB - 1.5 * BB, WB + BH, 1.5 * BB, BH, 
+				  "Save selection");
+    b2->callback(action_cb, (void*)"Save");
+    _browser = new Fl_Hold_Browser(WB, 2 * WB + 3 * BH, BB, 5 * BH);
+    _browser->callback(browse_cb, this);
+    for(unsigned int i = 0; i < _faces.size(); i++){
+      _browser->add("test");
+    }
+    Fl_Scroll *s = new Fl_Scroll(2 * WB + BB, 2 * WB + 3 * BH, width - 3 * WB - BB, 5 * BH);
+    for(int i = 0; i < 20; i++){
+      _input[i] = new Fl_Value_Input(2 * WB + BB, 2 * WB + (i + 3) * BH, BB, BH);
+      _input[i]->align(FL_ALIGN_RIGHT);
+      _input[i]->callback(update_cb, this);
+      _input[i]->minimum(0.);
+      _input[i]->maximum(1.);
+      _input[i]->step(0.01);
+      //_input[i]->hide();
+    }
+    s->end();
+    _uvPlot = new uvPlot(WB, 3 * WB + 8 * BH, width - 2 * WB, height - 5 * WB - 9 * BH);
+    _uvPlot->end();
+    Fl_Button *b3 = new Fl_Button(width - 2 * WB - 2 * BB, height - WB - BH, BB, BH, 
+				  "Compute");
+    b3->callback(compute_cb, this);
+    Fl_Button *b4 = new Fl_Button(width - WB - BB, height - WB - BH, BB, BH, 
+				  "Cancel");
+    b4->callback(close_cb, _window);
+    _window->end();
+    _window->hotspot(_window);
+    _window->resizable(_uvPlot);
+    _window->size_range(width, 0.75 * height);
+  }
+  void show()
+  {
+    _window->show();
+    for(int i = 0; i < 3; i++)
+      if(_select[i]->value()) _select[i]->do_callback();
+  }
+  uvPlot *uv() { return _uvPlot; }
+  projectionFace *face() { return 0; } // return current projection face
+  void reset(){} // reset all parameters to default values
+  void save(){} // save all parameters to a file
+  void load(){}  // load parameters from file
+};
+
+void browse_cb(Fl_Widget *w, void *data)
 {
-  struct _editor{
-    Fl_Window *window;
-    Fl_Value_Input *scale[3], *rotation[3], *translation[3];
-    Fl_Button *cancel, *print;
-  };
-  static _editor *editor = 0;
-
-  const int BH = 2 * GetFontSize() + 1;
-  const int BB = 7 * GetFontSize();
-  const int WB = 7;
-  
+  projectionEditor *e = (projectionEditor*)data;
   SBoundingBox3d bounds = GMODEL->bounds();
 
-  if(!editor){
-    editor = new _editor;
-    editor->window = new Dialog_Window(2 * BB + 3 * WB, 10 * BH + 3 * WB);
-
-    editor->scale[0] = new Fl_Value_Input(WB, WB, BB, BH, "Scale Factor X");
-    editor->scale[1] = new Fl_Value_Input(WB, WB + BH, BB, BH, "Scale Factor Y");
-    editor->scale[2] = new Fl_Value_Input(WB, WB + 2 * BH, BB, BH, "Scale Factor Z");
-    for(int i = 0; i < 3; i++){
-      editor->scale[i]->align(FL_ALIGN_RIGHT);
-      editor->scale[i]->maximum(CTX.lc * 10.);
-      editor->scale[i]->minimum(CTX.lc / 100.);
-      editor->scale[i]->step(CTX.lc / 100.);
-    }
+  // change active projection face 
+  //e->face()->setVisibility(true);
 
-    editor->rotation[0] = new Fl_Value_Input(WB, WB + 3*BH, BB, BH, "Rotation X");
-    editor->rotation[1] = new Fl_Value_Input(WB, WB + 4*BH, BB, BH, "Rotation Y");
-    editor->rotation[2] = new Fl_Value_Input(WB, WB + 5*BH, BB, BH, "Rotation Z");
-    for(int i = 0; i < 3; i++){
-      editor->rotation[i]->align(FL_ALIGN_RIGHT);
-      editor->rotation[i]->maximum(-180.);
-      editor->rotation[i]->minimum(180.);
-      editor->rotation[i]->step(0.1);
-    }
+  // set and show parameters accordingly
 
-    editor->translation[0] = new Fl_Value_Input(WB, WB + 6*BH, BB, BH, "Translation X");
-    editor->translation[1] = new Fl_Value_Input(WB, WB + 7*BH, BB, BH, "Translation Y");
-    editor->translation[2] = new Fl_Value_Input(WB, WB + 8*BH, BB, BH, "Translation Z");
-    for(int i = 0; i < 3; i++){
-      editor->translation[i]->align(FL_ALIGN_RIGHT);    
-      editor->translation[i]->maximum(bounds.max()[i] + CTX.lc);
-      editor->translation[i]->minimum(bounds.min()[i] - CTX.lc);
-      editor->translation[i]->step(CTX.lc / 100.);
-    }
+}
 
-    editor->print = new Fl_Button(WB , 2 * WB + 9 * BH, BB, BH, "Print");
-    editor->cancel = new Fl_Button(2 * WB + BB, 2 * WB + 9 * BH, BB, BH, "Cancel");
+void update_cb(Fl_Widget *w, void *data)
+{
+  projectionEditor *e = (projectionEditor*)data;
+
+  // get all parameters from GUI and modify projection surface accordingly
+
+  // project all selected points and update u,v display
+  std::vector<double> u, v;
+  for(int i = 0; i < 5000; i++){
+    u.push_back((double)rand() / (double)RAND_MAX);
+    v.push_back((double)rand() / (double)RAND_MAX);
+  }
+  e->uv()->fill(u, v);
+}
+
+void select_cb(Fl_Widget *w, void *data)
+{
+  char *str = (char*)data;
+  int what;
 
-    editor->window->end();
-    editor->window->hotspot(editor->window);
+  if(!strcmp(str, "Points")){
+    CTX.pick_elements = 0;
+    what = ENT_POINT;
   }
-  
-  editor->window->label(title);
-  for(int i = 0; i < 3; i++){
-    editor->scale[i]->value(CTX.lc);
-    editor->rotation[i]->value(0.);
-    editor->translation[i]->value(0.5 * (bounds.max()[i] + bounds.min()[i]));
+  else if(!strcmp(str, "Elements")){
+    CTX.pick_elements = 1;
+    what = ENT_ALL;
   }
+  else if(!strcmp(str, "Surfaces")){
+    CTX.pick_elements = 0;
+    what = ENT_SURFACE;
+  }
+  else
+    return;
 
-  p->setScale(SVector3(editor->scale[0]->value(),
-		       editor->scale[1]->value(),
-		       editor->scale[2]->value()));
-  p->setRotation(SVector3(editor->rotation[0]->value(),
-			  editor->rotation[1]->value(),
-			  editor->rotation[2]->value()));
-  p->setTranslation(SVector3(editor->translation[0]->value(),
-			     editor->translation[1]->value(),
-			     editor->translation[2]->value()));
-  Draw();
+  std::vector<GVertex*> vertices;
+  std::vector<GEdge*> edges;
+  std::vector<GFace*> faces;
+  std::vector<GRegion*> regions;
+  std::vector<MElement*> elements;
+
+  std::vector<MElement*> ele;
+  std::vector<GEntity*> ent;
 
-  editor->window->show();
+  while(1) {
+    CTX.mesh.changed = ENT_ALL;
+    Draw();
 
-  while(editor->window->shown()){
-    Fl::wait();
-    for (;;) {
-      Fl_Widget* o = Fl::readqueue();
-      if(!o) break;
-      if(o == editor->window || o == editor->cancel){
-	editor->window->hide();
-      	return 0;
+    if(ele.size() || ent.size())
+      Msg(ONSCREEN, "Select %s\n"
+	  "[Press 'e' to end selection, 'u' to undo last selection or 'q' to abort]", str);
+    else
+      Msg(ONSCREEN, "Select %s\n"
+	  "[Press 'e' to end selection or 'q' to abort]", str);
+
+    char ib = SelectEntity(what, vertices, edges, faces, regions, elements);
+    if(ib == 'l') {
+      if(CTX.pick_elements){
+	for(unsigned int i = 0; i < elements.size(); i++){
+	  if(elements[i]->getVisibility() != 2){
+	    elements[i]->setVisibility(2); ele.push_back(elements[i]);
+	  }
+	}
+      }
+      else{
+	for(unsigned int i = 0; i < vertices.size(); i++){
+	  if(vertices[i]->getSelection() != 1){
+	    vertices[i]->setSelection(1); ent.push_back(vertices[i]);
+	  }
+	}
+	for(unsigned int i = 0; i < faces.size(); i++){
+	  if(faces[i]->getSelection() != 1){
+	    faces[i]->setSelection(1); ent.push_back(faces[i]);
+	  }
+	}
       }
-      else if(o == editor->print){
-	for(int i = 0; i < 3; i++) printf("%g\n", editor->scale[i]->value());
-	for(int i = 0; i < 3; i++) printf("%g\n", editor->rotation[i]->value());
-	for(int i = 0; i < 3; i++) printf("%g\n", editor->translation[i]->value());
+    }
+    if(ib == 'r') {
+      if(CTX.pick_elements){
+	for(unsigned int i = 0; i < elements.size(); i++)
+	  elements[i]->setVisibility(1);
       }
       else{
-	p->setScale(SVector3(editor->scale[0]->value(),
-			     editor->scale[1]->value(),
-			     editor->scale[2]->value()));
-	p->setRotation(SVector3(editor->rotation[0]->value(),
-				editor->rotation[1]->value(),
-				editor->rotation[2]->value()));
-	p->setTranslation(SVector3(editor->translation[0]->value(),
-				   editor->translation[1]->value(),
-				   editor->translation[2]->value()));
+	for(unsigned int i = 0; i < vertices.size(); i++)
+	  vertices[i]->setSelection(0);
+	for(unsigned int i = 0; i < faces.size(); i++)
+	  faces[i]->setSelection(0);
+      }
+    }
+    if(ib == 'u') {
+      if(CTX.pick_elements){
+	if(ele.size()){
+	  ele[ele.size() - 1]->setVisibility(1);
+	  ele.pop_back();
+	}
       }
-      Draw();
+      else{
+	if(ent.size()){
+	  ent[ent.size() - 1]->setSelection(0);
+	  ent.pop_back();
+	}
+      }
+    }
+    if(ib == 'e') {
+      ZeroHighlight();
+      ele.clear();
+      ent.clear();
+    }
+    if(ib == 'q') {
+      ZeroHighlight();
+      break;
     }
   }
-  return 0;
+
+  CTX.mesh.changed = ENT_ALL;
+  CTX.pick_elements = 0;
+  Draw();  
+  Msg(ONSCREEN, "");
 }
 
-void mesh_parameterize_cb(Fl_Widget* w, void* data)
+void close_cb(Fl_Widget *w, void *data)
+{
+  if(data) ((Fl_Window *) data)->hide();
+}
+
+void action_cb(Fl_Widget *w, void *data)
 {
-  projectionFace *p = new parabolicCylinder(GMODEL, 10000);
+  char *str = (char*)data;
+  if(!strcmp(str, "Hide")){
+    CTX.hide_unselected = !CTX.hide_unselected;
+    CTX.mesh.changed = ENT_ALL;
+  }
+  else if(!strcmp(str, "Save")){
+    Msg(GERROR, "Save not implemented yet!");
+  }
+  Draw();
+}
 
-  GMODEL->add(p);
+void compute_cb(Fl_Widget *w, void *data)
+{
+  projectionEditor *e = (projectionEditor*)data;
 
-  CTX.mesh.changed = ENT_SURFACE;
+  Msg(GERROR, "Compute!");
+}
 
+void mesh_parameterize_cb(Fl_Widget* w, void* data)
+{
+  // display geometry surfaces
   opt_geometry_surfaces(0, GMSH_SET | GMSH_GUI, 1);
- 
-  Draw();
-  
-  projection_editor("Projection editor", p);
 
-  Msg(INFO, "Model added: %d faces", GMODEL->numFace());
+  // create one instance of each available projection surface
+  std::vector<projectionFace*> faces;
+  if(faces.empty()){
+    faces.push_back(new parabolicCylinder(GMODEL, 10000));
+    faces.push_back(new parabolicCylinder(GMODEL, 10001));
+  }
+
+  // make each projection surface invisible and 
+  for(unsigned int i = 0; i < faces.size(); i++){
+    faces[i]->setVisibility(false);
+    GMODEL->add(faces[i]);
+  }
+ 
+  // launch editor
+  static projectionEditor *editor = 0;
+  if(!editor) editor = new projectionEditor(faces);
+  editor->show();
 }
diff --git a/Geo/GEntity.cpp b/Geo/GEntity.cpp
index fab13016032702b83f9b0d0d4613fcf1c31b0835..0f9e9c2b3fe286eb87508bdf33cff54788524694 100644
--- a/Geo/GEntity.cpp
+++ b/Geo/GEntity.cpp
@@ -1,4 +1,4 @@
-// $Id: GEntity.cpp,v 1.10 2006-11-30 13:55:20 geuzaine Exp $
+// $Id: GEntity.cpp,v 1.11 2007-07-26 13:10:48 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -36,6 +36,12 @@ GEntity::~GEntity()
   if(meshRep) delete meshRep; 
 }
 
+char GEntity::getVisibility()
+{ 
+  if(CTX.hide_unselected && !CTX.pick_elements && !getSelection()) return false;
+  return _visible; 
+}
+
 bool GEntity::useColor()
 { 
   int r = CTX.UNPACK_RED(_color);
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 07a07ac5abd36b50c1558873bc60611d32e86a28..50ad8d0e8999d961b8f27e4d8816d04e338ae6f6 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -184,7 +184,7 @@ class GEntity {
   virtual SBoundingBox3d bounds() const { throw; }
 
   // Get the visibility flag
-  virtual char getVisibility(){ return _visible; }
+  virtual char getVisibility();
 
   // Set the visibility flag
   virtual void setVisibility(char val, bool recursive=false){ _visible = val; }
diff --git a/Geo/GModelIO_Fourier.cpp b/Geo/GModelIO_Fourier.cpp
deleted file mode 100644
index 91b4729dee5ee2a20ce6159fbfd7ded13ef68726..0000000000000000000000000000000000000000
--- a/Geo/GModelIO_Fourier.cpp
+++ /dev/null
@@ -1,1038 +0,0 @@
-#include "GModel.h"
-#include "GVertex.h"
-#include "GEdge.h"
-#include "GFace.h"
-#include "gmshFace.h"
-#include "Message.h"
-#include "Context.h"
-#include "Views.h"
-#include "Range.h"
-#include "model.h"
-#include "meshGFace.h"
-
-extern Context_T CTX;
-
-#if defined(HAVE_FOURIER_MODEL)
-
-static model *FM = 0;
-
-class fourierVertex : public GVertex {
- private:
-  MVertex *_v;
- public:
-  fourierVertex(GModel *m, MVertex *v) : GVertex(m, v->getNum()), _v(v)
-  {
-    mesh_vertices.push_back(v);
-  }
-  virtual ~fourierVertex() {}
-  virtual GPoint point() const { return GPoint(_v->x(), _v->y(), _v->z(), this); }
-  virtual double x() const { return _v->x(); }
-  virtual double y() const { return _v->y(); }
-  virtual double z() const { return _v->z(); }
-  virtual double prescribedMeshSizeAtVertex() const { return 0.1; }
-  ModelType getNativeType() const { return FourierModel; }
-};
-
-class fourierEdge : public GEdge {
- public:
-  fourierEdge(GModel *m, int num, GVertex *v1, GVertex *v2);
-  virtual ~fourierEdge() {}
-  double period() const { throw ; }
-  virtual bool periodic(int dim=0) const { return false; }
-  virtual Range<double> parBounds(int i) const { throw; }
-  virtual GeomType geomType() const { throw; }
-  virtual bool degenerate(int) const { return false; }
-  virtual bool continuous(int dim) const { return true; }
-  virtual GPoint point(double p) const { throw; }
-  virtual GPoint closestPoint(const SPoint3 & queryPoint) { throw; }
-  virtual int containsPoint(const SPoint3 &pt) const { throw; }
-  virtual int containsParam(double pt) const { throw; }
-  virtual SVector3 firstDer(double par) const { throw; }
-  virtual SPoint2 reparamOnFace(GFace * face, double epar, int dir) const { throw; }
-  virtual double parFromPoint(const SPoint3 &pt) const { throw; }
-  virtual int minimumMeshSegments () const { throw; }
-  virtual int minimumDrawSegments () const { throw; }
-  ModelType getNativeType() const { return FourierModel; }
-};
-
-class fourierFace : public GFace {
- private:
-  // flag to know if is the face is already meshed
-  int _discrete;
-  // floag to know if the face is just a plane
-  int _plane;
-  // vertices and edges purely local to the face (not shared with the model)
-  GVertex *_v[4];
-  GEdge *_e[4];
- public:
-  fourierFace(GModel *m, int num);
-  fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole);
-  virtual ~fourierFace();
-
-  void meshBoundary();
-
-  Range<double> parBounds(int i) const; 
-  virtual int paramDegeneracies(int dir, double *par) { return 0; }
-  
-  virtual GPoint point(double par1, double par2) const; 
-  virtual GPoint closestPoint(const SPoint3 & queryPoint) const ; 
-  
-  virtual int containsPoint(const SPoint3 &pt) const;  
-  virtual int containsParam(const SPoint2 &pt) const; 
-  
-  virtual SVector3 normal(const SPoint2 &param) const; 
-  virtual Pair<SVector3,SVector3> firstDer(const SPoint2 &param) const {throw;} 
-  
-  virtual GEntity::GeomType geomType() const; 
-  virtual int geomDirection() const { return 1; }
-  
-  virtual bool continuous(int dim) const { return true; }
-  virtual bool periodic(int dim) const { return false; }
-  virtual bool degenerate(int dim) const { return false; }
-  virtual double period(int dir) const {throw;}
-  ModelType getNativeType() const { return FourierModel; }
-  void * getNativePtr() const {throw;} 
-  virtual bool surfPeriodic(int dim) const {throw;}
-  virtual SPoint2 parFromPoint(const SPoint3 &) const;
-};
-
-void debugVertices(std::vector<MVertex*> &vertices, std::string file, 
-		   bool parametric, int num=0)
-{
-  char name[256];
-  sprintf(name, "%s_%d.pos", file.c_str(), num);
-  FILE *fp = fopen(name, "w");
-  fprintf(fp, "View \"debug\"{\n");
-  for(unsigned int i = 0; i < vertices.size(); i++){
-    double x, y, z;
-    if(parametric){
-      vertices[i]->getParameter(0, x); 
-      vertices[i]->getParameter(1, y); 
-      z = 0;
-    }
-    else{
-      x = vertices[i]->x(); 
-      y = vertices[i]->y(); 
-      z = vertices[i]->z();
-    }
-    fprintf(fp, "SP(%g,%g,%g){%d};\n", x, y, z, i);
-  }
-  fprintf(fp, "};\n");    
-  fclose(fp);
-}
-
-template<class T>
-void debugElements(std::vector<T*> &elements, std::string file, 
-		   bool parametric, int num=0)
-{
-  char name[256];
-  sprintf(name, "%s_%d.pos", file.c_str(), num);
-  FILE *fp = fopen(name, "w");
-  fprintf(fp, "View \"debug\"{\n");
-  for(unsigned int i = 0; i < elements.size(); i++){
-    fprintf(fp, "%s(", elements[i]->getStringForPOS());
-    for(int j = 0; j < elements[i]->getNumVertices(); j++){
-      MVertex *v = elements[i]->getVertex(j);
-      if(j) fprintf(fp, ",");
-      double x, y, z;
-      if(parametric){
-	v->getParameter(0, x);
-	v->getParameter(1, y);
-	z = 0;
-      }
-      else{
-	x = v->x(); y = v->y(); z = v->z();
-      }
-      fprintf(fp, "%g,%g,%g", x, y, z);
-    }
-    fprintf(fp, "){");
-    for(int j = 0; j < elements[i]->getNumVertices(); j++){
-      MVertex *v = elements[i]->getVertex(j);
-      if(j) fprintf(fp, ",");
-      if(v->getData()){
-	double pou = *(double*)v->getData();
-	fprintf(fp, "%g", pou);
-      }
-      else{
-	fprintf(fp, "%d", v->onWhat()->tag());
-      }
-    }
-    fprintf(fp, "};\n");
-  }
-  fprintf(fp, "};\n");
-  fclose(fp);
-}
-
-class meshCartesian{
-public:
-  void operator() (GFace *gf)
-  {  
-    int M = (int)(30. / CTX.mesh.lc_factor), N = (int)(30. / CTX.mesh.lc_factor);
-
-#if 1
-    switch(gf->tag()){
-    case 0: M = 30; N = 22; break; // falcon front
-    case 1: M = 30; N = 70; break; // falcon back
-    case 2: break; // wing front
-    case 3: M = 25; N = 30; break; // wing top (N along the fuselage)
-    case 4: M = 25; N = 30; break; // wing bottom
-    }
-    //M /= 2 ;
-    //N /= 2 ;
-#endif
-
-    for(int i = 0; i < M; i++){
-      for(int j = 0; j < N; j++){
-	double u = i/(double)(M - 1);
-	double v = j/(double)(N - 1);
-	if(gf->tag()==0){
-	  v = sin(M_PI/2.*v - M_PI/2.) + 1.; 
-	}
-
-	//if(gf->tag() == 2){
-	  // hack for sttr report 2 (ship model, top right patch)
-	  // XYZ_values_ship_Top(u,v=0) = 
-	  // XYZ_values_ship_Med(u*matching_const1_+matching_const2_,v=1)
-	  // where 0=<u=<1, 
-	  //       matching_const1_=0.523540136032613,
-	  //       matching_const2_=0.475747890863610. 
-	//u = u * 0.5268 + 0.475;
-	//}
-
-	GPoint p = gf->point(u, v);
-	double pou=1;
-	//FM->GetPou(gf->tag(), u, v, pou);
-	gf->mesh_vertices.push_back
-	  (new MDataFaceVertex<double>(p.x(), p.y(), p.z(), gf, u, v, pou));
-      }
-    }
-    for(int i = 0; i < M - 1; i++){
-      for(int j = 0; j < N - 1; j++){
-#if 1
-	MQuadrangle *q = new MQuadrangle(gf->mesh_vertices[i * N + j],
-					 gf->mesh_vertices[(i + 1) * N + j],
-					 gf->mesh_vertices[(i + 1) * N + (j + 1)],
-					 gf->mesh_vertices[i * N + (j + 1)]);
-	//if(FM->GetOrientation(gf->tag()) < 0) 
-	q->revert();
-	gf->quadrangles.push_back(q);
-#else
-	MVertex *v1 = gf->mesh_vertices[i * N + j];
-	MVertex *v2 = gf->mesh_vertices[(i + 1) * N + j];
-	MVertex *v3 = gf->mesh_vertices[(i + 1) * N + (j + 1)];
-	MVertex *v4 = gf->mesh_vertices[i * N + (j + 1)];
-
-	const double tol = 1.e-6;
-	//if(v1->distance(v2) > tol && v1->distance(v3) > tol && v2->distance(v3) > tol)
-	{
-	  MTriangle *t = new MTriangle(v1, v2, v3);
-	  t->revert();
-	  gf->triangles.push_back(t);
-	}
-	//if(v1->distance(v3) > tol && v1->distance(v4) > tol && v3->distance(v4) > tol)
-	{
-	  MTriangle *t = new MTriangle(v1, v3, v4);
-	  t->revert();
-	  gf->triangles.push_back(t);
-	}
-#endif
-      }
-    }
-  }
-};
-
-class meshGmsh{
-public:
-  void operator() (GFace *gf)
-  {  
-    fourierFace *ff = dynamic_cast<fourierFace*>(gf);
-    if(!ff) {
-      Msg(GERROR, "face %d is not Fourier", gf->tag());
-      return;
-    }
-    ff->meshBoundary();
-    meshGFace mesh; 
-    mesh(ff);
-  }
-};
-
-class computePartitionOfUnity{
-public:
-  void operator() (GFace *gf)
-  {  
-    // we only normalize the partition of unity amongst patches that
-    // overlap; we don't normalize across hard edges
-    std::vector<int> overlaps;
-    FM->GetNeighbor(gf->tag(), overlaps);
-    for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
-      MVertex *v = gf->mesh_vertices[i];
-      std::vector<std::pair<GFace*, SPoint2> > param;
-      for(unsigned int j = 0; j < overlaps.size(); j++){
-	GFace *gf2 = gf->model()->faceByTag(overlaps[j]);
-	SPoint2 uv2 = gf2->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
-	if(gf2->containsParam(uv2)){
-	  GPoint xyz2 = gf2->point(uv2);
-	  const double tol = 1.e-2; // Warning: tolerance
-	  if(fabs(xyz2.x() - v->x()) < tol && 
-	     fabs(xyz2.y() - v->y()) < tol && 
-	     fabs(xyz2.z() - v->z()) < tol)
-	    param.push_back(std::pair<GFace*, SPoint2>(gf2, uv2));
-	}
-      }
-      if(param.empty()){
-      	Msg(GERROR, "Vertex %d does not belong to any patch", v->getNum());
-      }
-      else{
-	double allPou = 0.;
-	for(unsigned int i = 0; i < param.size(); i++){
-	  //double u2 = param[i].second[0], v2 = param[i].second[1];
-	  double pou2=1;
-	  //FM->GetPou(param[i].first->tag(), u2, v2, pou2);
-	  allPou += pou2;
-	}
-	if(v->getData()){
-	  double *pou = (double*)v->getData();
-	  *pou = *pou / allPou;
-	}
-	else{
-	  Msg(GERROR, "Vertex %d has no POU data", v->getNum());
-	}
-      }
-    }
-
-    //debugElements(gf->quadrangles, "pou", false, gf->tag());
-  }
-};
-
-class createGroove{
-public:
-  void operator() (GFace *gf)
-  {  
-    // mark elements in the groove with '-1' visibility tag
-    for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
-      MElement *e = gf->quadrangles[i];
-      for(int j = 0; j < e->getNumVertices(); j++){
-	void *data = e->getVertex(j)->getData();
-	if(data){
-	  double pou = *(double*)data;
-	  if(pou < 0.51)
-	    e->setVisibility(-1);
-	}
-      }
-    }
-
-    // remove elements in the groove
-    std::vector<MQuadrangle*> newq;
-    for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
-      if(gf->quadrangles[i]->getVisibility() < 0) 
-	delete gf->quadrangles[i];
-      else
-	newq.push_back(gf->quadrangles[i]);
-    gf->quadrangles = newq;
-
-    // remove vertices in the groove
-    std::vector<MVertex*> newv;
-    for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
-      gf->mesh_vertices[i]->setVisibility(-1);
-    for(unsigned int i = 0; i < gf->quadrangles.size(); i++)
-      for(int j = 0; j < gf->quadrangles[i]->getNumVertices(); j++)
-	gf->quadrangles[i]->getVertex(j)->setVisibility(1);
-    for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
-      if(gf->mesh_vertices[i]->getVisibility() < 0) 
-	delete gf->mesh_vertices[i];
-      else
-	newv.push_back(gf->mesh_vertices[i]);
-    gf->mesh_vertices = newv;
-  }
-};
-
-void getOrderedBoundaryLoops(std::vector<MElement*> &elements, 
-			     std::vector<std::vector<MVertex*> > &loops)
-{
-  typedef std::pair<MVertex*, MVertex*> vpair;
-  std::map<vpair, vpair> edges;
-  for(unsigned int i = 0; i < elements.size(); i++){
-    for(int j = 0; j < elements[i]->getNumEdges(); j++){
-      MEdge e = elements[i]->getEdge(j);
-      vpair p(e.getMinVertex(), e.getMaxVertex());
-      if(edges.count(p)) edges.erase(p);
-      else edges[p] = vpair(e.getVertex(0), e.getVertex(1));
-    }
-  }
-
-  std::map<MVertex*, MVertex*> connect;
-  for(std::map<vpair, vpair>::iterator it = edges.begin(); it != edges.end(); it++)
-    connect[it->second.first] = it->second.second;
-
-  loops.resize(1);
-  while(connect.size()){
-    if(loops[loops.size() - 1].empty()) 
-      loops[loops.size() - 1].push_back(connect.begin()->first);
-    MVertex *prev = loops[loops.size() - 1][loops[loops.size() - 1].size() - 1];
-    MVertex *next = connect[prev];
-    connect.erase(prev);
-    loops[loops.size() - 1].push_back(next);
-    if(next == loops[loops.size() - 1][0])
-      loops.resize(loops.size() + 1);
-  }
-
-  if(loops[loops.size() - 1].empty())
-    loops.pop_back();
-  
-  Msg(INFO, "Found %d loop(s) in boundary", loops.size());
-}
-
-void classifyFaces(GFace *gf, std::set<GFace*> &connected, std::set<GFace*> &other)
-{
-  connected.insert(gf); 
-  std::set<MVertex*> verts;
-  for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++)
-    verts.insert(gf->mesh_vertices[i]);
-
-  for(int i = 0; i < gf->model()->numFace() - 1; i++){ // worst case
-    for(GModel::fiter it = gf->model()->firstFace(); it != gf->model()->lastFace(); it++){
-      if(connected.find(*it) == connected.end()){
-	for(unsigned int j = 0; j < (*it)->mesh_vertices.size(); j++){
-	  if(std::find(verts.begin(), verts.end(), (*it)->mesh_vertices[j]) != verts.end()){
-	    connected.insert(*it);
-	    for(unsigned int k = 0; k < (*it)->mesh_vertices.size(); k++)
-	      verts.insert((*it)->mesh_vertices[k]);
-	    break;
-	  }
-	}
-      }
-    }
-  }
-
-  for(GModel::fiter it = gf->model()->firstFace(); it != gf->model()->lastFace(); it++)
-    if(connected.find(*it) == connected.end()) 
-      other.insert(*it);
-
-  Msg(INFO, "Found %d face(s) connected to face %d", (int)connected.size(), gf->tag());
-  for(std::set<GFace*>::iterator it = connected.begin(); it != connected.end(); it++)
-    Msg(INFO, "   %d", (*it)->tag());
-}
-
-void getIntersectingBoundaryParts(GFace *gf, std::vector<MElement*> &elements,
-				  std::vector<std::vector<std::vector<MVertex*> > > &parts)
-{
-  std::vector<std::vector<MVertex*> > loops;
-  getOrderedBoundaryLoops(elements, loops);
-  parts.resize(loops.size());
-
-  for(unsigned int i = 0; i < loops.size(); i++){
-    // last vertex in loop is equal to the first vertex
-    bool newpart = true;
-    for(unsigned int j = 0; j < loops[i].size() - 1; j++){
-      MVertex *v = loops[i][j];
-      SPoint2 uv = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
-      if(gf->containsParam(uv)){
-	GPoint xyz = gf->point(uv);
-	const double tol = 1.e-2; // Warning: tolerance
-	if(fabs(xyz.x() - v->x()) < tol && 
-	   fabs(xyz.y() - v->y()) < tol && 
-	   fabs(xyz.z() - v->z()) < tol){
-	  if(newpart){
-	    parts[i].resize(parts[i].size() + 1);
-	    newpart = false;
-	  }
-	  v->setParameter(0, uv[0]);
-	  v->setParameter(1, uv[1]);
-	  parts[i][parts[i].size() - 1].push_back(v);
-	}
-	else{
-	  newpart = true;
-	}
-      }
-      else{
-	newpart = true;
-      }
-    }
-  }
-
-#if 0
-  static int nn = 0;
-  debugElements(elements, "elements", false, nn++);
-  for(unsigned int i = 0; i < loops.size(); i++)
-    debugVertices(loops[i], "loops", false, nn++);
-  for(unsigned int i = 0; i < parts.size(); i++)
-    for(unsigned int j = 0; j < parts[i].size(); j++)
-      debugVertices(parts[i][j], "parts", false, nn++);
-#endif
-}
-
-void meshGrout(GFace *gf, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole)
-{ 
-  debugVertices(hole, "hole", true, gf->tag());
-  debugVertices(loop, "loop", true, gf->tag());
-
-  fourierFace *grout = new fourierFace(gf, loop, hole);
-  meshGFace mesh; 
-  mesh(grout);
-  orientMeshGFace orient; 
-  orient(grout);
-  //debugElements(grout->triangles, "grout", true);
-  for(unsigned int i = 0; i < grout->triangles.size(); i++)
-    gf->triangles.push_back(grout->triangles[i]);
-  for(unsigned int i = 0; i < loop.size(); i++)
-    gf->mesh_vertices.push_back(loop[i]);
-  for(unsigned int i = 0; i < hole.size(); i++)
-    gf->mesh_vertices.push_back(hole[i]);
-  for(unsigned int i = 0; i < grout->mesh_vertices.size(); i++)
-    gf->mesh_vertices.push_back(grout->mesh_vertices[i]);
-  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)
-{
-  subloops.resize(1);
-  subloops[0].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
-      subloops.resize(subloops.size() + 1);
-    }
-    else{
-      subloops[subloops.size() - 1].push_back(loop[i]);
-    }
-  }
-  subloops[subloops.size() - 1].push_back(loop[loop.size() - 1]);
-
-  if(subloops.size() > 1){
-    for(unsigned int i = 0; i < subloops.size(); i++){
-      //if(subloops[i].size()) debugVertices(subloops[i], "x", false, i);
-    }
-    Msg(INFO, "HAVE SUBLOOPS!");
-    return true;
-  }
-  return false;
-}
-
-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)
-  {  
-    if(gf->tag() > 2) return;
-
-    Msg(INFO, "Processing grout for face %d", gf->tag());
-
-    std::set<GFace*> connected, other;
-    classifyFaces(gf, connected, other);
-
-    std::vector<MElement*> connectedElements;
-    for(std::set<GFace*>::iterator it = connected.begin(); it != connected.end(); it++){
-      for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
-	connectedElements.push_back((*it)->triangles[i]);
-      for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
-	connectedElements.push_back((*it)->quadrangles[i]);
-    }
-
-    std::vector<MElement*> otherElements;
-    for(std::set<GFace*>::iterator it = other.begin(); it != other.end(); it++){
-      for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
-	otherElements.push_back((*it)->triangles[i]);
-      for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++)
-	otherElements.push_back((*it)->quadrangles[i]);
-    }
-
-    std::vector<std::vector<std::vector<MVertex*> > > inside;
-    getIntersectingBoundaryParts(gf, connectedElements, inside);
-    Msg(INFO, "%d inside loop(s)", (int)inside.size());
-    for(unsigned int i = 0; i < inside.size(); i++)
-      Msg(INFO, "   inside loop %d intersect has %d part(s)", i, (int)inside[i].size());
-
-    std::vector<std::vector<std::vector<MVertex*> > > outside;
-    getIntersectingBoundaryParts(gf, otherElements, outside);
-    Msg(INFO, "%d outside loop(s)", (int)outside.size());
-    for(unsigned int i = 0; i < outside.size(); i++)
-      Msg(INFO, "   outside loop %d intersect has %d part(s)", i, (int)outside[i].size());
-
-    if(inside.size() == 1 && inside[0].size() == 1){
-      Msg(INFO, "*** CASE 1: grout has a single hole in one part ***");
-      meshGroutWithHole(gf, inside, outside);
-    }
-    else if(!outside.size()){
-      Msg(INFO, "*** CASE 2: grout has no outside contributions ***");
-      for(unsigned int i = 0; i < inside.size(); i++)
-	meshGroutWithoutHole(gf, inside[i]);
-    }
-    else{
-      Msg(WARNING, "*** TODO: UNKNOWN CASE ***");
-    }
-  }
-};
-
-fourierEdge::fourierEdge(GModel *model, int num, GVertex *v1, GVertex *v2)
-  : GEdge(model, num, v1, v2)
-{
-}
-
-fourierFace::fourierFace(GModel *m, int num)
-  : GFace(m, num)
-{
-  for(int i = 0; i < 4; i++){ _v[i] = 0; _e[i] = 0; }
-  _discrete = 1;
-  _plane = 0;
-}
-
-fourierFace::fourierFace(GFace *f, std::vector<MVertex*> &loop, std::vector<MVertex*> &hole)
-  : GFace(f->model(), f->tag())
-{
-  for(int i = 0; i < 4; i++){ _v[i] = 0; _e[i] = 0; }
-  _discrete = 0;
-  _plane = 0;
-
-  if(!loop.size()){
-    Msg(GERROR, "No vertices in exterior loop");
-    return; 
-  }
-  
-  _v[0] = new fourierVertex(f->model(), loop[0]);
-  _v[1] = new fourierVertex(f->model(), loop[loop.size() - 2]);
-  _e[0] = new fourierEdge(f->model(), 1, _v[0], _v[1]);
-  _e[0]->addFace(this);
-  _e[1] = new fourierEdge(f->model(), 2, _v[1], _v[0]);
-  _e[1]->addFace(this);
-  for(unsigned int i = 1; i < loop.size() - 2; i++)
-    _e[0]->mesh_vertices.push_back(loop[i]);
-
-  l_edges.push_back(_e[0]); l_dirs.push_back(-1);
-  l_edges.push_back(_e[1]); l_dirs.push_back(-1);
-
-  if(hole.size()){
-    _v[2] = new fourierVertex(f->model(), hole[0]);
-    _v[3] = new fourierVertex(f->model(), hole[hole.size() - 2]);
-    _e[2] = new fourierEdge(f->model(), 3, _v[2], _v[3]);
-    _e[2]->addFace(this);
-    _e[3] = new fourierEdge(f->model(), 4, _v[3], _v[2]);
-    _e[3]->addFace(this);
-    for(unsigned int i = 1; i < hole.size() - 2; i++)
-      _e[2]->mesh_vertices.push_back(hole[i]);
-
-    l_edges.push_back(_e[2]); l_dirs.push_back(1);
-    l_edges.push_back(_e[3]); l_dirs.push_back(1);
-  }
-}
-
-fourierFace::~fourierFace()
-{
-  if(!_discrete){
-    // this face was just used temporarily for meshing, so don't
-    // delete the mesh vertices or the mesh elements: they have been
-    // transferred elsewhere
-    for(int i = 0; i < 4; i++){ 
-      if(_e[i]){
-	_e[i]->mesh_vertices.clear();
-	delete _e[i];
-      }
-    }
-    for(int i = 0; i < 4; i++){ 
-      if(_v[i]){
-	_v[i]->mesh_vertices.clear();
-	delete _v[i];
-      }
-    }
-    triangles.clear();
-    quadrangles.clear();
-    mesh_vertices.clear();
-    l_edges.clear();
-  }
-}
-
-void fourierFace::meshBoundary()
-{
-  /*
-  const double tol = 1.e-6;
-
-  _discrete = 0;
-
-  int nu=0, nv=0;
-  FM->GetNum(tag(), nu, nv);
-
-  printf("patch=%d  nu=%d  nv=%d\n", tag(), nu, nv);
-  // youngae bug when not using chebyshev refinement:
-  //if(tag() == 6) nu = 36;
-  
-
-  std::vector<double> u, v;
-  FM->GetBoundary_Points(tag(), u, v);
-
-  if(2*nu+2*nv != (int)u.size()){
-    Msg(INFO, "Special planar patch from YoungAe: %d", tag());
-#if 1 // transfinite, by hand -- WARNING
-    _plane = 1; // to enable smoothing of transfinite meshes
-    // coarse:
-    //nu = 14; nv = 9;
-    // finer:
-    //nu = 24; nv = 14;
-    // fine:
-    nu = 52; nv = 28;
-
-    // remove duplicates
-    std::vector<MVertex*> verts;
-    for(unsigned int i = 0; i < u.size() - 1; i++){
-      if(!i || fabs(u[i]-u[i-1]) > tol || fabs(v[i]-v[i-1]) > tol){
-	GPoint p = point(u[i], v[i]);
-	verts.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
-      }
-    }
-    int corners[4] = {0, nu - 1, nu + nv - 2,  2 * nu + nv - 3};
-    for(int i = 0; i < 4; i++)
-      _v[i] = new fourierVertex(model(), verts[corners[i]]);
-    meshAttributes.Method = TRANSFINI;
-    meshAttributes.recombine = 1;
-    meshAttributes.transfiniteArrangement = 1;
-    meshAttributes.corners.clear();
-    for(int i = 0; i < 4; i++)
-      meshAttributes.corners.push_back(_v[i]);
-    _e[0] = new fourierEdge(model(), 1, _v[0], _v[1]);
-    _e[0]->addFace(this);
-    _e[1] = new fourierEdge(model(), 2, _v[1], _v[2]);
-    _e[1]->addFace(this);
-    _e[2] = new fourierEdge(model(), 3, _v[2], _v[3]);
-    _e[2]->addFace(this);
-    _e[3] = new fourierEdge(model(), 4, _v[3], _v[0]);
-    _e[3]->addFace(this);
-    for(int i = corners[0] + 1; i < corners[1]; i++)
-      _e[0]->mesh_vertices.push_back(verts[i]);
-    for(int i = corners[1] + 1; i < corners[2]; i++)
-      _e[1]->mesh_vertices.push_back(verts[i]);
-    for(int i = corners[2] + 1; i < corners[3]; i++)
-      _e[2]->mesh_vertices.push_back(verts[i]);
-    for(int i = corners[3] + 1; i < (int)verts.size(); i++)
-      _e[3]->mesh_vertices.push_back(verts[i]);
-    l_edges.push_back(_e[0]); l_dirs.push_back(1);
-    l_edges.push_back(_e[1]); l_dirs.push_back(1);
-    l_edges.push_back(_e[2]); l_dirs.push_back(1);
-    l_edges.push_back(_e[3]); l_dirs.push_back(1);
-#else // unstructured
-    GPoint p0 = point(u[0], v[0]);
-    GPoint p1 = point(u[1], v[1]);
-    MVertex *v0 = new MFaceVertex(p0.x(), p0.y(), p0.z(), this, u[0], v[0]);
-    MVertex *v1 = new MFaceVertex(p1.x(), p1.y(), p1.z(), this, u[1], v[1]);
-    _v[0] = new fourierVertex(model(), v0);
-    _v[1] = new fourierVertex(model(), v1);
-    _e[0] = new fourierEdge(model(), 1, _v[0], _v[1]);
-    _e[0]->addFace(this);
-    _e[1] = new fourierEdge(model(), 2, _v[1], _v[0]);
-    _e[1]->addFace(this);
-    for(unsigned int i = 2; i < u.size() - 1; i++){
-      if(fabs(u[i]-u[i-1]) > tol || fabs(v[i]-v[i-1]) > tol){
-	GPoint p = point(u[i], v[i]);
-	MVertex *vv = new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]);
-	_e[1]->mesh_vertices.push_back(vv);
-      }
-    }
-    l_edges.push_back(_e[0]); l_dirs.push_back(1);
-    l_edges.push_back(_e[1]); l_dirs.push_back(1);
-#endif
-    return;
-  }
-  
-  int corners[4] = {0, nu, nu + nv, 2 * nu + nv};
-  for(int i = 0; i < 4; i++){
-    double uu = u[corners[i]], vv = v[corners[i]];
-    GPoint p = point(uu, vv);
-    MVertex *newv = new MFaceVertex(p.x(), p.y(), p.z(), this, uu, vv);
-    _v[i] = new fourierVertex(model(), newv);
-  }
-  
-  meshAttributes.Method = TRANSFINI;
-  meshAttributes.recombine = 1;
-  meshAttributes.transfiniteArrangement = 1;
-  meshAttributes.corners.clear();
-  for(int i = 0; i < 4; i++)
-    meshAttributes.corners.push_back(_v[i]);
-      
-  _e[0] = new fourierEdge(model(), 1, _v[0], _v[1]);
-  _e[0]->addFace(this);
-  _e[1] = new fourierEdge(model(), 2, _v[1], _v[2]);
-  _e[1]->addFace(this);
-  _e[2] = new fourierEdge(model(), 3, _v[2], _v[3]);
-  _e[2]->addFace(this);
-  _e[3] = new fourierEdge(model(), 4, _v[3], _v[0]);
-  _e[3]->addFace(this);
-
-  for(int i = corners[0] + 1; i < corners[1] - 1; i++){
-    GPoint p = point(u[i], v[i]);
-    _e[0]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
-  }
-  for(int i = corners[1] + 1; i < corners[2] - 1; i++){
-    GPoint p = point(u[i], v[i]);
-    _e[1]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
-  }
-  for(int i = corners[2] + 1; i < corners[3] - 1; i++){
-    GPoint p = point(u[i], v[i]);
-    _e[2]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
-  }
-  for(int i = corners[3] + 1; i < (int)u.size() - 1; i++){
-    GPoint p = point(u[i], v[i]);
-    _e[3]->mesh_vertices.push_back(new MFaceVertex(p.x(), p.y(), p.z(), this, u[i], v[i]));
-  }
-
-  l_edges.push_back(_e[0]); l_dirs.push_back(1);
-  l_edges.push_back(_e[1]); l_dirs.push_back(1);
-  l_edges.push_back(_e[2]); l_dirs.push_back(1);
-  l_edges.push_back(_e[3]); l_dirs.push_back(1);
-  */
-}
-
-Range<double> fourierFace::parBounds(int i) const
-{
-  double min, max;
-  FM->GetParamBounds(tag(), i, min, max);
-  return Range<double>(min, max);
-}
-  
-GPoint fourierFace::point(double par1, double par2) const
-{
-  double x, y, z;
-  FM->GetPoint(tag(), par1, par2, x, y, z);
-  return GPoint(x, y, z);
-}
-
-GPoint fourierFace::closestPoint(const SPoint3 & queryPoint) const
-{
-  throw;
-}
-  
-int fourierFace::containsPoint(const SPoint3 &pt) const
-{
-  throw;
-}
-
-int fourierFace::containsParam(const SPoint2 &pt) const
-{
-  double umin, umax, vmin, vmax;
-  FM->GetParamBounds(tag(), 0, umin, umax);
-  FM->GetParamBounds(tag(), 1, vmin, vmax);
-  const double tol = 1.e-6;
-  if(pt[0] < umin - tol || pt[0] > umax + tol) return 0;
-  if(pt[1] < vmin - tol || pt[1] > vmax + tol) return 0;
-  return 1;
-}
-  
-SVector3 fourierFace::normal(const SPoint2 &param) const
-{
-  throw;
-}
-
-GEntity::GeomType fourierFace::geomType() const
-{
-  if(_discrete)
-    return  GEntity::DiscreteSurface;
-  else if(_plane)
-    return  GEntity::Plane;
-  else
-    return GEntity::Unknown;
-}
-
-SPoint2 fourierFace::parFromPoint(const SPoint3 &p) const
-{
-  double u, v;
-  FM->GetParamFromPoint(tag(), p.x(), p.y(), p.z(), u, v);
-  return SPoint2(u, v);
-}
-
-void cleanUpAndMergeAllFaces(GModel *m)
-{
-  GFace *newgf = new gmshFace(m, 1);
-  std::set<MVertex*, MVertexLessThanLexicographic> pos;
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
-    GFace *gf = *it;
-    pos.insert(gf->mesh_vertices.begin(), gf->mesh_vertices.end());
-    gf->mesh_vertices.clear();
-  }
-  std::set<MVertex*, MVertexLessThanLexicographic>::iterator itp;
-  for(itp = pos.begin(); itp != pos.end(); itp++)
-    newgf->mesh_vertices.push_back(*itp);
-
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
-    GFace *gf = *it;
-    for(unsigned int i = 0; i < gf->triangles.size(); i++){
-      std::vector<MVertex*> v(3);
-      for(int j = 0; j < 3; j++){
-	itp = pos.find(gf->triangles[i]->getVertex(j));
-	if(itp == pos.end())
-	  Msg(GERROR, "Could not find vertex");
-	else
-	  v[j] = *itp;
-      }
-      delete gf->triangles[i];
-      if(v[0] != v[1] && v[0] != v[2] && v[1] != v[2])
-	newgf->triangles.push_back(new MTriangle(v));
-    }
-    for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
-      std::vector<MVertex*> v(4);
-      for(int j = 0; j < 4; j++){
-	itp = pos.find(gf->quadrangles[i]->getVertex(j));
-	if(itp == pos.end())
-	  Msg(GERROR, "Could not find vertex");
-	else
-	  v[j] = *itp;
-      }
-      delete gf->quadrangles[i];
-      //if(v[0] != v[1] && v[0] != v[2] && v[0] != v[3] && 
-      // v[1] != v[2] && v[1] != v[3] && v[2] != v[3])
-	newgf->quadrangles.push_back(new MQuadrangle(v));
-    }
-    m->remove(gf);
-  }
-  m->add(newgf);
-}
-
-int GModel::readFourier(const std::string &name)
-{
-  FM = new model(name);
-
-  CTX.terminal = 1;
-  
-  Msg(INFO, "Fourier model created: %d patches", FM->GetNumPatches());
-
-  // create one face per patch
-  for(int i = 0; i < FM->GetNumPatches(); i++)
-    add(new fourierFace(this, i));
-
-  // mesh each face with quads
-  std::for_each(firstFace(), lastFace(), meshCartesian());
-
-  cleanUpAndMergeAllFaces(this);
-  return 1;
-
-  // mesh each face using the standard gmsh algorithms
-  //std::for_each(firstFace(), lastFace(), meshGmsh());
-  //return 1;
-
-  // compute partition of unity
-  //std::for_each(firstFace(), lastFace(), computePartitionOfUnity());
-  //return 1;
-
-  // create grooves
-  std::for_each(firstFace(), lastFace(), createGroove());
-
-  // create grout
-  std::for_each(firstFace(), lastFace(), createGrout());
-
-  // 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.terminal = 0;
-
-  CTX.mesh.changed = ENT_ALL;
-
-  return 1;
-}
-
-#else
-
-int GModel::readFourier(const std::string &name)
-{
-  Msg(GERROR, "Fourier model is not compiled in this version of Gmsh");
-  return 1;
-}
-
-#endif
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index b2b1a866f0ebacce76c1a826f29ec355bc89f44e..cc5245408f247ff455d1126b893393b846cdde8d 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.35 2007-04-12 08:47:24 remacle Exp $
+// $Id: MElement.cpp,v 1.36 2007-07-26 13:10:48 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -24,10 +24,19 @@
 #include "GEntity.h"
 #include "Numeric.h"
 #include "Message.h"
+#include "Context.h"
+
+extern Context_T CTX;
 
 int MElement::_globalNum = 0;
 double MElementLessThanLexicographic::tolerance = 1.e-6;
 
+char MElement::getVisibility()
+{
+  if(CTX.hide_unselected && _visible < 2) return false;
+  return _visible; 
+}
+
 double MElement::minEdge()
 {
   double m = 1.e25;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7cf90bcdf42e1c7f3e3a7d4cd1f05739352906f1..9e9a23bf9d096295a35fff841037c76ebd4988a9 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -70,7 +70,7 @@ class MElement
   virtual void setPartition(int num){ _partition = (short)num; }
 
   // get/set the visibility flag
-  virtual char getVisibility(){ return _visible; }
+  virtual char getVisibility();
   virtual void setVisibility(char val){ _visible = val; }
 
   // get the vertices
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 6eb2e0d8d422ed023e0434e7c3c683daee997a64..513369450f99af6959979b2031475c7c2ffd177b 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.198 2007-03-11 20:18:58 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.199 2007-07-26 13:10:48 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //