diff --git a/Common/Context.h b/Common/Context.h
index 0cde4d0bfb3db673780e36fb57533e01a5d84ae8..a775c1bb1bec672b005a124cfb7e4a7b71bbbbb8 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -150,7 +150,7 @@ public :
     int extrude_spline_points, old_newreg;
     double normals, tangents;
     double scaling_factor;
-    int auto_coherence;
+    int auto_coherence, highlight_orphans;
     double tolerance;
     double snap[3];
     int occ_fix_small_edges, occ_fix_small_faces, occ_sew_faces;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 214b50480e4c098d33cd31ec8fd0687d781b18a0..28cef443dbf6e3afcdcf665b00971605f74d85fd 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -781,6 +781,9 @@ StringXNumber GeometryOptions_Number[] = {
   { F|O, "ExtrudeSplinePoints" , opt_geometry_extrude_spline_points, 5. ,
     "Number of control points for splines created during extrusion" },
 
+  { F|O, "HighlightOrphans" , opt_geometry_highlight_orphans, 0. , 
+    "Highlight orphan entities (lines connected to a single surface, etc.)?" }, 
+
   { F|O, "Light" , opt_geometry_light , 1. , 
     "Enable lighting for the geometry" },
   { F|O, "Lines" , opt_geometry_lines , 1. , 
@@ -849,7 +852,7 @@ StringXNumber MeshOptions_Number[] = {
   { F|O, "Algorithm" , opt_mesh_algo2d , ALGO_2D_MESHADAPT ,
     "2D mesh algorithm (1=meshadapt, 2=delaunay)" }, 
   { F|O, "Algorithm3D" , opt_mesh_algo3d , ALGO_3D_NETGEN ,
-    "3D mesh algorithm (1=delaunay, 4=netgen)" }, 
+    "3D mesh algorithm (1=delaunay, 4=netgen, 5=tetgen)" }, 
   { F|O, "AngleSmoothNormals" , opt_mesh_angle_smooth_normals , 30.0 ,
     "Threshold angle below which normals are not smoothed" }, 
 
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 625e9ab12cc8498b40df65f5236c7852c885fb12..9f00485660bf2c024d5b77bd818ee31e44f96d50 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.329 2007-01-18 13:18:42 geuzaine Exp $
+// $Id: Options.cpp,v 1.330 2007-01-22 16:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -3775,6 +3775,17 @@ double opt_geometry_auto_coherence(OPT_ARGS_NUM)
   return CTX.geom.auto_coherence;
 }
 
+double opt_geometry_highlight_orphans(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.geom.highlight_orphans = (int)val;
+#if defined(HAVE_FLTK)
+  if(WID && (action & GMSH_GUI))
+    WID->geo_butt[10]->value(CTX.geom.highlight_orphans);
+#endif
+  return CTX.geom.highlight_orphans;
+}
+
 double opt_geometry_tolerance(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
@@ -4764,6 +4775,10 @@ double opt_mesh_second_order_incomplete(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
     CTX.mesh.second_order_incomplete = (int)val;
+#if defined(HAVE_FLTK)
+  if(WID && (action & GMSH_GUI))
+    WID->mesh_butt[4]->value(CTX.mesh.second_order_incomplete);
+#endif
   return CTX.mesh.second_order_incomplete;
 }
 
diff --git a/Common/Options.h b/Common/Options.h
index 8a82b7902e0905879af50a5a3011f6d794e91b6d..adeed0ef82a00636475ac3fda94e387cbc4baa2e 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -379,6 +379,7 @@ double opt_general_light51(OPT_ARGS_NUM);
 double opt_general_light52(OPT_ARGS_NUM);
 double opt_general_light53(OPT_ARGS_NUM);
 double opt_geometry_auto_coherence(OPT_ARGS_NUM);
+double opt_geometry_highlight_orphans(OPT_ARGS_NUM);
 double opt_geometry_tolerance(OPT_ARGS_NUM);
 double opt_geometry_normals(OPT_ARGS_NUM);
 double opt_geometry_tangents(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index d8222ccaca68c20bae2a5d11b8c057a5ed5197bf..3cee5aa787facd955380b2f885638d473a222610 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.505 2007-01-18 09:12:44 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.506 2007-01-22 16:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -1006,6 +1006,7 @@ void geometry_options_ok_cb(CALLBACK_ARGS)
   opt_geometry_surfaces_num(0, GMSH_SET, WID->geo_butt[6]->value());
   opt_geometry_volumes_num(0, GMSH_SET, WID->geo_butt[7]->value());
   opt_geometry_auto_coherence(0, GMSH_SET, WID->geo_butt[8]->value());
+  opt_geometry_highlight_orphans(0, GMSH_SET, WID->geo_butt[10]->value());
   opt_geometry_light(0, GMSH_SET, WID->geo_butt[9]->value());
 
   opt_geometry_normals(0, GMSH_SET, WID->geo_value[0]->value());
@@ -1046,8 +1047,11 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
     }
   }
 
+  opt_mesh_reverse_all_normals(0, GMSH_SET, WID->mesh_butt[0]->value());
+  opt_mesh_lc_from_curvature(0, GMSH_SET, WID->mesh_butt[1]->value());
   opt_mesh_optimize(0, GMSH_SET, WID->mesh_butt[2]->value());
   opt_mesh_order(0, GMSH_SET, WID->mesh_butt[3]->value()? 2 : 1);
+  opt_mesh_second_order_incomplete(0, GMSH_SET, WID->mesh_butt[4]->value());
   opt_mesh_constrained_bgmesh(0, GMSH_SET, WID->mesh_butt[5]->value());
   opt_mesh_points(0, GMSH_SET, WID->mesh_butt[6]->value());
   opt_mesh_lines(0, GMSH_SET, WID->mesh_butt[7]->value());
@@ -1070,10 +1074,8 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
   opt_mesh_cut_plane_only_volume(0, GMSH_SET, WID->mesh_butt[23]->value());
   opt_mesh_light(0, GMSH_SET, WID->mesh_butt[17]->value());
   opt_mesh_light_two_side(0, GMSH_SET, WID->mesh_butt[18]->value());
-  opt_mesh_reverse_all_normals(0, GMSH_SET, WID->mesh_butt[0]->value());
   opt_mesh_smooth_normals(0, GMSH_SET, WID->mesh_butt[19]->value());
   opt_mesh_light_lines(0, GMSH_SET, WID->mesh_butt[20]->value());
-  opt_mesh_lc_from_curvature(0, GMSH_SET, WID->mesh_butt[1]->value());
 
   opt_mesh_nb_smoothing(0, GMSH_SET, WID->mesh_value[0]->value());
   opt_mesh_lc_factor(0, GMSH_SET, WID->mesh_value[2]->value());
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 24a7ab3fe7e92856c8f2f50ad59dac1ba3115fe0..4062120690d6d64cdc801856b6b4de27c3cf5422 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.592 2007-01-18 09:12:45 geuzaine Exp $
+// $Id: GUI.cpp,v 1.593 2007-01-22 16:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -2210,10 +2210,14 @@ void GUI::create_option_window()
       Fl_Group *o = new Fl_Group(L + WB, WB + BH, width - 2 * WB, height - 2 * WB - BH, "Color");
       o->hide();
 
-      Fl_Scroll *s = new Fl_Scroll(L + 2 * WB, 2 * WB + 1 * BH, IW + 20, height - 4 * WB - 1 * BH);
+      geo_butt[10] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 1 * BH, BW, BH, "Highlight orphan entities");
+      geo_butt[10]->type(FL_TOGGLE_BUTTON);
+      geo_butt[10]->callback(geometry_options_ok_cb);
+
+      Fl_Scroll *s = new Fl_Scroll(L + 2 * WB, 2 * WB + 2 * BH, IW + 20, height - 4 * WB - 2 * BH);
       int i = 0;
       while(GeometryOptions_Color[i].str) {
-        geo_col[i] = new Fl_Button(L + 2 * WB, 2 * WB + (1 + i) * BH, IW, BH, GeometryOptions_Color[i].str);
+        geo_col[i] = new Fl_Button(L + 2 * WB, 2 * WB + (2 + i) * BH, IW, BH, GeometryOptions_Color[i].str);
         geo_col[i]->callback(color_cb, (void *)GeometryOptions_Color[i].function);
         i++;
       }
