diff --git a/Common/Context.h b/Common/Context.h
index c385df687097aa139fa57a3309cb482e5504a073..531cb0fcf13d45fb82b89c5f67fd75cb1548602b 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -215,6 +215,7 @@ public :
     double eps_line_width_factor, eps_point_size_factor;
     int jpeg_quality, jpeg_smoothing;
     int gif_dither, gif_sort, gif_interlace, gif_transparent;
+    int geo_labels;
     int text;
   } print;
 
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 84a5d86d0e712bf852a0b3886a4b65323507aa50..69aa38998a9e6b71684ba7229eb12a4137915bf3 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1365,6 +1365,8 @@ StringXNumber PrintOptions_Number[] = {
   { F|O, "Format" , opt_print_format , FORMAT_AUTO , 
     "File format (10=automatic)" }, 
 
+  { F|O, "GeoLabels" , opt_print_geo_labels , 1. ,
+    "Save labels in unrolled Gmsh geometries" },
   { F|O, "GifDither" , opt_print_gif_dither , 0. ,
     "Apply dithering to GIF output" },
   { F|O, "GifInterlace" , opt_print_gif_interlace , 0. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 0c6aad7bd07325bd198d9374017ab7beaadbb813..edc644fdeecc109d4005b8e37301abca8caae955 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.337 2007-03-13 16:11:42 remacle Exp $
+// $Id: Options.cpp,v 1.338 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -6972,6 +6972,13 @@ double opt_print_jpeg_smoothing(OPT_ARGS_NUM)
   return CTX.print.jpeg_smoothing;
 }
 
+double opt_print_geo_labels(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.print.geo_labels = (int)val;
+  return CTX.print.geo_labels;
+}
+
 double opt_print_gif_dither(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 6d083a38a92f7d255919d0fdfd2de534c50837ae..2ffe7e58119a77652a3e478aae4e73759539c888 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -635,6 +635,7 @@ double opt_print_eps_line_width_factor(OPT_ARGS_NUM);
 double opt_print_eps_point_size_factor(OPT_ARGS_NUM);
 double opt_print_jpeg_quality(OPT_ARGS_NUM);
 double opt_print_jpeg_smoothing(OPT_ARGS_NUM);
+double opt_print_geo_labels(OPT_ARGS_NUM);
 double opt_print_gif_dither(OPT_ARGS_NUM);
 double opt_print_gif_sort(OPT_ARGS_NUM);
 double opt_print_gif_interlace(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index b02f8eff678f884760f41ecfc7aa636d021d521a..9801ca047f63df58dc87a9f2236d7a03349ce92f 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.518 2007-03-13 16:11:42 remacle Exp $
+// $Id: Callbacks.cpp,v 1.519 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -649,7 +649,7 @@ void file_merge_cb(CALLBACK_ARGS)
 int _save_msh(char *name){ return msh_dialog(name); }
 int _save_pos(char *name){ return generic_mesh_dialog(name, "POS Options", FORMAT_POS); }
 int _save_options(char *name){ return options_dialog(name); }
-int _save_geo(char *name){ CreateOutputFile(name, FORMAT_GEO); return 1; }
+int _save_geo(char *name){ return geo_dialog(name); }
 int _save_cgns(char *name){ CreateOutputFile(name, FORMAT_CGNS); return 1; }
 int _save_unv(char *name){ return generic_mesh_dialog(name, "UNV Options", FORMAT_UNV); }
 int _save_med(char *name){ return generic_mesh_dialog(name, "MED Options", FORMAT_MED); }
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 984e715e08ea73425505a03388fda6466231eefd..9cd44e3dc712ab7686ec78bebd0106fc88d70833 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.601 2007-03-13 16:11:42 remacle Exp $
+// $Id: GUI.cpp,v 1.602 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -2253,7 +2253,7 @@ void GUI::create_option_window()
       mesh_value[0]->align(FL_ALIGN_RIGHT);
       mesh_value[0]->callback(mesh_options_ok_cb);
 
-      mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 5 * BH, IW, BH, "Element orders (1-5)");
+      mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 5 * BH, IW, BH, "Element order");
       mesh_value[3]->minimum(1);
       mesh_value[3]->maximum(5);
       mesh_value[3]->step(1);
@@ -2282,7 +2282,7 @@ void GUI::create_option_window()
 #endif
       mesh_butt[2]->callback(mesh_options_ok_cb);
 
-      mesh_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 10 * BH, BW, BH, "Optimize high order mesh (currently for 2D - plane)");
+      mesh_butt[3] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 10 * BH, BW, BH, "Optimize high order mesh (currently only for 2D/plane)");
       mesh_butt[3]->type(FL_TOGGLE_BUTTON);
       mesh_butt[3]->callback(mesh_options_ok_cb);
 
