diff --git a/Common/Options.cpp b/Common/Options.cpp
index 7d2cde1fd82632f6a82b004893c9d921636146ca..a6f2c9a18bececee17a3fd3ddeed8b773f2acf70 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.283 2006-08-06 22:58:47 geuzaine Exp $
+// $Id: Options.cpp,v 1.284 2006-08-07 00:08:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -147,7 +147,7 @@ void Init_Options(int num)
   CTX.mesh.histogram = 0;
   CTX.mesh.changed = 1;
   CTX.mesh.oldxtrude = CTX.mesh.oldxtrude_recombine = 0; // old extrusion mesh generator
-  CTX.mesh.check_duplicates = 0; // check for duplicate nodes in Read_Mesh
+  CTX.mesh.check_duplicates = 0; // check for duplicate nodes when loading a mesh
   CTX.mesh.bgmesh_type = WITHPOINTS;
   CTX.post.combine_time = 0; // try to combine_time views at startup
   CTX.post.plugin_draw_function = NULL;
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 23f6e550ad03384077491ef33968e66749e19f4c..0362a7835749ae48a8d9dc60e56ebda64781f9aa 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.422 2006-08-06 22:58:47 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.423 2006-08-07 00:08:07 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -649,24 +649,12 @@ int _save_lc(char *name)
   return 1;
 }
 
-int _save_gref(char *name)
-{
-  CreateOutputFile(name, CTX.mesh.format = FORMAT_GREF);
-  return 1;
-}
-
 int _save_unv(char *name)
 {
   CreateOutputFile(name, CTX.mesh.format = FORMAT_UNV);
   return 1;
 }
 
-int _save_p3d(char *name)
-{
-  CreateOutputFile(name, CTX.mesh.format = FORMAT_P3D);
-  return 1;
-}
-
 int _save_vrml(char *name)
 {
   CreateOutputFile(name, CTX.mesh.format = FORMAT_VRML);
@@ -789,11 +777,9 @@ void file_save_as_cb(CALLBACK_ARGS)
     {"Gmsh unrolled geometry (*.geo)", _save_geo},
     {"Gmsh mesh (*.msh)", _save_msh},
     {"Gmsh mesh statistics (*.pos)", _save_lc},
-    {"GREF mesh (*.gref)", _save_gref},
     {"I-DEAS universal mesh (*.unv)", _save_unv},
-    {"PLOT3D formatted ASCII mgrid (*.p3d)", _save_p3d},
     {"VRML surface mesh (*.wrl)", _save_vrml},
-    {"STL triangulation (*.stl)", _save_stl},
+    {"STL surface mesh (*.stl)", _save_stl},
     {"GIF (*.gif)", _save_gif},
 #if defined(HAVE_LIBJPEG)
     {"JPEG (*.jpg)", _save_jpeg},
diff --git a/Geo/Makefile b/Geo/Makefile
index 69ab0e8c9a94e6b359d30d467be685f47b1de7e3..078d396bf595c9a7e0c5adbc0e7dedb363e3adc5 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.82 2006-08-06 22:58:49 geuzaine Exp $
+# $Id: Makefile,v 1.83 2006-08-07 00:08:08 geuzaine Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -138,19 +138,20 @@ GVertex.o: GVertex.cpp GVertex.h GEntity.h Range.h SPoint3.h \
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 GEdge.o: GEdge.cpp GEdge.h GEntity.h Range.h SPoint3.h SBoundingBox3d.h \
   MVertex.h ../Common/GmshDefines.h GVertex.h GPoint.h SVector3.h \
-  SPoint2.h MElement.h
+  SPoint2.h MElement.h ../Numeric/Numeric.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 GFace.o: GFace.cpp GFace.h GPoint.h GEntity.h Range.h SPoint3.h \
-  SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h MElement.h SPoint2.h \
-  SVector3.h Pair.h GEdge.h GVertex.h
+  SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h MElement.h \
+  ../Numeric/Numeric.h SPoint2.h SVector3.h Pair.h GEdge.h GVertex.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 GRegion.o: GRegion.cpp GRegion.h GEntity.h Range.h SPoint3.h \
-  SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h MElement.h GFace.h \
-  GPoint.h SPoint2.h SVector3.h Pair.h
+  SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h MElement.h \
+  ../Numeric/Numeric.h GFace.h GPoint.h SPoint2.h SVector3.h Pair.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 GModel.o: GModel.cpp GModel.h GVertex.h GEntity.h Range.h SPoint3.h \
   SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h GPoint.h GEdge.h \
-  SVector3.h SPoint2.h MElement.h GFace.h Pair.h GRegion.h
+  SVector3.h SPoint2.h MElement.h ../Numeric/Numeric.h GFace.h Pair.h \
+  GRegion.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 GModelIO.o: GModelIO.cpp ../Common/Message.h ../Common/GmshDefines.h \
   gmshRegion.h ../Mesh/Mesh.h ../DataStr/List.h ../DataStr/Tree.h \
@@ -167,17 +168,17 @@ GModelIO.o: GModelIO.cpp ../Common/Message.h ../Common/GmshDefines.h \
 MVertex.o: MVertex.cpp MVertex.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
-  GEntity.h Range.h SPoint3.h SBoundingBox3d.h
+  ../Numeric/Numeric.h GEntity.h Range.h SPoint3.h SBoundingBox3d.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 gmshModel.o: gmshModel.cpp gmshModel.h GModel.h GVertex.h GEntity.h \
   Range.h SPoint3.h SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h \
-  GPoint.h GEdge.h SVector3.h SPoint2.h MElement.h GFace.h Pair.h \
-  GRegion.h ../Mesh/Mesh.h ../DataStr/List.h ../DataStr/Tree.h \
-  ../DataStr/avl.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Vertex.h \
-  ../Mesh/Simplex.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Face.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Edge.h ../Mesh/Vertex.h \
-  ../Mesh/Simplex.h ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
-  ../Common/SmoothNormals.h ../Numeric/Numeric.h ../Mesh/Metric.h \
+  GPoint.h GEdge.h SVector3.h SPoint2.h MElement.h ../Numeric/Numeric.h \
+  GFace.h Pair.h GRegion.h ../Mesh/Mesh.h ../DataStr/List.h \
+  ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Edge.h \
+  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Geo/ExtrudeParams.h \
+  ../Common/VertexArray.h ../Common/SmoothNormals.h ../Mesh/Metric.h \
   ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h \
   Geo.h ../Parser/OpenFile.h ../DataStr/Tools.h ../DataStr/List.h \
   ../DataStr/Tree.h ../Common/Message.h gmshVertex.h gmshFace.h \
@@ -185,43 +186,45 @@ gmshModel.o: gmshModel.cpp gmshModel.h GModel.h GVertex.h GEntity.h \
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 gmshEdge.o: gmshEdge.cpp gmshModel.h GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h GPoint.h \
-  GEdge.h SVector3.h SPoint2.h MElement.h GFace.h Pair.h GRegion.h \
-  gmshEdge.h gmshVertex.h ../Mesh/Mesh.h ../DataStr/List.h \
-  ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Edge.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Geo/ExtrudeParams.h \
-  ../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
-  ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
-  ../Mesh/Matrix.h ../Mesh/Interpolation.h ../Mesh/Vertex.h \
-  ../Mesh/Mesh.h CAD.h ExtrudeParams.h Geo.h ../Common/Context.h
+  GEdge.h SVector3.h SPoint2.h MElement.h ../Numeric/Numeric.h GFace.h \
+  Pair.h GRegion.h gmshEdge.h gmshVertex.h ../Mesh/Mesh.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
+  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Vertex.h \
+  ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h \
+  ../Mesh/Interpolation.h ../Mesh/Vertex.h ../Mesh/Mesh.h CAD.h \
+  ExtrudeParams.h Geo.h ../Common/Context.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 gmshFace.o: gmshFace.cpp gmshModel.h GModel.h GVertex.h GEntity.h Range.h \
   SPoint3.h SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h GPoint.h \
-  GEdge.h SVector3.h SPoint2.h MElement.h GFace.h Pair.h GRegion.h \
-  gmshEdge.h gmshVertex.h ../Mesh/Mesh.h ../DataStr/List.h \
-  ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Edge.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Geo/ExtrudeParams.h \
-  ../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
-  ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
-  ../Mesh/Matrix.h gmshFace.h ../Mesh/Interpolation.h ../Mesh/Vertex.h \
-  ../Mesh/Mesh.h CAD.h ExtrudeParams.h Geo.h ../Mesh/Utils.h \
-  ../Mesh/Vertex.h ../Mesh/Mesh.h
+  GEdge.h SVector3.h SPoint2.h MElement.h ../Numeric/Numeric.h GFace.h \
+  Pair.h GRegion.h gmshEdge.h gmshVertex.h ../Mesh/Mesh.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
+  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Vertex.h \
+  ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h gmshFace.h \
+  ../Mesh/Interpolation.h ../Mesh/Vertex.h ../Mesh/Mesh.h CAD.h \
+  ExtrudeParams.h Geo.h ../Mesh/Utils.h ../Mesh/Vertex.h ../Mesh/Mesh.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 gmshRegion.o: gmshRegion.cpp gmshModel.h GModel.h GVertex.h GEntity.h \
   Range.h SPoint3.h SBoundingBox3d.h MVertex.h ../Common/GmshDefines.h \
-  GPoint.h GEdge.h SVector3.h SPoint2.h MElement.h GFace.h Pair.h \
-  GRegion.h gmshEdge.h gmshVertex.h ../Mesh/Mesh.h ../DataStr/List.h \
-  ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Edge.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Geo/ExtrudeParams.h \
-  ../Common/VertexArray.h ../Common/SmoothNormals.h ../Numeric/Numeric.h \
-  ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
-  ../Mesh/Matrix.h gmshFace.h gmshRegion.h ../Mesh/Interpolation.h \
-  ../Mesh/Vertex.h ../Mesh/Mesh.h CAD.h ExtrudeParams.h Geo.h
+  GPoint.h GEdge.h SVector3.h SPoint2.h MElement.h ../Numeric/Numeric.h \
+  GFace.h Pair.h GRegion.h gmshEdge.h gmshVertex.h ../Mesh/Mesh.h \
+  ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h \
+  ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h \
+  ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
+  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Mesh/Metric.h ../Mesh/Vertex.h \
+  ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h gmshFace.h \
+  gmshRegion.h ../Mesh/Interpolation.h ../Mesh/Vertex.h ../Mesh/Mesh.h \
+  CAD.h ExtrudeParams.h Geo.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
 SVector3.o: SVector3.cpp SVector3.h SPoint3.h
 # 1 "/Users/geuzaine/.gmsh/Geo//"
diff --git a/Graphics/CreateFile.cpp b/Graphics/CreateFile.cpp
index 21aa08bac5973f3192d96feff293c73efdd579b1..b541d302c3dbbd9404c8e4c25ca682c2d8e1eb5e 100644
--- a/Graphics/CreateFile.cpp
+++ b/Graphics/CreateFile.cpp
@@ -1,4 +1,4 @@
-// $Id: CreateFile.cpp,v 1.86 2006-08-06 22:58:49 geuzaine Exp $
+// $Id: CreateFile.cpp,v 1.87 2006-08-07 00:08:08 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -55,8 +55,6 @@ int GuessFileFormatFromFileName(char *name)
   else if(!strcmp(ext, ".opt"))     return FORMAT_OPT;
   else if(!strcmp(ext, ".msh"))     return FORMAT_MSH;
   else if(!strcmp(ext, ".unv"))     return FORMAT_UNV;
-  else if(!strcmp(ext, ".p3d"))     return FORMAT_P3D;
-  else if(!strcmp(ext, ".dmg"))     return FORMAT_DMG;
   else if(!strcmp(ext, ".stl"))     return FORMAT_STL;
   else if(!strcmp(ext, ".pos"))     return FORMAT_LC;
   else if(!strcmp(ext, ".gif"))     return FORMAT_GIF;
@@ -73,32 +71,10 @@ int GuessFileFormatFromFileName(char *name)
   else if(!strcmp(ext, ".svg"))     return FORMAT_SVG;
   else if(!strcmp(ext, ".ppm"))     return FORMAT_PPM;
   else if(!strcmp(ext, ".yuv"))     return FORMAT_YUV;
-  else if(!strcmp(ext, ".gref"))    return FORMAT_GREF;
-  else if(!strcmp(ext, ".Gref"))    return FORMAT_GREF;
   else if(!strcmp(ext, ".wrl"))     return FORMAT_VRML;
   else                              return -1;
 }
 
-char *GetStringForFileFormat(int format)
-{
-  switch(format){
-  case FORMAT_PPM: return "PPM";
-  case FORMAT_YUV: return "YUV";
-  case FORMAT_GIF: return "GIF";
-  case FORMAT_JPEG: return "JPEG";
-  case FORMAT_JPEGTEX: return "JPEG";
-  case FORMAT_PNG: return "PNG";
-  case FORMAT_PNGTEX: return "PNG";
-  case FORMAT_PS: return "PS";
-  case FORMAT_EPS: return "EPS";
-  case FORMAT_EPSTEX: return "EPS";
-  case FORMAT_PDF: return "PDF";
-  case FORMAT_PDFTEX: return "PDF";
-  case FORMAT_SVG: return "SVG";
-  default: return "";
-  }
-}
-
 void GetDefaultFileName(int format, char *name)
 {
   char ext[10] = "";
@@ -107,10 +83,7 @@ void GetDefaultFileName(int format, char *name)
   case FORMAT_MSH:  strcpy(ext, ".msh"); break;
   case FORMAT_VRML: strcpy(ext, ".wrl"); break;
   case FORMAT_UNV:  strcpy(ext, ".unv"); break;
-  case FORMAT_GREF: strcpy(ext, ".Gref"); break;
-  case FORMAT_DMG:  strcpy(ext, ".dmg"); break;
   case FORMAT_STL:  strcpy(ext, ".stl"); break;
-  case FORMAT_P3D:  strcpy(ext, ".p3d"); break;
   default: break;
   }
   strcat(name, ext);
@@ -134,8 +107,7 @@ void CreateOutputFile(char *filename, int format)
   GLint height = viewport[3] - viewport[1];
 
   bool printEndMessage = true;
-  if(format != FORMAT_AUTO)
-    Msg(INFO, "Writing %s file '%s'", GetStringForFileFormat(format), name);
+  if(format != FORMAT_AUTO) Msg(INFO, "Writing '%s'", name);
 
   switch (format) {
 
@@ -169,12 +141,6 @@ void CreateOutputFile(char *filename, int format)
     GMODEL->writeUNV(name, CTX.mesh.scaling_factor);
     break;
 
-  case FORMAT_P3D:
-  case FORMAT_DMG:
-  case FORMAT_GREF:
-    Print_Mesh(name, format);
-    break;
-
   case FORMAT_LC:
     GMODEL->writePOS(name);
     break;
@@ -353,10 +319,7 @@ void CreateOutputFile(char *filename, int format)
     break;
   }
 
-  if(printEndMessage){
-    Msg(INFO, "Wrote %s file '%s'", GetStringForFileFormat(format), name);
-    Msg(STATUS2N, "Wrote '%s'", name);
-  }
+  if(printEndMessage) Msg(STATUS2, "Wrote '%s'", name);
 
   CTX.print.format = oldformat;
   Draw();
diff --git a/Graphics/Makefile b/Graphics/Makefile
index 26a8c430c18112229bc21d1b8d699bbf50c1166d..7b6d0f7b9bc1f68d820a66abafaa24ac1f800158 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.86 2006-08-04 14:28:02 geuzaine Exp $
+# $Id: Makefile,v 1.87 2006-08-07 00:08:08 geuzaine Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -100,11 +100,11 @@ Mesh.o: Mesh.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
   ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/GPoint.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/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h ../Geo/MElement.h \
-  ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h \
-  ../Geo/GEntity.h ../Geo/MElement.h Draw.h ../Common/Views.h \
-  ../Common/ColorTable.h ../Common/VertexArray.h \
-  ../Common/SmoothNormals.h ../Numeric/Numeric.h ../Common/GmshMatrix.h \
+  ../Numeric/Numeric.h ../Geo/GFace.h ../Geo/GPoint.h ../Geo/GEntity.h \
+  ../Geo/MElement.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+  ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/MElement.h Draw.h \
+  ../Common/Views.h ../Common/ColorTable.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h ../Common/GmshMatrix.h \
   ../Common/AdaptiveViews.h ../Common/GmshMatrix.h ../Common/Context.h \
   gl2ps.h
 # 1 "/Users/geuzaine/.gmsh/Graphics//"
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 6db452b814d09eed6c2a2256eb10dad2eedc1f15..5a859f62d25fd313e60b3781a6d0740693e95481 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.112 2006-08-04 14:28:03 geuzaine Exp $
+# $Id: Makefile,v 1.113 2006-08-07 00:08:08 geuzaine Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -60,8 +60,6 @@ SRC = 1D_Mesh.cpp \
       MeshQuality.cpp \
       Create.cpp \
       Generator.cpp \
-      Print_Mesh.cpp \
-      Read_Mesh.cpp \
       DiscreteSurface.cpp \
       SwapEdge.cpp \
       Utils.cpp \
@@ -394,31 +392,6 @@ Generator.o: Generator.cpp BDS.h ../Common/Views.h ../Common/ColorTable.h \
   ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
   ../Geo/MElement.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