@@ -2303,6 +2307,10 @@ void GUI::create_option_window()
       mesh_butt[3]->type(FL_TOGGLE_BUTTON);
       mesh_butt[3]->callback(mesh_options_ok_cb);
 
+      mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 10 * BH, BW, BH, "Use incomplete second order elements (8-node quads, etc.)");
+      mesh_butt[4]->type(FL_TOGGLE_BUTTON);
+      mesh_butt[4]->callback(mesh_options_ok_cb);
+
       o->end();
     }
 
diff --git a/Fltk/Main.cpp b/Fltk/Main.cpp
index 17c5a2f5c3a980dbd2e7dea6221489b46c955c9d..daaed384ef39db2d1f7c2b0b364c4544088e0713 100644
--- a/Fltk/Main.cpp
+++ b/Fltk/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.102 2006-12-16 15:44:28 geuzaine Exp $
+// $Id: Main.cpp,v 1.103 2007-01-22 16:31:42 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -138,6 +138,8 @@ int main(int argc, char *argv[])
       }
       else if(CTX.batch == -1)
         CreateOutputFile(CTX.output_filename, FORMAT_GEO);
+      else if(CTX.batch == -2)
+	GMODEL->checkMeshCoherence();
       exit(0);
     }
   }
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 512798b9c41f6b56d1522c5c892f55dbdb9b857d..a64f6f6a6166a143b26a80216b90a74a748b7984 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -98,11 +98,11 @@ class GEdge : public GEntity {
   virtual std::string getAdditionalInfoString();
 
   // tells if the edge is a 3D edge (in opposition with a trimmed curve on a surface)
-  virtual bool is3D () const {return true;}
+  virtual bool is3D() const {return true;}
 
   // the length of the model edge
-  inline double length () const {return _length;}
-  inline void   setLength (const double l) {_length = l;}
+  inline double length() const {return _length;}
+  inline void setLength(const double l) {_length = l;}
 
   // one can impose the mesh size at an edge
   virtual double prescribedMeshSizeAtVertex() const {return meshAttributes.meshSize;}
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index ba3225263939303a5cf75bd6765c7caf87319993..e3f31b1f4ce0b35b6a3d823fdef46b83bfd48a48 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.28 2007-01-18 13:18:42 geuzaine Exp $
+// $Id: GModel.cpp,v 1.29 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -296,3 +296,90 @@ void GModel::deleteMeshPartitions()
   meshPartitions.clear();
 }
 
+static int checkVertices(std::vector<MVertex*> &vertices,
+			 std::set<MVertex*, MVertexLessThanLexicographic> &pos)
+{
+  int num = 0;
+  for(unsigned int i = 0; i < vertices.size(); i++){
+    MVertex *v = vertices[i];
+    std::set<MVertex*, MVertexLessThanLexicographic>::iterator it = pos.find(v);
+    if(it == pos.end()){
+      pos.insert(v);
+    }
+    else{
+      Msg(INFO, "Vertices %d and %d have identical position (%g, %g, %g)", 
+	  (*it)->getNum(), v->getNum(), v->x(), v->y(), v->z());
+      num++;
+    }
+  }
+  return num;
+}
+
+template <class T>
+static int checkElements(std::vector<T*> &elements,
+			 std::set<MElement*, MElementLessThanLexicographic> &pos)
+{
+  int num = 0;
+  for(unsigned int i = 0; i < elements.size(); i++){
+    MElement *e = elements[i];
+    std::set<MElement*, MVertexLessThanLexicographic>::iterator it = pos.find(e);
+    if(it == pos.end()){
+      pos.insert(e);
+    }
+    else{
+      Msg(INFO, "Elements %d and %d have identical barycenter",
+	  (*it)->getNum(), e->getNum());
+      num++;
+    }
+  }
+  return num;
+}
+
+void GModel::checkMeshCoherence()
+{
+  int numEle = numElement();
+  if(!numEle) return;
+
+  Msg(INFO, "Checking mesh coherence (%d elements)", numEle);
+
+  // check for duplicate mesh vertices
+  {
+    double old_tol = MVertexLessThanLexicographic::tolerance; 
+    MVertexLessThanLexicographic::tolerance = CTX.geom.tolerance;
+    std::set<MVertex*, MVertexLessThanLexicographic> pos;
+    int num = 0;
+    for(viter it = firstVertex(); it != lastVertex(); ++it)
+      num += checkVertices((*it)->mesh_vertices, pos);
+    for(eiter it = firstEdge(); it != lastEdge(); ++it)
+      num += checkVertices((*it)->mesh_vertices, pos);
+    for(fiter it = firstFace(); it != lastFace(); ++it)
+      num += checkVertices((*it)->mesh_vertices, pos);
+    for(riter it = firstRegion(); it != lastRegion(); ++it)
+      num += checkVertices((*it)->mesh_vertices, pos);
+    if(num) Msg(WARNING, "%d duplicate vertices", num);
+    MVertexLessThanLexicographic::tolerance = old_tol;
+  }
+
+  // check for duplicate elements
+  {
+    double old_tol = MElementLessThanLexicographic::tolerance; 
+    MElementLessThanLexicographic::tolerance = CTX.geom.tolerance;
+    std::set<MElement*, MElementLessThanLexicographic> pos;
+    int num = 0;
+    for(eiter it = firstEdge(); it != lastEdge(); ++it)
+      num += checkElements((*it)->lines, pos);
+    for(fiter it = firstFace(); it != lastFace(); ++it){
+      num += checkElements((*it)->triangles, pos);
+      num += checkElements((*it)->quadrangles, pos);
+    }
+    for(riter it = firstRegion(); it != lastRegion(); ++it){
+      num += checkElements((*it)->tetrahedra, pos);
+      num += checkElements((*it)->hexahedra, pos);
+      num += checkElements((*it)->prisms, pos);
+      num += checkElements((*it)->pyramids, pos);
+    }
+    if(num) Msg(WARNING, "%d duplicate elements", num);
+    MElementLessThanLexicographic::tolerance = old_tol;
+  }
+
+}
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 20365702e15c1d36c8d9727ea4de5b5bd7e5ee38..99d25366e206fc0fc69dbdca40c96a5b8813d6bb 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -145,6 +145,9 @@ class GModel
   // Deletes all the partitions
   virtual void deleteMeshPartitions();
 