diff --git a/Fltk/GUI_Extras.cpp b/Fltk/GUI_Extras.cpp
index 1a039607541e68deaa63858c315948b83b4ac69f..da8a881a73a8f798187604e85054b9c89c9f4fb7 100644
--- a/Fltk/GUI_Extras.cpp
+++ b/Fltk/GUI_Extras.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI_Extras.cpp,v 1.32 2006-12-18 21:19:55 geuzaine Exp $
+// $Id: GUI_Extras.cpp,v 1.33 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -592,6 +592,59 @@ int options_dialog(char *name)
   return 0;
 }
 
+// geo dialog
+
+int geo_dialog(char *name)
+{
+  struct _geo_dialog{
+    Fl_Window *window;
+    Fl_Check_Button *b;
+    Fl_Button *ok, *cancel;
+  };
+  static _geo_dialog *dialog = NULL;
+
+  const int BH = 2 * CTX.fontsize + 1;
+  const int BB = 7 * CTX.fontsize + 9;
+  const int WB = 7;
+
+  if(!dialog){
+    dialog = new _geo_dialog;
+    int h = 3 * WB + 2 * BH, w = 2 * BB + 3 * WB, y = WB;
+    // not a "Dialog_Window" since it is modal 
+    dialog->window = new Fl_Double_Window(w, h, "GEO Options");
+    dialog->window->box(GMSH_WINDOW_BOX);
+    dialog->b = new Fl_Check_Button(WB, y, 2 * BB + WB, BH, "Save physical group labels"); y += BH;
+    dialog->b->type(FL_TOGGLE_BUTTON);
+    dialog->ok = new Fl_Return_Button(WB, y + WB, BB, BH, "OK");
+    dialog->cancel = new Fl_Button(2 * WB + BB, y + WB, BB, BH, "Cancel");
+    dialog->window->set_modal();
+    dialog->window->end();
+    dialog->window->hotspot(dialog->window);
+  }
+  
+  dialog->b->value(CTX.print.geo_labels ? 1 : 0);
+  dialog->window->show();
+
+  while(dialog->window->shown()){
+    Fl::wait();
+    for (;;) {
+      Fl_Widget* o = Fl::readqueue();
+      if (!o) break;
+      if (o == dialog->ok) {
+	opt_print_geo_labels(0, GMSH_SET | GMSH_GUI, dialog->b->value() ? 1 : 0);
+	CreateOutputFile(name, FORMAT_GEO);
+	dialog->window->hide();
+	return 1;
+      }
+      if (o == dialog->window || o == dialog->cancel){
+	dialog->window->hide();
+	return 0;
+      }
+    }
+  }
+  return 0;
+}
+
 // Generic save mesh dialog
 
 int generic_mesh_dialog(char *name, char *title, int format)
diff --git a/Fltk/GUI_Extras.h b/Fltk/GUI_Extras.h
index c5b70db991ebe4e087221c4df8e03d415583a958..584a9af38a57bffb3b0eda183bf78c95834cbdea 100644
--- a/Fltk/GUI_Extras.h
+++ b/Fltk/GUI_Extras.h
@@ -31,6 +31,7 @@ int perspective_editor();
 
 int jpeg_dialog(char *filename);
 int gif_dialog(char *filename);
+int geo_dialog(char *filename);
 int generic_bitmap_dialog(char *filename, char *title, int format);
 int generic_mesh_dialog(char *filename, char *title, int format);
 int gl2ps_dialog(char *filename, char *title, int format);
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 8a5410c50214362810f4fa046764449e35c52146..4140af7a5a62764164ddaef40016e7e82911f7a6 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -161,7 +161,7 @@ class GModel
   // Gmsh native CAD format
   int importTHEM();
   int readGEO(const std::string &name);