-Print_Mesh.o: Print_Mesh.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 ../Geo/Geo.h ../Geo/CAD.h ../Mesh/Mesh.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
-  ../Mesh/Vertex.h ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
-  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
-  ../Common/SmoothNormals.h ../Common/GmshDefines.h ../Mesh/Metric.h \
-  ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h ../Mesh/Matrix.h \
-  ../Geo/ExtrudeParams.h Mesh.h Create.h Vertex.h ../Common/Context.h
-# 1 "/Users/geuzaine/.gmsh/Mesh//"
-Read_Mesh.o: Read_Mesh.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 \
-  ../Geo/Geo.h ../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Vertex.h \
-  ../Mesh/Element.h ../Mesh/Face.h ../Mesh/Vertex.h ../Mesh/Element.h \
-  ../Mesh/Edge.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
-  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
-  ../Common/SmoothNormals.h ../Numeric/Numeric.h ../Common/GmshDefines.h \
-  ../Mesh/Metric.h ../Mesh/Vertex.h ../Mesh/Simplex.h ../Mesh/Mesh.h \
-  ../Mesh/Matrix.h ../Geo/ExtrudeParams.h Mesh.h 3D_Mesh.h Create.h \
-  Vertex.h ../Common/Context.h PartitionMesh.h
-# 1 "/Users/geuzaine/.gmsh/Mesh//"
 DiscreteSurface.o: DiscreteSurface.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 \
@@ -470,13 +443,12 @@ meshGEdge.o: meshGEdge.cpp meshGEdge.h ../Geo/GEdge.h ../Geo/GEntity.h \
   ../Geo/SPoint3.h ../Geo/MVertex.h ../Common/GmshDefines.h \
   ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/MVertex.h ../Geo/GPoint.h \
   ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
-  ../Geo/MElement.h ../Geo/MVertex.h ../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 \
-  Utils.h Vertex.h Mesh.h Element.h Simplex.h Face.h Edge.h \
-  ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
-  ../Common/SmoothNormals.h ../Numeric/Numeric.h Metric.h Matrix.h \
-  ../Common/Context.h
+  ../Geo/MElement.h ../Geo/MVertex.h ../Numeric/Numeric.h \
+  ../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 Utils.h Vertex.h Mesh.h Element.h \
+  Simplex.h Face.h Edge.h ../Geo/ExtrudeParams.h ../Common/VertexArray.h \
+  ../Common/SmoothNormals.h Metric.h Matrix.h ../Common/Context.h
 # 1 "/Users/geuzaine/.gmsh/Mesh//"
 meshGFace.o: meshGFace.cpp meshGFace.h ../Geo/SPoint2.h ../Geo/GVertex.h \
   ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
diff --git a/Mesh/Mesh.h b/Mesh/Mesh.h
index 58b039877b50918acbffdcaea8abb49c0d4eeb51..4ae621d1e374400a6da4c711e15c854f5c255e89 100644
--- a/Mesh/Mesh.h
+++ b/Mesh/Mesh.h
@@ -362,8 +362,6 @@ void mai3d(int Asked);
 void Init_Mesh0();
 void Init_Mesh();
 void Print_Geo(char *c);
-void Print_Mesh(char *c, int Type);
-void Read_Mesh(Mesh *M, FILE *fp, char *filename, int Type);
 void GetStatistics(double s[50]);
 void GetDefaultMeshFileName(int Type, char *name);
 
