diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 87174d8eca44c9191c1f1495e9d24ebe3c15a0bc..a5f26b76741cf0cca8afe882c42dfe38a891df1e 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.77 2006-08-16 22:43:55 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.78 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -292,8 +292,9 @@ void Get_Options(int argc, char *argv[])
             WriteView(*(Post_View **) List_Pointer(CTX.post.list, j),
                       argv[i + 1], 1, j ? 1 : 0);
 	  // convert any mesh to the latest format
-	  if(GMODEL->numVertex()){
+	  if(GMODEL->getMeshStatus() > 0){
 	    CTX.mesh.msh_file_version = 2.0;
+	    CTX.mesh.msh_binary = 1;
 	    CreateOutputFile(argv[i + 1], FORMAT_MSH);
 	  }
         }
diff --git a/Common/Context.h b/Common/Context.h
index e3f9716bd6c25a20a35e118167b5a5e7c5d4699b..54fb93bddb3a46685ff309ccdfce4cda31426b41 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -182,7 +182,7 @@ public :
       return val;
     }
     int oldxtrude, oldxtrude_recombine;
-    int allow_degenerated_extrude, save_all, stl_binary;
+    int allow_degenerated_extrude, save_all, stl_binary, msh_binary;
     char *triangle_options;
     int smooth_normals, reverse_all_normals;
     double angle_smooth_normals;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 4e0b623307f6c55fbe41064cbf8837d268799ae6..e2689e356d51c40a35dce9db34d7a1f1d6d8ecd1 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -917,8 +917,10 @@ StringXNumber MeshOptions_Number[] = {
     "Minimum number of points used to mesh a circle" },
   { F|O, "MinimumElementSizeFact" , opt_mesh_min_elem_size_fact, 500. ,
     "Minimum element size factor in the Remesher" },
+  { F|O, "MshBinary" , opt_mesh_msh_binary , 0. , 
+    "Write MSH files in binary format?" },
   { F|O, "MshFileVersion" , opt_mesh_msh_file_version , 1.0 , 
-    "Version of the `msh' file format to use" },
+    "Version of the MSH file format to use" },
 
   { F|O, "NbElemsPerRadiusOfCurv" , opt_mesh_nb_elem_per_rc, 5. ,
     "Number of elements per radius of curvature in the remesher" },
diff --git a/Common/Options.cpp b/Common/Options.cpp
index d327fdc2cca8734ec6c67fcb060ece184d989439..0f4ccad4e509da099dcbcbbe29f1514e0a75043a 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.300 2006-08-19 04:24:02 geuzaine Exp $
+// $Id: Options.cpp,v 1.301 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -4626,6 +4626,13 @@ double opt_mesh_msh_file_version(OPT_ARGS_NUM)
   return CTX.mesh.msh_file_version;
 }
 
+double opt_mesh_msh_binary(OPT_ARGS_NUM)
+{
+  if(action & GMSH_SET)
+    CTX.mesh.msh_binary = (int)val;
+  return CTX.mesh.msh_binary;
+}
+
 double opt_mesh_stl_binary(OPT_ARGS_NUM)
 {
   if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 0cad7a65065914a6bedb4a0ce88b1b4f93b55311..cee5f195f72fa70558e5e79657a14427602dd3a8 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -447,6 +447,7 @@ double opt_mesh_light_lines(OPT_ARGS_NUM);
 double opt_mesh_light_two_side(OPT_ARGS_NUM);
 double opt_mesh_format(OPT_ARGS_NUM);
 double opt_mesh_msh_file_version(OPT_ARGS_NUM);
+double opt_mesh_msh_binary(OPT_ARGS_NUM);
 double opt_mesh_stl_binary(OPT_ARGS_NUM);
 double opt_mesh_nb_smoothing(OPT_ARGS_NUM);
 double opt_mesh_stl_distance_tol(OPT_ARGS_NUM);
diff --git a/Fltk/GUI_Extras.cpp b/Fltk/GUI_Extras.cpp
index ba4b10c38c4827fc22989beec573e5698732ae7d..5870e53324b0708da13fcde08a1317530dd56180 100644
--- a/Fltk/GUI_Extras.cpp
+++ b/Fltk/GUI_Extras.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI_Extras.cpp,v 1.19 2006-08-19 04:24:03 geuzaine Exp $
+// $Id: GUI_Extras.cpp,v 1.20 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -531,9 +531,10 @@ int msh_dialog(char *name)
   };
   static _msh_dialog *dialog = NULL;
 
-  static Fl_Menu_Item versionmenu[] = {
+  static Fl_Menu_Item formatmenu[] = {
     {"Version 1.0", 0, 0, 0},
-    {"Version 2.0", 0, 0, 0},
+    {"Version 2.0 ASCII", 0, 0, 0},
+    {"Version 2.0 Binary", 0, 0, 0},
     {0}
   };
 
@@ -544,7 +545,7 @@ int msh_dialog(char *name)
     dialog->window = new Fl_Double_Window(200, h, "MSH Options"); y = 10;
     dialog->window->box(GMSH_WINDOW_BOX);
     dialog->c = new Fl_Choice(10, y, 130, 25, "Format"); y+= 25;
-    dialog->c->menu(versionmenu);
+    dialog->c->menu(formatmenu);
     dialog->c->align(FL_ALIGN_RIGHT);
     dialog->b = new Fl_Check_Button(10, y, 180, 25, "Save all (ignore physicals)"); y += 25;
     dialog->b->type(FL_TOGGLE_BUTTON);
@@ -557,7 +558,8 @@ int msh_dialog(char *name)
     dialog->window->hotspot(dialog->window);
   }
   
-  dialog->c->value((CTX.mesh.msh_file_version==1.0) ? 0 : 1);
+  dialog->c->value((CTX.mesh.msh_file_version == 1.0) ? 0 : 
+		   CTX.mesh.msh_binary ? 2 : 1);
   dialog->b->value(CTX.mesh.save_all);
   dialog->window->show();
 