-  int writeGEO(const std::string &name);
+  int writeGEO(const std::string &name, bool printLabels=true);
 
   // Fourier model
   int readFourier(const std::string &name);
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index f2a6fab26d6d7be5f6faef2bef556513ced35026..fcf48f119bbed7586894a7fe3134b7b1bdcb769d 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Geo.cpp,v 1.10 2007-03-13 09:25:50 geuzaine Exp $
+// $Id: GModelIO_Geo.cpp,v 1.11 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -149,11 +149,11 @@ class writeGVertexGEO {
     if(gv->getNativeType() == GEntity::GmshModel){
       Vertex *v = (Vertex*)gv->getNativePtr();
       if(!v) return;
-      fprintf(geo, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+      fprintf(geo, "Point (%d) = {%.16g, %.16g, %.16g, %.16g};\n",
 	      v->Num, v->Pos.X, v->Pos.Y, v->Pos.Z, v->lc);
     }
     else{
-      fprintf(geo, "Point(%d) = {%.16g, %.16g, %.16g, %.16g};\n",
+      fprintf(geo, "Point (%d) = {%.16g, %.16g, %.16g, %.16g};\n",
 	      gv->tag(), gv->x(), gv->y(), gv->z(), 
 	      gv->prescribedMeshSizeAtVertex());
     }
@@ -249,7 +249,7 @@ class writeGEdgeGEO {
 	  for(int i = 1; i < ge->minimumDrawSegments(); i++){
 	    double u = umin + (double)i / ge->minimumDrawSegments() * (umax - umin);
 	    GPoint p = ge->point(u);
-	    fprintf(geo, "Point(p%d + %d) = {%.16g, %.16g, %.16g, 1.e+22};\n", 
+	    fprintf(geo, "Point (p%d + %d) = {%.16g, %.16g, %.16g, 1.e+22};\n", 
 		    ge->tag(), i, p.x(), p.y(), p.z());
 	  }
 	  fprintf(geo, "CatmullRom (%d) = {%d", ge->tag(), ge->getBeginVertex()->tag());
@@ -328,24 +328,27 @@ class writePhysicalGroupGEO {
  private :
   FILE *geo;
   int dim;
+  bool printLabels;
   std::map<int, std::string> &oldLabels, &newLabels;
  public :
-  writePhysicalGroupGEO(FILE *fp, int i, 
+  writePhysicalGroupGEO(FILE *fp, int i, bool labels,
 			std::map<int, std::string> &o,
 			std::map<int, std::string> &n)
-    : dim(i), oldLabels(o), newLabels(n)
+    : dim(i), printLabels(labels), oldLabels(o), newLabels(n)
   { 
     geo = fp ? fp : stdout; 
   }
   void operator () (std::pair<const int, std::vector<GEntity *> > &g)
   {
     std::string oldName, newName;
-    if(oldLabels.count(g.first)) {
-      oldName = oldLabels[g.first];
-      fprintf(geo, "%s = %d;\n", oldName.c_str(), g.first);
-    }
-    else if(newLabels.count(g.first)) {
-      newName = newLabels[g.first];
+    if(printLabels){
+      if(oldLabels.count(g.first)) {
+	oldName = oldLabels[g.first];
+	fprintf(geo, "%s = %d;\n", oldName.c_str(), g.first);
+      }
+      else if(newLabels.count(g.first)) {
+	newName = newLabels[g.first];
+      }
     }
 
     switch (dim) {
@@ -369,7 +372,7 @@ class writePhysicalGroupGEO {
   }
 };
 
-int GModel::writeGEO(const std::string &name)
+int GModel::writeGEO(const std::string &name, bool printLabels)
 {
   FILE *fp = fopen(name.c_str(), "w");
 
@@ -395,7 +398,7 @@ int GModel::writeGEO(const std::string &name)
   getPhysicalGroups(groups);
   for(int i = 0; i < 4; i++)
     std::for_each(groups[i].begin(), groups[i].end(), 
-		  writePhysicalGroupGEO(fp, i, labels, physicalNames));
+		  writePhysicalGroupGEO(fp, i, printLabels, labels, physicalNames));
 
   if(fp) fclose(fp);
   return 1;
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 4f2fecbcd33685559e859681439d474b26f4854c..d9f857fecd033aeab1e2ac65583731dc8239bff6 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.11 2007-03-13 09:25:50 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.12 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -2023,24 +2023,24 @@ int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor
     return 0;
   }
 
-  fprintf(fp, "%d\n", faces.size() + regions.size());
+  fprintf(fp, "%d\n", (int)(faces.size() + regions.size()));
   
   for(unsigned int i = 0; i < faces.size(); i++)
     fprintf(fp, "%d %d 1\n", 
-	    faces[i]->transfinite_vertices.size(),
-	    faces[i]->transfinite_vertices[0].size());
+	    (int)faces[i]->transfinite_vertices.size(),
+	    (int)faces[i]->transfinite_vertices[0].size());
 
   for(unsigned int i = 0; i < regions.size(); i++)
     fprintf(fp, "%d %d %d\n", 
-	    regions[i]->transfinite_vertices.size(),
-	    regions[i]->transfinite_vertices[0].size(),
-	    regions[i]->transfinite_vertices[0][0].size());
+	    (int)regions[i]->transfinite_vertices.size(),
+	    (int)regions[i]->transfinite_vertices[0].size(),
+	    (int)regions[i]->transfinite_vertices[0][0].size());
   
   for(unsigned int i = 0; i < faces.size(); i++){
     GFace *gf = faces[i];
     for(int coord = 0; coord < 3; coord++){
-      for(int j = 0; j < gf->transfinite_vertices.size(); j++){
-	for(int k = 0; k < gf->transfinite_vertices[j].size(); k++){
+      for(unsigned int j = 0; j < gf->transfinite_vertices.size(); j++){
+	for(unsigned int k = 0; k < gf->transfinite_vertices[j].size(); k++){
 	  MVertex *v = gf->transfinite_vertices[j][k];
 	  double d = (coord == 0) ? v->x() : (coord == 1) ? v->y() : v->z();
 	  fprintf(fp, "%g ", d * scalingFactor);
@@ -2053,9 +2053,9 @@ int GModel::writeP3D(const std::string &name, bool saveAll, double scalingFactor
   for(unsigned int i = 0; i < regions.size(); i++){
     GRegion *gr = regions[i];
     for(int coord = 0; coord < 3; coord++){
-      for(int j = 0; j < gr->transfinite_vertices.size(); j++){
-	for(int k = 0; k < gr->transfinite_vertices[j].size(); k++){
-	  for(int l = 0; l < gr->transfinite_vertices[j][k].size(); l++){
+      for(unsigned int j = 0; j < gr->transfinite_vertices.size(); j++){
+	for(unsigned int k = 0; k < gr->transfinite_vertices[j].size(); k++){
+	  for(unsigned int l = 0; l < gr->transfinite_vertices[j][k].size(); l++){
 	    MVertex *v = gr->transfinite_vertices[j][k][l];
 	    double d = (coord == 0) ? v->x() : (coord == 1) ? v->y() : v->z();
 	    fprintf(fp, "%g ", d * scalingFactor);
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 99756033e6a9cd6ea01942b19f536271e189a8e3..0989b2848a8a664198c7cd132d4f5a831f3e96a7 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.11 2007-03-11 20:18:58 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.12 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -75,7 +75,7 @@ void MVertex::writeUNV(FILE *fp, double scalingFactor)
   char tmp[128];
   sprintf(tmp, "%25.16E%25.16E%25.16E\n", x() * scalingFactor, 
 	  y() * scalingFactor, z() * scalingFactor);
-  for(int i = 0; i < strlen(tmp); i++) if(tmp[i] == 'E') tmp[i] = 'D';
+  for(unsigned int i = 0; i < strlen(tmp); i++) if(tmp[i] == 'E') tmp[i] = 'D';
   fprintf(fp, tmp);
 }
 
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index f215c3a021f142ce9c9a8111835fd184078ef264..9c5ff830ba92cfe4b3151fd7f0db6bcec62162b1 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.10 2007-03-16 10:03:40 remacle Exp $
+// $Id: HighOrder.cpp,v 1.11 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -601,55 +601,43 @@ void SetOrder1(GModel *m)
     removeHighOrderVertices(*it);
 }
 
-bool straightLine ( std::vector<MVertex *> &l , MVertex *n1, MVertex *n2)
+bool straightLine(std::vector<MVertex*> &l, MVertex *n1, MVertex *n2)
 {
   // x = a * t + b
   // x1 = b
   // x2 = a + b
-  for (int i=0;i<l.size();i++)
-    {
-      MVertex *v = l[i];
-      double b = n1->x();
-      double a = n2->x() - b;
-      double t = (v->x() - b)/a;
-      double by = n1->y();
-      double ay = n2->y() - by;
-      double y = ay * t + by;
-      if (fabs(y-v->y()) > 1.e-11 * CTX.lc)
-	{
-	  //	  printf("coucou %g %g\n",y,v->y());
-	  return false;      
-	}
+  for(unsigned int i = 0; i < l.size(); i++){
+    MVertex *v = l[i];
+    double b = n1->x();
+    double a = n2->x() - b;
+    double t = (v->x() - b) / a;
+    double by = n1->y();
+    double ay = n2->y() - by;
+    double y = ay * t + by;
+    if(fabs(y-v->y()) > 1.e-11 * CTX.lc){
+      return false;      
     }
-
+  }
   return true;
-  
 }
 
-void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ) 
+void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
 {
-  double u[7] = {0,1,0,.5,0,.5,.3333};
-  double v[7] = {0,0,1,.5,.5,0,.3333};
   double j[2][2];  
-
   int n = 3;
-
-  for (int i=0;i<n;i++)
-    {
-      for (int k=0;k<n-i;k++)
-	{
-	  t->jac ( (double)i/(n-1), (double)k/(n-1) , j );
-	  double det = det2x2 ( j );
-	  //	  printf("det = %12.5E\n",det);
-	  minJ = std::min ( det , minJ );
-	  maxJ = std::max ( det , maxJ );
-	}
+  for(int i = 0; i < n; i++){
+    for(int k = 0; k < n - i; k++){
+      t->jac((double)i / (n - 1), (double)k / (n - 1), j);
+      double det = det2x2(j);
+      minJ = std::min(det, minJ);
+      maxJ = std::max(det, maxJ);
     }
+  }
 }
 
-void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
+void smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
 {
-  typedef std::map<std::pair<MVertex*, MVertex*> , std::vector<MElement*> > edge2tris;
+  typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
   edge2tris e2t;
   for(unsigned int i = 0; i < gf->triangles.size(); i++){
     MTriangle *t = gf->triangles[i];
@@ -660,137 +648,104 @@ void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
     }
   }
 
-  for (  edge2tris::iterator it = e2t.begin() ; it != e2t.end() ; ++it)
-    {
-      std::pair<MVertex*, MVertex*> edge = it->first;
-      std::vector<MVertex*> e1,e2,e3,e4,e;
-      std::vector<MElement*> triangles = it->second;
-
-      if (triangles.size() == 2)
-	{
-	  MVertex *n2 = edge.first; 
-	  MVertex *n4 = edge.second;
-	  MTriangle *t1 = (MTriangle*)triangles[0];
-	  MTriangle *t2 = (MTriangle*)triangles[1];
-
-	  MVertex *n1 = t1->getOtherVertex (n2,n4);
-	  MVertex *n3 = t2->getOtherVertex (n2,n4);
-	  if (n1<n2)
-	    e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1,n2)];
-	  else
-	    e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2,n1)];
-	  if (n2<n3)
-	    e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2,n3)];
-	  else
-	    e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3,n2)];
-	  if (n3<n4)
-	    e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3,n4)];
-	  else
-	    e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4,n3)];
-	  if (n4<n1)
-	    e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4,n1)];
-	  else
-	    e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1,n4)];
-	  if (n2<n4)
-	    e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2,n4)];
-	  else
-	    e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4,n2)];
-	  
-	  if ((!straightLine ( e1,n1,n2) ||
-	       !straightLine ( e2,n2,n3) ||
-	       !straightLine ( e3,n3,n4) ||
-	       !straightLine ( e4,n4,n1) ))
-	      //	      &&  straightLine ( e,n2,n4))
-	    {
-	      for (int i=0 ; i < e.size() ; i++)
-		{
-		  double v = (double)(i+1) / (e.size()+1);
-		  double u = 1.-v;
-		  
-// 		  printf("%g %g\n",u,v);
- 
- 		  MVertex *vert  = (n2<n4)?e[i]:e[e.size()-i-1];
- 		  MVertex *vert1 = (n1<n2)?e1[e1.size()-i-1]:e1[i];
- 		  MVertex *vert3 = (n3<n4)?e3[i]:e3[e3.size()-i-1];
-		  MVertex *vert4 = (n4<n1)?e4[e4.size()-i-1]:e4[i];
- 		  MVertex *vert2 = (n2<n3)?e2[i]:e2[e2.size()-i-1];
-		  
-		  
-// 		  printf("n1 %g %g\n",n1->x(),n1->y());
-// 		  printf("n2 %g %g\n",n2->x(),n2->y());
-// 		  printf("n3 %g %g\n",n3->x(),n3->y());
-// 		  printf("n4 %g %g\n",n4->x(),n4->y());
-		  
-// 		  printf("vert1 %g %g\n",vert1->x(),vert1->y());
-// 		  printf("vert2 %g %g\n",vert2->x(),vert2->y());
-// 		  printf("vert3 %g %g\n",vert3->x(),vert3->y());
-// 		  printf("vert4 %g %g\n",vert4->x(),vert4->y());
-		  
-// 		  printf("vert %g %g\n",vert->x(),vert->y());
-		  
-		  
-		  vert->x() = vert->x() + 0.05 * ( 
-		    (1.-u) * vert4->x() + u * vert2->x() +
-		    (1.-v) * vert1->x() + v * vert3->x() -
-		    ( (1.-u)*(1.-v) * n1->x() 
-		      + u * (1.-v) * n2->x() 
-		      + u*v*n3->x() 
-		      + (1.-u) * v * n4->x()) - vert->x());
-		    vert->y() = vert->y() + 0.05 * (
-		    (1.-u) * vert4->y() + u * vert2->y() +
-		    (1.-v) * vert1->y() + v * vert3->y() -
-		    ( (1.-u)*(1.-v) * n1->y() 
-		      + u * (1.-v) * n2->y() 
-		      + u*v*n3->y() 
-		      + (1.-u) * v * n4->y()) - vert->y());
-		  
-		    //		  printf("vert %g %g\n\n",vert->x(),vert->y());
-		}
-	      
-	    }
+  for(edge2tris::iterator it = e2t.begin(); it != e2t.end(); ++it){
+    std::pair<MVertex*, MVertex*> edge = it->first;
+    std::vector<MVertex*> e1, e2, e3, e4, e;
+    std::vector<MElement*> triangles = it->second;
+    if(triangles.size() == 2){
+      MVertex *n2 = edge.first; 
+      MVertex *n4 = edge.second;
+      MTriangle *t1 = (MTriangle*)triangles[0];
+      MTriangle *t2 = (MTriangle*)triangles[1];
+      MVertex *n1 = t1->getOtherVertex(n2, n4);
+      MVertex *n3 = t2->getOtherVertex(n2, n4);
+      if(n1 < n2)
+	e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1, n2)];
+      else
+	e1 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n1)];
+      if(n2 < n3)
+	e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n3)];
+      else
+	e2 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3, n2)];
+      if(n3 < n4)
+	e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n3, n4)];
+      else
+	e3 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n3)];
+      if(n4 < n1)
+	e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n1)];
+      else
+	e4 = edgeVertices[std::make_pair<MVertex*, MVertex*> (n1, n4)];
+      if(n2 < n4)
+	e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n4)];
+      else
+	e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n2)];
+      
+      if((!straightLine(e1, n1, n2) || !straightLine(e2, n2, n3) ||
+	  !straightLine(e3, n3, n4) || !straightLine(e4, n4, n1))){
+	for(unsigned int i = 0; i < e.size(); i++){
+	  double v = (double)(i + 1) / (e.size() + 1);
+	  double u = 1. - v;
+	  MVertex *vert  = (n2 < n4) ? e[i] : e[e.size() - i - 1];
+	  MVertex *vert1 = (n1 < n2) ? e1[e1.size() - i - 1] : e1[i];
+	  MVertex *vert3 = (n3 < n4) ? e3[i] : e3[e3.size() - i - 1];
+	  MVertex *vert4 = (n4 < n1) ? e4[e4.size() - i - 1] : e4[i];
+	  MVertex *vert2 = (n2 < n3) ? e2[i] : e2[e2.size() - i - 1];
+	  vert->x() = vert->x() + 0.05 * ( (1.-u) * vert4->x() + u * vert2->x() +
+					   (1.-v) * vert1->x() + v * vert3->x() -
+					   ( (1.-u)*(1.-v) * n1->x() 
+					     + u * (1.-v) * n2->x() 
+					     + u*v*n3->x() 
+					     + (1.-u) * v * n4->x()) - vert->x());
+	  vert->y() = vert->y() + 0.05 * ( (1.-u) * vert4->y() + u * vert2->y() +
+					   (1.-v) * vert1->y() + v * vert3->y() -
+					   ( (1.-u)*(1.-v) * n1->y() 
+					     + u * (1.-v) * n2->y() 
+					     + u*v*n3->y() 
+					     + (1.-u) * v * n4->y()) - vert->y());
 	}
-    }    
+      }
+    }
+  }    
 }
 