diff --git a/Mesh/Print_Mesh.cpp b/Mesh/Print_Mesh.cpp
deleted file mode 100644
index 00872cffe2dd90ec9bc8a357b2096d527a28ead3..0000000000000000000000000000000000000000
--- a/Mesh/Print_Mesh.cpp
+++ /dev/null
@@ -1,1947 +0,0 @@
-// $Id: Print_Mesh.cpp,v 1.77 2006-08-05 10:05:45 geuzaine Exp $
-//
-// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-// 
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-#include "Gmsh.h"
-#include "Numeric.h"
-#include "Geo.h"
-#include "CAD.h"
-#include "Mesh.h"
-#include "Create.h"
-#include "Context.h"
-#include <list>
-#include <map>
-
-extern Context_T CTX;
-extern Mesh *THEM;
-
-// Write mesh in native MSH format
-
-#define LINE            1
-#define TRIANGLE        2
-#define QUADRANGLE      3
-#define TETRAHEDRON     4
-#define HEXAHEDRON      5
-#define PRISM           6
-#define PYRAMID         7
-#define LINE_2          8
-#define TRIANGLE_2      9
-#define QUADRANGLE_2   10
-#define TETRAHEDRON_2  11
-#define HEXAHEDRON_2   12
-#define PRISM_2        13
-#define PYRAMID_2      14
-#define POINT          15
-
-static int MSH_VOL_NUM, MSH_SUR_NUM, MSH_LIN_NUM;
-static int MSH_NODE_NUM, MSH_ELEMENT_NUM, MSH_ADD;
-static int MSH_PHYSICAL_NUM, MSH_PHYSICAL_ORI;
-static FILE *MSHFILE;
-
-static void _msh_print_node(void *a, void *b)
-{
-  Vertex *V = *(Vertex **) a;
-
-  MSH_NODE_NUM++;
-  if(CTX.mesh.renumber_nodes_continuous)
-    V->Num = MSH_NODE_NUM;
-
-  fprintf(MSHFILE, "%d %.16g %.16g %.16g\n",
-          V->Num,
-          V->Pos.X * CTX.mesh.scaling_factor,
-          V->Pos.Y * CTX.mesh.scaling_factor,
-          V->Pos.Z * CTX.mesh.scaling_factor);
-}
-
-static void _msh_process_nodes()
-{
-  int i, j, Num;
-  PhysicalGroup *p;
-  Vertex *pv, **ppv, v;
-
-  for(i = 0; i < List_Nbr(THEM->PhysicalGroups); i++) {
-    List_Read(THEM->PhysicalGroups, i, &p);
-    if(p->Typ == MSH_PHYSICAL_POINT) {
-      for(j = 0; j < List_Nbr(p->Entities); j++) {
-        List_Read(p->Entities, j, &Num);
-        pv = &v;
-        pv->Num = abs(Num);
-        if(!Tree_Search(THEM->Vertices, &pv)) {
-          if((ppv = (Vertex **) Tree_PQuery(THEM->Points, &pv)))
-            Tree_Add(THEM->Vertices, ppv);
-        }
-      }
-    }
-  }
-
-  if(CTX.mesh.msh_file_version == 2.0)
-    fprintf(MSHFILE, "$Nodes\n");
-  else
-    fprintf(MSHFILE, "$NOD\n");
-  fprintf(MSHFILE, "%d\n", Tree_Nbr(THEM->Vertices));
-
-  MSH_NODE_NUM = 0;
-  Tree_Action(THEM->Vertices, _msh_print_node);
-  if(CTX.mesh.msh_file_version == 2.0)
-    fprintf(MSHFILE, "$EndNodes\n");
-  else
-    fprintf(MSHFILE, "$ENDNOD\n");
-}
-
-int _get_partition(int index)
-{
-  MeshPartition **part = (MeshPartition**)List_Pointer_Test(THEM->Partitions, index);
-  if(!part)
-    return -1; // OK, no partitions available
-  else
-    return (*part)->Num; // partition number
-}
-
-static void _msh_print_simplex(void *a, void *b)
-{
-  int i, type, nbn, nbs = 0;
-
-  SimplexBase *s = *(SimplexBase **) a;
-
-  if(MSH_VOL_NUM && (MSH_VOL_NUM != s->iEnt))
-    return;
-
-  if(MSH_SUR_NUM && (MSH_SUR_NUM != s->iEnt))
-    return;
-
-  if(MSH_LIN_NUM && (MSH_LIN_NUM != s->iEnt))
-    return;
-
-  if(!MSH_ADD) {
-    MSH_ELEMENT_NUM++;
-    return;
-  }
-
-  if(!s->V[2]) {
-    nbn = 2;
-    if(s->VSUP) {
-      type = LINE_2;
-      nbs = 1;
-    }
-    else
-      type = LINE;
-  }
-  else if(!s->V[3]) {
-    nbn = 3;
-    if(s->VSUP) {
-      type = TRIANGLE_2;
-      nbs = 3;
-    }
-    else
-      type = TRIANGLE;
-  }
-  else {
-    nbn = 4;
-    if(s->VSUP) {
-      type = TETRAHEDRON_2;
-      nbs = 6;
-      if(s->Volume_Simplexe() < 0) {
-	Vertex *temp;
-	temp = s->V[0];	s->V[0] = s->V[1]; s->V[1] = temp;
-	temp = s->VSUP[1]; s->VSUP[1] = s->VSUP[2]; s->VSUP[2] = temp;
-	temp = s->VSUP[5]; s->VSUP[5] = s->VSUP[3]; s->VSUP[3] = temp;
-      }
-    }
-    else{
-      type = TETRAHEDRON;
-      if(s->Volume_Simplexe() < 0) {
-	Vertex *temp;
-	temp = s->V[0]; s->V[0] = s->V[1]; s->V[1] = temp;
-      }
-    }
-  }
-
-  if(CTX.mesh.msh_file_version == 2.0){
-    int part = _get_partition(s->iPart);
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : s->iEnt, s->iEnt);
-    if(part >= 0) fprintf(MSHFILE, " %d", part);
-  }
-  else
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type,
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : s->iEnt, s->iEnt,
-	    nbn + nbs);
-
-  if(MSH_PHYSICAL_ORI > 0) {
-    for(i = 0; i < nbn; i++)
-      fprintf(MSHFILE, " %d", s->V[i]->Num);
-    for(i = 0; i < nbs; i++)
-      fprintf(MSHFILE, " %d", s->VSUP[i]->Num);
-  }
-  else {
-    for(i = 0; i < nbn; i++)
-      fprintf(MSHFILE, " %d", s->V[nbn - i - 1]->Num);
-    for(i = 0; i < nbs; i++)
-      fprintf(MSHFILE, " %d", s->VSUP[nbs - i - 1]->Num);
-  }
-
-  fprintf(MSHFILE, "\n");
-}
-
-static void _msh_print_quadrangle(void *a, void *b)
-{
-  int i, type, nbn, nbs = 0;
-
-  Quadrangle *q = *(Quadrangle **) a;
-
-  if(MSH_SUR_NUM && (MSH_SUR_NUM != q->iEnt))
-    return;
-
-  if(!MSH_ADD) {
-    MSH_ELEMENT_NUM++;
-    return;
-  }
-
-  nbn = 4;
-  if(q->VSUP) {
-    type = QUADRANGLE_2;
-    nbs = 4 + 1;
-  }
-  else
-    type = QUADRANGLE;
-
-  if(CTX.mesh.msh_file_version == 2.0){
-    int part = _get_partition(q->iPart);
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : q->iEnt, q->iEnt);
-    if(part >= 0) fprintf(MSHFILE, " %d", part);
-  }
-  else
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type,
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : q->iEnt, q->iEnt,
-	    nbn + nbs);
-
-  if(MSH_PHYSICAL_ORI > 0) {
-    for(i = 0; i < nbn; i++)
-      fprintf(MSHFILE, " %d", q->V[i]->Num);
-    for(i = 0; i < nbs; i++)
-      fprintf(MSHFILE, " %d", q->VSUP[i]->Num);
-  }
-  else {
-    for(i = 0; i < nbn; i++)
-      fprintf(MSHFILE, " %d", q->V[nbn - i - 1]->Num);
-    if(nbs){
-      for(i = 0; i < nbs - 1; i++)
-	fprintf(MSHFILE, " %d", q->VSUP[nbs - i - 2]->Num);
-      fprintf(MSHFILE, " %d", q->VSUP[nbs - 1]->Num);
-    }
-  }
-
-  fprintf(MSHFILE, "\n");
-}
-
-static void _msh_print_hexahedron(void *a, void *b)
-{
-  int i, type, nbn, nbs = 0;
-
-  Hexahedron *h = *(Hexahedron **) a;
-
-  if(MSH_VOL_NUM && (MSH_VOL_NUM != h->iEnt))
-    return;
-
-  if(!MSH_ADD) {
-    MSH_ELEMENT_NUM++;
-    return;
-  }
-
-  nbn = 8;
-  if(h->VSUP) {
-    type = HEXAHEDRON_2;
-    nbs = 12 + 6 + 1;
-    if(h->Orientation() < 0) {
-      Vertex *temp;
-      temp = h->V[0]; h->V[0] = h->V[2]; h->V[2] = temp;
-      temp = h->V[4]; h->V[4] = h->V[6]; h->V[6] = temp;
-      Vertex *old[12];
-      for(int i = 0; i < 12; i++) old[i] = h->VSUP[i];
-      h->VSUP[0] = old[3]; h->VSUP[1] = old[5]; h->VSUP[2] = old[6];
-      h->VSUP[3] = old[0]; h->VSUP[4] = old[4]; h->VSUP[5] = old[1];
-      h->VSUP[6] = old[2]; h->VSUP[7] = old[7]; h->VSUP[8] = old[10];
-      h->VSUP[9] = old[11]; h->VSUP[10] = old[8]; h->VSUP[11] = old[9];
-    }
-  }
-  else {
-    type = HEXAHEDRON;
-    if(h->Orientation() < 0) {
-      Vertex *temp;
-      temp = h->V[0]; h->V[0] = h->V[2]; h->V[2] = temp;
-      temp = h->V[4]; h->V[4] = h->V[6]; h->V[6] = temp;
-    }
-  }
-
-  if(CTX.mesh.msh_file_version == 2.0){
-    int part = _get_partition(h->iPart);
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : h->iEnt, h->iEnt);
-    if(part >= 0) fprintf(MSHFILE, " %d", part);
-  }
-  else
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type,
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : h->iEnt, h->iEnt,
-	    nbn + nbs);
-
-  for(i = 0; i < nbn; i++)
-    fprintf(MSHFILE, " %d", h->V[i]->Num);
-  for(i = 0; i < nbs; i++)
-    fprintf(MSHFILE, " %d", h->VSUP[i]->Num);
-
-  fprintf(MSHFILE, "\n");
-}
-
-static void _msh_print_prism(void *a, void *b)
-{
-  int i, type, nbn, nbs = 0;
-
-  Prism *p = *(Prism **) a;
-
-  if(MSH_VOL_NUM && (MSH_VOL_NUM != p->iEnt))
-    return;
-
-  if(!MSH_ADD) {
-    MSH_ELEMENT_NUM++;
-    return;
-  }
-
-  nbn = 6;
-  if(p->VSUP) {
-    type = PRISM_2;
-    nbs = 9 + 3;
-    if(p->Orientation() < 0) {
-      Vertex *temp;
-      temp = p->V[0]; p->V[0] = p->V[1]; p->V[1] = temp;
-      temp = p->V[3]; p->V[3] = p->V[4]; p->V[4] = temp;
-      temp = p->VSUP[1]; p->VSUP[1] = p->VSUP[3]; p->VSUP[3] = temp;
-      temp = p->VSUP[2]; p->VSUP[2] = p->VSUP[4]; p->VSUP[4] = temp;
-      temp = p->VSUP[7]; p->VSUP[7] = p->VSUP[8]; p->VSUP[8] = temp;
-    }
-  }
-  else {
-    type = PRISM;
-    if(p->Orientation() < 0) {
-      Vertex *temp;
-      temp = p->V[0]; p->V[0] = p->V[1]; p->V[1] = temp;
-      temp = p->V[3]; p->V[3] = p->V[4]; p->V[4] = temp;
-    }
-  }
-
-  if(CTX.mesh.msh_file_version == 2.0){
-    int part = _get_partition(p->iPart);
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt);
-    if(part >= 0) fprintf(MSHFILE, " %d", part);
-  }
-  else
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type,
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt,
-	    nbn + nbs);
-
-  for(i = 0; i < nbn; i++)
-    fprintf(MSHFILE, " %d", p->V[i]->Num);
-  for(i = 0; i < nbs; i++)
-    fprintf(MSHFILE, " %d", p->VSUP[i]->Num);
-
-  fprintf(MSHFILE, "\n");
-}
-
-static void _msh_print_pyramid(void *a, void *b)
-{
-  int i, type, nbn, nbs = 0;
-
-  Pyramid *p = *(Pyramid **) a;
-
-  if(MSH_VOL_NUM && (MSH_VOL_NUM != p->iEnt))
-    return;
-
-  if(!MSH_ADD) {
-    MSH_ELEMENT_NUM++;
-    return;
-  }
-
-  nbn = 5;
-  if(p->VSUP) {
-    type = PYRAMID_2;
-    nbs = 8 + 1;
-    if(p->Orientation() < 0) {
-      Vertex *temp;
-      temp = p->V[0]; p->V[0] = p->V[2]; p->V[2] = temp;
-      temp = p->VSUP[0]; p->VSUP[0] = p->VSUP[3]; p->VSUP[3] = temp;
-      temp = p->VSUP[1]; p->VSUP[1] = p->VSUP[5]; p->VSUP[5] = temp;
-      temp = p->VSUP[2]; p->VSUP[2] = p->VSUP[6]; p->VSUP[6] = temp;
-    }
-  }
-  else {
-    type = PYRAMID;
-    if(p->Orientation() < 0) {
-      Vertex *temp;
-      temp = p->V[0]; p->V[0] = p->V[2]; p->V[2] = temp;
-    }
-  }
-
-  if(CTX.mesh.msh_file_version == 2.0){
-    int part = _get_partition(p->iPart);
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type, (part >= 0) ? 3 : 2, 
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt);
-    if(part >= 0) fprintf(MSHFILE, " %d", part);
-  }
-  else
-    fprintf(MSHFILE, "%d %d %d %d %d",
-	    MSH_ELEMENT_NUM++, type,
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : p->iEnt, p->iEnt,
-	    nbn + nbs);
-
-  for(i = 0; i < nbn; i++)
-    fprintf(MSHFILE, " %d", p->V[i]->Num);
-  for(i = 0; i < nbs; i++)
-    fprintf(MSHFILE, " %d", p->VSUP[i]->Num);
-
-  fprintf(MSHFILE, "\n");
-}
-
-static void _msh_print_point(Vertex *V)
-{
-  if(!MSH_ADD) {
-    MSH_ELEMENT_NUM++;
-    return;
-  }
-
-  if(CTX.mesh.msh_file_version == 2.0)
-    fprintf(MSHFILE, "%d %d 2 %d %d %d\n",
-	    MSH_ELEMENT_NUM++, POINT, MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : V->Num, V->Num,
-	    V->Num);
-  else
-    fprintf(MSHFILE, "%d %d %d %d 1 %d\n",
-	    MSH_ELEMENT_NUM++, POINT,
-	    MSH_PHYSICAL_NUM ? MSH_PHYSICAL_NUM : V->Num, V->Num, V->Num);
-}
-
-static void _msh_print_elements()
-{
-  int i, j, k, Num;
-
-  PhysicalGroup *p;
-  Volume *pV;
-  Surface *ps;
-  Curve *pc;
-  Vertex *pv, v;
-
-  List_T *ListCurves = Tree2List(THEM->Curves);
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-  List_T *ListVolumes = Tree2List(THEM->Volumes);
-
-  for(i = 0; i < List_Nbr(THEM->PhysicalGroups); i++) {
-    List_Read(THEM->PhysicalGroups, i, &p);
-    MSH_PHYSICAL_NUM = p->Num;
-    MSH_VOL_NUM = MSH_SUR_NUM = MSH_LIN_NUM = 0;
-
-    switch (p->Typ) {
-
-    case MSH_PHYSICAL_POINT:
-      for(j = 0; j < List_Nbr(p->Entities); j++) {
-        pv = &v;
-        List_Read(p->Entities, j, &Num);
-        pv->Num = abs(Num);
-        MSH_PHYSICAL_ORI = sign(Num);
-        if(Tree_Query(THEM->Vertices, &pv))
-          _msh_print_point(pv);
-      }
-      break;
-
-    case MSH_PHYSICAL_LINE:
-      if(CTX.mesh.oldxtrude) {  //for old extrusion mesh generator
-        for(k = 0; k < List_Nbr(ListVolumes); k++) {
-          List_Read(ListVolumes, k, &pV);
-          for(j = 0; j < List_Nbr(p->Entities); j++) {
-            List_Read(p->Entities, j, &Num);
-            MSH_LIN_NUM = abs(Num);
-            MSH_PHYSICAL_ORI = sign(Num);
-            Tree_Action(pV->Lin_Surf, _msh_print_simplex);
-          }
-        }
-      }
-      else{
-	for(k = 0; k < List_Nbr(ListCurves); k++) {
-	  List_Read(ListCurves, k, &pc);
-	  for(j = 0; j < List_Nbr(p->Entities); j++) {
-	    List_Read(p->Entities, j, &Num);
-	    MSH_LIN_NUM = abs(Num);
-	    MSH_PHYSICAL_ORI = sign(Num);
-	    Tree_Action(pc->Simplexes, _msh_print_simplex);
-	    Tree_Action(pc->SimplexesBase, _msh_print_simplex);
-	  }
-	}
-      }
-      break;
-
-    case MSH_PHYSICAL_SURFACE:
-      if(CTX.mesh.oldxtrude) {  //for old extrusion mesh generator
-        for(k = 0; k < List_Nbr(ListVolumes); k++) {
-          List_Read(ListVolumes, k, &pV);
-          for(j = 0; j < List_Nbr(p->Entities); j++) {
-            List_Read(p->Entities, j, &Num);
-            MSH_SUR_NUM = abs(Num);
-            MSH_PHYSICAL_ORI = sign(Num);
-            Tree_Action(pV->Simp_Surf, _msh_print_simplex);
-            Tree_Action(pV->Quad_Surf, _msh_print_quadrangle);
-          }
-        }
-      }
-      else{
-	for(k = 0; k < List_Nbr(ListSurfaces); k++) {
-	  List_Read(ListSurfaces, k, &ps);
-	  for(j = 0; j < List_Nbr(p->Entities); j++) {
-	    List_Read(p->Entities, j, &Num);
-	    MSH_SUR_NUM = abs(Num);
-	    MSH_PHYSICAL_ORI = sign(Num);
-	    Tree_Action(ps->Simplexes, _msh_print_simplex);
-	    Tree_Action(ps->SimplexesBase, _msh_print_simplex);
-	    Tree_Action(ps->Quadrangles, _msh_print_quadrangle);
-	  }
-	}
-      }
-      break;
-
-    case MSH_PHYSICAL_VOLUME:
-      for(k = 0; k < List_Nbr(ListVolumes); k++) {
-        List_Read(ListVolumes, k, &pV);
-        for(j = 0; j < List_Nbr(p->Entities); j++) {
-          List_Read(p->Entities, j, &Num);
-          MSH_VOL_NUM = abs(Num);
-          MSH_PHYSICAL_ORI = sign(Num);
-          Tree_Action(pV->Simplexes, _msh_print_simplex);
-          Tree_Action(pV->SimplexesBase, _msh_print_simplex);
-          Tree_Action(pV->Hexahedra, _msh_print_hexahedron);
-          Tree_Action(pV->Prisms, _msh_print_prism);
-          Tree_Action(pV->Pyramids, _msh_print_pyramid);
-        }
-      }
-      break;
-
-    default:
-      Msg(GERROR, "Unknown type of physical group");
-      break;
-    }
-
-  }
-
-  List_Delete(ListCurves);
-  List_Delete(ListSurfaces);
-  List_Delete(ListVolumes);
-}
-
-
-static void _get_all_model_points ( std::list<Vertex*> &mp )
-{
-  List_T *curves = Tree2List(THEM->Curves);
-  std::set<Vertex*> points;
-  for(int i = 0; i < List_Nbr(curves); i++){
-    Curve *c;
-    List_Read(curves, i, &c);
-    if (c->Num >=0)
-      {
-	if (c->beg && points.find(c->beg) == points.end())
-	  {
-	    points.insert(c->beg);
-	    mp.push_back(c->beg); 
-	  }
-	if (c->end && points.find(c->end) == points.end())
-	  {
-	    points.insert(c->end);
-	    mp.push_back(c->end); 
-	  }
-      }
-  }
-}
-
-static void _msh_print_all_modelpoints()
-{
-  std::list<Vertex*> mp;
-  _get_all_model_points ( mp );
-  std::list<Vertex*>::iterator it = mp.begin();
-  while (it != mp.end())
-    {
-      Vertex *v = (*it);
-      _msh_print_point(v);
-      it++;
-    }
-}
-
-static void _msh_print_all_curves(void *a, void *b)
-{
-  Curve *c = *(Curve **) a; 
-
-  Tree_Action(c->Simplexes, _msh_print_simplex);
-  Tree_Action(c->SimplexesBase, _msh_print_simplex);
-}
-
-static void _msh_print_all_surfaces(void *a, void *b)
-{
-  Surface *s = *(Surface **) a;
-  Tree_Action(s->Simplexes, _msh_print_simplex);
-  Tree_Action(s->SimplexesBase, _msh_print_simplex);
-  Tree_Action(s->Quadrangles, _msh_print_quadrangle);
-}
-
-static void _msh_print_all_simpsurf(void *a, void *b)
-{
-  Volume *v = *(Volume **) a;
-  Tree_Action(v->Simp_Surf, _msh_print_simplex);
-  Tree_Action(v->Quad_Surf, _msh_print_quadrangle);
-}
-
-static void _msh_print_all_linsurf(void *a, void *b)
-{
-  Volume *v = *(Volume **) a;
-  Tree_Action(v->Lin_Surf, _msh_print_simplex);
-}
-
-static void _msh_print_all_volumes(void *a, void *b)
-{
-  Volume *v = *(Volume **) a;
-  Tree_Action(v->Simplexes, _msh_print_simplex);
-  Tree_Action(v->SimplexesBase, _msh_print_simplex);
-  Tree_Action(v->Hexahedra, _msh_print_hexahedron);
-  Tree_Action(v->Prisms, _msh_print_prism);
-  Tree_Action(v->Pyramids, _msh_print_pyramid);
-}
-
-static void _msh_print_all_elements()
-{
-  MSH_PHYSICAL_NUM = 0;
-  MSH_PHYSICAL_ORI = 1;
-  MSH_LIN_NUM = MSH_SUR_NUM = MSH_VOL_NUM = 0;
-
-
-  _msh_print_all_modelpoints();
-
-  if(CTX.mesh.oldxtrude) {
-    Tree_Action(THEM->Volumes, _msh_print_all_simpsurf);
-    Tree_Action(THEM->Volumes, _msh_print_all_linsurf);
-  }
-  else {
-    Tree_Action(THEM->Curves, _msh_print_all_curves);
-    Tree_Action(THEM->Surfaces, _msh_print_all_surfaces);
-  }
-
-  Tree_Action(THEM->Volumes, _msh_print_all_volumes);
-}
-
-static void _msh_process_elements()
-{
-  MSH_ADD = 0;
-  MSH_ELEMENT_NUM = 1;
-
-  if(!List_Nbr(THEM->PhysicalGroups) || CTX.mesh.save_all) {
-    Msg(INFO, "Saving all elements (discarding physical groups)");
-    _msh_print_all_elements();
-  }
-  else
-    _msh_print_elements();
-
-  if(CTX.mesh.msh_file_version == 2.0)
-    fprintf(MSHFILE, "$Elements\n");
-  else
-    fprintf(MSHFILE, "$ELM\n");
-
-  fprintf(MSHFILE, "%d\n", MSH_ELEMENT_NUM - 1);
-
-  if(MSH_ELEMENT_NUM == 1)
-    Msg(WARNING, "No elements to save");
-
-  MSH_ADD = 1;
-  MSH_ELEMENT_NUM = 1;
-  if(!List_Nbr(THEM->PhysicalGroups) || CTX.mesh.save_all)
-    _msh_print_all_elements();
-  else
-    _msh_print_elements();
-
-  if(CTX.mesh.msh_file_version == 2.0)
-    fprintf(MSHFILE, "$EndElements\n");
-  else
-    fprintf(MSHFILE, "$ENDELM\n");
-}
-
-void Print_Mesh_MSH(FILE *fp)
-{
-  MSHFILE = fp;
-  if(CTX.mesh.msh_file_version == 1.0){
-    // OK, no header
-  }
-  else if(CTX.mesh.msh_file_version == 2.0){
-    fprintf(MSHFILE, "$MeshFormat\n");
-    fprintf(MSHFILE, "%g %d %d\n", CTX.mesh.msh_file_version,
-	    LIST_FORMAT_ASCII, (int)sizeof(double));
-    fprintf(MSHFILE, "$EndMeshFormat\n");
-  }
-  else{
-    Msg(GERROR, "Unknown MSH file version to generate (%g)", 
-	CTX.mesh.msh_file_version);
-    return;
-  }
-  _msh_process_nodes();
-  _msh_process_elements();
-  Msg(INFO, "%d nodes", MSH_NODE_NUM);
-  Msg(INFO, "%d elements", MSH_ELEMENT_NUM - 1);
-}
-
-// Write mesh in UNV format
-
-// IDEAS records
-#define HEADER       151
-#define UNITS        164
-#define NODES        2411
-#define ELEMENTS     2412
-#define RESNODE      55
-#define RESELEM      56
-#define RESVECT      57
-#define GROUPOFNODES 790
-
-// IDEAS elements
-#define BEAM         21
-#define BEAM2        24
-#define THINSHLL     91
-#define THINSHLL2    92
-#define QUAD         94
-#define QUAD2        95 // Ca c'est une impro !!!
-#define SOLIDFEM     111
-#define WEDGE        112
-#define BRICK        115
-#define SOLIDFEM2    118
-
-static int ELEMENT_ID;
-static Tree_T *tree;
-static int UNV_VOL_NUM;
-static FILE *UNVFILE;
-
-static void _unv_process_nodes()
-{
-  Vertex *v;
-  List_T *Nodes = Tree2List(THEM->Vertices);
-
-  fprintf(UNVFILE, "%6d\n", -1);
-  fprintf(UNVFILE, "%6d\n", NODES);
-  int nbnod = List_Nbr(Nodes);
-
-  for(int i = 0; i < nbnod; i++) {
-    List_Read(Nodes, i, &v);
-    int idnod = v->Num;
-    double x = v->Pos.X * CTX.mesh.scaling_factor;
-    double y = v->Pos.Y * CTX.mesh.scaling_factor;
-    double z = v->Pos.Z * CTX.mesh.scaling_factor;
-    fprintf(UNVFILE, "%10d%10d%10d%10d\n", idnod, 1, 1, 11);
-    char tmp[128];
-    // ugly hack to print number with 'D+XX' exponents
-    sprintf(tmp, "%25.16E%25.16E%25.16E\n", x, y, z);
-    tmp[21] = 'D';
-    tmp[46] = 'D';
-    tmp[71] = 'D';
-    fprintf(UNVFILE, tmp);    
-  }
-
-  List_Delete(Nodes);
-  fprintf(UNVFILE, "%6d\n", -1);
-}
-
-static void _unv_print_record(int num, int fetyp, int geo, int n, int nsup, 
-			      Vertex **v, Vertex **vsup)
-{
-  fprintf(UNVFILE, "%10d%10d%10d%10d%10d%10d\n",
-	  num, fetyp, geo, geo, 7, n + nsup);
-  if(fetyp == BEAM || fetyp == BEAM2)
-    fprintf(UNVFILE, "%10d%10d%10d\n", 0, 0, 0);
-  int ntot = 0;
-  for(int k = 0; k < n; k++) {
-    fprintf(UNVFILE, "%10d", v[k]->Num);
-    if(ntot % 8 == 7)
-      fprintf(UNVFILE, "\n");
-    ntot++;
-  }
-  for(int k = 0; k < nsup; k++) {
-    fprintf(UNVFILE, "%10d", vsup[k]->Num);
-    if(ntot % 8 == 7)
-      fprintf(UNVFILE, "\n");
-    ntot++;
-  }
-  if(ntot - 1 % 8 != 7)
-    fprintf(UNVFILE, "\n");
-}
-
-static void _unv_process_1D_elements()
-{
-  List_T *ListCurves = Tree2List(THEM->Curves);
-  List_T *Elements;
-  SimplexBase *sx;
-  Curve *c;
-
-  for(int i = 0; i < List_Nbr(ListCurves); i++) {
-    List_Read(ListCurves, i, &c);
-    for(int simtype = 0; simtype < 2; simtype ++){
-      Elements = (!simtype) ? Tree2List(c->Simplexes) : Tree2List(c->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Elements); j++) {
-	List_Read(Elements, j, &sx);
-	if(sx->VSUP)
-	  _unv_print_record(sx->Num, BEAM2, c->Num, 2, 1, &sx->V[0], sx->VSUP);
-	else 
-	  _unv_print_record(sx->Num, BEAM, c->Num, 2, 0, &sx->V[0], NULL);
-      }
-      List_Delete(Elements);
-    }
-  }
-
-  List_Delete(ListCurves);
-}
-
-static void _unv_process_2D_elements()
-{
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-  List_T *Elements;
-  Surface *s;
-  SimplexBase *sx;
-  Quadrangle *qx;
-
-  for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-      
-    // triangles
-    for(int simtype = 0; simtype < 2; simtype++){
-      Elements = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Elements); j++) {
-	List_Read(Elements, j, &sx);
-	if(sx->VSUP)
-	  _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 3, &sx->V[0], sx->VSUP);
-	else
-	  _unv_print_record(abs(sx->Num), THINSHLL, s->Num, 3, 0, &sx->V[0], NULL);
-      }
-      List_Delete(Elements);
-    }
-    
-    // quadrangles
-    Elements = Tree2List(s->Quadrangles);
-    for(int j = 0; j < List_Nbr(Elements); j++) {
-      List_Read(Elements, j, &qx);
-      if(qx->VSUP)
-	_unv_print_record(abs(qx->Num), QUAD, s->Num, 4, 4+1, &qx->V[0], qx->VSUP);
-      else
-	_unv_print_record(abs(qx->Num), QUAD2, s->Num, 4, 0, &qx->V[0], NULL);
-    }
-    List_Delete(Elements);
-  }
-  List_Delete(ListSurfaces);
-}
-
-static void _unv_process_3D_elements()
-{
-  List_T *ListVolumes = Tree2List(THEM->Volumes);
-  List_T *Elements;
-  SimplexBase *sx;
-  Hexahedron *hx;
-  Prism *px;
-  Volume *v;
-
-  for(int i = 0; i < List_Nbr(ListVolumes); i++) {
-    List_Read(ListVolumes, i, &v);
-
-    // Tets
-    for(int simtype = 0; simtype < 2; simtype++){
-      Elements = (!simtype) ? Tree2List(v->Simplexes) : Tree2List(v->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Elements); j++) {
-	List_Read(Elements, j, &sx);
-	if(sx->Volume_Simplexe() < 0) {
-	  Vertex *temp;
-	  temp = sx->V[0];
-	  sx->V[0] = sx->V[1];
-	  sx->V[1] = temp;
-	  if(sx->VSUP){
-	    temp = sx->VSUP[1];
-	    sx->VSUP[1] = sx->VSUP[2];
-	    sx->VSUP[2] = temp;
-	    temp = sx->VSUP[5];
-	    sx->VSUP[5] = sx->VSUP[3];
-	    sx->VSUP[3] = temp;
-	  }
-	}
-	if(sx->VSUP)
-	  _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 6, &sx->V[0], sx->VSUP);
-	else
-	  _unv_print_record(ELEMENT_ID++, SOLIDFEM, v->Num, 4, 0, &sx->V[0], NULL);
-      }
-      List_Delete(Elements);
-    }
-
-    // Prisms
-    Elements = Tree2List(v->Prisms);
-    for(int j = 0; j < List_Nbr(Elements); j++) {
-      List_Read(Elements, j, &px);
-      if(px->VSUP)
-	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 9+3, &px->V[0], px->VSUP);
-      else
-	_unv_print_record(ELEMENT_ID++, WEDGE, v->Num, 6, 0, &px->V[0], NULL);
-    }
-    List_Delete(Elements);
-
-    // Hexas
-    Elements = Tree2List(v->Hexahedra);
-    for(int j = 0; j < List_Nbr(Elements); j++) {
-      List_Read(Elements, j, &hx);
-      if(hx->VSUP)
-	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 12+6+1, &hx->V[0], hx->VSUP);
-      else
-	_unv_print_record(ELEMENT_ID++, BRICK, v->Num, 8, 0, &hx->V[0], NULL);
-    }
-    List_Delete(Elements);
-  }
-  List_Delete(ListVolumes);
-}
-
-static void _unv_add_vertex(void *a, void *b)
-{
-  Vertex *v = *(Vertex **) a;
-  if(Tree_Search(tree, &v->Num))
-    return;
-  Tree_Add(tree, &v->Num);
-  fprintf(UNVFILE, "%10d%10d%2d%2d%2d%2d%2d%2d\n", v->Num, 1, 0, 1, 0, 0, 0, 0);
-  fprintf(UNVFILE, "   0.0000000000000000D+00   1.0000000000000000D+00   0.0000000000000000D+00\n");
-  fprintf(UNVFILE, "   0.0000000000000000D+00   0.0000000000000000D+00   0.0000000000000000D+00\n");
-  fprintf(UNVFILE, "%10d%10d%10d%10d%10d%10d\n", 0, 0, 0, 0, 0, 0);
-}
-
-static void _unv_add_simplex_vertices(void *a, void *b)
-{
-  SimplexBase *s = *(SimplexBase **) a;
-  if(s->iEnt != UNV_VOL_NUM)
-    return;
-  for(int i = 0; i < 4; i++)
-    _unv_add_vertex(&s->V[i], NULL);
-}
-
-static void _unv_add_hexahedron_vertices(void *a, void *b)
-{
-  Hexahedron *h = *(Hexahedron **) a;
-  if(h->iEnt != UNV_VOL_NUM)
-    return;
-  for(int i = 0; i < 8; i++)
-    _unv_add_vertex(&h->V[i], NULL);
-}
-
-static void _unv_add_prism_vertices(void *a, void *b)
-{
-  Prism *p = *(Prism **) a;
-  if(p->iEnt != UNV_VOL_NUM)
-    return;
-  for(int i = 0; i < 6; i++)
-    _unv_add_vertex(&p->V[i], NULL);
-}
-
-static void _unv_add_pyramid_vertices(void *a, void *b)
-{
-  Pyramid *p = *(Pyramid **) a;
-  if(p->iEnt != UNV_VOL_NUM)
-    return;
-  for(int i = 0; i < 5; i++)
-    _unv_add_vertex(&p->V[i], NULL);
-}
-
-static void _unv_process_groups()
-{
-  Volume *pV;
-  Surface *ps, s;
-  Curve *pc, c;
-  Vertex *pv, v;
-  PhysicalGroup *p;
-  List_T *ListVolumes;
-
-  for(int i = 0; i < List_Nbr(THEM->PhysicalGroups); i++) {
-
-    List_Read(THEM->PhysicalGroups, i, &p);
-
-    fprintf(UNVFILE, "%6d\n", -1);
-    fprintf(UNVFILE, "%6d\n", GROUPOFNODES);
-    fprintf(UNVFILE, "%10d%10d\n", p->Num, 1);
-    fprintf(UNVFILE, "LOAD SET %2d\n", 1);
-
-    tree = Tree_Create(sizeof(int), fcmp_absint);
-
-    switch (p->Typ) {
-    case MSH_PHYSICAL_VOLUME:
-      ListVolumes = Tree2List(THEM->Volumes);
-      for(int k = 0; k < List_Nbr(ListVolumes); k++) {
-        List_Read(ListVolumes, k, &pV);
-        for(int j = 0; j < List_Nbr(p->Entities); j++) {
-          List_Read(p->Entities, j, &UNV_VOL_NUM);
-          Tree_Action(pV->Simplexes, _unv_add_simplex_vertices);
-          Tree_Action(pV->SimplexesBase, _unv_add_simplex_vertices);
-          Tree_Action(pV->Hexahedra, _unv_add_hexahedron_vertices);
-          Tree_Action(pV->Prisms, _unv_add_prism_vertices);
-          Tree_Action(pV->Pyramids, _unv_add_pyramid_vertices);
-        }
-      }
-      List_Delete(ListVolumes);
-      break;
-    case MSH_PHYSICAL_SURFACE:
-      for(int j = 0; j < List_Nbr(p->Entities); j++) {
-        ps = &s;
-        List_Read(p->Entities, j, &ps->Num);
-        if(Tree_Query(THEM->Surfaces, &ps))
-          Tree_Action(ps->Vertices, _unv_add_vertex);
-      }
-      break;
-    case MSH_PHYSICAL_LINE:
-      for(int j = 0; j < List_Nbr(p->Entities); j++) {
-        pc = &c;
-        List_Read(p->Entities, j, &pc->Num);
-        if(Tree_Query(THEM->Curves, &pc))
-          for(int k = 0; k < List_Nbr(pc->Vertices); k++)
-            _unv_add_vertex(List_Pointer(pc->Vertices, k), NULL);
-      }
-      break;
-    case MSH_PHYSICAL_POINT:
-      for(int j = 0; j < List_Nbr(p->Entities); j++) {
-        pv = &v;
-        List_Read(p->Entities, j, &pv->Num);
-        if(Tree_Query(THEM->Vertices, &pv))
-          _unv_add_vertex(&pv, NULL);
-      }
-      break;
-    }
-
-    Tree_Delete(tree);
-
-    fprintf(UNVFILE, "%6d\n", -1);
-  }
-}
-
-void Print_Mesh_UNV(FILE *fp)
-{
-  UNVFILE = fp;
-  _unv_process_nodes();
-  fprintf(UNVFILE, "%6d\n", -1);
-  fprintf(UNVFILE, "%6d\n", ELEMENTS);
-  ELEMENT_ID = 1;
-  _unv_process_3D_elements();
-  _unv_process_2D_elements();
-  _unv_process_1D_elements();
-  fprintf(UNVFILE, "%6d\n", -1);
-  _unv_process_groups();
-}
-
-// Write mesh in Gref format
-
-static FILE *GREFFILE;
-
-static void _gref_consecutive_nodes(Tree_T *ConsecutiveNTree,
-				    Tree_T *ConsecutiveETree)
-{
-  SimplexBase *sx;
-  Quadrangle *qx;
-  Surface *s;
-  int nbnod, nbedges, nbdof;
-
-  int newnum = 0;
-
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-  for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-    // triangles
-    for(int simtype = 0; simtype < 2; simtype++){
-      List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Triangles); j++) {
-	List_Read(Triangles, j, &sx);
-	for(int k = 0; k < 3; k++) {
-	  if(sx->V[k]->Frozen >= 0) {
-	    sx->V[k]->Frozen = --newnum;
-	    Tree_Insert(ConsecutiveNTree, &(sx->V[k]));
-	  }
-	}
-      }
-      List_Delete(Triangles);
-    }
-    // quads
-    List_T *Quadrangles = Tree2List(s->Quadrangles);
-    for(int j = 0; j < List_Nbr(Quadrangles); j++) {
-      List_Read(Quadrangles, j, &qx);
-      for(int k = 0; k < 4; k++) {
-        if(qx->V[k]->Frozen >= 0) {
-          qx->V[k]->Frozen = --newnum;
-          Tree_Insert(ConsecutiveNTree, &(qx->V[k]));
-        }
-      }
-    }
-    List_Delete(Quadrangles);
-  }
-
-  nbnod = -newnum;
-
-  for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-    // triangles
-    for(int simtype = 0; simtype < 2; simtype++){
-      List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Triangles); j++) {
-	List_Read(Triangles, j, &sx);
-	for(int k = 0; k < 3; k++) {
-	  if(sx->VSUP[k]->Frozen >= 0) {
-	    sx->VSUP[k]->Frozen = --newnum;
-	    Tree_Insert(ConsecutiveETree, &(sx->VSUP[k]));
-	  }
-	}
-      }
-      List_Delete(Triangles);
-    }
-    // quads
-    List_T *Quadrangles = Tree2List(s->Quadrangles);
-    for(int j = 0; j < List_Nbr(Quadrangles); j++) {
-      List_Read(Quadrangles, j, &qx);
-      for(int k = 0; k < 4; k++) {
-	if(qx->VSUP[k]->Frozen >= 0) {
-	  qx->VSUP[k]->Frozen = --newnum;
-	  Tree_Insert(ConsecutiveETree, &(qx->VSUP[k]));
-	}
-      }
-    }
-    List_Delete(Quadrangles);
-  }
-
-  List_Delete(ListSurfaces);
-
-  nbedges = -newnum - nbnod;
-  nbdof = nbnod + nbedges;
-  Msg(INFO, "%d Dofs", nbdof);
-}
-
-static void _gref_end_consecutive_nodes()
-{
-  SimplexBase *sx;
-  Quadrangle *qx;
-  Surface *s;
-
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-  for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-    // triangles
-    for(int simtype = 0; simtype < 2; simtype++){
-      List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Triangles); j++) {
-	List_Read(Triangles, j, &sx);
-	for(int k = 0; k < 3; k++)
-	  sx->V[k]->Frozen = 0;
-	for(int k = 0; k < 3; k++)
-	  sx->VSUP[k]->Frozen = 0;
-      }
-      List_Delete(Triangles);
-    }
-    // quads
-    List_T *Quadrangles = Tree2List(s->Quadrangles);
-    for(int j = 0; j < List_Nbr(Quadrangles); j++) {
-      List_Read(Quadrangles, j, &qx);
-      for(int k = 0; k < 4; k++)
-        qx->V[k]->Frozen = 0;
-      for(int k = 0; k < 4; k++)
-	qx->VSUP[k]->Frozen = 0;
-    }
-    List_Delete(Quadrangles);
-  }
-  List_Delete(ListSurfaces);
-}
-
-static int _gref_compare_frozen(const void *a, const void *b)
-{
-  Vertex *q, *w;
-  q = *(Vertex **) a;
-  w = *(Vertex **) b;
-  return w->Frozen - q->Frozen;
-}
-
-static int _gref_process_nodes(Tree_T *ConsecutiveNTree, 
-			       Tree_T *ConsecutiveETree)
-{
-  int i, nbtriqua;
-  Vertex *v;
-  Surface *s;
-  List_T *Nodes;
-
-  Tree_Action(THEM->Curves, Degre2_Curve);
-  Tree_Action(THEM->Surfaces, Degre2_Surface);
- 
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-  nbtriqua = 0;
-  for(i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-    nbtriqua += Tree_Nbr(s->Simplexes) + Tree_Nbr(s->SimplexesBase);
-    nbtriqua += Tree_Nbr(s->Quadrangles);
-  }
-  List_Delete(ListSurfaces);
-
-  _gref_consecutive_nodes(ConsecutiveNTree, ConsecutiveETree);
-
-  fprintf(GREFFILE, "%d %d %d\n", nbtriqua, Tree_Nbr(ConsecutiveNTree),
-          Tree_Nbr(ConsecutiveNTree) + Tree_Nbr(ConsecutiveETree));
-
-  Nodes = Tree2List(ConsecutiveNTree);
-  for(i = 0; i < List_Nbr(Nodes); i++) {
-    List_Read(Nodes, i, &v);
-    fprintf(GREFFILE, "%21.16e ", v->Pos.X * CTX.mesh.scaling_factor);
-    if(i % 3 == 2)
-      fprintf(GREFFILE, "\n");
-  }
-  if((List_Nbr(Nodes) - 1) % 3 != 2)
-    fprintf(GREFFILE, "\n");
-  for(i = 0; i < List_Nbr(Nodes); i++) {
-    List_Read(Nodes, i, &v);
-    fprintf(GREFFILE, "%21.16e ", v->Pos.Y * CTX.mesh.scaling_factor);
-    if(i % 3 == 2)
-      fprintf(GREFFILE, "\n");
-  }
-  if((List_Nbr(Nodes) - 1) % 3 != 2)
-    fprintf(GREFFILE, "\n");
-  i = Tree_Nbr(ConsecutiveNTree);
-  List_Delete(Nodes);
-  return i;
-}
-
-static int _gref_find_physical(Vertex *v)
-{
-  PhysicalGroup *p;
-  Curve *c;
-  int i, j;
-  for(i = 0; i < List_Nbr(THEM->PhysicalGroups); i++) {
-    List_Read(THEM->PhysicalGroups, i, &p);
-    if(p->Typ == MSH_PHYSICAL_POINT) {
-      if(List_Search(p->Entities, &v->Num, fcmp_absint)) {
-        return p->Num;
-      }
-    }
-  }
-
-  if(v->ListCurves) {
-    for(i = 0; i < List_Nbr(THEM->PhysicalGroups); i++) {
-      List_Read(THEM->PhysicalGroups, i, &p);
-      if(p->Typ == MSH_PHYSICAL_LINE) {
-        for(j = 0; j < List_Nbr(v->ListCurves); j++) {
-          List_Read(v->ListCurves, j, &c);
-          if(List_Search(p->Entities, &c->Num, fcmp_absint)) {
-            return p->Num;
-          }
-        }
-      }
-    }
-  }
-  return 0;
-}
-
-static void _gref_process_boundary_conditions(Tree_T *TRN, Tree_T *TRE)
-{
-  int i, ent;
-  Vertex *v;
-
-  List_T *Nodes = Tree2List(TRN);
-  for(i = 0; i < List_Nbr(Nodes); i++) {
-    List_Read(Nodes, i, &v);
-    ent = _gref_find_physical(v);
-    fprintf(GREFFILE, "%d %d ", ent, ent);
-    if(i % 3 == 2)
-      fprintf(GREFFILE, "\n");
-  }
-  if((List_Nbr(Nodes) - 1) % 3 != 2)
-    fprintf(GREFFILE, "\n");
-  List_Delete(Nodes);
-
-  Nodes = Tree2List(TRE);
-  for(i = 0; i < List_Nbr(Nodes); i++) {
-    List_Read(Nodes, i, &v);
-    ent = _gref_find_physical(v);
-    fprintf(GREFFILE, "%d %d ", ent, ent);
-    if(i % 3 == 2)
-      fprintf(GREFFILE, "\n");
-  }
-  if((List_Nbr(Nodes) - 1) % 3 != 2)
-    fprintf(GREFFILE, "\n");
-  List_Delete(Nodes);
-}
-
-static void _gref_process_elements(int nn)
-{
-  SimplexBase *sx;
-  Quadrangle *qx;
-  Surface *s;
-
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-
-  for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-    // triangles
-    for(int simtype = 0; simtype < 2; simtype++){
-      List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Triangles); j++) {
-	List_Read(Triangles, j, &sx);
-	fprintf(GREFFILE, "%d %d %d\n", -sx->V[0]->Frozen,
-		-sx->V[1]->Frozen, -sx->V[2]->Frozen);
-      }
-      List_Delete(Triangles);
-    }
-    // quads
-    List_T *Quadrangles = Tree2List(s->Quadrangles);
-    for(int j = 0; j < List_Nbr(Quadrangles); j++) {
-      List_Read(Quadrangles, j, &qx);
-      fprintf(GREFFILE, "%d %d %d %d\n", -qx->V[0]->Frozen,
-	      -qx->V[1]->Frozen, -qx->V[2]->Frozen, -qx->V[3]->Frozen);
-    }
-    List_Delete(Quadrangles);
-  }
-
-  // Degres de Liberte
-  for(int i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &s);
-    // triangles
-    for(int simtype = 0; simtype < 2; simtype++){
-      List_T *Triangles = (!simtype) ? Tree2List(s->Simplexes) : Tree2List(s->SimplexesBase);
-      for(int j = 0; j < List_Nbr(Triangles); j++) {
-	List_Read(Triangles, j, &sx);
-	fprintf(GREFFILE, "%d %d %d %d %d %d %d %d %d %d %d %d\n",
-		-2 * sx->V[0]->Frozen - 1,
-		-2 * sx->V[0]->Frozen,
-		-2 * sx->VSUP[0]->Frozen - 1,
-		-2 * sx->VSUP[0]->Frozen,
-		-2 * sx->V[1]->Frozen - 1,
-		-2 * sx->V[1]->Frozen,
-		-2 * sx->VSUP[1]->Frozen - 1,
-		-2 * sx->VSUP[1]->Frozen,
-		-2 * sx->V[2]->Frozen - 1,
-		-2 * sx->V[2]->Frozen,
-		-2 * sx->VSUP[2]->Frozen - 1, -2 * sx->VSUP[2]->Frozen);
-      }
-      List_Delete(Triangles);
-    }
-    // quads
-    List_T *Quadrangles = Tree2List(s->Quadrangles);
-    for(int j = 0; j < List_Nbr(Quadrangles); j++) {
-      List_Read(Quadrangles, j, &qx);
-      fprintf(GREFFILE, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d\n",
-	      -2 * qx->V[0]->Frozen - 1,
-	      -2 * qx->V[0]->Frozen,
-	      -2 * qx->VSUP[0]->Frozen - 1,
-	      -2 * qx->VSUP[0]->Frozen,
-	      -2 * qx->V[1]->Frozen - 1,
-	      -2 * qx->V[1]->Frozen,
-	      -2 * qx->VSUP[1]->Frozen - 1,
-	      -2 * qx->VSUP[1]->Frozen,
-	      -2 * qx->V[2]->Frozen - 1,
-	      -2 * qx->V[2]->Frozen,
-	      -2 * qx->VSUP[2]->Frozen - 1,
-	      -2 * qx->VSUP[2]->Frozen,
-	      -2 * qx->V[3]->Frozen - 1,
-	      -2 * qx->V[3]->Frozen,
-	      -2 * qx->VSUP[3]->Frozen - 1, -2 * qx->VSUP[3]->Frozen);
-    }
-    List_Delete(Quadrangles);
-  }
-  List_Delete(ListSurfaces);
-}
-
-void Print_Mesh_GREF(FILE *fp)
-{
-  GREFFILE = fp;
-  Tree_T *TRN = Tree_Create(sizeof(Vertex *), _gref_compare_frozen);
-  Tree_T *TRE = Tree_Create(sizeof(Vertex *), _gref_compare_frozen);
-  _gref_process_nodes(TRN, TRE);
-  _gref_process_elements(Tree_Nbr(TRN));
-  _gref_process_boundary_conditions(TRN, TRE);
-  Tree_Delete(TRN);
-  Tree_Delete(TRE);
-  _gref_end_consecutive_nodes();
-}
-
-// Write mesh in VRML 1 format
-
-static List_T *wrlnodes = NULL;
-static FILE *WRLFILE;
-
-static void _wrl_print_node(void *a, void *b)
-{
-  Vertex *V = *(Vertex **) a;
-  fprintf(WRLFILE, "%.16g %.16g %.16g,\n",
-          V->Pos.X * CTX.mesh.scaling_factor,
-          V->Pos.Y * CTX.mesh.scaling_factor,
-          V->Pos.Z * CTX.mesh.scaling_factor);
-  List_Add(wrlnodes, &V->Num);
-}
-
-static void _wrl_process_nodes()
-{
-  if(!wrlnodes)
-    wrlnodes = List_Create(Tree_Size(THEM->Vertices), 100, sizeof(int));
-  else
-    List_Reset(wrlnodes);
-  fprintf(WRLFILE, "#VRML V1.0 ascii\n");
-  fprintf(WRLFILE, "#created by Gmsh\n");
-  fprintf(WRLFILE, "Coordinate3 {\n");
-  fprintf(WRLFILE, "  point [\n");
-  Tree_Action(THEM->Vertices, _wrl_print_node);
-  fprintf(WRLFILE, "  ]\n");
-  fprintf(WRLFILE, "}\n");
-}
-
-static void _wrl_print_line(void *a, void *b)
-{
-  SimplexBase *s = *(SimplexBase **) a;
-  for(int i = 0; i < 2; i++){
-    if(s->V[i]){
-      int j = List_ISearch(wrlnodes, &s->V[i]->Num, fcmp_int);
-      if(j < 0)
-	Msg(GERROR, "Unknown node %d in line %d", s->V[i]->Num, s->Num);
-      else
-	fprintf(WRLFILE, "%d,", j);
-    }
-  }
-  fprintf(WRLFILE, "-1,\n");
-}
-
-static void _wrl_print_triangle(void *a, void *b)
-{
-  SimplexBase *s = *(SimplexBase **) a;
-  for(int i = 0; i < 3; i++){
-    if(s->V[i]){
-      int j = List_ISearch(wrlnodes, &s->V[i]->Num, fcmp_int);
-      if(j < 0)
-	Msg(GERROR, "Unknown node %d in triangle %d", s->V[i]->Num, s->Num);
-      else
-	fprintf(WRLFILE, "%d,", j);
-    }
-  }
-  fprintf(WRLFILE, "-1,\n");
-}
-
-static void _wrl_print_quadrangle(void *a, void *b)
-{
-  Quadrangle *q = *(Quadrangle **) a;
-  for(int i = 0; i < 4; i++){
-    if(q->V[i]){
-      int j = List_ISearch(wrlnodes, &q->V[i]->Num, fcmp_int);
-      if(j < 0)
-	Msg(GERROR, "Unknown node %d in quadrangle %d", q->V[i]->Num, q->Num);
-      else
-	fprintf(WRLFILE, "%d,", j);
-    }
-  }
-  fprintf(WRLFILE, "-1,\n");
-}
-
-static void _wrl_print_all_curves(void *a, void *b)
-{
-  Curve *c = *(Curve **) a;
-  if(c->Num < 0)
-    return;
-  fprintf(WRLFILE, "DEF Curve%d IndexedLineSet {\n", c->Num);
-  fprintf(WRLFILE, "  coordIndex [\n");
-  Tree_Action(c->Simplexes, _wrl_print_line);
-  Tree_Action(c->SimplexesBase, _wrl_print_line);
-  fprintf(WRLFILE, "  ]\n");
-  fprintf(WRLFILE, "}\n");
-}
-
-static void _wrl_print_all_surfaces(void *a, void *b)
-{
-  Surface *s = *(Surface **) a;
-  fprintf(WRLFILE, "DEF Surface%d IndexedFaceSet {\n", s->Num);
-  fprintf(WRLFILE, "  coordIndex [\n");
-  Tree_Action(s->Simplexes, _wrl_print_triangle);
-  Tree_Action(s->SimplexesBase, _wrl_print_triangle);
-  Tree_Action(s->Quadrangles, _wrl_print_quadrangle);
-  fprintf(WRLFILE, "  ]\n");
-  fprintf(WRLFILE, "}\n");
-}
-
-static void _wrl_process_elements()
-{
-  if(!wrlnodes)
-    Msg(GERROR, "VRML node list does not exist");
-  else {
-    Tree_Action(THEM->Curves, _wrl_print_all_curves);
-    Tree_Action(THEM->Surfaces, _wrl_print_all_surfaces);
-  }
-}
-
-void Print_Mesh_WRL(FILE *fp)
-{
-  WRLFILE = fp;
-  _wrl_process_nodes();
-  _wrl_process_elements();
-}
-
-// Write surface mesh in STL format
-
-static FILE *STLFILE;
-
-static void _stl_print_triangle(void *a, void *b)
-{
-  SimplexBase *s = *(SimplexBase **) a;
-  
-  if(!s->V[2]) return;
-  
-  double n[3];
-  normal3points(s->V[0]->Pos.X, s->V[0]->Pos.Y, s->V[0]->Pos.Z, 
-		s->V[1]->Pos.X, s->V[1]->Pos.Y, s->V[1]->Pos.Z, 
-		s->V[2]->Pos.X, s->V[2]->Pos.Y, s->V[2]->Pos.Z, n);
-
-  fprintf(STLFILE, "facet normal %g %g %g\n", n[0], n[1], n[2]);
-  fprintf(STLFILE, "  outer loop\n");
-  fprintf(STLFILE, "    vertex %g %g %g\n", s->V[0]->Pos.X, s->V[0]->Pos.Y, s->V[0]->Pos.Z);
-  fprintf(STLFILE, "    vertex %g %g %g\n", s->V[1]->Pos.X, s->V[1]->Pos.Y, s->V[1]->Pos.Z);
-  fprintf(STLFILE, "    vertex %g %g %g\n", s->V[2]->Pos.X, s->V[2]->Pos.Y, s->V[2]->Pos.Z);
-  fprintf(STLFILE, "  endloop\n");
-  fprintf(STLFILE, "endfacet\n");
-}
-
-static void _stl_print_quadrangle(void *a, void *b)
-{
-  Quadrangle *q = *(Quadrangle **) a;
-  
-  double n[3];
-  normal3points(q->V[0]->Pos.X, q->V[0]->Pos.Y, q->V[0]->Pos.Z, 
-		q->V[1]->Pos.X, q->V[1]->Pos.Y, q->V[1]->Pos.Z, 
-		q->V[2]->Pos.X, q->V[2]->Pos.Y, q->V[2]->Pos.Z, n);
-
-  fprintf(STLFILE, "facet normal %g %g %g\n", n[0], n[1], n[2]);
-  fprintf(STLFILE, "  outer loop\n");
-  fprintf(STLFILE, "    vertex %g %g %g\n", q->V[0]->Pos.X, q->V[0]->Pos.Y, q->V[0]->Pos.Z);
-  fprintf(STLFILE, "    vertex %g %g %g\n", q->V[1]->Pos.X, q->V[1]->Pos.Y, q->V[1]->Pos.Z);
-  fprintf(STLFILE, "    vertex %g %g %g\n", q->V[2]->Pos.X, q->V[2]->Pos.Y, q->V[2]->Pos.Z);
-  fprintf(STLFILE, "  endloop\n");
-  fprintf(STLFILE, "endfacet\n");
-  fprintf(STLFILE, "facet normal %g %g %g\n", n[0], n[1], n[2]);
-  fprintf(STLFILE, "  outer loop\n");
-  fprintf(STLFILE, "    vertex %g %g %g\n", q->V[0]->Pos.X, q->V[0]->Pos.Y, q->V[0]->Pos.Z);
-  fprintf(STLFILE, "    vertex %g %g %g\n", q->V[2]->Pos.X, q->V[2]->Pos.Y, q->V[2]->Pos.Z);
-  fprintf(STLFILE, "    vertex %g %g %g\n", q->V[3]->Pos.X, q->V[3]->Pos.Y, q->V[3]->Pos.Z);
-  fprintf(STLFILE, "  endloop\n");
-  fprintf(STLFILE, "endfacet\n");
-}
-
-static void _stl_print_all_surfaces(void *a, void *b)
-{
-  Surface *s = *(Surface **) a;
-  Tree_Action(s->Simplexes, _stl_print_triangle);
-  Tree_Action(s->SimplexesBase, _stl_print_triangle);
-  Tree_Action(s->Quadrangles, _stl_print_quadrangle);
-}
-
-void Print_Mesh_STL(FILE *fp)
-{
-  STLFILE = fp;
-  fprintf(STLFILE, "solid Created by Gmsh\n");
-  Tree_Action(THEM->Surfaces, _stl_print_all_surfaces);
-  fprintf(STLFILE, "endsolid Created by Gmsh\n");
-}
-
-// Write mesh in DMG format
-
-static int _dmg_is_topologic(Vertex *v, List_T *curves)
-{
-  Curve *c;
-  for(int i = 0; i < List_Nbr(curves); i++) {
-    List_Read(curves, i, &c);
-    if(!compareVertex(&v, &c->beg))
-      return 1;
-  }
-  return 0;
-}
-
-void Print_Mesh_DMG(FILE *fp)
-{
-  int i, j;
-  List_T *ll, *l;
-  Vertex *v;
-  Curve *c;
-  Surface *s;
-  int k;
-
-  l = Tree2List(THEM->Points);
-  ll = Tree2List(THEM->Curves);
-
-  k = 0;
-  for(i = 0; i < List_Nbr(l); i++) {
-    List_Read(l, i, &v);
-    if(_dmg_is_topologic(v, ll)) {
-      k++;
-    }
-  }
-
-  // write first the global infos 
-
-  fprintf(fp, "%d %d %d %d \n", Tree_Nbr(THEM->Volumes),
-          Tree_Nbr(THEM->Surfaces),
-          Tree_Nbr(THEM->Curves) / 2,     // the 2 is for the reverse curves
-          k);
-
-  // then write the bounding box
-
-  fprintf(fp, "%12.5E %12.5E %12.5E \n", 
-	  CTX.min[0], CTX.min[1], CTX.min[2]);
-  fprintf(fp, "%12.5E %12.5E %12.5E \n", 
-	  CTX.max[0], CTX.max[1], CTX.max[2]);
-
-  // write the points
-  k = 0;
-  for(i = 0; i < List_Nbr(l); i++) {
-    List_Read(l, i, &v);
-    if(_dmg_is_topologic(v, ll)) {
-      v->Frozen = k++;
-      fprintf(fp, "%d %12.5E %12.5E %12.5E \n", 
-	      v->Frozen, v->Pos.X, v->Pos.Y, v->Pos.Z);
-    }
-  }
-  List_Delete(l);
-
-  // write the curves
-  l = ll;
-  k = 0;
-  for(i = 0; i < List_Nbr(l); i++) {
-    List_Read(l, i, &c);
-    if(c->Num > 0) {
-      c->ipar[3] = k;
-      Curve *cinv = FindCurve(-c->Num, THEM);
-      cinv->ipar[3] = k++;
-      fprintf(fp, "%d %d %d \n", 
-	      c->ipar[3], c->beg->Frozen, c->end->Frozen);
-    }
-  }
-
-  List_Delete(l);
-
-  // write the surfaces
-  l = Tree2List(THEM->Surfaces);
-
-  for(i = 0; i < List_Nbr(l); i++) {
-    List_Read(l, i, &s);
-
-    int numEdgeLoop[2000], iLoop = 0;
-    Vertex *beg = NULL;
-    numEdgeLoop[iLoop] = 0;
-    int deb = 1;
-    for(j = 0; j < List_Nbr(s->Generatrices); j++) {
-      List_Read(s->Generatrices, j, &c);
-      if(deb) {
-        beg = c->beg;
-        deb = 0;
-      }
-      Msg(INFO, "beg->%d end->%d", c->beg->Num, c->end->Num);
-      (numEdgeLoop[iLoop])++;
-      if(c->end == beg) {
-        iLoop++;
-        numEdgeLoop[iLoop] = 0;
-        deb = 1;
-      }
-    }
-    s->ipar[3] = i;
-    fprintf(fp, "%d %d\n", i, iLoop);
-    int iEdge = 0;
-    for(k = 0; k < iLoop; k++) {
-      fprintf(fp, "%d ", numEdgeLoop[k]);
-      for(j = 0; j < numEdgeLoop[k]; j++) {
-        List_Read(s->Generatrices, iEdge++, &c);
-        fprintf(fp, "%d %d ", abs(c->ipar[3]), (c->Num > 0) ? 1 : -1);
-      }
-      fprintf(fp, "\n");
-    }
-  }
-  List_Delete(l);
-
-  // write the volumes (2 b continued...)
-}
-
-// Write mesh in Plot3D format, ASCII structured, multi-zone
-// FIXME: still need to implement this for extruded grids
-
-static FILE *P3DFILE;
-
-int _p3d_cmp_entities(const void *a, const void *b)
-{
-  Element **e1 = (Element **) a;
-  Element **e2 = (Element **) b;
-
-  return (*e1)->iEnt - (*e2)->iEnt;
-}
-
-int _p3d_cmp_surf_num(const void *a, const void *b)
-{
-  Surface **e1 = (Surface **) a;
-  Surface **e2 = (Surface **) b;
-
-  return (*e1)->Num - (*e2)->Num;
-}
-
-int _p3d_cmp_vol_num(const void *a, const void *b)
-{
-  Volume **e1 = (Volume **) a;
-  Volume **e2 = (Volume **) b;
-
-  return (*e1)->Num - (*e2)->Num;
-}
-
-void _p3d_print_quads(List_T *ListQuads, int Nu, int Nv)
-{
-  int i = 0, j = 0, c = 1, curQuad = 0;
-  double coord;
-  Quadrangle *pQ;
-
-  for (c = 0; c < 3; c++) {
-    for (j = 0; j < (Nu - 1); j++) {
-      for (i = 0; i < (Nv - 1); i++) {
-        curQuad = i + (j * (Nv - 1));
-        List_Read(ListQuads, curQuad, &pQ);
-        coord = (c==0) ? pQ->V[0]->Pos.X : ((c==1) ? pQ->V[0]->Pos.Y : pQ->V[0]->Pos.Z);
-        fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
-      }
-      coord = (c==0) ? pQ->V[3]->Pos.X : ((c==1) ? pQ->V[3]->Pos.Y : pQ->V[3]->Pos.Z);
-      fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
-    } j--;
-    for (i = 0; i < (Nv - 1); i++) {
-      curQuad = i + (j * (Nv - 1));
-      List_Read(ListQuads, curQuad, &pQ);
-      coord = (c==0) ? pQ->V[1]->Pos.X : ((c==1) ? pQ->V[1]->Pos.Y : pQ->V[1]->Pos.Z);
-      fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
-    }
-    coord = (c==0) ? pQ->V[2]->Pos.X : ((c==1) ? pQ->V[2]->Pos.Y : pQ->V[2]->Pos.Z);
-    fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
-  }
-}
-
-void _p3d_print_hex(List_T *ListHex, int Nu, int Nv, int Nw)
-{
-  int i = 0, j = 0, k = 0, c = 0, curHex = 0;
-  double coord;
-  Hexahedron *pH;
-
-  for (c = 0; c < 3; c++) {
-    for (k = 0; k < (Nu - 1); k++) {
-      for (j = 0; j < (Nv - 1); j++) {
-        for (i = 0; i < (Nw - 1); i++) {
-          curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
-          List_Read(ListHex, curHex, &pH);
-          coord = (c==0) ? pH->V[0]->Pos.X : ((c==1) ? pH->V[0]->Pos.Y : pH->V[0]->Pos.Z);
-          fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
-        }
-       coord = (c==0) ? pH->V[4]->Pos.X : ((c==1) ? pH->V[4]->Pos.Y : pH->V[4]->Pos.Z);
-       fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
-      } j--;
-      for (i = 0; i < (Nw - 1); i++) {
-        curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
-        List_Read(ListHex, curHex, &pH);
-        coord = (c==0) ? pH->V[3]->Pos.X : ((c==1) ? pH->V[3]->Pos.Y : pH->V[3]->Pos.Z);
-        fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
-      }
-      coord = (c==0) ? pH->V[7]->Pos.X : ((c==1) ? pH->V[7]->Pos.Y : pH->V[7]->Pos.Z);
-      fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
-    } k--;
-    for (j = 0; j < (Nv - 1); j++) {
-      for (i = 0; i < (Nw - 1); i++) {
-        curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
-        List_Read(ListHex, curHex, &pH);
-        coord = (c==0) ? pH->V[1]->Pos.X : ((c==1) ? pH->V[1]->Pos.Y : pH->V[1]->Pos.Z);
-        fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
-      }
-      coord = (c==0) ? pH->V[5]->Pos.X : ((c==1) ? pH->V[5]->Pos.Y : pH->V[5]->Pos.Z);
-      fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
-    } j--;
-    for (i = 0; i < (Nw - 1); i++) {
-      curHex = i + (j * (Nw - 1)) + (k * (Nw - 1) * (Nv - 1));
-      List_Read(ListHex, curHex, &pH);
-      coord = (c==0) ? pH->V[2]->Pos.X : ((c==1) ? pH->V[2]->Pos.Y : pH->V[2]->Pos.Z);
-      fprintf(P3DFILE, "%g ", coord * CTX.mesh.scaling_factor);
-    }
-    coord = (c==0) ? pH->V[6]->Pos.X : ((c==1) ? pH->V[6]->Pos.Y : pH->V[6]->Pos.Z);
-    fprintf(P3DFILE, "%g\n", coord * CTX.mesh.scaling_factor);
-  }
-}
-
-void _p3d_print_all_elements()
-{
-  int i, Nv, Nu, ElemCnt = 0;
-  Volume *pV;
-  Surface *pS;
-  List_T *ListQuads;
-  List_T *ListHex;
-
-  List_T *ListSurfaces = Tree2List(THEM->Surfaces);
-  List_T *ListVolumes = Tree2List(THEM->Volumes);
-
-  List_T *ListStructSurf = List_Create(1,1,sizeof(Surface *));
-  List_T *ListStructVol = List_Create(1,1,sizeof(Volume *));
-
-  List_Sort(ListSurfaces, _p3d_cmp_surf_num);
-  List_Sort(ListVolumes, _p3d_cmp_vol_num);
-
-  for(i = 0; i < List_Nbr(ListSurfaces); i++) {
-    List_Read(ListSurfaces, i, &pS);
-    if (pS->Nu &&
-        pS->Nv &&
-        (Tree_Nbr(pS->Quadrangles) == (pS->Nu - 1) * (pS->Nv - 1))){
-        ElemCnt++;
-        List_Add(ListStructSurf, &pS);
-    }
-  }
-
-  for(i = 0; i < List_Nbr(ListVolumes); i++) {
-    List_Read(ListVolumes, i, &pV);
-    List_Read(pV->Surfaces, 0, &pS);
-    Nu = pS->Nu;
-    Nv = pS->Nv;
-    List_Read(pV->Surfaces, 1, &pS);
-    if (Nu &&
-        Nv &&
-        pS->Nv &&
-        (Tree_Nbr(pV->Hexahedra) == (Nu - 1) * (Nv - 1) * (pS->Nv - 1))){
-        ElemCnt++;
-        List_Add(ListStructVol, &pV);
-    }
-  }
-
-  fprintf(P3DFILE, "%d\n", ElemCnt);
-
-  for(i = 0; i < List_Nbr(ListStructSurf); i++) {
-    List_Read(ListStructSurf, i, &pS);
-    fprintf(P3DFILE, "%d %d 1\n", pS->Nv, pS->Nu);
-  }
-
-  for(i = 0; i < List_Nbr(ListStructVol); i++) {
-    List_Read(ListStructVol, i, &pV);
-    List_Read(pV->Surfaces, 0, &pS);
-    Nu = pS->Nu;
-    Nv = pS->Nv;
-    List_Read(pV->Surfaces, 1, &pS);
-    fprintf(P3DFILE, "%d %d %d\n", pS->Nv, Nv, Nu);
-  }
-
-  for(i = 0; i < List_Nbr(ListStructSurf); i++) {
-    List_Read(ListStructSurf, i, &pS);
-    ListQuads = Tree2List(pS->Quadrangles);
-    List_Sort(ListQuads, _p3d_cmp_entities);
-    _p3d_print_quads(ListQuads, pS->Nu, pS->Nv);
-    List_Delete(ListQuads);
-  }
-
-  for(i = 0; i < List_Nbr(ListStructVol); i++) {
-    List_Read(ListStructVol, i, &pV);
-    List_Read(pV->Surfaces, 0, &pS);
-    Nu = pS->Nu;
-    Nv = pS->Nv;
-    List_Read(pV->Surfaces, 1, &pS);
-    ListHex = Tree2List(pV->Hexahedra);
-    List_Sort(ListHex, _p3d_cmp_entities);
-    _p3d_print_hex(ListHex, Nu, Nv, pS->Nv);
-    List_Delete(ListHex);
-  }
-
-  List_Delete(ListSurfaces);
-  List_Delete(ListVolumes);
-  List_Delete(ListStructSurf);
-  List_Delete(ListStructVol);
-}
-
-void _p3d_print_elements()
-{
-  // FIXME: to do
-}
-
-void Print_Mesh_P3D(FILE *fp)
-{
-  P3DFILE = fp;
-
-  if(!List_Nbr(THEM->PhysicalGroups) || CTX.mesh.save_all) {
-    Msg(INFO, "Saving all elements (discarding physical groups)");
-    _p3d_print_all_elements();
-  }
-  else{
-    Msg(WARNING, "No Plot3d format of physicals yet, saving all elements!");
-    //_p3d_print_elements();
-    _p3d_print_all_elements();
-  }
-}
-
-// Public Print_Mesh routine
-
-void GetDefaultMeshFileName(int Type, char *name)
-{
-  char ext[10] = "";
-  strcpy(name, THEM->name);
-  switch(Type){
-  case FORMAT_MSH:  strcpy(ext, ".msh"); break;
-  case FORMAT_VRML: strcpy(ext, ".wrl"); break;
-  case FORMAT_UNV:  strcpy(ext, ".unv"); break;
-  case FORMAT_GREF: strcpy(ext, ".Gref"); break;
-  case FORMAT_DMG:  strcpy(ext, ".dmg"); break;
-  case FORMAT_STL:  strcpy(ext, ".stl"); break;
-  case FORMAT_P3D:  strcpy(ext, ".p3d"); break;
-  default: Msg(GERROR, "Unknown mesh file format %d", Type); break;
-  }
-  strcat(name, ext);
-}
-
-void Print_Mesh(char *filename, int Type)
-{
-  char name[256];
-
-  if(CTX.threads_lock) {
-    Msg(INFO, "I'm busy! Ask me that later...");
-    return;
-  }
-
-  CTX.threads_lock = 1;
-
-  if(!filename)
-    GetDefaultMeshFileName(Type, name);
-  else
-    strcpy(name, filename);
-
-  FILE *fp = fopen(name, "w");
-  if(!fp) {
-    Msg(GERROR, "Unable to open file '%s'", name);
-    CTX.threads_lock = 0;
-    return;
-  }
-
-  switch(Type){
-  case FORMAT_MSH:  Print_Mesh_MSH(fp); break;
-  case FORMAT_VRML: Print_Mesh_WRL(fp); break;
-  case FORMAT_UNV:  Print_Mesh_UNV(fp); break;
-  case FORMAT_GREF: Print_Mesh_GREF(fp); break;
-  case FORMAT_DMG:  Print_Mesh_DMG(fp); break;
-  case FORMAT_STL:  Print_Mesh_STL(fp); break;
-  case FORMAT_P3D:  Print_Mesh_P3D(fp); break;
-  default:
-    break;
-  }
-
-  fclose(fp);
-
-  CTX.threads_lock = 0;
-}
diff --git a/Mesh/Read_Mesh.cpp b/Mesh/Read_Mesh.cpp
deleted file mode 100644
index 8892e98254a60c4a49d93b4a435fc446baa26806..0000000000000000000000000000000000000000
--- a/Mesh/Read_Mesh.cpp
+++ /dev/null
@@ -1,1070 +0,0 @@
-// $Id: Read_Mesh.cpp,v 1.106 2006-07-24 14:05:50 geuzaine Exp $
-//
-// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
-//
-// This program is free software; you can redistribute it and/or modify
-// it under the terms of the GNU General Public License as published by
-// the Free Software Foundation; either version 2 of the License, or
-// (at your option) any later version.
-//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with this program; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
-// USA.
-// 
-// Please report all bugs and problems to <gmsh@geuz.org>.
-
-#include <map>
-#include <vector>
-#include "Gmsh.h"
-#include "Geo.h"
-#include "CAD.h"
-#include "Mesh.h"
-#include "3D_Mesh.h"
-#include "Create.h"
-#include "Numeric.h"
-#include "Context.h"
-#include "PartitionMesh.h"
-
-extern Context_T CTX;
-
-// Read mesh in native MSH format
-
-#define LGN1 1
-#define TRI1 2
-#define QUA1 3
-#define TET1 4
-#define HEX1 5
-#define PRI1 6
-#define PYR1 7
-#define LGN2 8
-#define TRI2 9
-#define QUA2 10
-#define TET2 11
-#define HEX2 12
-#define PRI2 13
-#define PYR2 14
-#define PNT  15
-
-#define NB_NOD_MAX_ELM 30
-
-// If a "normal" elementary entity does not exist, we create a
-// "discrete" entity, i.e., an entity entirely defined by its mesh
-
-Curve *addElementaryCurve(Mesh * M, int Num)
-{
-  Curve *c;
-  if(!(c = FindCurve(Num, M))) {
-    c = Create_Curve(Num, MSH_SEGM_DISCRETE, 0, NULL, NULL, -1, -1, 0., 1.);
-    Tree_Add(M->Curves, &c);
-    CreateReversedCurve(M, c);
-  }
-  return c;
-}
-
-class nodex{
- public:
-  Vertex *v;
-  int left, right; // index left/right segment in simplex list
-  nodex(Vertex *vv) : v(vv), left(-1), right(-1) {}
-};
-
-void addVerticesInCurve(void *a, void *b)
-{
-  Curve *c = *(Curve**)a;
-
-  // check if already done or if nothing to do
-  if(List_Nbr(c->Vertices) || !Tree_Nbr(c->SimplexesBase)) 
-    return;
-
-  List_T *elm = Tree2List(c->SimplexesBase);
-
-  // construct node x segment connectivity
-  std::map<int, nodex*> nodes;
-  for(int i = 0; i < List_Nbr(elm); i++){
-    SimplexBase *s = *(SimplexBase**)List_Pointer(elm, i);
-    for(int j = 0; j < 2; j++){
-      if(!nodes.count(s->V[j]->Num))
-	nodes[s->V[j]->Num] = new nodex(s->V[j]);
-      if(!j)
-	nodes[s->V[j]->Num]->right = i;
-      else
-	nodes[s->V[j]->Num]->left = i;
-    }
-  }
-
-  // find starting element (or use the first one for a closed curve)
-  std::map<int, nodex*>::const_iterator beg = nodes.begin();
-  std::map<int, nodex*>::const_iterator end = nodes.end();
-  int start = 0;
-  while(beg != end){
-    if((*beg).second->left < 0){ // no element to the left
-      start = (*beg).second->right;
-      break;
-    }
-    beg++;
-  }
-
-  // add nodes by following the line segments
-  if(start < 0){
-    Msg(GERROR, "Something is wrong in mesh of curve %d", c->Num);
-  }
-  else{
-    int N = List_Nbr(elm);
-    SimplexBase *s = *(SimplexBase**)List_Pointer(elm, start);
-    List_Add(c->Vertices, &s->V[0]);
-    while(N > 0){
-      if(s->VSUP) List_Add(c->Vertices, &s->VSUP[0]);
-      List_Add(c->Vertices, &s->V[1]);
-      nodex *n = nodes[s->V[1]->Num];
-      if(n->left != start){
-	Msg(GERROR, "Wrong orientation of element %d in curve %d", s->Num, c->Num);
-	break;
-      }
-      start = n->right;
-      if(start >= 0){
-	s = *(SimplexBase**)List_Pointer(elm, start);
-      }
-      else if(N != 1){
-	Msg(GERROR, "Something is wrong in mesh of curve %d", c->Num);
-	break;
-      }
-      N--;
-    }
-  }
-  
-  List_Delete(elm);
-}
-
-Surface *addElementarySurface(Mesh * M, int Num)
-{
-  Surface *s;
-  if(!(s = FindSurface(Num, M))) {
-    s = Create_Surface(Num, MSH_SURF_DISCRETE);
-    Tree_Add(M->Surfaces, &s);
-  }
-  return s;
-}
-
-Volume *addElementaryVolume(Mesh * M, int Num)
-{
-  Volume *v;
-  if(!(v = FindVolume(Num, M))) {
-    v = Create_Volume(Num, MSH_VOLUME_DISCRETE);
-    Tree_Add(M->Volumes, &v);
-  }
-  return v;
-}
-
-void addPhysicalGroup(Tree_T * groups, int Type, int Physical, int Elementary)
-{
-  // we add in a temporary group tree for performance reasons; the
-  // tree is converted back into a list in the mesh at the end of read_mesh
-  PhysicalGroup g, *pg;
-  pg = &g;
-  pg->Num = Physical;
-  pg->Typ = Type;
-  if(Tree_Query(groups, &pg)) {
-    List_Insert(pg->Entities, &Elementary, fcmp_int);
-  }
-  else {
-    List_T *tmp = List_Create(1, 1, sizeof(int));
-    List_Add(tmp, &Elementary);
-    pg = Create_PhysicalGroup(Physical, Type, tmp);
-    Tree_Add(groups, &pg);
-    List_Delete(tmp);
-  }
-}
-
-int addMeshPartition(int Num, Mesh * M)
-{
-  MeshPartition P, *p, **pp;
-  p = &P;
-  p->Num = Num;
-  if((pp = (MeshPartition**)List_PQuery(M->Partitions, &p, compareMeshPartitionNum))){
-    return (*pp)->Index;
-  }
-  else{
-    p = (MeshPartition*)Malloc(sizeof(MeshPartition));
-    p->Num = Num;
-    p->Visible = VIS_GEOM | VIS_MESH;
-    p->Index = List_Nbr(M->Partitions);
-    List_Add(M->Partitions, &p);
-    return p->Index;
-  }
-  return 0;
-}
-
-int getNbrNodes(int Type)
-{
-  switch (Type) {
-  case PNT : return 1;
-  case LGN1: return 2;
-  case LGN2: return 2 + 1;
-  case TRI1: return 3;
-  case TRI2: return 3 + 3;
-  case QUA1: return 4;
-  case QUA2: return 4 + 4 + 1;
-  case TET1: return 4;
-  case TET2: return 4 + 6;
-  case HEX1: return 8;
-  case HEX2: return 8 + 12 + 6 + 1;
-  case PRI1: return 6;
-  case PRI2: return 6 + 9 + 3;
-  case PYR1: return 5;
-  case PYR2: return 5 + 8 + 1;
-  default: return 0;
-  }
-}
-
-double SetLC(Vertex *v1, Vertex *v2, Vertex *v3 = 0, Vertex *v4 = 0)
-{ 
-  double lc1 = sqrt((v1->Pos.X - v2->Pos.X) * (v1->Pos.X - v2->Pos.X) +
-		    (v1->Pos.Y - v2->Pos.Y) * (v1->Pos.Y - v2->Pos.Y) +
-		    (v1->Pos.Z - v2->Pos.Z) * (v1->Pos.Z - v2->Pos.Z));
-  if(!v3){
-    double lc = lc1 * CTX.mesh.lc_factor;
-    v1->lc = v2->lc = lc;
-    return lc;
-  }
-  else{
-    double lc2 = sqrt((v1->Pos.X - v3->Pos.X) * (v1->Pos.X - v3->Pos.X) +
-		      (v1->Pos.Y - v3->Pos.Y) * (v1->Pos.Y - v3->Pos.Y) +
-		      (v1->Pos.Z - v3->Pos.Z) * (v1->Pos.Z - v3->Pos.Z));
-    double lc3 = sqrt((v2->Pos.X - v3->Pos.X) * (v2->Pos.X - v3->Pos.X) +
-		      (v2->Pos.Y - v3->Pos.Y) * (v2->Pos.Y - v3->Pos.Y) +
-		      (v2->Pos.Z - v3->Pos.Z) * (v2->Pos.Z - v3->Pos.Z));
-    double lc = DMAX(lc1, DMAX(lc2, lc3)) * CTX.mesh.lc_factor;
-    v1->lc = v2->lc = v3->lc = lc;
-    if(v4) v4->lc = lc;
-    return lc;
-  }
-}
-
-int getPartition(const std::multimap<int,int> &nod2proc, int nbNod, Vertex *verts)
-{
-  std::map<int,int> proc2nbnod;
-  for(int i = 0; i < nbNod; i++){
-    std::multimap<int,int>::const_iterator beg = nod2proc.lower_bound(verts[i].Num);      
-    std::multimap<int,int>::const_iterator end = nod2proc.upper_bound(verts[i].Num);
-    while(beg != end){
-      int iProc = (*beg).second;	  
-      if (proc2nbnod.find(iProc) == proc2nbnod.end()) proc2nbnod[iProc] = 1;
-      else proc2nbnod[iProc] = proc2nbnod[iProc]+1;
-      ++beg;
-    }
-  }
-  {
-    std::map<int,int>::const_iterator beg = proc2nbnod.begin();      
-    std::map<int,int>::const_iterator end = proc2nbnod.end();
-    while(beg!=end){
-      if((*beg).second == nbNod) return (*beg).first;
-      beg++;
-    }    
-  }
-  return 0;  
-}
-
-void Read_Mesh_MSH(Mesh * M, FILE * fp)
-{
-  char String[256];
-  double version = 1.0;
-  int format = LIST_FORMAT_ASCII, size = sizeof(double);
-  int Nbr_Nodes, Nbr_Elements, i_Node, i_Element;
-  int Num, Type, Physical, Elementary, Partition, i, j;
-  int NbTags, Tag;
-  double x, y, z, lc1, lc2;
-  Vertex *vert, verts[NB_NOD_MAX_ELM], *vertsp[NB_NOD_MAX_ELM], **vertspp;
-  SimplexBase *simp;
-  Quadrangle *quad;
-  Hexahedron *hex;
-  Prism *pri;
-  Pyramid *pyr;
-  Curve *c;
-  Surface *s;
-  Volume *v;
-  Tree_T *Duplicates = NULL;
-  Tree_T *groups = List2Tree(M->PhysicalGroups, comparePhysicalGroup);
-  std::multimap<int,int> nod2proc;
-
-  while(1) {
-    do {
-      if(!fgets(String, sizeof(String), fp))
-	break;
-      if(feof(fp))
-        break;
-    } while(String[0] != '$');
-
-    if(feof(fp))
-      break;
-
-    /*  F o r m a t  */
-
-    if(!strncmp(&String[1], "MeshFormat", 10)) {
-      fscanf(fp, "%lf %d %d\n", &version, &format, &size);
-
-      if(version == 2.0){
-	Msg(INFO, "Detected mesh file format %g", version);
-      }
-      else{
-	Msg(GERROR, "Unknown MSH file version to read (%g)", version);
-	return;
-      }
-
-      if(format == 0)
-        format = LIST_FORMAT_ASCII;
-      else if(format == 1)
-        format = LIST_FORMAT_BINARY;
-      else {
-        Msg(GERROR, "Unknown data format for mesh");
-        return;
-      }
-    }
-
-    /*  POINTS -- this field is deprecated, and will eventually disappear */
-
-    else if(!strncmp(&String[1], "PTS", 3) ||
-	    !strncmp(&String[1], "Points", 6)) {
-
-      fscanf(fp, "%d", &Nbr_Nodes);
-      Msg(INFO, "%d points", Nbr_Nodes);
-
-      for(i_Node = 0; i_Node < Nbr_Nodes; i_Node++) {
-        fscanf(fp, "%d %lf %lf %lf %lf %lf", &Num, &x, &y, &z, &lc1,
-               &lc2);
-        vert = Create_Vertex(Num, x, y, z, lc1, lc2);
-        if(!Tree_Insert(M->Points, &vert)){
-	  Msg(GERROR, "Point %d already exists", vert->Num);
-	  Free_Vertex(&vert, 0);
-	}
-      }
-    }
-
-    /*  NODE'S PROCESSORS */
-
-    else if(!strncmp(&String[1], "PARA", 4)){
-      fscanf(fp, "%d", &Nbr_Nodes);
-      Msg(INFO, "%d parallel nodes", Nbr_Nodes);
-      for(i_Node = 0; i_Node < Nbr_Nodes; i_Node++) {
-	int nbProc;
-	fscanf(fp, "%d %d", &Num, &nbProc);
-	for (int iProc = 0; iProc < nbProc; iProc++){
-	  int iProcNum;
-	  fscanf(fp, "%d", &iProcNum);
-	  nod2proc.insert(std::pair<int, int>(Num, iProcNum));
-	}
-      }
-    }
-    
-    /*  NODES  */
-
-    else if(!strncmp(&String[1], "NOD", 3) ||
-	    !strncmp(&String[1], "NOE", 3) ||
-	    !strncmp(&String[1], "Nodes", 5)) {
-
-      fscanf(fp, "%d", &Nbr_Nodes);
-      Msg(INFO, "%d nodes", Nbr_Nodes);
-
-      if(CTX.mesh.check_duplicates)
-        Duplicates = Tree_Create(sizeof(Vertex *), comparePosition);
-
-      int NN = (Nbr_Nodes > 100000) ? Nbr_Nodes/50 : Nbr_Nodes/10;
-
-      for(i_Node = 0; i_Node < Nbr_Nodes; i_Node++) {
-        fscanf(fp, "%d %lf %lf %lf", &Num, &x, &y, &z);
-        vert = Create_Vertex(Num, x, y, z, 1.0, 0.0);
-        if(!Tree_Insert(M->Vertices, &vert)){
-	  Msg(GERROR, "Node %d already exists", vert->Num);
-	  Free_Vertex(&vert, 0);
-	}
-	else if(CTX.mesh.check_duplicates) {
-          if((vertspp = (Vertex **) Tree_PQuery(Duplicates, &vert)))
-            Msg(GERROR, "Nodes %d and %d have identical coordinates (%g, %g, %g)",
-                Num, (*vertspp)->Num, x, y, z);
-          else
-            Tree_Add(Duplicates, &vert);
-        }
-	if(NN && (i_Node % NN == NN - 1))
-	  Msg(PROGRESS, "Read %d nodes", i_Node + 1);
-      }
-      Msg(PROGRESS, "");
-
-      if(CTX.mesh.check_duplicates)
-        Tree_Delete(Duplicates);
-    }
-
-    /* ELEMENTS */
-
-    else if(!strncmp(&String[1], "ELM", 3) ||
-	    !strncmp(&String[1], "Elements", 8)) {
-
-      fscanf(fp, "%d", &Nbr_Elements);
-      Msg(INFO, "%d elements", Nbr_Elements);
-
-      if(CTX.mesh.check_duplicates)
-        Duplicates = Tree_Create(sizeof(Vertex *), comparePosition);
-
-      int NN = (Nbr_Elements > 100000) ? Nbr_Elements/50 : Nbr_Elements/10;
-
-      for(i_Element = 0; i_Element < Nbr_Elements; i_Element++) {
-	
-	if(version <= 1.0){
-	  fscanf(fp, "%d %d %d %d %d",
-		 &Num, &Type, &Physical, &Elementary, &Nbr_Nodes);
-	  Partition = 1;
-	  int Nbr_Nodes_Check = getNbrNodes(Type);
-	  if(!Nbr_Nodes_Check){
-	    Msg(GERROR, "Unknown type for element %d", Num); 
-	    return;
-	  }
-	  if(Nbr_Nodes != Nbr_Nodes_Check){
-	    Msg(GERROR, "Wrong number of nodes for element %d", Num);
-	    return;
-	  }
-	}
-	else{
-	  fscanf(fp, "%d %d %d", &Num, &Type, &NbTags);
-	  Elementary = Physical = Partition = 1;
-	  for(j = 0; j < NbTags; j++){
-	    fscanf(fp, "%d", &Tag);	    
-	    if(j == 0)
-	      Physical = Tag;
-	    else if(j == 1)
-	      Elementary = Tag;
-	    else if(j == 2)
-	      Partition = Tag;
-	    // ignore any other tags for now
-	  }
-	  Nbr_Nodes = getNbrNodes(Type);
-	  if(!Nbr_Nodes){
-	    Msg(GERROR, "Unknown type for element %d", Num); 
-	    return;
-	  }
-	}
-
-        for(j = 0; j < Nbr_Nodes; j++){
-	  fscanf(fp, "%d", &verts[j].Num);
-	}
-	
-	if (nod2proc.size()){
-	  Partition = getPartition ( nod2proc ,Nbr_Nodes , verts );
-	}
-	
-        for(i = 0; i < Nbr_Nodes; i++) {
-          vertsp[i] = &verts[i];
-          if(!(vertspp = (Vertex **) Tree_PQuery(M->Vertices, &vertsp[i])))
-            Msg(GERROR, "Unknown vertex %d in element %d", verts[i].Num, Num);
-          else
-            vertsp[i] = *vertspp;
-        }
-
-        if(CTX.mesh.check_duplicates) {
-          vert = Create_Vertex(Num, 0., 0., 0., 1.0, 0.0);
-          for(i = 0; i < Nbr_Nodes; i++) {
-            vert->Pos.X += vertsp[i]->Pos.X;
-            vert->Pos.Y += vertsp[i]->Pos.Y;
-            vert->Pos.Z += vertsp[i]->Pos.Z;
-          }
-          vert->Pos.X /= (double)Nbr_Nodes;
-          vert->Pos.Y /= (double)Nbr_Nodes;
-          vert->Pos.Z /= (double)Nbr_Nodes;
-          if((vertspp = (Vertex **) Tree_PQuery(Duplicates, &vert)))
-            Msg(WARNING, "Elements %d and %d have identical barycenters",
-                Num, (*vertspp)->Num);
-          else
-            Tree_Add(Duplicates, &vert);
-        }
-
-        switch (Type) {
-        case LGN1:
-        case LGN2:
-	  c = addElementaryCurve(M, abs(Elementary));
-	  addPhysicalGroup(groups, MSH_PHYSICAL_LINE, Physical, abs(Elementary));
-          simp = Create_SimplexBase(vertsp[0], vertsp[1], NULL, NULL);
-          simp->Num = Num;
-          simp->iEnt = Elementary;
-          simp->iPart = addMeshPartition(Partition, M);
-	  SetLC(vertsp[0], vertsp[1]);
-	  if(Type == LGN2){
-	    simp->VSUP = (Vertex **) Malloc(1 * sizeof(Vertex *));
-	    simp->VSUP[0] = vertsp[2];
-	    simp->VSUP[0]->Degree = 2;
-	  }
-          if(!Tree_Insert(c->SimplexesBase, &simp)){
-	    Msg(GERROR, "Line element %d already exists", simp->Num);
-	    Free_SimplexBase(&simp, 0);
-	  }
-	  // we don't insert the vertices in the list of vertices at
-	  // this point, since we need the list to be ordered (and
-	  // consistent!), and simply doing a List_Insert and sorting
-	  // according to node numbers is not enough (in addition to
-	  // being slow): see addVerticesInCurve() below.
-          break;
-        case TRI1:
-        case TRI2:
-	  s = addElementarySurface(M, Elementary);
-	  addPhysicalGroup(groups, MSH_PHYSICAL_SURFACE, Physical, Elementary);
-          simp = Create_SimplexBase(vertsp[0], vertsp[1], vertsp[2], NULL);
-          simp->Num = Num;
-          simp->iEnt = Elementary;
-          simp->iPart = addMeshPartition(Partition, M);
-	  SetLC(vertsp[0], vertsp[1], vertsp[2]);
-	  if(Type == TRI2){
-	    simp->VSUP = (Vertex **) Malloc(3 * sizeof(Vertex *));
-	    for(i = 0; i < 3; i++){
-	      simp->VSUP[i] = vertsp[i+3];
-	      simp->VSUP[i]->Degree = 2;
-	    }
-	  }
-          if(!Tree_Insert(s->SimplexesBase, &simp)){
-	    Msg(GERROR, "Triangle %d already exists", simp->Num);
-	    Free_SimplexBase(&simp, 0);
-	  }
-	  else{
-	    for(i = 0; i < Nbr_Nodes; i++)
-	      Tree_Insert(s->Vertices, &vertsp[i]);
-	  }
-          break;
-        case QUA1:
-        case QUA2:
-	  s = addElementarySurface(M, Elementary);
-	  addPhysicalGroup(groups, MSH_PHYSICAL_SURFACE, Physical, Elementary);
-          quad = Create_Quadrangle(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
-          quad->Num = Num;
-          quad->iEnt = Elementary;
-          quad->iPart = addMeshPartition(Partition, M);
-	  SetLC(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
-	  if(Type == QUA2){
-	    quad->VSUP = (Vertex **) Malloc((4 + 1) * sizeof(Vertex *));
-	    for(i = 0; i < 4 + 1; i++){
-	      quad->VSUP[i] = vertsp[i+4];
-	      quad->VSUP[i]->Degree = 2;
-	    }
-	  }
-          if(!Tree_Insert(s->Quadrangles, &quad)){
-	    Msg(GERROR, "Quadrangle %d already exists", quad->Num);
-	    Free_SimplexBase(&quad, 0);
-	  }
-	  else{
-	    for(i = 0; i < Nbr_Nodes; i++)
-	      Tree_Insert(s->Vertices, &vertsp[i]);
-	  }
-          break;
-        case TET1:
-        case TET2:
-	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
-          simp = Create_SimplexBase(vertsp[0], vertsp[1], vertsp[2], vertsp[3]);
-          simp->Num = Num;
-          simp->iEnt = Elementary;
-          simp->iPart = addMeshPartition(Partition, M);
-	  if(Type == TET2){
-	    simp->VSUP = (Vertex **) Malloc(6 * sizeof(Vertex *));
-	    for(i = 0; i < 6; i++){
-	      simp->VSUP[i] = vertsp[i+4];
-	      simp->VSUP[i]->Degree = 2;
-	    }
-	  }
-          if(!Tree_Insert(v->SimplexesBase, &simp)){
-	    Msg(GERROR, "Tetrahedron %d already exists", simp->Num);
-	    Free_SimplexBase(&simp, 0);
-	  }
-#if 0 // removed to speed things up (not used at the moment anyway)
-	  else{
-	    for(i = 0; i < Nbr_Nodes; i++)
-	      Tree_Insert(v->Vertices, &vertsp[i]);
-	  }
-#endif
-          break;
-        case HEX1:
-        case HEX2:
-	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
-          hex = Create_Hexahedron(vertsp[0], vertsp[1], vertsp[2], vertsp[3],
-                                  vertsp[4], vertsp[5], vertsp[6], vertsp[7]);
-          hex->Num = Num;
-          hex->iEnt = Elementary;
-          hex->iPart = addMeshPartition(Partition, M);
-	  if(Type == HEX2){
-	    hex->VSUP = (Vertex **) Malloc((12 + 6 + 1) * sizeof(Vertex *));
-	    for(i = 0; i < 12 + 6 + 1; i++){
-	      hex->VSUP[i] = vertsp[i+8];
-	      hex->VSUP[i]->Degree = 2;
-	    }
-	  }
-          if(!Tree_Insert(v->Hexahedra, &hex)){
-	    Msg(GERROR, "Hexahedron %d already exists", hex->Num);
-	    Free_Hexahedron(&hex, 0);
-	  }
-#if 0 // removed to speed things up (not used at the moment anyway)
-	  else{
-	    for(i = 0; i < Nbr_Nodes; i++)
-	      Tree_Insert(v->Vertices, &vertsp[i]);
-	  }
-#endif
-          break;
-        case PRI1:
-        case PRI2:
-	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
-          pri = Create_Prism(vertsp[0], vertsp[1], vertsp[2],
-                             vertsp[3], vertsp[4], vertsp[5]);
-          pri->Num = Num;
-          pri->iEnt = Elementary;
-          pri->iPart = addMeshPartition(Partition, M);
-	  if(Type == PRI2){
-	    pri->VSUP = (Vertex **) Malloc((9 + 3) * sizeof(Vertex *));
-	    for(i = 0; i < 9 + 3; i++){
-	      pri->VSUP[i] = vertsp[i+6];
-	      pri->VSUP[i]->Degree = 2;
-	    }
-	  }
-          if(!Tree_Insert(v->Prisms, &pri)){
-	    Msg(GERROR, "Prism %d already exists", pri->Num);
-	    Free_Prism(&pri, 0);
-	  }
-#if 0 // removed to speed things up (not used at the moment anyway)
-	  else{
-	    for(i = 0; i < Nbr_Nodes; i++)
-	      Tree_Insert(v->Vertices, &vertsp[i]);
-	  }
-#endif
-          break;
-        case PYR1:
-        case PYR2:
-	  v = addElementaryVolume(M, Elementary);
-	  addPhysicalGroup(groups, MSH_PHYSICAL_VOLUME, Physical, Elementary);
-          pyr = Create_Pyramid(vertsp[0], vertsp[1], vertsp[2],
-                               vertsp[3], vertsp[4]);
-          pyr->Num = Num;
-          pyr->iEnt = Elementary;
-          pyr->iPart = addMeshPartition(Partition, M);
-	  if(Type == PYR2){
-	    pyr->VSUP = (Vertex **) Malloc((8 + 1) * sizeof(Vertex *));
-	    for(i = 0; i < 8 + 1; i++){
-	      pyr->VSUP[i] = vertsp[i+5];
-	      pyr->VSUP[i]->Degree = 2;
-	    }
-	  }
-          if(!Tree_Insert(v->Pyramids, &pyr)){
-	    Msg(GERROR, "Pyramid %d already exists", pyr->Num);
-	    Free_Pyramid(&pyr, 0);
-	  }
-#if 0 // removed to speed things up (not used at the moment anyway)
-	  else{
-	    for(i = 0; i < Nbr_Nodes; i++)
-	      Tree_Insert(v->Vertices, &vertsp[i]);
-	  }
-#endif
-          break;
-        case PNT:
-	  addPhysicalGroup(groups, MSH_PHYSICAL_POINT, Physical, Elementary);
-	  // we need to make a new one: vertices in M->Vertices and
-	  // M->Points should never point to the same memory location
-	  vert = Create_Vertex(vertsp[0]->Num, vertsp[0]->Pos.X, vertsp[0]->Pos.Y, 
-			       vertsp[0]->Pos.Z, vertsp[0]->lc, vertsp[0]->w);
-          if(!Tree_Insert(M->Points, &vert)){
-	    Msg(WARNING, "Point %d already exists", vert->Num);
-	    Free_Vertex(&vert, 0);
-	  }
-          break;
-        default:
-	  Msg(GERROR, "Unknown type for element %d", Num); 
-          break;
-        }
-
-	if(NN && (i_Element % NN == NN - 1))
-	  Msg(PROGRESS, "Read %d elements", i_Element + 1);
-      }
-      Msg(PROGRESS, "");
-
-      if(CTX.mesh.check_duplicates) {
-        Tree_Action(Duplicates, Free_Vertex);
-        Tree_Delete(Duplicates);
-      }
-
-    }
-
-    do {
-      if(!fgets(String, sizeof(String), fp))
-	Msg(GERROR, "Prematured end of mesh file");
-      if(feof(fp))
-        Msg(GERROR, "Prematured end of mesh file");
-    } while(String[0] != '$');
-
-  }
-
-  if(Tree_Nbr(M->Volumes))
-    M->status = 3;
-  else if(Tree_Nbr(M->Surfaces))
-    M->status = 2;
-  else if(Tree_Nbr(M->Curves))
-    M->status = 1;
-  else if(Tree_Nbr(M->Points))
-    M->status = 0;
-  else
-    M->status = -1;
-
-  // add vertices in curves
-  Tree_Action(M->Curves, addVerticesInCurve);
-
-  // For efficiency reasons, we store the partition index (and not the
-  // partition number) in the various mesh elements. We need to
-  // re-sort the list according to these indices to allow direct
-  // access through List_Pointer & co.
-  List_Sort(M->Partitions, compareMeshPartitionIndex);
-
-  // Transfer the temp group tree back into the mesh
-  List_Delete(M->PhysicalGroups);
-  M->PhysicalGroups = Tree2List(groups);
-  Tree_Delete(groups);
-}
-
-// Read mesh in VTK format
-
-void Read_Mesh_VTK(Mesh * m, FILE *fp)
-{
-  char line[256], dumline1[256], dumline2[256];
-  int i;
-  int NbFaces, NbVertices, Vertex1, Vertex2, Vertex3, NbVerticesOnFace;
-  double x, y, z;
-  Vertex *v1, *v2, *v3;
-
-  fgets(line, 255, fp);
-  fgets(line, 255, fp);
-  fgets(line, 255, fp);
-  fgets(line, 255, fp);
-  fgets(line, 255, fp);
-
-  sscanf(line, "%s %d %s", dumline1, &NbVertices, dumline2);
-  Surface *surf = Create_Surface(1, MSH_SURF_DISCRETE);
-  Tree_Add(m->Surfaces, &surf);
-  for(i = 0; i < NbVertices; i++) {
-    fscanf(fp, "%le %le %le", &x, &y, &z);
-    Vertex *v = Create_Vertex(i, x, y, z, 1.0, 1.0);
-    Tree_Add(m->Vertices, &v);
-    Tree_Add(surf->Vertices, &v);
-    v->ListSurf = List_Create(1, 1, sizeof(Surface *));
-    List_Add(v->ListSurf, &surf);
-  }
-
-  fscanf(fp, "%s %d %d", dumline1, &NbFaces, &i);
-  for(int i = 0; i < NbFaces; i++) {
-    fscanf(fp, "%d", &NbVerticesOnFace);
-    Simplex *s;
-    if(NbVerticesOnFace == 3) {
-      fscanf(fp, "%d %d %d", &Vertex1, &Vertex2, &Vertex3);
-      v1 = FindVertex(Vertex1, m);
-      v2 = FindVertex(Vertex2, m);
-      v3 = FindVertex(Vertex3, m);
-      if(!v1 || !v2 || !v3){
-	Msg(GERROR, "Bad vertex reference in VTK file: aborting");
-	return;
-      }
-      else{
-	s = Create_Simplex(v1, v2, v3, NULL);
-	s->Num = i;
-	s->iEnt = 1;
-      }
-    }
-    if(!(surf = FindSurface(1, m))) {
-      surf = Create_Surface(1, MSH_SURF_DISCRETE);
-      Tree_Add(m->Surfaces, &surf);
-    }
-    Tree_Add(surf->Simplexes, &s);
-  }
-
-  if(NbFaces)
-    m->status = 2;
-
-  Volume *vol = Create_Volume(1, MSH_VOLUME_DISCRETE);
-  vol->Surfaces = List_Create(1, 1, sizeof(Surface *));
-  List_Add(vol->Surfaces, &surf);
-  Tree_Add(m->Volumes, &vol);
-}
-
-// Read mesh in SMS format
-
-#define ENTITY_VERTEX 0
-#define ENTITY_EDGE   1
-#define ENTITY_FACE   2
-#define ENTITY_REGION 3
-#define ENTITY_NONE   4
-
-void Read_Mesh_SMS(Mesh * m, FILE * in)
-{
-  char line[1023];
-  int i, patch, nbPts;
-  int NbRegions, NbFaces, NbEdges, NbVertices, NbPoints,
-    GEntityType, GEntityId, EntityNbConnections, Dummy,
-    Edge1, Edge2, Edge3, Edge4, Face1, Face2, Face3, Face4;
-  int VertexId1, VertexId2, NbEdgesOnFace, NbFacesOnRegion;
-  double x, y, z, u, v;
-  List_T *AllEdges, *AllFaces;
-  Vertex *v1 = NULL, *v2 = NULL, *v3 = NULL, *v4 = NULL;
-
-  fscanf(in, "%s %d", line, &Dummy);
-  fscanf(in, "%d %d %d %d %d", &NbRegions, &NbFaces, &NbEdges, &NbVertices,
-         &NbPoints);
-
-  Msg(INFO, "Reading a mesh in scorec format");
-  Msg(INFO, "%d Vertices", NbVertices);
-
-  for(i = 0; i < NbVertices; i++) {
-    fscanf(in, "%d", &GEntityId);
-    if(GEntityId) {
-      fscanf(in, "%d %d %lf %lf %lf", &GEntityType, &EntityNbConnections, &x,
-             &y, &z);
-      Vertex *vert = Create_Vertex(i, x, y, z, 1.0, 1.0);
-      Tree_Add(m->Vertices, &vert);
-      switch (GEntityType) {
-      case 0:
-	{
-	  // we need to make a new one: vertices in m->Vertices and
-	  // m->Points should never point to the same memory location
-	  Vertex *pnt = Create_Vertex(i, x, y, z, 1.0, 1.0);
-	  Tree_Add(m->Points, &pnt);
-	}
-        break;
-      case 1:
-        fscanf(in, "%le", &u);
-        break;
-      case 2:
-        fscanf(in, "%le %le %d", &u, &v, &patch);
-        break;
-      case 3:
-        break;
-      }
-    }
-  }
-
-  Msg(INFO, "%d edges", NbEdges);
-  AllEdges = List_Create(NbEdges, 1, sizeof(Edge));
-  Edge e;
-
-  for(int i = 0; i < NbEdges; i++) {
-    fscanf(in, "%d", &GEntityId);
-
-    if(GEntityId) {
-      fscanf(in, "%d %d %d %d %d", &GEntityType, &VertexId1, &VertexId2,
-             &EntityNbConnections, &nbPts);
-      for(int j = 0; j < nbPts; j++) {
-        switch (GEntityType) {
-        case 0:
-          break;
-        case 1:
-          fscanf(in, "%le", &u);
-          break;
-        case 2:
-          fscanf(in, "%le %le %d", &u, &v, &patch);
-          break;
-        case 3:
-          break;
-        }
-      }
-      e.Points = NULL;
-      Vertex *v1 = FindVertex(VertexId1 - 1, m);
-      Vertex *v2 = FindVertex(VertexId2 - 1, m);
-      e.V[0] = v1;
-      e.V[1] = v2;
-      List_Add(AllEdges, &e);
-      switch (GEntityType) {
-      case ENTITY_EDGE:
-        Simplex * s = Create_Simplex(v1, v2, NULL, NULL);
-        Curve *c;
-        if(!(c = FindCurve(GEntityId, m))) {
-          c = Create_Curve(GEntityId, MSH_SEGM_DISCRETE, 1, NULL, NULL, -1, -1, 0, 1);
-          Tree_Add(m->Curves, &c);
-        }
-        s->iEnt = GEntityId;
-        s->Num = i;
-        Tree_Add(c->Simplexes, &s);
-      }
-    }
-  }
-
-  AllFaces = List_Create(NbFaces, 1, sizeof(Simplex *));
-
-  Volume *vol = Create_Volume(1, MSH_VOLUME_DISCRETE);
-  vol->Surfaces = List_Create(1, 1, sizeof(Surface *));
-  Tree_Add(m->Volumes, &vol);
-
-  Msg(INFO, "%d faces", NbFaces);
-  for(int i = 0; i < NbFaces; i++) {
-    fscanf(in, "%d", &GEntityId);
-    if(GEntityId) {
-      fscanf(in, "%d %d", &GEntityType, &NbEdgesOnFace);
-
-      if(NbEdgesOnFace == 3) {
-        fscanf(in, "%d %d %d %d", &Edge1, &Edge2, &Edge3, &nbPts);
-        List_Read(AllEdges, abs(Edge1) - 1, &e);
-        if(Edge1 > 0)
-          v1 = e.V[0];
-        else
-          v1 = e.V[1];
-        List_Read(AllEdges, abs(Edge2) - 1, &e);
-        if(Edge2 > 0)
-          v2 = e.V[0];
-        else
-          v2 = e.V[1];
-        List_Read(AllEdges, abs(Edge3) - 1, &e);
-        if(Edge3 > 0)
-          v3 = e.V[0];
-        else
-          v3 = e.V[1];
-        v4 = NULL;
-      }
-      else if(NbEdgesOnFace == 4) {
-        fscanf(in, "%d %d %d %d %d", &Edge1, &Edge2, &Edge3, &Edge4, &nbPts);
-        List_Read(AllEdges, abs(Edge1) - 1, &e);
-        if(Edge1 > 0)
-          v1 = e.V[0];
-        else
-          v1 = e.V[1];
-        List_Read(AllEdges, abs(Edge2) - 1, &e);
-        if(Edge2 > 0)
-          v2 = e.V[0];
-        else
-          v2 = e.V[1];
-        List_Read(AllEdges, abs(Edge3) - 1, &e);
-        if(Edge3 > 0)
-          v3 = e.V[0];
-        else
-          v3 = e.V[1];
-        List_Read(AllEdges, abs(Edge4) - 1, &e);
-        if(Edge4 > 0)
-          v4 = e.V[0];
-        else
-          v4 = e.V[1];
-      }
-      else {
-        Msg(GERROR, "Wrong number pf edges on face (%d)", NbEdgesOnFace);
-      }
-      for(int j = 0; j < nbPts; j++) {
-        switch (GEntityType) {
-        case 0:
-          break;
-        case 1:
-          fscanf(in, "%le", &u);
-          break;
-        case 2:
-          fscanf(in, "%le %le %d", &u, &v, &patch);
-          break;
-        case 3:
-          break;
-        }
-      }
-
-      Simplex *s = Create_Simplex(v1, v2, v3, v4);
-      s->Num = i + 1;
-      s->iEnt = GEntityId + 10000;
-
-      Surface *surf;
-      List_Add(AllFaces, &s);
-
-      switch (GEntityType) {
-      case ENTITY_REGION:
-        break;
-      case ENTITY_FACE:
-        if(!(surf = FindSurface(GEntityId + 10000, m))) {
-          surf = Create_Surface(GEntityId + 10000, MSH_SURF_DISCRETE);
-          if(!NbRegions)
-            List_Add(vol->Surfaces, &surf);
-          Tree_Add(m->Surfaces, &surf);
-        }
-        Tree_Add(surf->Vertices, &s->V[0]);
-        Tree_Add(surf->Vertices, &s->V[1]);
-        Tree_Add(surf->Vertices, &s->V[2]);
-        Tree_Add(surf->Simplexes, &s);
-      }
-    }
-  }
-
-  Msg(INFO, "%d region", NbRegions);
-
-  for(int i = 0; i < NbRegions; i++) {
-    fscanf(in, "%d", &GEntityId);
-    if(GEntityId) {
-      fscanf(in, "%d", &NbFacesOnRegion);
-      Simplex *myS1, *myS2;
-      if(NbFacesOnRegion == 4) {
-        fscanf(in, "%d %d %d %d %d", &Face1, &Face2, &Face3, &Face4, &Dummy);
-        List_Read(AllFaces, abs(Face1) - 1, &myS1);
-        List_Read(AllFaces, abs(Face2) - 1, &myS2);
-        v1 = myS1->V[0];
-        v2 = myS1->V[1];
-        v3 = myS1->V[2];
-        for(int hh = 0; hh < 3; hh++)
-          if(compareVertex(&v1, &myS2->V[hh]) &&
-             compareVertex(&v2, &myS2->V[hh]) &&
-             compareVertex(&v3, &myS2->V[hh]))
-            v4 = myS2->V[hh];
-      }
-      if(!v1 || !v2 || !v3 || !v4) {
-        Msg(GERROR, "%d", NbFacesOnRegion);
-        Msg(GERROR, "%p %p %p %p", v1, v2, v3, v4);
-	Msg(GERROR, "%p %p %p", myS1->V[0], myS1->V[1], myS1->V[2]);
-        Msg(GERROR, "%p %p %p", myS2->V[0], myS2->V[1], myS2->V[2]);
-        return;
-      }
-      Simplex *s = Create_Simplex(v1, v2, v3, v4);
-
-      if(!(vol = FindVolume(GEntityId, m))) {
-        vol = Create_Volume(GEntityId, MSH_VOLUME_DISCRETE);
-        Tree_Add(m->Volumes, &vol);
-      }
-      s->iEnt = GEntityId;
-      Tree_Insert(vol->Simplexes, &s);
-      Tree_Insert(m->Simplexes, &s);
-    }
-  }
-
-  List_Delete(AllEdges);
-  List_Delete(AllFaces);
-
-  if(Tree_Nbr(m->Volumes)) {
-    m->status = 3;
-  }
-  else if(Tree_Nbr(m->Surfaces)) {
-    m->status = 2;
-  }
-  else if(Tree_Nbr(m->Curves)) {
-    m->status = 1;
-  }
-  else if(Tree_Nbr(m->Points))
-    m->status = 0;
-  else
-    m->status = -1;
-}
-
-// Public Read_Mesh routine
-
-void Read_Mesh(Mesh * M, FILE * fp, char *filename, int type)
-{
-  if(filename)
-    Msg(INFO, "Reading mesh file '%s'", filename);
-
-  switch (type) {
-  case FORMAT_MSH: Read_Mesh_MSH(M, fp); break;
-  case FORMAT_SMS: Read_Mesh_SMS(M, fp); break;
-  case FORMAT_VTK: Read_Mesh_VTK(M, fp); break;
-  default:
-    Msg(GERROR, "Unkown mesh file format");
-    return;
-  }
-
-  CTX.mesh.changed = 1;
-
-  if(filename){
-    Msg(INFO, "Read mesh file '%s'", filename);
-    Msg(STATUS2N, "Read '%s'", filename);
-  }
-
-  if(CTX.mesh.nbPartitions != 1)
-    PartitionMesh(M, CTX.mesh.nbPartitions);
-}
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index b81392b7ce4cdb268f91006f273673c8a28a22ca..88e53394c90d76ff83c4e3039ce8394b8abf3be5 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.99 2006-08-04 14:28:03 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.100 2006-08-07 00:08:09 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -344,13 +344,6 @@ int MergeProblem(char *name, int warn_if_missing)
       SetBoundingBox();
       status = THEM->status;
     }
-    else if(!strncmp(tmp, "sms", 3)) {
-      if(THEM->status < 0)
-	mai3d(0);
-      Read_Mesh(THEM, fp, name, FORMAT_SMS);
-      SetBoundingBox();
-      status = THEM->status;
-    }
     else if(!strncmp(tmp, "$PostFormat", 11) ||
 	    !strncmp(tmp, "$View", 5)) {
       ReadView(fp, name);
diff --git a/Plugin/StructuralSolver.cpp b/Plugin/StructuralSolver.cpp
index 08fc9492d30dd4695086659c1888e2f0e70b7b65..187af7f4882c5d4834fc3010e2732131ed5950a9 100644
--- a/Plugin/StructuralSolver.cpp
+++ b/Plugin/StructuralSolver.cpp
@@ -42,7 +42,8 @@ Structural_BeamSection:: Structural_BeamSection( const char *direct, std::string
   //  printf("%s/%s\n", direct,name.c_str());
   // read the section (msh format)
   FILE *f = fopen (temp,"r");  
-  Read_Mesh (&m, f, temp,FORMAT_MSH);
+  //need to rewrite this whole plugin using GModel
+  //Read_Mesh (&m, f, temp,FORMAT_MSH); 
   fclose(f);
   // get rid of the extension
   name.erase(name.find("."));