+  // Performs various coherence tests on the mesh
+  virtual void checkMeshCoherence();
+
   // A container for smooth normals
   smooth_normals *normals;
 
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 826be7b6165f948309bb6c8c510973ae5dff3b16..3558bf4939683920b485720fac1f2aab59a93885 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.28 2007-01-18 09:12:45 geuzaine Exp $
+// $Id: MElement.cpp,v 1.29 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -118,6 +118,7 @@ int quadfaces_pyramid[1][4] = {
 };
 
 int MElement::_globalNum = 0;
+double MElementLessThanLexicographic::tolerance = 1.e-6;
 
 double MElement::minEdge()
 {
diff --git a/Geo/MElement.h b/Geo/MElement.h
index ade08451b7652b87a308d3b1f0d731ac2b7805c9..349c850ed71110eba88cb5e941e2396e4b342f0b 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -160,6 +160,22 @@ class MElement
   virtual char *getStringForBDF() = 0;
 };
 
+class MElementLessThanLexicographic{
+ public:
+  static double tolerance;
+  bool operator()(MElement *e1, MElement *e2) const
+  {
+    SPoint3 b1 = e1->barycenter();
+    SPoint3 b2 = e2->barycenter();
+    if(b1.x() - b2.x() >  tolerance) return true;
+    if(b1.x() - b2.x() < -tolerance) return false;
+    if(b1.y() - b2.y() >  tolerance) return true;
+    if(b1.y() - b2.y() < -tolerance) return false;
+    if(b1.z() - b2.z() >  tolerance) return true;
+    return false;
+  }
+};
+
 class MLine : public MElement {
  protected:
   MVertex *_v[2];
@@ -348,48 +364,43 @@ class MTriangle6 : public MTriangle {
 };
 
 template <class T> 
-void sort3 ( T * t [3])
+void sort3(T *t[3])
 {
   T *temp;
-  if (t[0] > t[1])
-    {
-      temp = t[1];
-      t[1] = t[0];
-      t[0] = temp;
-    }
-  if (t[1] > t[2])
-    {
-      temp = t[2];
-      t[2] = t[1];
-      t[1] = temp;
-    }
-  if (t[0] > t[1])
-    {
-      temp = t[1];
-      t[1] = t[0];
-      t[0] = temp;
-    }
+  if(t[0] > t[1]){
+    temp = t[1];
+    t[1] = t[0];
+    t[0] = temp;
+  }
+  if(t[1] > t[2]){
+    temp = t[2];
+    t[2] = t[1];
+    t[1] = temp;
+  }
+  if(t[0] > t[1]){
+    temp = t[1];
+    t[1] = t[0];
+    t[0] = temp;
+  }
 }
 
 struct compareMTriangleLexicographic
 {
-  bool operator () ( MTriangle *t1 , MTriangle *t2 ) const
+  bool operator () (MTriangle *t1, MTriangle *t2) const
   {
-    MVertex *_v1[3] = {t1->getVertex(0),t1->getVertex(1),t1->getVertex(2)};
-    MVertex *_v2[3] = {t2->getVertex(0),t2->getVertex(1),t2->getVertex(2)};
-
+    MVertex *_v1[3] = {t1->getVertex(0), t1->getVertex(1), t1->getVertex(2)};
+    MVertex *_v2[3] = {t2->getVertex(0), t2->getVertex(1), t2->getVertex(2)};
     sort3(_v1);
     sort3(_v2);
-    if (_v1[0] < _v2[0]) return true;
-    if (_v1[0] > _v2[0]) return false;
-    if (_v1[1] < _v2[1]) return true;
-    if (_v1[1] > _v2[1]) return false;
-    if (_v1[2] < _v2[2]) return true;
+    if(_v1[0] < _v2[0]) return true;
+    if(_v1[0] > _v2[0]) return false;
+    if(_v1[1] < _v2[1]) return true;
+    if(_v1[1] > _v2[1]) return false;
+    if(_v1[2] < _v2[2]) return true;
     return false;
   }
 };
 
-
 class MQuadrangle : public MElement {
  protected:
   MVertex *_v[4];
diff --git a/Geo/Makefile b/Geo/Makefile
index a4b8193d30bd45d948fcb3a29a64bbf434eada40..daf028ef6b301750041c15418008acf88d7f07fa 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.123 2006-12-21 17:10:15 geuzaine Exp $
+# $Id: Makefile,v 1.124 2007-01-22 16:31:43 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -109,13 +109,13 @@ gmshVertex.o: gmshVertex.cpp GFace.h GPoint.h GEntity.h Range.h SPoint3.h \
   MVertex.h SPoint2.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \
   MFace.h ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
   ExtrudeParams.h Pair.h gmshVertex.h Geo.h ../DataStr/Tree.h \
-  ../DataStr/avl.h GeoInterpolation.h
+  ../DataStr/avl.h GeoInterpolation.h ../Common/Message.h
 gmshEdge.o: gmshEdge.cpp GFace.h GPoint.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h ../Common/GmshDefines.h GEdgeLoop.h GEdge.h GVertex.h \
   MVertex.h SPoint2.h SVector3.h MElement.h MEdge.h ../Common/Hash.h \
   MFace.h ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
   ExtrudeParams.h Pair.h gmshEdge.h Geo.h ../DataStr/Tree.h \
-  ../DataStr/avl.h gmshVertex.h GeoInterpolation.h
+  ../DataStr/avl.h gmshVertex.h GeoInterpolation.h ../Common/Message.h
 gmshFace.o: gmshFace.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/Geo/OCCVertex.cpp b/Geo/OCCVertex.cpp
index fe92a1eb88005e696f0aac4e6fe53b61b499d2ba..004104904963f8521596d6c426062eb2c22517b4 100644
--- a/Geo/OCCVertex.cpp
+++ b/Geo/OCCVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCVertex.cpp,v 1.10 2007-01-16 11:31:41 geuzaine Exp $
+// $Id: OCCVertex.cpp,v 1.11 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -85,6 +85,9 @@ SPoint2 OCCVertex::reparamOnFace(GFace *gf, int dir) const
     }
     ++it;
   }