@@ -567,7 +569,10 @@ int msh_dialog(char *name)
       Fl_Widget* o = Fl::readqueue();
       if (!o) break;
       if (o == dialog->ok) {
-	opt_mesh_msh_file_version(0, GMSH_SET | GMSH_GUI, dialog->c->value() + 1);
+	opt_mesh_msh_file_version(0, GMSH_SET | GMSH_GUI, 
+				  (dialog->c->value() == 0) ? 1. : 2.);
+	opt_mesh_msh_binary(0, GMSH_SET | GMSH_GUI, 
+			    (dialog->c->value() == 2) ? 1 : 0);
 	opt_mesh_save_all(0, GMSH_SET | GMSH_GUI, dialog->b->value());
 	CreateOutputFile(name, FORMAT_MSH);
 	dialog->window->hide();
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 0164c298a9477946137b80eade8856e81d2df06a..aff08b32a84aeeadfeba40b5f0c0a34a13af020c 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.14 2006-08-19 08:26:47 remacle Exp $
+// $Id: GModel.cpp,v 1.15 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -167,11 +167,9 @@ SBoundingBox3d GModel::bounds()
 
 int GModel::getMeshStatus()
 {
-
-  // VERSION CG WAS WRONG
-  // A SURFACE CAN BE MESHESD WITHOUD INTERNAL VERTICES
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    if((*it)->tetrahedra.size() ||(*it)->hexahedra.size() || (*it)->prisms.size() || (*it)->pyramids.size()) return 3;
+    if((*it)->tetrahedra.size() ||(*it)->hexahedra.size() || 
+       (*it)->prisms.size() || (*it)->pyramids.size()) return 3;
   for(fiter it = firstFace(); it != lastFace(); ++it)
     if((*it)->triangles.size() || (*it)->quadrangles.size()) return 2;
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
diff --git a/Geo/GModel.h b/Geo/GModel.h
index b30e6d73edc016a3364f781e204d18677a301ce2..ca52a4e19d567f6f5e70195eaeb14d041318e774 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -119,8 +119,8 @@ class GModel
 
   // IO routines
   int readMSH(const std::string &name);
-  int writeMSH(const std::string &name, double version=1.0, bool saveAll=false,
-	       double scalingFactor=1.0);
+  int writeMSH(const std::string &name, double version=1.0, bool binary=false,
+	       bool saveAll=false, double scalingFactor=1.0);
   int writePOS(const std::string &name, double scalingFactor=1.0);
   int readSTL(const std::string &name, double tolerance=1.e-3);
   int writeSTL(const std::string &name, bool binary=false, double scalingFactor=1.0);
diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index 2854347247b9d517d7a3b98c863373be5986288b..c0614545e057f24be134ece2115f7fa809d2e611 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO.cpp,v 1.26 2006-08-19 04:24:03 geuzaine Exp $
+// $Id: GModelIO.cpp,v 1.27 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -30,8 +30,20 @@
 #include "MElement.h"
 #include "SBoundingBox3d.h"
 
+static void swapBytes(char *array, int size, int n)
+{
+  char *x = new char[size];
+  for(int i = 0; i < n; i++) {
+    char *a = &array[i * size];
+    memcpy(x, a, size);
+    for(int c = 0; c < size; c++)
+      a[size - 1 - c] = x[c];
+  }
+  delete [] x;
+}
+
 template<class T>
-void addElements(std::vector<T*> &dst, const std::vector<MElement*> &src)
+static void addElements(std::vector<T*> &dst, const std::vector<MElement*> &src)
 {
   for(unsigned int i = 0; i < src.size(); i++) dst.push_back((T*)src[i]);
 }
@@ -40,8 +52,7 @@ static void storeElementsInEntities(GModel *m,
 				    std::map<int, std::vector<MElement*> > &map)
 {
   std::map<int, std::vector<MElement*> >::const_iterator it = map.begin();
-  std::map<int, std::vector<MElement*> >::const_iterator ite = map.end();
-  for(; it != ite; ++it){
+  for(; it != map.end(); ++it){
     if(!it->second.size()) continue;
     int numEdges = it->second[0]->getNumEdges();
     switch(numEdges){
@@ -118,8 +129,7 @@ static void associateEntityWithVertices(GModel *m)
 static void storeVerticesInEntities(std::map<int, MVertex*> &vertices)
 {
   std::map<int, MVertex*>::const_iterator it = vertices.begin();
-  std::map<int, MVertex*>::const_iterator ite = vertices.end();
-  for(; it != ite; ++it){
+  for(; it != vertices.end(); ++it){
     MVertex *v = it->second;
     GEntity *ge = v->onWhat();
     if(ge) 
@@ -133,11 +143,13 @@ static void storeVerticesInEntities(std::vector<MVertex*> &vertices)
 {
   for(unsigned int i = 0; i < vertices.size(); i++){
     MVertex *v = vertices[i];
-    GEntity *ge = v->onWhat();
-    if(ge) 
-      ge->mesh_vertices.push_back(v);
-    else
-      delete v; // we delete all unused vertices
+    if(v){ // the vector can have null entries (first or last element)
+      GEntity *ge = v->onWhat();
+      if(ge) 
+	ge->mesh_vertices.push_back(v);
+      else
+	delete v; // we delete all unused vertices
+    }
   }
 }
 
@@ -145,8 +157,7 @@ static void storePhysicalTagsInEntities(GModel *m, int dim,
 					std::map<int, std::map<int, std::string> > &map)
 {
   std::map<int, std::map<int, std::string> >::const_iterator it = map.begin();
-  std::map<int, std::map<int, std::string> >::const_iterator ite = map.end();
-  for(; it != ite; ++it){
+  for(; it != map.end(); ++it){
     GEntity *ge = 0;
     switch(dim){
     case 0: ge = m->vertexByTag(it->first); break;
@@ -156,8 +167,7 @@ static void storePhysicalTagsInEntities(GModel *m, int dim,
     }
     if(ge){
       std::map<int, std::string>::const_iterator it2 = it->second.begin();
-      std::map<int, std::string>::const_iterator ite2 = it->second.end();
-      for(; it2 != ite2; ++it2)
+      for(; it2 != it->second.end(); ++it2)
 	ge->physicals.push_back(it2->first);
     }
   }
@@ -181,8 +191,120 @@ static int getNumVerticesForElementTypeMSH(int type)
   case PRI2: return 6 + 9 + 3;
   case PYR1: return 5;
   case PYR2: return 5 + 8 + 1;
-  default: return 0;
+  default: 
+    Msg(GERROR, "Unknown type of element for MSH format");
+    return 0;
+  }
+}
+
+template<class T>
+static void createElementMSH(GModel *m, int num, int type, int physical, 
+			     int elementary, int partition, int n[30], T &v, 
+			     std::map<int, std::vector<MVertex*> > &points,
+			     std::map<int, std::vector<MElement*> > elements[7],
+			     std::map<int, std::map<int, std::string> > physicals[4])
+{
+  int dim = 0;
+
+  switch (type) {
+  case PNT:
+    points[elementary].push_back(v[n[0]]);
+    dim = 0;
+    break;
+  case LGN1:
+    elements[0][elementary].push_back
+      (new MLine(v[n[0]], v[n[1]], num, partition));
+    dim = 1;
+    break;
+  case LGN2:
+    elements[0][elementary].push_back
+      (new MLine2(v[n[0]], v[n[1]], v[n[2]], num, partition));
+    dim = 1;
+    break;
+  case TRI1:
+    elements[1][elementary].push_back
+      (new MTriangle(v[n[0]], v[n[1]], v[n[2]], num, partition));
+    dim = 2;
+    break;
+  case TRI2:
+    elements[1][elementary].push_back
+      (new MTriangle2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		      num, partition));
+    dim = 2;
+    break;
+  case QUA1:
+    elements[2][elementary].push_back
+      (new MQuadrangle(v[n[0]], v[n[1]], v[n[2]], v[n[3]], num, partition));
+    dim = 2;
+    break;
+  case QUA2:
+    elements[2][elementary].push_back
+      (new MQuadrangle2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+			v[n[6]], v[n[7]], v[n[8]], num, partition));
+    dim = 2;
+    break;
+  case TET1:
+    elements[3][elementary].push_back
+      (new MTetrahedron(v[n[0]], v[n[1]], v[n[2]], v[n[3]], num, partition));
+    dim = 3; 
+    break;
+  case TET2:
+    elements[3][elementary].push_back
+      (new MTetrahedron2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+			 v[n[6]], v[n[7]], v[n[8]], v[n[9]], num, partition));
+    dim = 3; 
+    break;
+  case HEX1:
+    elements[4][elementary].push_back
+      (new MHexahedron(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		       v[n[6]], v[n[7]], num, partition));
+    dim = 3; 
+    break;
+  case HEX2:
+    elements[4][elementary].push_back
+      (new MHexahedron2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+			v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
+			v[n[12]], v[n[13]], v[n[14]], v[n[15]], v[n[16]], v[n[17]], 
+			v[n[18]], v[n[19]], v[n[20]], v[n[21]], v[n[22]], v[n[23]], 
+			v[n[24]], v[n[25]], v[n[26]], num, partition));
+    dim = 3; 
+    break;
+  case PRI1: 
+    elements[5][elementary].push_back
+      (new MPrism(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		  num, partition));
+    dim = 3; 
+    break;
+  case PRI2: 
+    elements[5][elementary].push_back
+      (new MPrism2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		   v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
+		   v[n[12]], v[n[13]], v[n[14]], v[n[15]], v[n[16]], v[n[17]], 
+		   num, partition));
+    dim = 3; 
+    break;
+  case PYR1: 
+    elements[6][elementary].push_back
+      (new MPyramid(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], num, partition));
+    dim = 3; 
+    break;
+  case PYR2: 
+    elements[6][elementary].push_back
+      (new MPyramid2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
+		     v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
+		     v[n[12]], v[n[13]], num, partition));
+    dim = 3; 
+    break;
+  default:
+    Msg(GERROR, "Unknown type (%d) for element %d", type, num); 
+    break;
   }
+  
+  if(physical && (!physicals[dim].count(elementary) || 
+		  !physicals[dim][elementary].count(physical)))
+    physicals[dim][elementary][physical] = "unnamed";
+  
+  if(partition) m->getMeshPartitions().insert(partition);
 }
 
 int GModel::readMSH(const std::string &name)
@@ -194,8 +316,10 @@ int GModel::readMSH(const std::string &name)
   }
 
   double version = 1.0;
+  bool binary = false, swap = false;
   char str[256];
-  std::map<int, MVertex*> vertices;
+  std::map<int, MVertex*> vertexMap;
+  std::vector<MVertex*> vertexVector;
   std::map<int, std::vector<MVertex*> > points;
   std::map<int, std::vector<MElement*> > elements[7];
   std::map<int, std::map<int, std::string> > physicals[4];
@@ -212,177 +336,145 @@ int GModel::readMSH(const std::string &name)
 
     if(!strncmp(&str[1], "MeshFormat", 10)) {
 
+      if(!fgets(str, sizeof(str), fp)) return 0;
       int format, size;
-      fscanf(fp, "%lf %d %d\n", &version, &format, &size);
+      if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) return 0;
+
+      if(format){
+	binary = true;
+	Msg(INFO, "Reading binary MSH file");
+	int one;
+	if(fread(&one, sizeof(int), 1, fp) != 1) return 0;
+	if(one != 1){
+	  swap = true;
+	  Msg(INFO, "Swapping bytes in binary file");
+	}
+      }
 
     }
     else if(!strncmp(&str[1], "NO", 2) || !strncmp(&str[1], "Nodes", 5)) {
 
+      if(!fgets(str, sizeof(str), fp)) return 0;
       int numVertices;
-      fscanf(fp, "%d", &numVertices);
+      if(sscanf(str, "%d", &numVertices) != 1) return 0;
       Msg(INFO, "%d vertices", numVertices);
 
       int progress = (numVertices > 100000) ? numVertices / 25 : 0;
+      int minVertex = numVertices + 1, maxVertex = -1;
       for(int i = 0; i < numVertices; i++) {
 	int num;
-	double x, y, z;
-        fscanf(fp, "%d %lf %lf %lf", &num, &x, &y, &z);
-	if(vertices.count(num))
+	double xyz[3];
+	if(!binary){
+	  if(fscanf(fp, "%d %lf %lf %lf", &num, &xyz[0], &xyz[1], &xyz[2]) != 4) return 0;
+	}
+	else{
+	  if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
+	  if(swap) swapBytes((char*)&num, sizeof(int), 1);
+	  if(fread(xyz, sizeof(double), 3, fp) != 3) return 0;
+	  if(swap) swapBytes((char*)&xyz, sizeof(double), 3);
+	}
+	minVertex = std::min(minVertex, num);
+	maxVertex = std::max(maxVertex, num);
+	if(vertexMap.count(num))
 	  Msg(WARNING, "Skipping duplicate vertex %d", num);
 	else
-	  vertices[num] = new MVertex(x, y, z);
+	  vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2]);
 	if(progress && (i % progress == progress - 1))
 	  Msg(PROGRESS, "Read %d vertices", i + 1);
       }