-void checkHighOrderTriangles ( GModel *m )
+void checkHighOrderTriangles(GModel *m)
 {
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
-    {
-      bool twod = true;
-//       for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
-//  	MTriangle *t = (*it)->triangles[i];
-//  	if(t->getVertex(0)->z() != 0.0 || t->getVertex(1)->z() != 0.0 ||t->getVertex(2)->z() != 0.0)twod = false;
-//       }      
-      if (twod)
-	{      
-	  double minJ = 1.e22;
-	  double maxJ = -1.e22;
-	  for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
-	    double minJloc = 1.e22;
-	    double maxJloc = -1.e22;	    
-	    MTriangle *t = (*it)->triangles[i];
-	    if (t->getPolynomialOrder() > 1 && t->getPolynomialOrder() < 6)
-	      {
-		getMinMaxJac (t,minJloc,maxJloc);
-		minJ = std::min(minJ,minJloc);
-		maxJ = std::max(maxJ,maxJloc);
-		if (minJloc*maxJloc < 0)
-		  Msg(GERROR,"Triangle %d %d %d has negative jacobians in GFace %d",t->getVertex(0)->getNum(),
-		      t->getVertex(1)->getNum(),t->getVertex(2)->getNum(),(*it)->tag());
-	      }
-	  }
-	  if (minJ != 1.e22)
-	    {
-	      if (minJ*maxJ < 0)
-		Msg(GERROR,"There exists triangles with negative jacobians in GFace %d",(*it)->tag());
-	      else 
-		Msg(INFO,"No negative jacobians detected on model face %d range = (%g,%g)",(*it)->tag(),minJ,maxJ);	  
-	    }
-	}  
-    }
+  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+    bool twod = true;
+    for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+      MTriangle *t = (*it)->triangles[i];
+      if(t->getVertex(0)->z() != 0.0 || 
+	 t->getVertex(1)->z() != 0.0 ||
+	 t->getVertex(2)->z() != 0.0)
+	twod = false;
+    }      
+    if(twod){ // only perform the test for 2D/plane faces for now
+      double minJ = 1.e22;
+      double maxJ = -1.e22;
+      for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+	double minJloc = 1.e22;
+	double maxJloc = -1.e22;	    
+	MTriangle *t = (*it)->triangles[i];
+	if(t->getPolynomialOrder() > 1 && t->getPolynomialOrder() < 6){
+	  getMinMaxJac (t, minJloc, maxJloc);
+	  minJ = std::min(minJ, minJloc);
+	  maxJ = std::max(maxJ, maxJloc);
+	  if(minJloc * maxJloc < 0)
+	    Msg(WARNING, "Triangle %d (%d %d %d) has negative Jacobian", t->getNum(),
+		t->getVertex(0)->getNum(), t->getVertex(1)->getNum(), t->getVertex(2)->getNum());
+	}
+      }
+      if(minJ != 1.e22){
+	if(minJ * maxJ < 0)
+	  Msg(GERROR, "Some triangles have negative Jacobians in surface %d", (*it)->tag());
+	else 
+	  Msg(INFO, "No negative Jacobians detected on model face %d: range = (%g,%g)",
+	      (*it)->tag(), minJ, maxJ);
+      }
+    }  
+  }
 }  
 