+
+  // normally never here
+  return GVertex::reparamOnFace(gf, dir);
 }
 
 #endif
diff --git a/Graphics/Geom.cpp b/Graphics/Geom.cpp
index 82c139d761b154742404cef188459c800f245963..f6749ae19230ba1ff175c6c13851e66cab521f64 100644
--- a/Graphics/Geom.cpp
+++ b/Graphics/Geom.cpp
@@ -1,4 +1,4 @@
-// $Id: Geom.cpp,v 1.127 2006-12-15 03:15:32 geuzaine Exp $
+// $Id: Geom.cpp,v 1.128 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -52,6 +52,14 @@ class drawGVertex {
       glColor4ubv((GLubyte *) & CTX.color.geom.point);
     }
     
+    if(CTX.geom.highlight_orphans){
+      std::list<GEdge*> edges = v->edges();
+      if(edges.size() == 0)
+	glColor3d(1., 0., 0.);
+      else if(edges.size() == 1)
+	glColor3d(1., 0.6, 0.);
+    }
+
     if(CTX.geom.points) {
       if(CTX.geom.point_type > 0) {
 	if(v->getSelection())
@@ -109,6 +117,14 @@ class drawGEdge {
       glColor4ubv((GLubyte *) & CTX.color.geom.line);
     }
     
+    if(CTX.geom.highlight_orphans){
+      std::list<GFace*> faces = e->faces();
+      if(faces.size() == 0)
+	glColor3d(1., 0., 0.);
+      else if(faces.size() == 1)
+	glColor3d(1., 0.6, 0.);
+    }
+
     Range<double> t_bounds = e->parBounds(0);
     double t_min = t_bounds.low();
     double t_max = t_bounds.high();
@@ -386,7 +402,7 @@ public :
       gl2psLineWidth(CTX.geom.line_width / 2. * CTX.print.eps_line_width_factor);
       glColor4ubv((GLubyte *) & CTX.color.geom.surface);
     }
-    
+
     if(f->geomType() == GEntity::Plane)
       _drawPlaneGFace(f);
     else if(f->geomType() == GEntity::ProjectionSurface)
diff --git a/Graphics/Makefile b/Graphics/Makefile
index 5819d155e1a346a05f616e9b85f3fc47b8e1741f..b9eb25515dd8cab06fa8bbdcb051ca17b9a308f0 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.109 2006-11-27 22:22:15 geuzaine Exp $
+# $Id: Makefile,v 1.110 2007-01-22 16:31:43 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -107,7 +107,7 @@ Mesh.o: Mesh.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
   ../Common/GmshMatrix.h ../Geo/MRep.h ../Geo/GEdge.h ../Geo/GFace.h \
   ../Geo/GRegion.h ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MElement.h \
-  ../Common/OS.h gl2ps.h tc.h
+  ../Common/OS.h gl2ps.h
 Geom.o: Geom.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
   ../DataStr/List.h ../DataStr/Tree.h ../Common/GmshUI.h Draw.h \
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 0afb581af20def482a8a1875e3241d67fd416464..91c846d8ab08e6bfbea9caeac5438c84dbcc5702 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.191 2006-12-05 18:34:58 geuzaine Exp $
+// $Id: Mesh.cpp,v 1.192 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -27,7 +27,6 @@
 #include "MRep.h"
 #include "OS.h"
 #include "gl2ps.h"
-#include "tc.h"
 
 extern GModel *GMODEL;
 extern Context_T CTX;
@@ -178,14 +177,7 @@ static void drawElementLabels(GEntity *e, std::vector<T*> &elements,
 	sprintf(str, "%d", e->tag());
       else
 	sprintf(str, "%d", ele->getNum());
-      if(e->dim() == 1){
-	double offset = 0.3 * (CTX.mesh.line_width + CTX.gl_fontsize) * CTX.pixel_equiv_x;
-	glRasterPos3d(pc.x() + offset / CTX.s[0], 
-		      pc.y() + offset / CTX.s[1], 
-		      pc.z() + offset / CTX.s[2]);
-      }
-      else
-	glRasterPos3d(pc.x(), pc.y(), pc.z());
+      glRasterPos3d(pc.x(), pc.y(), pc.z());
       Draw_String(str);
     }
   }
@@ -244,10 +236,7 @@ static void drawVertexLabel(GEntity *e, MVertex *v, int partition=-1)
     sprintf(str, "%d", e->tag());
   else
     sprintf(str, "%d", v->getNum());
-  double offset = 0.3 * (CTX.mesh.point_size + CTX.gl_fontsize) * CTX.pixel_equiv_x;
-  glRasterPos3d(v->x() + offset / CTX.s[0],
-		v->y() + offset / CTX.s[1],
-		v->z() + offset / CTX.s[2]);
+  glRasterPos3d(v->x(), v->y(), v->z());
   Draw_String(str);
 }
 
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 6bc744d029610775147303d9a4dc82fd618f648d..cfd021039d83774ceeab77409fc43d7967f6b701 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.159 2007-01-12 13:16:59 remacle Exp $
+# $Id: Makefile,v 1.160 2007-01-22 16:31:43 geuzaine Exp $
 #
 # Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 #
@@ -122,25 +122,26 @@ meshGEdgeExtruded.o: meshGEdgeExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
   ../Common/SmoothNormals.h ../Common/Message.h