-      Msg(PROGRESS, "");
+      if(progress) Msg(PROGRESS, "");
+      
+      // If the vertex numbering is dense, tranfer the map into a
+      // vector to speed up element creation
+      if((int)vertexMap.size() == numVertices && 
+	 ((minVertex == 1 && maxVertex == numVertices) ||
+	  (minVertex == 0 && maxVertex == numVertices - 1))){
+	Msg(INFO, "Vertex numbering is dense");
+	vertexVector.resize(vertexMap.size() + 1);
+	if(minVertex == 1) 
+	  vertexVector[0] = 0;
+	else
+	  vertexVector[numVertices] = 0;
+	std::map<int, MVertex*>::const_iterator it = vertexMap.begin();
+	for(; it != vertexMap.end(); ++it)
+	  vertexVector[it->first] = it->second;
+	vertexMap.clear();
+      }
 
     }
     else if(!strncmp(&str[1], "ELM", 3) || !strncmp(&str[1], "Elements", 8)) {
 
+      if(!fgets(str, sizeof(str), fp)) return 0;
       int numElements;
-      fscanf(fp, "%d", &numElements);
+      sscanf(str, "%d", &numElements);
       Msg(INFO, "%d elements", numElements);
 
       int progress = (numElements > 100000) ? numElements / 25 : 0;
-      for(int i = 0; i < numElements; i++) {
-	int num, type, physical = 0, elementary = 0, partition = 0, numVertices;
-	if(version <= 1.0){
-	  fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices);
-	  int check = getNumVerticesForElementTypeMSH(type);
-	  if(!check){
-	    Msg(GERROR, "Unknown type for element %d", num); 
-	    continue;
+      if(!binary){
+	for(int i = 0; i < numElements; i++) {
+	  int num, type, physical = 0, elementary = 0, partition = 0, numVertices;
+	  if(version <= 1.0){
+	    fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices);
+	    if(numVertices != getNumVerticesForElementTypeMSH(type)) return 0;
 	  }
-	  if(numVertices != check){
-	    Msg(GERROR, "Wrong number of vertices (%d) for element %d", numVertices, num);
-	    continue;
+	  else{
+	    int numTags;
+	    fscanf(fp, "%d %d %d", &num, &type, &numTags);
+	    for(int j = 0; j < numTags; j++){
+	      int tag;
+	      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
+	    }
+	    if(!(numVertices = getNumVerticesForElementTypeMSH(type))) return 0;
 	  }
+	  int indices[30];
+	  for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &indices[j]);
+	  if(vertexVector.size())
+	    createElementMSH(this, num, type, physical, elementary, partition, indices,
+			     vertexVector, points, elements, physicals);
+	  else
+	    createElementMSH(this, num, type, physical, elementary, partition, indices,
+			     vertexMap, points, elements, physicals);
+	  if(progress && (i % progress == progress - 1))
+	    Msg(PROGRESS, "Read %d elements", i + 1);
 	}
-	else{
-	  int numTags;
-	  fscanf(fp, "%d %d %d", &num, &type, &numTags);
-	  for(int j = 0; j < numTags; j++){
-	    int tag;
-	    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
-	  }
-	  numVertices = getNumVerticesForElementTypeMSH(type);
-	  if(!numVertices){
-	    Msg(GERROR, "Unknown type (%d) for element %d", type, num); 
-	    continue;
+      }
+      else{
+	int numElementsPartial = 0;
+	while(numElementsPartial < numElements){
+	  int tags[3];
+	  if(fread(tags, sizeof(int), 3, fp) != 3) return 0;
+	  if(swap) swapBytes((char*)&tags, sizeof(int), 3);
+	  int type = tags[0];
+	  int numElms = tags[1];
+	  int numTags = tags[2];
+	  unsigned int n = 1 + numTags + getNumVerticesForElementTypeMSH(type);
+	  int *data = new int[n];
+	  for(int i = 0; i < numElms; i++) {
+	    if(fread(data, sizeof(int), n, fp) != n) return 0;
+	    if(swap) swapBytes((char*)&data, sizeof(int), n);
+	    int num = data[0];
+	    int physical = (numTags > 0) ? data[4 - numTags] : 0;
+	    int elementary = (numTags > 1) ? data[4 - numTags + 1] : 0;
+	    int partition = (numTags > 2) ? data[4 - numTags + 2] : 0;
+	    int *verts = &data[numTags + 1];
+	    if(vertexVector.size())
+	      createElementMSH(this, num, type, physical, elementary, partition,
+			       verts, vertexVector, points, elements, physicals);
+	    else
+	      createElementMSH(this, num, type, physical, elementary, partition, 
+			       verts, vertexMap, points, elements, physicals);
+	    if(progress && ((numElementsPartial + i) % progress == progress - 1))
+	      Msg(PROGRESS, "Read %d elements", i + 1);
 	  }
+	  delete [] data;
+	  numElementsPartial += numElms;
 	}
-	int n[30];
-        for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &n[j]);
-	int dim = 0;
-	std::map<int, MVertex*> &v(vertices);
-        switch (type) {
-        case PNT:
-	  points[elementary].push_back(v[n[0]]);
-	  dim = 0;
-          break;
-        case LGN1:
-	  elements[0][elementary].push_back
-	    (new MLine(v[n[0]], v[n[1]], num, partition));
-	  dim = 1;
-          break;
-        case LGN2:
-	  elements[0][elementary].push_back
-	    (new MLine2(v[n[0]], v[n[1]], v[n[2]], num, partition));
-	  dim = 1;
-          break;
-        case TRI1:
-	  elements[1][elementary].push_back
-	    (new MTriangle(v[n[0]], v[n[1]], v[n[2]], num, partition));
-	  dim = 2;
-          break;
-        case TRI2:
-	  elements[1][elementary].push_back
-	    (new MTriangle2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			    num, partition));
-	  dim = 2;
-          break;
-        case QUA1:
-	  elements[2][elementary].push_back
-	    (new MQuadrangle(v[n[0]], v[n[1]], v[n[2]], v[n[3]], num, partition));
-	  dim = 2;
-          break;
-        case QUA2:
-	  elements[2][elementary].push_back
-	    (new MQuadrangle2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			      v[n[6]], v[n[7]], v[n[8]], num, partition));
-	  dim = 2;
-          break;
-	case TET1:
-	  elements[3][elementary].push_back
-	    (new MTetrahedron(v[n[0]], v[n[1]], v[n[2]], v[n[3]], num, partition));
-	  dim = 3; 
-	  break;
-	case TET2:
-	  elements[3][elementary].push_back
-	    (new MTetrahedron2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			       v[n[6]], v[n[7]], v[n[8]], v[n[9]], num, partition));
-	  dim = 3; 
-	  break;
-	case HEX1:
-	  elements[4][elementary].push_back
-	    (new MHexahedron(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			     v[n[6]], v[n[7]], num, partition));
-	  dim = 3; 
-	  break;
-	case HEX2:
-	  elements[4][elementary].push_back
-	    (new MHexahedron2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			      v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
-			      v[n[12]], v[n[13]], v[n[14]], v[n[15]], v[n[16]], v[n[17]], 
-			      v[n[18]], v[n[19]], v[n[20]], v[n[21]], v[n[22]], v[n[23]], 
-			      v[n[24]], v[n[25]], v[n[26]], num, partition));
-	  dim = 3; 
-	  break;
-	case PRI1: 
-	  elements[5][elementary].push_back
-	    (new MPrism(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			num, partition));
-	  dim = 3; 
-	  break;
-	case PRI2: 
-	  elements[5][elementary].push_back
-	    (new MPrism2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			 v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
-			 v[n[12]], v[n[13]], v[n[14]], v[n[15]], v[n[16]], v[n[17]], 
-			 num, partition));
-	  dim = 3; 
-	  break;
-	case PYR1: 
-	  elements[6][elementary].push_back
-	    (new MPyramid(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], num, partition));
-	  dim = 3; 
-	  break;
-	case PYR2: 
-	  elements[6][elementary].push_back
-	    (new MPyramid2(v[n[0]], v[n[1]], v[n[2]], v[n[3]], v[n[4]], v[n[5]], 
-			   v[n[6]], v[n[7]], v[n[8]], v[n[9]], v[n[10]], v[n[11]], 
-			   v[n[12]], v[n[13]], num, partition));
-	  dim = 3; 
-	  break;
-        default:
-	  Msg(GERROR, "Unknown type (%d) for element %d", type, num); 
-          break;
-        }
-
-	if(physical && (!physicals[dim].count(elementary) || 
-			!physicals[dim][elementary].count(physical)))
-	    physicals[dim][elementary][physical] = "unnamed";
-	
-	if(partition) meshPartitions.insert(partition);
-	
-	if(progress && (i % progress == progress - 1))
-	  Msg(PROGRESS, "Read %d elements", i + 1);
       }