-
 void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
 {
   // replace all the elements in the mesh with second order elements
@@ -826,16 +781,13 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
   for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
     setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
 
-
-  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
-    {
-      if (CTX.mesh.smooth_internal_edges)
-	for (int i=0;i<10;i++)
-	  smoothInternalEdges ( *it , edgeVertices);
-    }
-
-  checkHighOrderTriangles ( m );
-
+  if(CTX.mesh.smooth_internal_edges){
+    for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+      for (int i = 0; i < 10; i++) smoothInternalEdges(*it, edgeVertices);
+  }
+  
+  checkHighOrderTriangles(m);
+  
   double t2 = Cpu();
   Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
   Msg(STATUS1, "Mesh");
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index 0acc05dd76c2425cb3072b4f2ef115e123ede7ef..f69c5c6e30adec4585466138d22b675ffd1646d4 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionExtruded.cpp,v 1.12 2007-03-11 20:18:58 geuzaine Exp $
+// $Id: meshGRegionExtruded.cpp,v 1.13 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -199,7 +199,7 @@ void carveHole(std::vector<T*> &elements, double distance, ANNkd_tree *kdtree)
   ANNdistArray dist = new ANNdist[1];
   std::vector<T*> temp;
   for(unsigned int i = 0; i < elements.size(); i++){
-    for(unsigned int j = 0; j < elements[i]->getNumVertices(); j++){
+    for(int j = 0; j < elements[i]->getNumVertices(); j++){
       MVertex *v = elements[i]->getVertex(j);
       double xyz[3] = {v->x(), v->y(), v->z()};
       kdtree->annkSearch(xyz, 1, index, dist);
diff --git a/Parser/CreateFile.cpp b/Parser/CreateFile.cpp
index 520b58c5495ee3c984f112c932cbf597a0891f86..1f00105eca0d9a99b572f279c40451e51cf18b8a 100644
--- a/Parser/CreateFile.cpp
+++ b/Parser/CreateFile.cpp
@@ -1,4 +1,4 @@
-// $Id: CreateFile.cpp,v 1.15 2007-03-11 20:18:58 geuzaine Exp $
+// $Id: CreateFile.cpp,v 1.16 2007-03-18 12:05:16 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -191,7 +191,7 @@ void CreateOutputFile(char *filename, int format)
     break;
 
   case FORMAT_GEO:
-    GMODEL->writeGEO(name);
+    GMODEL->writeGEO(name, CTX.print.geo_labels);
     break;
 
 #if defined(HAVE_FLTK)
diff --git a/Parser/Makefile b/Parser/Makefile
index 73bd9a10b972b51835a8db788bf9a105bd2eee90..15e421393b2403f1d828dbfacaa514e38fbdea28 100644
--- a/Parser/Makefile
+++ b/Parser/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.116 2007-02-26 08:25:46 geuzaine Exp $
+# $Id: Makefile,v 1.117 2007-03-18 12:05:16 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -131,9 +131,9 @@ OpenFile.o: OpenFile.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../Geo/SBoundingBox3d.h Parser.h OpenFile.h ../Common/CommandLine.h \
   ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
   ../Common/SmoothData.h ../Common/AdaptiveViews.h ../Common/GmshMatrix.h \
-  ../Graphics/ReadImg.h ../Common/OS.h ../Common/GmshUI.h \
-  ../Graphics/Draw.h ../Graphics/SelectBuffer.h ../Fltk/GUI.h \
-  ../Fltk/Opengl_Window.h ../Fltk/Colorbar_Window.h \
+  ../Graphics/ReadImg.h ../Common/OS.h ../Mesh/HighOrder.h \
+  ../Common/GmshUI.h ../Graphics/Draw.h ../Graphics/SelectBuffer.h \
+  ../Fltk/GUI.h ../Fltk/Opengl_Window.h ../Fltk/Colorbar_Window.h \
   ../Fltk/Popup_Button.h ../Fltk/SpherePosition_Widget.h
 CreateFile.o: CreateFile.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/configure b/configure
index a3dcbfb1e22fb6a453cd0313982b42bb211b5662..84e4a07c2027cf4d8089b96704575f6a0a3f5921 100755
--- a/configure
+++ b/configure
@@ -2788,7 +2788,7 @@ case "$UNAME" in
       UNAME="${UNAME}-no-cygwin"
       CC="${CC} -mno-cygwin"
       CXX="${CXX} -mno-cygwin"
-      LINKER="${LINKER} -mno-cygwin"
+                  LINKER="${LINKER} -mno-cygwin -Wl,--stack,10000000"
     fi
     ;;
 esac
diff --git a/configure.in b/configure.in
index 6d7c61656eb22cc83c3cbf3da7d49998ed5834a0..f10bc1b023425a653c547c0d3b3817cbec5f9a25 100644
--- a/configure.in
+++ b/configure.in
@@ -1,4 +1,4 @@
-dnl $Id: configure.in,v 1.121 2007-02-02 23:50:33 geuzaine Exp $
+dnl $Id: configure.in,v 1.122 2007-03-18 12:05:16 geuzaine Exp $
 dnl
 dnl Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 dnl
@@ -145,7 +145,9 @@ case "$UNAME" in
       UNAME="${UNAME}-no-cygwin"
       CC="${CC} -mno-cygwin"
       CXX="${CXX} -mno-cygwin"
-      LINKER="${LINKER} -mno-cygwin"
+      dnl need to increase default stack size passed to ld for
+      dnl recursive tet classification in new 3D Delaunay algo
+      LINKER="${LINKER} -mno-cygwin -Wl,--stack,10000000"
     fi
     ;;
 esac
diff --git a/doc/texinfo/opt_mesh.texi b/doc/texinfo/opt_mesh.texi
index 8ba288821e55c6dce29721cbb40599849231ab72..536f5c78fafd6dc8b3da08759a53df1ff96bacfb 100644
--- a/doc/texinfo/opt_mesh.texi
+++ b/doc/texinfo/opt_mesh.texi
@@ -95,10 +95,15 @@ Default value: @code{0}@*
 Saved in: @code{General.OptionsFileName}
 
 @item Mesh.ElementOrder
-Element order (1=linear elements, 2=quadratic elements)@*
+Element order (1=linear elements, N (<6) = elements of higher order)@*
 Default value: @code{1}@*
 Saved in: @code{General.OptionsFileName}
 
+@item Mesh.SmoothInternalEdges
+Number of smoothing steps of internal edges for high order meshes@*
+Default value: @code{0}@*
+Saved in: @code{General.OptionsFileName}
+
 @item Mesh.Explode
 Element shrinking factor (between 0 and 1)@*
 Default value: @code{1}@*
diff --git a/doc/texinfo/opt_print.texi b/doc/texinfo/opt_print.texi
index a5045b7f7a082bdadd8a926a37555edf123a136c..32a7898460dfa6e05ac57a542c81ce8cedc7d30f 100644
--- a/doc/texinfo/opt_print.texi
+++ b/doc/texinfo/opt_print.texi
@@ -49,6 +49,11 @@ File format (10=automatic)@*
 Default value: @code{10}@*
 Saved in: @code{General.OptionsFileName}
 
+@item Print.GeoLabels
+Save labels in unrolled Gmsh geometries@*
+Default value: @code{1}@*
+Saved in: @code{General.OptionsFileName}
+
 @item Print.GifDither
 Apply dithering to GIF output@*
 Default value: @code{0}@*
diff --git a/doc/texinfo/opt_view.texi b/doc/texinfo/opt_view.texi
index 24ab0a3388988119210fb44df1f30726da7f57f8..1ab71ff8f2694dc4925690de336fa113aee3ae15 100644
--- a/doc/texinfo/opt_view.texi
+++ b/doc/texinfo/opt_view.texi
@@ -396,42 +396,42 @@ Saved in: @code{General.OptionsFileName}
 
 @item View.Max
 Maximum value in the view (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.MaxX
 Maximum view coordinate along the X-axis (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.MaxY
 Maximum view coordinate along the Y-axis (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.MaxZ
 Maximum view coordinate along the Z-axis (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.Min
 Minimum value in the view (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.MinX
 Minimum view coordinate along the X-axis (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.MinY
 Minimum view coordinate along the Y-axis (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.MinZ
 Minimum view coordinate along the Z-axis (read-only)@*
-Default value: @code{0}@*
+Default value: @code{-2.22195e-146}@*
 Saved in: @code{-}
 
 @item View.NbIso