-meshGFace.o: meshGFace.cpp meshGFace.h DivideAndConquer.h \
-  BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
-  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
-  ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
+meshGFace.o: meshGFace.cpp meshGFace.h meshGFaceDelaunayInsertion.h \
+  ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
+  ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h \
+  ../Geo/SVector3.h ../Numeric/Numeric.h ../Common/Context.h \
+  ../DataStr/List.h DivideAndConquer.h BackgroundMesh.h ../Geo/GVertex.h \
+  ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/MVertex.h \
   ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
-  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h \
-  ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/MVertex.h ../Geo/MEdge.h \
-  ../Geo/MVertex.h ../Geo/SVector3.h ../Common/Hash.h ../Geo/MFace.h \
-  ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/Numeric.h \
-  ../Common/Context.h ../DataStr/List.h ../Geo/ExtrudeParams.h \
-  ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/GEdgeLoop.h \
-  ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
-  ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/MRep.h ../Geo/GEdge.h \
-  ../Geo/GFace.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h \
-  ../Geo/ExtrudeParams.h ../Geo/MVertex.h ../Geo/MEdge.h \
-  ../Geo/MElement.h ../Common/VertexArray.h ../Common/Message.h \
-  ../Common/OS.h BDS.h ../Common/Views.h ../Common/ColorTable.h \
-  ../Common/VertexArray.h ../Common/SmoothNormals.h \
-  ../Common/AdaptiveViews.h ../Common/GmshMatrix.h
+  ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+  ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/GFace.h ../Geo/GPoint.h \
+  ../Geo/GEntity.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/MElement.h \
+  ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h \
+  ../Geo/MRep.h ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h \
+  ../Geo/GEntity.h ../Geo/MElement.h ../Geo/ExtrudeParams.h \
+  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MElement.h \
+  ../Common/VertexArray.h ../Common/Message.h ../Common/OS.h BDS.h \
+  ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
+  ../Common/GmshMatrix.h
 meshGFaceTransfinite.o: meshGFaceTransfinite.cpp meshGFace.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
   ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshDefines.h \
@@ -169,6 +170,22 @@ meshGFaceExtruded.o: meshGFaceExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
   ../Common/SmoothNormals.h ../Common/Message.h
+meshGFaceDelaunayInsertion.o: meshGFaceDelaunayInsertion.cpp BDS.h \
+  ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+  ../Common/GmshDefines.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/MVertex.h \
+  ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/MElement.h \
+  ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+  ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
+  ../Geo/ExtrudeParams.h ../Geo/MElement.h ../Geo/SPoint2.h \
+  ../Geo/SVector3.h ../Geo/Pair.h ../Geo/ExtrudeParams.h \
+  ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Common/AdaptiveViews.h \
+  ../Common/GmshMatrix.h BackgroundMesh.h meshGFaceDelaunayInsertion.h \
+  ../Common/Message.h
 meshGRegion.o: meshGRegion.cpp meshGRegion.h \
   meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
   ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
@@ -184,7 +201,9 @@ meshGRegion.o: meshGRegion.cpp meshGRegion.h \
   ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
-  ../Common/SmoothNormals.h ../Geo/MRep.h ../Geo/GEdge.h ../Geo/GFace.h \
+  ../Common/SmoothNormals.h ../Geo/gmshRegion.h ../Geo/Geo.h \
+  ../DataStr/Tree.h ../DataStr/avl.h ../Geo/ExtrudeParams.h \
+  ../Geo/GRegion.h ../Geo/MRep.h ../Geo/GEdge.h ../Geo/GFace.h \
   ../Geo/GRegion.h ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MElement.h \
   ../Common/VertexArray.h ../Common/Message.h ../Common/OS.h BDS.h \
   ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
@@ -196,9 +215,16 @@ meshGRegionDelaunayInsertion.o: meshGRegionDelaunayInsertion.cpp \
   ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
   ../Common/Hash.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
   ../Numeric/Numeric.h ../Common/Context.h ../DataStr/List.h \
-  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
-  ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/MElement.h \
-  ../Geo/ExtrudeParams.h ../Common/Message.h
+  ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+  ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+  ../Geo/MVertex.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
+  ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+  ../Geo/SPoint2.h ../Geo/MElement.h ../Geo/ExtrudeParams.h \
+  ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/GEdgeLoop.h \
+  ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+  ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
+  ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
+  ../Common/SmoothNormals.h ../Common/Message.h
 meshGRegionTransfinite.o: meshGRegionTransfinite.cpp meshGFace.h \
   ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/Range.h \
   ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
@@ -227,12 +253,11 @@ meshGRegionExtruded.o: meshGRegionExtruded.cpp ../Geo/ExtrudeParams.h \
   ../Geo/GEdge.h ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h \
   ../Geo/Pair.h ../Geo/ExtrudeParams.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h ../Geo/ExtrudeParams.h ../Geo/SBoundingBox3d.h \
-  ../Common/SmoothNormals.h ../Common/Message.h
+  ../Common/SmoothNormals.h meshGFace.h meshGRegion.h ../Common/Message.h
 DivideAndConquer.o: DivideAndConquer.cpp ../Common/Gmsh.h \
   ../Common/Message.h ../DataStr/Malloc.h ../DataStr/List.h \
   ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h \
-  ../DataStr/Tree.h ../Numeric/Numeric.h DivideAndConquer.h \
-  ../Common/Context.h
+  ../DataStr/Tree.h ../Numeric/Numeric.h DivideAndConquer.h
 BackgroundMesh.o: BackgroundMesh.cpp ../Common/Gmsh.h ../Common/Message.h \
   ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
   ../DataStr/avl.h ../DataStr/Tools.h ../DataStr/List.h ../DataStr/Tree.h \
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index 73e3af890e984198d4cd2503467020ea4caf0362..d9f2bc452bdecba17c0945290efc450b45409551 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: SecondOrder.cpp,v 1.48 2006-11-27 22:22:17 geuzaine Exp $
+// $Id: SecondOrder.cpp,v 1.49 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -417,9 +417,12 @@ void Degre2(bool linear, bool incomplete)
   Msg(STATUS1, "Meshing second order...");
   double t1 = Cpu();
 
+  // first, make sure to remove any existsing second order vertices/elements
+  Degre1();
+
+  // then create new second order vertices/elements
   std::map<std::pair<MVertex*,MVertex*>, MVertex* > edgeVertices;
   std::map<std::vector<MVertex*>, MVertex* > faceVertices;
-
   for(GModel::eiter it = GMODEL->firstEdge(); it != GMODEL->lastEdge(); ++it)
     setSecondOrder(*it, edgeVertices, linear);
   for(GModel::fiter it = GMODEL->firstFace(); it != GMODEL->lastFace(); ++it)
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index c103dd28000499d2bf020e2eb4cee99815653105..af7d7b27ede06198fff3cb7127085ddb5b5caf89 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.26 2007-01-12 13:16:59 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.27 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -25,68 +25,55 @@
 #include "GFace.h"
 #include "MRep.h"
 #include "BackgroundMesh.h"