-      Msg(PROGRESS, "");
+      if(progress) Msg(PROGRESS, "");
 
     }
 
@@ -400,8 +492,7 @@ int GModel::readMSH(const std::string &name)
   // treat points separately
   {
     std::map<int, std::vector<MVertex*> >::const_iterator it = points.begin();
-    std::map<int, std::vector<MVertex*> >::const_iterator ite = points.end();
-    for(; it != ite; ++it){
+    for(; it != points.end(); ++it){
       GVertex *v = vertexByTag(it->first);
       if(!v){
 	v = new gmshVertex(this, it->first);
@@ -422,7 +513,10 @@ int GModel::readMSH(const std::string &name)
     (*it)->mesh_vertices.clear();
 
   // store the vertices in their associated geometrical entity
-  storeVerticesInEntities(vertices);
+  if(vertexVector.size())
+    storeVerticesInEntities(vertexVector);
+  else
+    storeVerticesInEntities(vertexMap);
 
   // store the physical tags
   for(int i = 0; i < 4; i++)  
@@ -432,21 +526,27 @@ int GModel::readMSH(const std::string &name)
   return 1;
 }
 
+static void writeTagMSH(FILE *fp, int type, int num, int numTags)
+{
+  int data[3] = {type, num, numTags};
+  fwrite(data, sizeof(int), 3, fp);
+}
+
 template<class T>
 static void writeElementsMSH(FILE *fp, const std::vector<T*> &ele, int saveAll, 
-			     double version, int &num, int elementary, 
+			     double version, bool binary, int &num, int elementary, 
 			     std::vector<int> &physicals)
 {
   for(unsigned int i = 0; i < ele.size(); i++)
     if(saveAll)
-      ele[i]->writeMSH(fp, version, ++num, elementary, 0);
+      ele[i]->writeMSH(fp, version, binary, ++num, elementary, 0);
     else
       for(unsigned int j = 0; j < physicals.size(); j++)
-	ele[i]->writeMSH(fp, version, ++num, elementary, physicals[j]);
+	ele[i]->writeMSH(fp, version, binary, ++num, elementary, physicals[j]);
 }
 
-int GModel::writeMSH(const std::string &name, double version, bool saveAll, 
-		     double scalingFactor)
+int GModel::writeMSH(const std::string &name, double version, bool binary, 
+		     bool saveAll, double scalingFactor)
 {
   FILE *fp = fopen(name.c_str(), "w");
   if(!fp){
@@ -454,6 +554,9 @@ int GModel::writeMSH(const std::string &name, double version, bool saveAll,
     return 0;
   }
 
+  // binary format exists only in version 2
+  if(binary) version = 2.0;
+
   // if there are no physicals we save all the elements
   if(noPhysicalGroups()) saveAll = true;
 
@@ -461,46 +564,75 @@ int GModel::writeMSH(const std::string &name, double version, bool saveAll,
   // continuous sequence
   int numVertices = renumberMeshVertices();
   
-  // get the number of elements
+  // get the number of elements (we assume that all the elements in a
+  // list have the same type, i.e., they are all of the same
+  // polynomial order)
+  std::map<int,int> elements;
+  for(viter it = firstVertex(); it != lastVertex(); ++it){
+    int p = (saveAll ? 1 : (*it)->physicals.size());
+    int n = p * (*it)->mesh_vertices.size();
+    if(n) elements[PNT] += n;
+  }
+  for(eiter it = firstEdge(); it != lastEdge(); ++it){
+    int p = (saveAll ? 1 : (*it)->physicals.size());
+    int n = p * (*it)->lines.size();
+    if(n) elements[(*it)->lines[0]->getTypeForMSH()] += n;
+  }
+  for(fiter it = firstFace(); it != lastFace(); ++it){
+    int p = (saveAll ? 1 : (*it)->physicals.size());
+    int n = p * (*it)->triangles.size();
+    if(n) elements[(*it)->triangles[0]->getTypeForMSH()] += n;
+    n = p * (*it)->quadrangles.size();
+    if(n) elements[(*it)->quadrangles[0]->getTypeForMSH()] += n;
+  }
+  for(riter it = firstRegion(); it != lastRegion(); ++it){
+    int p = (saveAll ? 1 : (*it)->physicals.size());
+    int n = p * (*it)->tetrahedra.size();
+    if(n) elements[(*it)->tetrahedra[0]->getTypeForMSH()] += n;
+    n = p * (*it)->hexahedra.size();
+    if(n) elements[(*it)->hexahedra[0]->getTypeForMSH()] += n;
+    n = p * (*it)->prisms.size();
+    if(n) elements[(*it)->prisms[0]->getTypeForMSH()] += n;
+    n = p * (*it)->pyramids.size();
+    if(n) elements[(*it)->pyramids[0]->getTypeForMSH()] += n;
+  }
+
   int numElements = 0;
-  for(viter it = firstVertex(); it != lastVertex(); ++it)
-    numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
-      (*it)->mesh_vertices.size();
-  for(eiter it = firstEdge(); it != lastEdge(); ++it)
-    numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
-      (*it)->lines.size();
-  for(fiter it = firstFace(); it != lastFace(); ++it)
-    numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
-      ((*it)->triangles.size() + (*it)->quadrangles.size());
-  for(riter it = firstRegion(); it != lastRegion(); ++it)
-    numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
-      ((*it)->tetrahedra.size() + (*it)->hexahedra.size() +
-       (*it)->prisms.size() + (*it)->pyramids.size());
+  std::map<int,int>::const_iterator it = elements.begin();
+  for(; it != elements.end(); ++it)
+    numElements += it->second;
 
-  if(version > 2.0){
+  if(version >= 2.0){
     fprintf(fp, "$MeshFormat\n");
-    fprintf(fp, "%g %d %d\n", version, 0, (int)sizeof(double));
+    fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
+    if(binary){
+      int one = 1;
+      fwrite(&one, sizeof(int), 1, fp);
+      fprintf(fp, "\n");
+    }
     fprintf(fp, "$EndMeshFormat\n");
     fprintf(fp, "$Nodes\n");
   }
   else
     fprintf(fp, "$NOD\n");
- 
+
   fprintf(fp, "%d\n", numVertices);
   for(viter it = firstVertex(); it != lastVertex(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++)
-      (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
   for(fiter it = firstFace(); it != lastFace(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
   for(riter it = firstRegion(); it != lastRegion(); ++it)
     for(unsigned int i = 0; i < (*it)->mesh_vertices.size(); i++) 
-      (*it)->mesh_vertices[i]->writeMSH(fp, scalingFactor);
+      (*it)->mesh_vertices[i]->writeMSH(fp, binary, scalingFactor);
 
-  if(version > 2.0){
+  if(binary) fprintf(fp, "\n");
+
+  if(version >= 2.0){
     fprintf(fp, "$EndNodes\n");
     fprintf(fp, "$Elements\n");
   }
@@ -510,34 +642,58 @@ int GModel::writeMSH(const std::string &name, double version, bool saveAll,
   }
 
   fprintf(fp, "%d\n", numElements);
-  int num = 0;
+  int num = 0, numTags = 3;
 
-  for(viter it = firstVertex(); it != lastVertex(); ++it){
-    writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, num,
+  if(binary && elements.count(PNT)) writeTagMSH(fp, PNT, elements[PNT], numTags);
+  for(viter it = firstVertex(); it != lastVertex(); ++it)
+    writeElementsMSH(fp, (*it)->mesh_vertices, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-  }
-  for(eiter it = firstEdge(); it != lastEdge(); ++it){
-    writeElementsMSH(fp, (*it)->lines, saveAll, version, num,
+
+  if(binary && elements.count(LGN1)) writeTagMSH(fp, LGN1, elements[LGN1], numTags);
+  else if(binary && elements.count(LGN2)) writeTagMSH(fp, LGN2, elements[LGN2], numTags);
+  for(eiter it = firstEdge(); it != lastEdge(); ++it)
+    writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-  }
-  for(fiter it = firstFace(); it != lastFace(); ++it){
-    writeElementsMSH(fp, (*it)->triangles, saveAll, version, num,
+
+  if(binary && elements.count(TRI1)) writeTagMSH(fp, TRI1, elements[TRI1], numTags);
+  else if(binary && elements.count(TRI2)) writeTagMSH(fp, TRI2, elements[TRI2], numTags);
+  for(fiter it = firstFace(); it != lastFace(); ++it)
+    writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-    writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, num,
+
+  if(binary && elements.count(QUA1)) writeTagMSH(fp, QUA1, elements[QUA1], numTags);
+  else if(binary && elements.count(QUA2)) writeTagMSH(fp, QUA2, elements[QUA2], numTags);
+  for(fiter it = firstFace(); it != lastFace(); ++it)
+    writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-  }
-  for(riter it = firstRegion(); it != lastRegion(); ++it){
-    writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, num,
+
+  if(binary && elements.count(TET1)) writeTagMSH(fp, TET1, elements[TET1], numTags);
+  else if(binary && elements.count(TET2)) writeTagMSH(fp, TET2, elements[TET2], numTags);
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-    writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, num,
+
+  if(binary && elements.count(HEX1)) writeTagMSH(fp, HEX1, elements[HEX1], numTags);
+  else if(binary && elements.count(HEX2)) writeTagMSH(fp, HEX2, elements[HEX2], numTags);
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-    writeElementsMSH(fp, (*it)->prisms, saveAll, version, num,
+
+  if(binary && elements.count(PRI1)) writeTagMSH(fp, PRI1, elements[PRI1], numTags);
+  else if(binary && elements.count(PRI2)) writeTagMSH(fp, PRI2, elements[PRI2], numTags);
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    writeElementsMSH(fp, (*it)->prisms, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-    writeElementsMSH(fp, (*it)->pyramids, saveAll, version, num,
+
+  if(binary && elements.count(PYR1)) writeTagMSH(fp, PYR1, elements[PYR1], numTags);
+  else if(binary && elements.count(PYR2)) writeTagMSH(fp, PYR2, elements[PYR2], numTags);
+  for(riter it = firstRegion(); it != lastRegion(); ++it)
+    writeElementsMSH(fp, (*it)->pyramids, saveAll, version, binary, num,
 		     (*it)->tag(), (*it)->physicals);
-  }
   
-  if(version > 2.0){
+  if(binary) fprintf(fp, "\n");
+
+  if(version >= 2.0){
     fprintf(fp, "$EndElements\n");
   }
   else{
@@ -601,18 +757,6 @@ int GModel::writePOS(const std::string &name, double scalingFactor)
   return 1;
 }
 
-static void swapBytes(char *array, int size, int n)
-{
-  char *x = new char[size];
-  for(int i = 0; i < n; i++) {
-    char *a = &array[i * size];
-    memcpy(x, a, size);
-    for(int c = 0; c < size; c++)
-      a[size - 1 - c] = x[c];
-  }
-  delete [] x;
-}
-
 int GModel::readSTL(const std::string &name, double tolerance)
 {
   FILE *fp = fopen(name.c_str(), "rb");
@@ -1035,8 +1179,7 @@ int GModel::writeUNV(const std::string &name, double scalingFactor)
 
   for(int dim = 0; dim < 4; dim++){
     std::map<int, std::vector<GEntity*> >::const_iterator it = physicals[dim].begin();
-    std::map<int, std::vector<GEntity*> >::const_iterator ite = physicals[dim].end();
-    for(; it != ite; ++it){
+    for(; it != physicals[dim].end(); ++it){
       // IDEAS GROUPOFNODES record
       fprintf(fp, "%6d\n", -1);
       fprintf(fp, "%6d\n", 790);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 23e1130251a39bb5682b19fdcd913a358ab67ddb..2ad65c42823797e2792cdc728e644d1baa05aead 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.9 2006-08-19 04:24:03 geuzaine Exp $
+// $Id: MElement.cpp,v 1.10 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -107,8 +107,8 @@ SPoint3 MElement::barycenter()
   return p;
 }
 
-void MElement::writeMSH(FILE *fp, double version, int num, int elementary, 
-			int physical)
+void MElement::writeMSH(FILE *fp, double version, bool binary, int num, 
+			int elementary, int physical)
 {
   // if necessary, change the ordering of the vertices to get positive
   // volume
@@ -117,32 +117,48 @@ void MElement::writeMSH(FILE *fp, double version, int num, int elementary,
   int n = getNumVertices();
   int type = getTypeForMSH();
 
-  fprintf(fp, "%d %d", num ? num : _num, type);
-  if(version < 2.0)
-    fprintf(fp, " %d %d %d", physical, elementary, n);
-  else
-    fprintf(fp, " 3 %d %d %d", physical, elementary, _partition);
-  
+  if(!binary){
+    fprintf(fp, "%d %d", num ? num : _num, type);
+    if(version < 2.0)
+      fprintf(fp, " %d %d %d", physical, elementary, n);
+    else
+      fprintf(fp, " 3 %d %d %d", physical, elementary, _partition);
+  }
+  else{
+    int tags[4] = {num ? num : _num, physical, elementary, _partition};
+    fwrite(tags, sizeof(int), 4, fp);
+  }
+
+  int verts[30];
+
   if(physical >= 0){
     for(int i = 0; i < n; i++)
-      fprintf(fp, " %d", getVertex(i)->getNum());
+      verts[i] = getVertex(i)->getNum();
   }
   else{
     int nn = n - getNumEdgeVertices() - getNumFaceVertices() - getNumVolumeVertices();
+    int j = 0;
     for(int i = 0; i < nn; i++)
-      fprintf(fp, " %d", getVertex(nn - i - 1)->getNum());
+      verts[j++] = getVertex(nn - i - 1)->getNum();
     int ne = getNumEdgeVertices();
     for(int i = 0; i < ne; i++)
-      fprintf(fp, " %d", getVertex(nn + ne - i - 1)->getNum());
+      verts[j++] = getVertex(nn + ne - i - 1)->getNum();
     int nf = getNumFaceVertices();
     for(int i = 0; i < nf; i++)
-      fprintf(fp, " %d", getVertex(nn + ne + nf - i - 1)->getNum());
+      verts[j++] = getVertex(nn + ne + nf - i - 1)->getNum();
     int nv = getNumVolumeVertices();
     for(int i = 0; i < nv; i++)
-      fprintf(fp, " %d", getVertex(n - i - 1)->getNum());
+      verts[j++] = getVertex(n - i - 1)->getNum();
   }
 
-  fprintf(fp, "\n");
+  if(!binary){
+    for(int i = 0; i < n; i++)
+      fprintf(fp, " %d", verts[i]);
+    fprintf(fp, "\n");
+  }
+  else{
+    fwrite(verts, sizeof(int), n, fp);
+  }
 }
 
 void MElement::writePOS(FILE *fp, double scalingFactor, int elementary)
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 004784d7cac0c20c404b4af48eddb5c704d3d411..f26241218c9a50f95b742644c4175e3b478f4658 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -127,8 +127,8 @@ class MElement
   virtual void setVolumePositive(){}
 
   // IO routines
-  virtual void writeMSH(FILE *fp, double version=1.0, int num=0, 
-			int elementary=1, int physical=1);
+  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, 
+			int num=0, int elementary=1, int physical=1);
   virtual void writePOS(FILE *fp, double scalingFactor=1.0,
 			int elementary=1);
   virtual void writeSTL(FILE *fp, bool binary=false, double scalingFactor=1.0);
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index fd582d3b874fe6b7b8e76cd08d2ed5bfb7379e4d..39a4db192f44ae63b7277e776d7ffef72c076ad0 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: MVertex.cpp,v 1.5 2006-08-15 06:26:52 geuzaine Exp $
+// $Id: MVertex.cpp,v 1.6 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -25,21 +25,38 @@
 int MVertex::_globalNum = 0;
 double MVertexLessThanLexicographic::tolerance = 1.e-6;
 
-void MVertex::writeMSH(FILE *fp, double scalingFactor)
+void MVertex::writeMSH(FILE *fp, bool binary, double scalingFactor)
 {
-  fprintf(fp, "%d %.16g %.16g %.16g\n", _num, x() * scalingFactor, 
-	  y() * scalingFactor, z() * scalingFactor);
+  if(!binary){
+    fprintf(fp, "%d %.16g %.16g %.16g\n", _num, 
+	    x() * scalingFactor, 
+	    y() * scalingFactor,
+	    z() * scalingFactor);
+  }
+  else{
+    fwrite(&_num, sizeof(int), 1, fp);
+    double data[3] = {x() * scalingFactor, y() * scalingFactor, z() * scalingFactor};
+    fwrite(data, sizeof(double), 3, fp);
+  }
 }
 
-void MVertex::writeMSH(FILE *fp, double version, int num, 
+void MVertex::writeMSH(FILE *fp, double version, bool binary, int num, 
 		       int elementary, int physical)
 {
-  fprintf(fp, "%d 15", num);
-  if(version < 2.0)
-    fprintf(fp, " %d %d 1", physical, elementary);
-  else
-    fprintf(fp, " 2 %d %d", physical, elementary);
-  fprintf(fp, " %d\n", _num);
+  if(!binary){
+    fprintf(fp, "%d 15", num);
+    if(version < 2.0)
+      fprintf(fp, " %d %d 1", physical, elementary);
+    else
+      fprintf(fp, " 2 %d %d", physical, elementary);
+    fprintf(fp, " %d\n", _num);
+  }
+  else{
+    int tags[4] = {num, physical, elementary, 0};
+    fwrite(tags, sizeof(int), 4, fp);
+    int verts[1] = {_num};
+    fwrite(verts, sizeof(int), 1, fp);
+  }
 }
 
 void MVertex::writeVRML(FILE *fp, double scalingFactor)
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 93d68783eb72e248576494a9cd823251b372c51a..788bbef907140fc3dd07501f5cade4852eef6c2c 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -73,8 +73,9 @@ class MVertex{
   inline void setNum(int num) { _num = num; }
 
   // IO routines
-  void writeMSH(FILE *fp, double scalingFactor=1.0);
-  void writeMSH(FILE *fp, double version, int num, int elementary, int physical);
+  void writeMSH(FILE *fp, bool binary=false, double scalingFactor=1.0);
+  void writeMSH(FILE *fp, double version, bool binary, int num, 
+		int elementary, int physical);
   void writeVRML(FILE *fp, double scalingFactor=1.0);
   void writeUNV(FILE *fp, double scalingFactor=1.0);
   void writeMESH(FILE *fp, double scalingFactor=1.0);
diff --git a/Geo/Makefile b/Geo/Makefile
index dcd32a7e76d55d25ec340318a0f92016c179b0e9..e5ab20422a6037470927f1a58283120a6e97835b 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.96 2006-08-19 08:26:47 remacle Exp $
+# $Id: Makefile,v 1.97 2006-08-19 18:48:06 geuzaine Exp $
 #
 # Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 #
@@ -26,7 +26,7 @@ LIB     = ../lib/libGmshGeo.a
 INCLUDE = -I../Common -I../DataStr -I../Geo -I../Mesh -I../Numeric\
           -I../Parser -I../Fltk\
           -I../contrib/NR
-CFLAGS  = -pg ${OPTIM} ${FLAGS} ${INCLUDE} 
+CFLAGS  = ${OPTIM} ${FLAGS} ${INCLUDE} 
 
 SRC = CAD.cpp \
       ExtrudeParams.cpp \
diff --git a/Graphics/Mesh.cpp b/Graphics/Mesh.cpp
index 7ca2bd049b35782bb742edadc0e739c536c7770f..dbcc15c5ff2213f3a222f000d6e7e1d53ec11536 100644
--- a/Graphics/Mesh.cpp
+++ b/Graphics/Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: Mesh.cpp,v 1.179 2006-08-19 08:26:47 remacle Exp $
+// $Id: Mesh.cpp,v 1.180 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -794,8 +794,6 @@ class drawMeshGRegion {
 
 void Draw_Mesh()
 {
-
-
   if(!CTX.mesh.draw) return;
   
   glPointSize(CTX.mesh.point_size);
diff --git a/Parser/CreateFile.cpp b/Parser/CreateFile.cpp
index 20e18f3b229eae9b419a8dbdcd3ee3f6d4ba6c89..2c7858a06a5d2ca845ae2cde6dfde9eec5df4383 100644
--- a/Parser/CreateFile.cpp
+++ b/Parser/CreateFile.cpp
@@ -1,4 +1,4 @@
-// $Id: CreateFile.cpp,v 1.2 2006-08-19 04:24:03 geuzaine Exp $
+// $Id: CreateFile.cpp,v 1.3 2006-08-19 18:48:06 geuzaine Exp $
 //
 // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 //
@@ -145,8 +145,8 @@ void CreateOutputFile(char *filename, int format)
     break;
 
   case FORMAT_MSH:
-    GMODEL->writeMSH(name, CTX.mesh.msh_file_version, CTX.mesh.save_all,
-		     CTX.mesh.scaling_factor);
+    GMODEL->writeMSH(name, CTX.mesh.msh_file_version, CTX.mesh.msh_binary, 
+		     CTX.mesh.save_all, CTX.mesh.scaling_factor);
     break;
 
   case FORMAT_STL:
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index f5ff203569a7d2ebc8580504f7c845b75022cbfd..aa2e554d021754368c0cb37e2bb152e6088c88af 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -1,5 +1,5 @@
 \input texinfo.tex @c -*-texinfo-*-
-@c $Id: gmsh.texi,v 1.211 2006-08-18 15:41:58 geuzaine Exp $
+@c $Id: gmsh.texi,v 1.212 2006-08-19 18:48:06 geuzaine Exp $
 @c
 @c Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
 @c
@@ -237,7 +237,8 @@ File formats
 Gmsh mesh file formats
 
 * Version 1.0::                 
-* Version 2.0::                 
+* Version 2.0 ASCII::           
+* Version 2.0 Binary::          
 
 Gmsh post-processing file formats
 
@@ -2870,18 +2871,19 @@ program is also a good example on how to read and write files in the
 
 @menu
 * Version 1.0::                 
-* Version 2.0::                 
+* Version 2.0 ASCII::           
+* Version 2.0 Binary::          
 @end menu
 
 @c .........................................................................
 @c Version 1.0
 @c .........................................................................
 
-@node Version 1.0, Version 2.0, Gmsh mesh file formats, Gmsh mesh file formats
+@node Version 1.0, Version 2.0 ASCII, Gmsh mesh file formats, Gmsh mesh file formats
 @subsection Version 1.0
 
 The @file{.msh} file format, version 1.0, is Gmsh's old native mesh file
-format, now superseded by the format described in @ref{Version 2.0}.
+format, now superseded by the format described in @ref{Version 2.0 ASCII}.
 
 In the @file{.msh} file format, version 1.0, the file is divided in two
 sections, defining the nodes (@code{$NOD}-@code{$ENDNOD}) and the elements
@@ -2987,16 +2989,17 @@ of the edges and quadrangular faces given in @ref{Gmsh node ordering}.
 @end table
 
 @c .........................................................................
-@c Version 2.0
+@c Version 2.0 ASCII
 @c .........................................................................
 
-@node Version 2.0,  , Version 1.0, Gmsh mesh file formats
-@subsection Version 2.0
+@node Version 2.0 ASCII, Version 2.0 Binary, Version 1.0, Gmsh mesh file formats
+@subsection Version 2.0 ASCII
 
-Version 2.0 of the @file{.msh} file format is Gmsh's new native mesh file
-format. It is very similar to the old one (@pxref{Version 1.0}), but is more
-general: it contains information about itself and allows to associate an
-arbitrary number of integer tags with each element.
+The version 2.0 of the @file{.msh} file format is Gmsh's new native mesh
+file format. It is very similar to the old one (@pxref{Version 1.0}),
+but is more general: it contains information about itself and allows to
+associate an arbitrary number of integer tags with each element. Ialso
+exists in both ASCII and binary form.
 
 The @file{.msh} file format, version 2.0, is divided in three sections,
 defining the file format (@code{$MeshFormat}-@code{$EndMeshFormat}), the
@@ -3027,7 +3030,7 @@ is an integer equal to 0 in the ASCII file format.
 
 @item @var{data-size}
 is an integer equal to the size of the floating point numbers used in the
-file (usually, @var{data-size} = sizeof(double)).
+file (currently only @var{data-size} = sizeof(double) is supported).
 
 @item @var{number-of-nodes}
 is the number of nodes in the mesh.
@@ -3114,6 +3117,107 @@ ordering of these additional nodes follows the ordering of the edges and
 quadrangular faces given in @ref{Gmsh node ordering}.
 @end table
 
+@c .........................................................................
+@c Version 2.0 BINARY
+@c .........................................................................
+
+@node Version 2.0 Binary,  , Version 2.0 ASCII, Gmsh mesh file formats
+@subsection Version 2.0 Binary
+
+The binary file format is similar to the ASCII format described in
+@ref{Version 2.0 ASCII}:
+
+@example
+$MeshFormat
+2.0 @var{file-type} @var{data-size}
+@var{one-binary}
+$EndMeshFormat
+$Nodes
+@var{number-of-nodes}
+@var{nodes-binary}
+$EndNodes
+$Elements
+@var{number-of-elements}
+@var{element-header-binary}
+@var{elements-binary}
+@var{element-header-binary}
+@var{elements-binary}
+@dots{}
+$EndElements
+@end example
+
+@noindent
+where
+@table @code
+@item @var{file-type}
+is an integer equal to 1.
+@item @var{data-size}
+has the same meaning as in the ASCII file format
+@item @var{one-binary}
+is an integer of value 1 written in binary form. This integer is used
+for detecting if the computer on which the binary file was written and
+the computer on which the file is read are of the same type (little or
+big endian).
+
+Here is a pseudo C code to write @var{one-binary}:
+@example
+int one = 1;
+fwrite(&one, sizeof(int), 1, file);
+@end example
+
+@item @var{number-of-nodes}
+has the same meaning as in the ASCII file format
+
+@item @var{nodes-binary}
+is the list of nodes in binary form, i.e., a array of
+@var{number-of-nodes} * (4 + 3 * @var{data-size}) bytes. For each node,
+the first 4 bytes contain the node number, the next (3 *
+@var{data-size}) bytes contain the three floating point coordinates
+
+Here is a pseudo C to write @var{nodes-binary}:
+@example
+for(i = 0; i < number_of_nodes; i++)@{
+  fwrite(&num, sizeof(int), 1, file);
+  double xyz[3] = @{node_i_x, node_i_y, node_i_z@};
+  fwrite(&xyz, sizeof(double), 3, file);
+@}
+@end example
+
+@item @var{number-of-elements}
+has the same meaning as in the ASCII file format
+
+@item @var{element-header-binary}
+is a list of 3 integers in binary form, i.e., an array of (3 * 4) bytes:
+the first four contain the type of the elements that follow (same as
+@var{elm-type} in the ASCII format), the next four contain the number of
+elements that follow, and the last four contain the number of tags per
+element (same as @var{number-of-tags} in the ASCII format).
+
+Here is a pseudo C code to write @var{element-header-binary}:
+@example
+int header[3] = @{elm_type, num_elm_follow, num_tags@};
+fwrite(&header, sizeof(int), 3, file);
+@end example
+
+@item @var{elements-binary}
+is a list of elements in binary form, i.e., an array of ``number of
+elements that follow'' * (4 + @var{number-of-tags} * 4 +
+#@var{node-number-list} * 4) bytes. For each element, the first four
+bytes contain the element number, the next (@var{number-of-tags} * 4)
+contain the tags, and the last (#@var{node-number-list} * 4) contain the
+node indices.
+
+Here is a pseudo C to write @var{elements-binary} for triangles with 3
+tags:
+@example
+for(i = 0; i < number_of_triangles; i++)@{
+  int data[7] = @{num, physical, elementary, partition, 
+                 node_1, node_2, node_3@};
+  fwrite(data, sizeof(int), 7, file);
+@}
+@end example
+@end table
+
 @c -------------------------------------------------------------------------
 @c Gmsh post-processing file formats
 @c -------------------------------------------------------------------------