-#include "Context.h"
 #include "Message.h"
 
-extern Context_T CTX;
-
-static GEdge * _myGEdge;
-static double  _myGEdgeLength, t_begin, t_end, lc_begin, lc_end;
-static Range<double> _myGEdgeBounds;
-
-double F_LC_ANALY (double xx, double yy, double zz)
-{
-  //  return 0.005 + 0.05*fabs (sin(5*xx) + sin(15*yy) + sin(15*zz));
-  //  return 0.02;
-  //  return 0.002 + 0.04*fabs (sin(6*xx) + sin(6*yy) + sin(6*zz));
-  return 0.003 + 0.05*fabs(sin(8*xx) + sin(8*yy) + sin(8*zz));
-  return 0.02 + 0.1*fabs(sin(3*xx) + sin(3*yy) + sin(3*zz));
-  return 0.01 + 0.1*fabs(sin((xx*xx+(zz-0.7)*(zz-0.7)-.25))); 
-  return 0.05 + 0.1*fabs(xx*yy);
-}
-
-double F_Lc(double t)
+double F_Lc(GEdge *ge, double t)
 {
-  GPoint p = _myGEdge -> point (t);
+  GPoint p = ge->point(t);
   double lc_here;
-  if (t == t_begin)
-    lc_here    = BGM_MeshSize(_myGEdge->getBeginVertex(), t , 0 , p.x(),p.y(),p.z());
-  else if (t == t_end)
-    lc_here    = BGM_MeshSize(_myGEdge->getEndVertex(), t , 0 , p.x(),p.y(),p.z());
+
+  Range<double> bounds = ge->parBounds(0);
+  double t_begin = bounds.low();
+  double t_end = bounds.high();
+
+  if(t == t_begin)
+    lc_here = BGM_MeshSize(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
+  else if(t == t_end)
+    lc_here = BGM_MeshSize(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
   else
-    lc_here    = BGM_MeshSize(_myGEdge, t , 0 , p.x(),p.y(),p.z());
-  SVector3 der      = _myGEdge -> firstDer(t) ;
+    lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
 
-  const double d    = norm(der);  
-  return d/lc_here;
+  SVector3 der = ge->firstDer(t);
+  const double d = norm(der);
+  return d / lc_here;
 }
 
-double F_Transfinite(double t)
+double F_Transfinite(GEdge *ge, double t)
 {
   double val, r;
 
-  SVector3 der = _myGEdge->firstDer(t) ;
+  SVector3 der = ge->firstDer(t) ;
   double d = norm(der);
 
-  double coef = _myGEdge->meshAttributes.coeffTransfinite;
-  int type = _myGEdge->meshAttributes.typeTransfinite;
-  int nbpt = _myGEdge->meshAttributes.nbPointsTransfinite;
+  double coef = ge->meshAttributes.coeffTransfinite;
+  int type = ge->meshAttributes.typeTransfinite;
+  int nbpt = ge->meshAttributes.nbPointsTransfinite;
 
   if(coef <= 0.0 || coef == 1.0) {
     // coef < 0 should never happen
-    val = d * coef / _myGEdgeLength;
+    val = d * coef / ge->length();
   }
   else {
     switch (abs(type)) {
 
-    case 1: // Geometric progression ar^i; Sum of n terms = THEC->l = a (r^n-1)/(r-1)
+    case 1: // Geometric progression ar^i; Sum of n terms = length = a (r^n-1)/(r-1)
       {
 	if(sign(type) >= 0)
 	  r = coef;
 	else
 	  r = 1. / coef;
-	double a = _myGEdgeLength * (r - 1.) / (pow(r, nbpt - 1.) - 1.);
-	int i = (int)(log(t * _myGEdgeLength / a * (r - 1.) + 1.) / log(r));
+	double a = ge->length() * (r - 1.) / (pow(r, nbpt - 1.) - 1.);
+	int i = (int)(log(t * ge->length() / a * (r - 1.) + 1.) / log(r));
 	val = d / (a * pow(r, (double)i));
       }
       break;
@@ -97,30 +84,31 @@ double F_Transfinite(double t)
 	if(coef > 1.0) {
 	  a = -4. * sqrt(coef - 1.) *
 	    atan2(1., sqrt(coef - 1.)) /
-	    ((double)nbpt *  _myGEdgeLength);
+	    ((double)nbpt *  ge->length());
 	}
 	else {
 	  a = 2. * sqrt(1. - coef) *
 	    log(fabs((1. + 1. / sqrt(1. - coef))
 		     / (1. - 1. / sqrt(1. - coef))))
-	    / ((double)nbpt * _myGEdgeLength);
+	    / ((double)nbpt * ge->length());
 	}
-	double b = -a * _myGEdgeLength * _myGEdgeLength / (4. * (coef - 1.));
-	val = d / (-a * DSQR(t * _myGEdgeLength - (_myGEdgeLength) * 0.5) + b);
+	double b = -a * ge->length() * ge->length() / (4. * (coef - 1.));
+	val = d / (-a * DSQR(t * ge->length() - (ge->length()) * 0.5) + b);
       }
       break;
       
     default:
       Msg(WARNING, "Unknown case in Transfinite Line mesh");
       val = 1.;
+      break;
     }
   }
   return val;
 }
 
-double F_One(double t)
+double F_One(GEdge *ge, double t)
 {
-  SVector3 der = _myGEdge->firstDer(t) ;
+  SVector3 der = ge->firstDer(t) ;
   return norm(der);
 }
 
@@ -134,24 +122,22 @@ double trapezoidal(IntPoint * P1, IntPoint * P2)
   return (0.5 * (P1->lc + P2->lc) * (P2->t - P1->t));
 }
 
-void RecursiveIntegration(IntPoint * from, IntPoint * to,
-                          double (*f) (double X), List_T * pPoints,
+void RecursiveIntegration(GEdge *ge, IntPoint * from, IntPoint * to,
+                          double (*f) (GEdge *e, double X), List_T * pPoints,
                           double Prec, int *depth)
 {
   IntPoint P, p1;
-  double err, val1, val2, val3;
 
   (*depth)++;
 
   P.t = 0.5 * (from->t + to->t);
-  P.lc = f(P.t);
+  P.lc = f(ge, P.t);
 
-  val1 = trapezoidal(from, to);
-  val2 = trapezoidal(from, &P);
-  val3 = trapezoidal(&P, to);
+  double val1 = trapezoidal(from, to);
+  double val2 = trapezoidal(from, &P);
+  double val3 = trapezoidal(&P, to);
+  double err = fabs(val1 - val2 - val3);
 
-  err = fabs(val1 - val2 - val3);
-  //  Msg(INFO,"Int %22.15 E %22.15 E %22.15 E\n", val1,val2,val3);
   if(((err < Prec) && (*depth > 1)) || (*depth > 25)) {
     List_Read(pPoints, List_Nbr(pPoints) - 1, &p1);
     P.p = p1.p + val2;
@@ -162,43 +148,46 @@ void RecursiveIntegration(IntPoint * from, IntPoint * to,
     List_Add(pPoints, to);
   }
   else {
-    RecursiveIntegration(from, &P, f, pPoints, Prec, depth);
-    RecursiveIntegration(&P, to, f, pPoints, Prec, depth);
+    RecursiveIntegration(ge, from, &P, f, pPoints, Prec, depth);
+    RecursiveIntegration(ge, &P, to, f, pPoints, Prec, depth);
   }
+
   (*depth)--;
 }
 
-double Integration(double t1, double t2, double (*f) (double X),
+double Integration(GEdge *ge, double t1, double t2, 
+		   double (*f) (GEdge *e, double X),
                    List_T * pPoints, double Prec)
 {
-  int depth;
   IntPoint from, to;
 
-  depth = 0;
+  int depth = 0;
 
   from.t = t1;
-  from.lc = f(from.t);
+  from.lc = f(ge, from.t);
   from.p = 0.0;
   List_Add(pPoints, &from);
 
   to.t = t2;
-  to.lc = f(to.t);
-  RecursiveIntegration(&from, &to, f, pPoints, Prec, &depth);
+  to.lc = f(ge, to.t);
+  RecursiveIntegration(ge, &from, &to, f, pPoints, Prec, &depth);
 
   List_Read(pPoints, List_Nbr(pPoints) - 1, &to);
-  return (to.p);
+  return to.p;
 }
 
-void deMeshGEdge :: operator() (GEdge *ge) 
+void deMeshGEdge::operator() (GEdge *ge) 
 {
-  for (unsigned int i=0;i<ge->mesh_vertices.size();i++) delete ge->mesh_vertices[i];
+  for (unsigned int i = 0; i < ge->mesh_vertices.size(); i++) 
+    delete ge->mesh_vertices[i];
   ge->mesh_vertices.clear();
-  for (unsigned int i=0;i<ge->lines.size();i++) delete ge->lines[i];
+  for (unsigned int i = 0; i < ge->lines.size(); i++) 
+    delete ge->lines[i];
   ge->lines.clear();
   if(ge->meshRep) ge->meshRep->destroy();
 }
 
-void meshGEdge :: operator() (GEdge *ge) 
+void meshGEdge::operator() (GEdge *ge) 
 {  
   if(ge->geomType() == GEntity::DiscreteCurve) return;
 
@@ -213,37 +202,26 @@ void meshGEdge :: operator() (GEdge *ge)
   // Create a list of integration points
   List_T *Points = List_Create(10, 10, sizeof(IntPoint));
 
-  // For avoiding the global variable :
-  // We have to change the Integration function in order
-  // to pass an extra argument... 
-  _myGEdge = ge;
-    
   // compute bounds
-  _myGEdgeBounds = ge->parBounds(0) ;
-  t_begin = _myGEdgeBounds.low();
-  t_end = _myGEdgeBounds.high();
+  Range<double> bounds = ge->parBounds(0);
+  double t_begin = bounds.low();
+  double t_end = bounds.high();
   
   // first compute the length of the curve by integrating one
-  _myGEdgeLength = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), 
-			       F_One, Points, 1.e-8);
-  ge->setLength (_myGEdgeLength);
+  double length = Integration(ge, t_begin, t_end, F_One, Points, 1.e-8);
+  ge->setLength(length);
 
   List_Reset(Points);
-
-  lc_begin = _myGEdge->getBeginVertex()->prescribedMeshSizeAtVertex();
-  lc_end = _myGEdge->getEndVertex()->prescribedMeshSizeAtVertex();
     
   // Integrate detJ/lc du 
   double a;
   int N;
   if(ge->meshAttributes.Method == TRANSFINI){
-    a = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), 
-		    F_Transfinite, Points, 1.e-8);
+    a = Integration(ge, t_begin, t_end, F_Transfinite, Points, 1.e-8);
     N = ge->meshAttributes.nbPointsTransfinite;
   }
   else{
-    a = Integration(_myGEdgeBounds.low(), _myGEdgeBounds.high(), 
-		    F_Lc, Points, 1.e-8);
+    a = Integration(ge, t_begin, t_end, F_Lc, Points, 1.e-8);
     N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
   }
   const double b = a / (double)(N - 1);
@@ -251,9 +229,8 @@ void meshGEdge :: operator() (GEdge *ge)
   int count = 1, NUMP = 1;
   IntPoint P1, P2;
 
-  // do not consider the first and the last vertex 
-  // those are not classified on this mesh edge
-
+  // do not consider the first and the last vertex (those are not
+  // classified on this mesh edge)
   if(N > 2){
     ge->mesh_vertices.resize(N - 2);
     while(NUMP < N - 1) {
diff --git a/Mesh/meshGEdge.h b/Mesh/meshGEdge.h
index 48eb8d843be2e1409506872231749b222829607a..aacf097068c64b64fc4a3c0a24fdec34da3f636d 100644
--- a/Mesh/meshGEdge.h
+++ b/Mesh/meshGEdge.h
@@ -21,16 +21,15 @@
 // Please report all bugs and problems to <gmsh@geuz.org>.
 
 class GEdge;
+
 // Create the mesh of the edge
-class meshGEdge 
-{
+class meshGEdge {
  public :
   void operator () ( GEdge * );
 };
 
 // destroy the mesh of the edge
-class deMeshGEdge 
-{
+class deMeshGEdge {
  public :
   void operator () ( GEdge * );
 };
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index b39ad0026fe5eb57cee6419d2bfa57fa4d170e08..f59147fbfbd46f228bccfa9b9ed9a0c0af3c35b5 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.53 2007-01-20 14:06:37 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.54 2007-01-22 16:31:43 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -84,7 +84,16 @@ void computeEdgeLoops(const GFace *gf,
   }
 }
 
-extern double F_LC_ANALY (double xx, double yy, double zz);
+double F_LC_ANALY(double xx, double yy, double zz)
+{
+  //  return 0.005 + 0.05*fabs (sin(5*xx) + sin(15*yy) + sin(15*zz));
+  //  return 0.02;
+  //  return 0.002 + 0.04*fabs (sin(6*xx) + sin(6*yy) + sin(6*zz));
+  return 0.003 + 0.05*fabs(sin(8*xx) + sin(8*yy) + sin(8*zz));
+  return 0.02 + 0.1*fabs(sin(3*xx) + sin(3*yy) + sin(3*zz));
+  return 0.01 + 0.1*fabs(sin((xx*xx+(zz-0.7)*(zz-0.7)-.25))); 
+  return 0.05 + 0.1*fabs(xx*yy);
+}
 
 double NewGetLc(BDS_Point *p)
 {
diff --git a/Mesh/meshGFaceExtruded.cpp b/Mesh/meshGFaceExtruded.cpp
index 91d0fdaeda091a3386c465d37badb35a7d2b4ab1..b04c979ad76c056c0bf6fde40514e86851a49bb3 100644
--- a/Mesh/meshGFaceExtruded.cpp
+++ b/Mesh/meshGFaceExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceExtruded.cpp,v 1.15 2007-01-16 11:31:42 geuzaine Exp $
+// $Id: meshGFaceExtruded.cpp,v 1.16 2007-01-22 16:31:44 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -52,7 +52,7 @@ void createQuaTri(std::vector<MVertex*> &v, GFace *to,
 	to->triangles.push_back(new MTriangle(v[2], v[1], v[0]));
 	to->triangles.push_back(new MTriangle(v[2], v[3], v[1]));
       }
-      else {
+      else{
 	to->triangles.push_back(new MTriangle(v[2], v[3], v[0]));
 	to->triangles.push_back(new MTriangle(v[0], v[3], v[1]));
       }
diff --git a/Mesh/meshGRegionExtruded.cpp b/Mesh/meshGRegionExtruded.cpp
index de516335fbc15e4ad9908c391205ef5c9176d12c..1c30fd029b15612aed3c25e58c1e84bd3be737c7 100644
--- a/Mesh/meshGRegionExtruded.cpp
+++ b/Mesh/meshGRegionExtruded.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionExtruded.cpp,v 1.9 2007-01-16 14:15:18 geuzaine Exp $
+// $Id: meshGRegionExtruded.cpp,v 1.10 2007-01-22 16:31:44 geuzaine Exp $
 //
 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
 //
@@ -288,8 +288,8 @@ void phase2(GRegion *gr,
 	    if(!edgeExists(v[3], v[1], edges_swap)) {
 	      deleteEdge(v[3], v[1], edges);
 	      createEdge(v[0], v[4], edges);
-	      createEdge(v[0], v[4], edges_swap);
 	      createEdge(v[3], v[1], edges_swap);
+	      createEdge(v[0], v[4], edges_swap);
 	    }
 	    else if(!edgeExists(v[4], v[2], edges_swap)) {
 	      deleteEdge(v[4], v[2], edges);
@@ -317,14 +317,14 @@ void phase2(GRegion *gr,
 	    else if(!edgeExists(v[1], v[5], edges_swap)) {
 	      deleteEdge(v[1], v[5], edges);
 	      createEdge(v[4], v[2], edges);
-	      createEdge(v[4], v[2], edges_swap);
 	      createEdge(v[1], v[5], edges_swap);
+	      createEdge(v[4], v[2], edges_swap);
 	    }
 	    else if(!edgeExists(v[3], v[2], edges_swap)) {
 	      deleteEdge(v[3], v[2], edges);
 	      createEdge(v[0], v[5], edges);
-	      createEdge(v[0], v[5], edges_swap);
 	      createEdge(v[3], v[2], edges_swap);
+	      createEdge(v[0], v[5], edges_swap);
 	    }
 	  }
 	}
@@ -474,5 +474,5 @@ int SubdivideExtrudedMesh(GModel *m)
   }
 
   MVertexLessThanLexicographic::tolerance = old_tol;   
-  return 0;
+  return 1;
 }
diff --git a/doc/LICENSE b/doc/LICENSE
index 80a5c5b74f75e9e9937ab12792cbd3f6e7573c5b..b53ee8e63dcfe9e8464eb0bf7be072c0f9757038 100644
--- a/doc/LICENSE
+++ b/doc/LICENSE
@@ -2,13 +2,13 @@ Gmsh is provided under the terms of the GNU General Public License
 (GPL) with the following exception:
 
   The copyright holders of Gmsh give you permission to combine Gmsh
-  with with code included in the standard release of Triangle (written
-  by Jonathan Shewchuk) and TetGen (written by Hang Si) under their
-  respective licenses. You may copy and distribute such a system
-  following the terms of the GNU GPL for Gmsh and the licenses of the
-  other code concerned, provided that you include the source code of
-  that other code when and as the GNU GPL requires distribution of
-  source code.
+  with code included in the standard release of Triangle (written by
+  Jonathan Shewchuk), TetGen (written by Hang Si) and Netgen (written
+  by Joachim Sch"oberl) under their respective licenses. You may copy
+  and distribute such a system following the terms of the GNU GPL for
+  Gmsh and the licenses of the other code concerned, provided that you
+  include the source code of that other code when and as the GNU GPL
+  requires distribution of source code.
 
   Note that people who make modified versions of Gmsh are not
   obligated to grant this special exception for their modified
diff --git a/doc/TODO b/doc/TODO
index 3f3fda17a8a6377a8f0b2e6b942414af0d666db3..0d0d02fc76c63d1d8a8e284c7f0df5acfbbf4a5d 100644
--- a/doc/TODO
+++ b/doc/TODO
@@ -1,4 +1,4 @@
-$Id: TODO,v 1.38 2007-01-19 15:34:05 geuzaine Exp $
+$Id: TODO,v 1.39 2007-01-22 16:31:44 geuzaine Exp $
 
 ********************************************************************
 
@@ -6,12 +6,11 @@ introduce Right/Left/Alternate for extruded meshes
 
 ********************************************************************
 
-color edges depending on number of adjacent surfaces (to check face
-sewing w/ tolerance)
+fix second order mesh for periodic surfaces
 
 ********************************************************************
 
-fix second order mesh for periodic surfaces
+recode full-quad recombine algo
 
 ********************************************************************
 
@@ -20,6 +19,10 @@ cases
 
 ********************************************************************
 
+recode BDS::extract edges
+
+********************************************************************
+
 physical groups->add->line/surface: pressing '-' after/while selecting 
 a surface should add it with reverse orientation?