diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index b0ade62bafd5fd1288b1158e80c4314b052283ff..ae9e4b86c40ed71cbdae59eb6978ee2b2b1ac84b 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1284,7 +1284,7 @@ StringXNumber PostProcessingOptions_Number[] = {
 
   { F|O, "Format" , opt_post_file_format , 0. ,
     "Default file format for post-processing views (0=ASCII view, 1=binary "
-    "view, 2=parsed view, 3=STL triangulation, 4=text, 5=mesh)" },
+    "view, 2=parsed view, 3=STL triangulation, 4=raw text, 5=Gmsh mesh, 6=MED file)" },
 
   { F|O, "HorizontalScales" , opt_post_horizontal_scales , 1. , 
     "Display value scales horizontally" },
diff --git a/Common/DummyBindings.h b/Common/DummyBindings.h
index 731c61a7626ede07869c00c598741f303f49faaa..890137a59f8b3c07aa959bbbf61e2a49817bc2dc 100644
--- a/Common/DummyBindings.h
+++ b/Common/DummyBindings.h
@@ -1,31 +1,29 @@
 #ifndef _DUMMY_BINDINGS_H_
 #define _DUMMY_BINDINGS_H_
+
 class classBinding {
-public:
+ public:
   void setDescription(std::string description){};
   template<typename parentType>
   void setParentClass(){}
   template <typename cb>
-  methodBinding *addMethod(std::string n, cb f){
-    return new methodBinding();
-  }
-  template <typename tObj, typename t0, typename t1, typename t2, typename t3 >
-  methodBinding *setConstructor(){}
-  template <typename tObj, typename t0, typename t1, typename t2 >
-  methodBinding *setConstructor(){}
+  methodBinding *addMethod(std::string n, cb f){ return new methodBinding(); }
+  template <typename tObj, typename t0, typename t1, typename t2, typename t3>
+  methodBinding *setConstructor(){ return 0; }
+  template <typename tObj, typename t0, typename t1, typename t2>
+  methodBinding *setConstructor(){ return 0; }
   template <typename tObj, typename t0, typename t1>
-  methodBinding *setConstructor(){return 0;}
+  methodBinding *setConstructor(){ return 0; }
   template <typename tObj, typename t0>
-  methodBinding *setConstructor(){return 0;}
+  methodBinding *setConstructor(){ return 0; }
   template<typename tObj>
-  methodBinding *setConstructor(){return 0;}
+  methodBinding *setConstructor(){ return 0; }
 };
 
 class binding {
-  public:
+ public:
   template<class t>
-  classBinding *addClass(std::string className){
-    return new classBinding();
-  }
+  classBinding *addClass(std::string className){ return new classBinding(); }
 };
+
 #endif
diff --git a/Common/LuaBindings.h b/Common/LuaBindings.h
index 290387e0bc576083568332d65d602741b8ef282c..b9d63d26a9f864c75c8c9e20b2aacf25e75fa0d5 100644
--- a/Common/LuaBindings.h
+++ b/Common/LuaBindings.h
@@ -1,6 +1,8 @@
 #ifndef _LUA_BINDINGS_H_
 #define _LUA_BINDINGS_H_
+
 #include "Bindings.h"
+
 #ifdef HAVE_LUA
 #include "BindingsDocTemplates.h"
 
diff --git a/Fltk/fileDialogs.cpp b/Fltk/fileDialogs.cpp
index ded35a1915cbf6c8c3942624604063297acdf951..266916a6261a9a9565143e792c12307070035e55 100644
--- a/Fltk/fileDialogs.cpp
+++ b/Fltk/fileDialogs.cpp
@@ -773,7 +773,7 @@ int mshFileDialog(const char *name)
       if (!o) break;
       if (o == dialog->ok) {
         opt_mesh_msh_file_version(0, GMSH_SET | GMSH_GUI, 
-                                  (dialog->c->value() == 0) ? 1.0 : 2.1);
+                                  (dialog->c->value() == 0) ? 1.0 : 2.2);
         opt_mesh_binary(0, GMSH_SET | GMSH_GUI, (dialog->c->value() == 2) ? 1 : 0);
         opt_mesh_save_all(0, GMSH_SET | GMSH_GUI, dialog->b->value() ? 1 : 0);
         opt_mesh_save_parametric(0, GMSH_SET | GMSH_GUI, dialog->p->value() ? 1 : 0);
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index bed591cf43fb66d6658d0bf0d651b3531268dc8d..0a9565d5eb8a41d74b3b98d3b0a61e76f825fa80 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -316,7 +316,7 @@ int GModel::readMSH(const std::string &name)
         if(numVertices > 100000) 
           Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
       }
-      // If the vertex numbering is dense, tranfer the map into a
+      // 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) ||
@@ -344,7 +344,8 @@ int GModel::readMSH(const std::string &name)
       Msg::ResetProgressMeter();
       if(!binary){
         for(int i = 0; i < numElements; i++) {
-          int num, type, physical = 0, elementary = 0, partition = 0, parentN = 0, numVertices;
+          int num, type, physical = 0, elementary = 0, partition = 0, parent = 0, numVertices;
+          std::vector<short> ghosts;
           if(version <= 1.0){
             if(fscanf(fp, "%d %d %d %d %d", &num, &type, &physical, &elementary, &numVertices) != 5)
               return 0;
@@ -353,14 +354,18 @@ int GModel::readMSH(const std::string &name)
           else{
             int numTags;
             if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
+            std::vector<int> tags(numTags);
+            int numPartitions = 0;
             for(int j = 0; j < numTags; j++){
               int tag;
               if(fscanf(fp, "%d", &tag) != 1) return 0;
-              if(j == 0)      physical = tag;
+              if(j == 0) physical = tag;
               else if(j == 1) elementary = tag;
-              else if(j == 2) partition = tag;
-              else if(j == 3) parentN = tag;
-              // ignore any other tags for now
+              else if(version < 2.2 && j == 2) partition = tag;
+              else if(version >= 2.2 && j == 2) numPartitions = tag;
+              else if(version >= 2.2 && j == 3) partition = tag;
+              else if(j >= 4 && j < 4 + numPartitions) ghosts.push_back(-tag);
+              else if(j == numTags - 1) parent = tag;
             }
             if(!(numVertices = MElement::getInfoMSH(type))) {
               if(type != MSH_POLYG_ && type != MSH_POLYH_) return 0;
@@ -379,15 +384,17 @@ int GModel::readMSH(const std::string &name)
           }
           MElement *e = createElementMSH(this, num, type, physical, elementary,
                                          partition, vertices, elements, physicals);
-          if(parentN) {
-            MElement *parent = NULL;
-            if(parents.find(parentN) == parents.end()){
-              parent = getParent(parentN, e->getDim(), elements);
-              if(parent) parents[parentN] = parent;
-              else Msg::Error("Parent element %d not found!",parentN);
+          for(unsigned int j = 0; j < ghosts.size(); j++)
+            _ghostCells.insert(std::pair<MElement*, short>(e, ghosts[j]));
+          if(parent) {
+            MElement *p = 0;
+            if(parents.find(parent) == parents.end()){
+              p = getParent(parent, e->getDim(), elements);
+              if(p) parents[parent] = p;
+              else Msg::Error("Parent element %d not found", parent);
             }
-            else parent = parents.find(parentN)->second;
-            e->setParent(parent);
+            else p = parents.find(parent)->second;
+            e->setParent(p);
           }
           if(numElements > 100000) 
             Msg::ProgressMeter(i + 1, numElements, "Reading elements");
@@ -410,9 +417,11 @@ int GModel::readMSH(const std::string &name)
             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 physical = (numTags > 0) ? data[1] : 0;
+            int elementary = (numTags > 1) ? data[2] : 0;
+            int numPartitions = (version >= 2.2 && numTags > 2) ? data[3] : 0;
+            int partition = (version < 2.2 && numTags > 2) ? data[3] : 
+              (version >= 2.2 && numTags > 3) ? data[4] : 0;
             int *indices = &data[numTags + 1];
             std::vector<MVertex*> vertices;
             if(vertexVector.size()){
@@ -421,8 +430,11 @@ int GModel::readMSH(const std::string &name)
             else{
               if(!getVertices(numVertices, indices, vertexMap, vertices)) return 0;
             }
-            createElementMSH(this, num, type, physical, elementary, partition, 
-                             vertices, elements, physicals);
+            MElement *e = createElementMSH(this, num, type, physical, elementary,
+                                           partition, vertices, elements, physicals);
+            if(numPartitions > 1)
+              for(int j = 0; j < numPartitions - 1; j++)
+                _ghostCells.insert(std::pair<MElement*, short>(e, -data[5 + j]));
             if(numElements > 100000)
               Msg::ProgressMeter(numElementsPartial + i + 1, numElements, 
                                  "Reading elements");
@@ -480,85 +492,50 @@ int GModel::readMSH(const std::string &name)
   return postpro ? 2 : 1;
 }
 
-static void writeElementHeaderMSH(bool binary, FILE *fp, std::map<int,int> &elements,
-                                  int t1, int t2=0, int t3=0, int t4=0, 
-                                  int t5=0, int t6=0, int t7=0, int t8=0, int t9=0)
-{
-  if(!binary) return;
-
-  int numTags = 3;
-  int data[3];
-  if(elements.count(t1)){
-    data[0] = t1;  data[1] = elements[t1];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t2 && elements.count(t2)){
-    data[0] = t2;  data[1] = elements[t2];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t3 && elements.count(t3)){
-    data[0] = t3;  data[1] = elements[t3];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t4 && elements.count(t4)){
-    data[0] = t4;  data[1] = elements[t4];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t5 && elements.count(t5)){
-    data[0] = t5;  data[1] = elements[t5];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t6 && elements.count(t6)){
-    data[0] = t6;  data[1] = elements[t6];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t7 && elements.count(t7)){
-    data[0] = t7;  data[1] = elements[t7];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t8 && elements.count(t8)){
-    data[0] = t8;  data[1] = elements[t8];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-  else if(t9 && elements.count(t9)){
-    data[0] = t9;  data[1] = elements[t9];  data[2] = numTags;
-    fwrite(data, sizeof(int), 3, fp);
-  }
-}
-
 template<class T>
-static void writeElementMSH(FILE *fp, T *ele, bool saveAll, 
+static void writeElementMSH(FILE *fp, GModel *model, T *ele, bool saveAll, 
                             double version, bool binary, int &num, int elementary, 
                             std::vector<int> &physicals, int parentNum)
 {
+  std::vector<short> ghosts;
+  if(model->getGhostCells().size()){
+    std::pair<std::multimap<MElement*, short>::iterator,
+              std::multimap<MElement*, short>::iterator> itp =
+      model->getGhostCells().equal_range(ele);
+    for(std::multimap<MElement*, short>::iterator it = itp.first;
+        it != itp.second; it++)
+      ghosts.push_back(it->second);
+  }
+  
   if(saveAll)
-    ele->writeMSH(fp, version, binary, ++num, elementary, 0, parentNum);
+    ele->writeMSH(fp, version, binary, ++num, elementary, 0, parentNum, &ghosts);
   else
     for(unsigned int j = 0; j < physicals.size(); j++)
-      ele->writeMSH(fp, version, binary, ++num, elementary, physicals[j], parentNum);
+      ele->writeMSH(fp, version, binary, ++num, elementary, physicals[j],
+                    parentNum, &ghosts);
 }
 
 template<class T>
-static void writeElementsMSH(FILE *fp, std::vector<T*> &ele, bool saveAll,
-                             double version, bool binary, int &num, int elementary,
-                             std::vector<int> &physicals)
+static void writeElementsMSH(FILE *fp, GModel *model, std::vector<T*> &ele, 
+                             bool saveAll, double version, bool binary, int &num, 
+                             int elementary, std::vector<int> &physicals)
 {
   for(unsigned int i = 0; i < ele.size(); i++)
-    writeElementMSH(fp, ele[i], saveAll, version, binary, num,
+    writeElementMSH(fp, model, ele[i], saveAll, version, binary, num,
                     elementary, physicals, 0);
 }
 
 template<class T>
-static void writeElementsMSH(FILE *fp, std::vector<T*> &ele, bool saveAll,
-                             double version, bool binary, int &num, int elementary,
-                             std::vector<int> &physicals,
+static void writeElementsMSH(FILE *fp, GModel *model, std::vector<T*> &ele,
+                             bool saveAll, double version, bool binary, int &num,
+                             int elementary, std::vector<int> &physicals,
                              std::map<MElement*, int> &parentsNum)
 {
   for(unsigned int i = 0; i < ele.size(); i++) {
     int parentNum = (ele[i]->getParent() &&
                      parentsNum.find(ele[i]->getParent()) != parentsNum.end()) ?
       parentsNum.find(ele[i]->getParent())->second : 0;
-    writeElementMSH(fp, ele[i], saveAll, version, binary, num,
+    writeElementMSH(fp, model, ele[i], saveAll, version, binary, num,
                     elementary, physicals, parentNum);
   }
 }
@@ -577,70 +554,43 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
 
   // get the number of vertices and index the vertices in a continuous
   // sequence
-  bool renumber = true; // FIXME
-  int numVertices = renumber ? indexMeshVertices(saveAll) : getNumMeshVertices();
+  int numVertices = indexMeshVertices(saveAll);
   
   // binary format exists only in version 2
   if(version > 1 || binary) 
-    version = 2.1;
+    version = 2.2;
   else
     version = 1.0;
-
-  // 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;
-  std::map<MElement *, GEntity *> parents[2];
-  int numParents = 0;
-  for(viter it = firstVertex(); it != lastVertex(); ++it){
-    int p = (saveAll ? 1 : (*it)->physicals.size());
-    int n = p * (*it)->points.size();
-    if(n) elements[(*it)->points[0]->getTypeForMSH()] += 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;
-  }
+  
+  // get the number of elements we need to save
+  int numElements = 0;
+  std::map<MElement*, GEntity*> parents[2];
+  for(viter it = firstVertex(); it != lastVertex(); ++it)
+    numElements += (saveAll ? 1 : (*it)->physicals.size()) * (*it)->points.size();
+  for(eiter it = firstEdge(); it != lastEdge(); ++it)
+    numElements += (saveAll ? 1 : (*it)->physicals.size()) * (*it)->lines.size();
   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;
-    n = p * (*it)->polygons.size();
-    if(n) elements[(*it)->polygons[0]->getTypeForMSH()] += n;
-    for(unsigned int i = 0; i < (*it)->polygons.size(); i++)
+    numElements += (saveAll ? 1 : (*it)->physicals.size()) * 
+      ((*it)->triangles.size() + (*it)->quadrangles.size() + (*it)->polygons.size());
+    for(unsigned int i = 0; i < (*it)->polygons.size(); i++){
       if((*it)->polygons[i]->ownsParent()) {
         parents[0][(*it)->polygons[i]->getParent()] = (*it);
-        numParents += p;
+        numElements += (saveAll ? 1 : (*it)->physicals.size());
       }
+    }
   }
   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;
-    n = p * (*it)->polyhedra.size();
-    if(n) elements[(*it)->polyhedra[0]->getTypeForMSH()] += n;
-    for(unsigned int i = 0; i < (*it)->polyhedra.size(); i++)
+    numElements += (saveAll ? 1 : (*it)->physicals.size()) *
+      ((*it)->tetrahedra.size() + (*it)->hexahedra.size() + (*it)->prisms.size() +
+       (*it)->pyramids.size() + (*it)->polyhedra.size());
+    for(unsigned int i = 0; i < (*it)->polyhedra.size(); i++){
       if((*it)->polyhedra[i]->ownsParent()) {
         parents[1][(*it)->polyhedra[i]->getParent()] = (*it);
-        numParents += p;
+        numElements += (saveAll ? 1 : (*it)->physicals.size());
       }
+    }
   }
 
-  int numElements = 0;
-  std::map<int, int>::const_iterator it = elements.begin();
-  for(; it != elements.end(); ++it)
-    numElements += it->second;
-  numElements += numParents;
-
   if(version >= 2.0){
     fprintf(fp, "$MeshFormat\n");
     fprintf(fp, "%g %d %d\n", version, binary ? 1 : 0, (int)sizeof(double));
@@ -676,9 +626,9 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
     for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
       entities[i]->mesh_vertices[j]->writeMSH(fp, binary, saveParametric, 
                                               scalingFactor);
-
+  
   if(binary) fprintf(fp, "\n");
-
+  
   if(version >= 2.0){
     if(saveParametric)
       fprintf(fp, "$EndParametricNodes\n");
@@ -693,121 +643,112 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
 
   fprintf(fp, "%d\n", numElements);
   int num = 0;
-
   std::map<MElement*, int> parentsNum;
-  std::map<MElement*, GEntity*>::const_iterator itP;
 
   // points
-  writeElementHeaderMSH(binary, fp, elements, MSH_PNT);
   for(viter it = firstVertex(); it != lastVertex(); ++it)
-    writeElementsMSH(fp, (*it)->points, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->points, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
 
   // lines
-  writeElementHeaderMSH(binary, fp, elements, MSH_LIN_2, MSH_LIN_3, MSH_LIN_4,
-                        MSH_LIN_5);
   for(eiter it = firstEdge(); it != lastEdge(); ++it)
-   writeElementsMSH(fp, (*it)->lines, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->lines, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
-
+  
   // triangles
-  writeElementHeaderMSH(binary, fp, elements, MSH_TRI_3, MSH_TRI_6, MSH_TRI_9,
-                        MSH_TRI_10, MSH_TRI_12, MSH_TRI_15, MSH_TRI_15I, MSH_TRI_21);
-  for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
-    if(itP->first->getType() == TYPE_TRI) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin(); 
+      it != parents[0].end(); it++)
+    if(it->first->getType() == TYPE_TRI) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(fiter it = firstFace(); it != lastFace(); ++it)
-    writeElementsMSH(fp, (*it)->triangles, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->triangles, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
 
   // quads
-  writeElementHeaderMSH(binary, fp, elements, MSH_QUA_4, MSH_QUA_9, MSH_QUA_8, 
-                        MSH_QUA_16, MSH_QUA_25, MSH_QUA_36, MSH_QUA_12, MSH_QUA_16I,
-                        MSH_QUA_20);
-  for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
-    if(itP->first->getType() == TYPE_QUA) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin(); 
+      it != parents[0].end(); it++)
+    if(it->first->getType() == TYPE_QUA) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(fiter it = firstFace(); it != lastFace(); ++it)
-    writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->quadrangles, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
 
   // polygons
-  writeElementHeaderMSH(binary, fp, elements, MSH_POLYG_);
-  for(itP = parents[0].begin(); itP != parents[0].end(); itP++)
-    if(itP->first->getType() == TYPE_POLYG) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[0].begin();
+      it != parents[0].end(); it++)
+    if(it->first->getType() == TYPE_POLYG) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(fiter it = firstFace(); it != lastFace(); it++)
-    writeElementsMSH(fp, (*it)->polygons, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->polygons, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals, parentsNum);
 
   // tets
-  writeElementHeaderMSH(binary, fp, elements, MSH_TET_4, MSH_TET_10, MSH_TET_20, 
-                        MSH_TET_35, MSH_TET_56, MSH_TET_52);
-  for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
-    if(itP->first->getType() == TYPE_TET) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+      it != parents[1].end(); it++)
+    if(it->first->getType() == TYPE_TET) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->tetrahedra, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
 
   // hexas
-  writeElementHeaderMSH(binary, fp, elements, MSH_HEX_8, MSH_HEX_27, MSH_HEX_20);
-  for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
-    if(itP->first->getType() == TYPE_HEX) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+      it != parents[1].end(); it++)
+    if(it->first->getType() == TYPE_HEX) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    writeElementsMSH(fp, (*it)->hexahedra, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->hexahedra, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
 
   // prisms
-  writeElementHeaderMSH(binary, fp, elements, MSH_PRI_6, MSH_PRI_18, MSH_PRI_15);
-  for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
-    if(itP->first->getType() == TYPE_PRI) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+      it != parents[1].end(); it++)
+    if(it->first->getType() == TYPE_PRI) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    writeElementsMSH(fp, (*it)->prisms, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->prisms, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
-
+  
   // pyramids
-  writeElementHeaderMSH(binary, fp, elements, MSH_PYR_5, MSH_PYR_14, MSH_PYR_13);
-  for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
-    if(itP->first->getType() == TYPE_PYR) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+      it != parents[1].end(); it++)
+    if(it->first->getType() == TYPE_PYR) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    writeElementsMSH(fp, (*it)->pyramids, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->pyramids, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals);
 
   // polyhedra
-  writeElementHeaderMSH(binary, fp, elements, MSH_POLYH_);
-  for(itP = parents[1].begin(); itP != parents[1].end(); itP++)
-    if(itP->first->getType() == TYPE_POLYH) {
-      writeElementMSH(fp, itP->first, saveAll, version, binary, num,
-                      itP->second->tag(), itP->second->physicals, 0);
-      parentsNum[itP->first] = num;
+  for(std::map<MElement*, GEntity*>::const_iterator it = parents[1].begin();
+      it != parents[1].end(); it++)
+    if(it->first->getType() == TYPE_POLYH) {
+      writeElementMSH(fp, this, it->first, saveAll, version, binary, num,
+                      it->second->tag(), it->second->physicals, 0);
+      parentsNum[it->first] = num;
     }
   for(riter it = firstRegion(); it != lastRegion(); ++it)
-    writeElementsMSH(fp, (*it)->polyhedra, saveAll, version, binary, num,
+    writeElementsMSH(fp, this, (*it)->polyhedra, saveAll, version, binary, num,
                      (*it)->tag(), (*it)->physicals, parentsNum);
   
   if(binary) fprintf(fp, "\n");
@@ -820,6 +761,20 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
   }
 
   fclose(fp);
+
+#if 1
+  if(_ghostCells.size()){
+    Msg::Info("Wrinting ghost cells in debug file");
+    fp = fopen("ghosts.pos", "w");
+    fprintf(fp, "View \"ghosts\"{\n");
+    for(std::multimap<MElement*, short>::iterator it = _ghostCells.begin();
+        it != _ghostCells.end(); it++)
+      it->first->writePOS(fp, false, true, false, false, false, false);
+    fprintf(fp, "};\n");
+    fclose(fp);
+  }
+#endif
+
   return 1;
 }
 
@@ -882,7 +837,7 @@ int GModel::writeDistanceMSH(const std::string &name, double version, bool binar
   for(std::map<MVertex*, double>::iterator it = distance.begin(); it != distance.end(); it++) {
       fprintf(f,"SP(%g,%g,%g){%g};\n",
 	      it->first->x(), it->first->y(), it->first->z(),
-	      it->second);      
+	      it->second);     
   }
   fprintf(f,"};\n");
   fclose(f);
@@ -929,7 +884,7 @@ int GModel::writePOS(const std::string &name, bool printElementary,
   bool f[6] = {printElementary, printElementNumber, printGamma, printEta, printRho,
                printDisto};
 
-  bool first = true;  
+  bool first = true; 
   std::string names;
   if(f[0]){
     if(first) first = false; else names += ",";
@@ -1412,7 +1367,7 @@ int GModel::readUNV(const std::string &name)
   while(!feof(fp)) {
     if(!fgets(buffer, sizeof(buffer), fp)) break;
     if(!strncmp(buffer, "    -1", 6)){
-      if(!fgets(buffer, sizeof(buffer), fp)) break;      
+      if(!fgets(buffer, sizeof(buffer), fp)) break;     
       if(!strncmp(buffer, "    -1", 6))
         if(!fgets(buffer, sizeof(buffer), fp)) break;
       int record = 0;
@@ -1566,7 +1521,7 @@ int GModel::writeUNV(const std::string &name, bool saveAll, bool saveGroupsOfNod
   for(unsigned int i = 0; i < entities.size(); i++)
     for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++) 
       entities[i]->mesh_vertices[j]->writeUNV(fp, scalingFactor);
-  fprintf(fp, "%6d\n", -1);  
+  fprintf(fp, "%6d\n", -1); 
 
   // elements
   fprintf(fp, "%6d\n", -1);
@@ -1873,7 +1828,7 @@ int GModel::writeFEA(const std::string &name, int elementTagType,
 	fprintf(fp,"%d %g %g %g\n", entities[i]->mesh_vertices[j]->getIndex(),
 		entities[i]->mesh_vertices[j]->x() * scalingFactor,
 		entities[i]->mesh_vertices[j]->y() * scalingFactor,
-		entities[i]->mesh_vertices[j]->z() * scalingFactor);  
+		entities[i]->mesh_vertices[j]->z() * scalingFactor);
 
   int iElement = 1;
   for(fiter it = firstFace(); it != lastFace(); ++it){
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 2d475aaff14bb80c5dfd3c3b1ad7ad0d167bb5b1..4dad3a651f437a0651e61f7a3d4a81e00e26f3fd 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -372,8 +372,8 @@ void MElement::interpolateCurl(double val[], double u, double v, double w, doubl
   f[2] = fy[0] - fx[1];
 }
 
-double MElement::interpolateDiv(double val[], double u, double v, double w, int stride,
-                                int order)
+double MElement::interpolateDiv(double val[], double u, double v, double w, 
+                                int stride, int order)
 {
   double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
   getJacobian(u, v, w, jac);
@@ -385,7 +385,8 @@ double MElement::interpolateDiv(double val[], double u, double v, double w, int
 }
 
 void MElement::writeMSH(FILE *fp, double version, bool binary, int num, 
-                        int elementary, int physical, int parentNum)
+                        int elementary, int physical, int parentNum,
+                        std::vector<short> *ghosts)
 {
   int type = getTypeForMSH();
 
@@ -400,12 +401,36 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
     fprintf(fp, "%d %d", num ? num : _num, type);
     if(version < 2.0)
       fprintf(fp, " %d %d %d", abs(physical), elementary, n);
-    else
-      fprintf(fp, " 3 %d %d %d", abs(physical), elementary, _partition);
+    else if(!_partition)
+      fprintf(fp, " 2 %d %d", abs(physical), elementary);
+    else if(!ghosts)
+      fprintf(fp, " 4 %d %d 1 %d", abs(physical), elementary, _partition);
+    else{
+      int numGhosts = ghosts->size();
+      fprintf(fp, " %d %d %d %d %d", 4 + numGhosts, abs(physical), elementary, 
+              1 + numGhosts, _partition);
+      for(unsigned int i = 0; i < ghosts->size(); i++)
+        fprintf(fp, " %d", -(*ghosts)[i]);
+    }
   }
   else{
-    int tags[4] = {num ? num : _num, abs(physical), elementary, _partition};
-    fwrite(tags, sizeof(int), 4, fp);
+    int numTags, numGhosts = 0;
+    if(!_partition) numTags = 2;
+    else if(!ghosts) numTags = 4;
+    else{
+      numGhosts = ghosts->size();
+      numTags = 4 + numGhosts;
+    }
+    // we write elements in blobs of single elements; this will lead
+    // to suboptimal reads, but it's much simpler when the number of
+    // tags change from element to element (third-party codes can
+    // still write MSH file optimized for reading speed, by grouping
+    // elements with the same number of tags in blobs)
+    int blob[60] = {type, 1, numTags, num ? num : _num, abs(physical), elementary, 
+                    1 + numGhosts, _partition};
+    if(ghosts)
+      for(unsigned int i = 0; i < numGhosts; i++) blob[8 + i] = -(*ghosts)[i];
+    fwrite(blob, sizeof(int), 4 + numTags, fp);
   }
 
   if(physical < 0) revert();
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 45f3ffdca3eda10045fe12e73e6bfab4573b0423..2b6c803f7d8a5534f16ad33776848f52b8ecadb9 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -226,8 +226,9 @@ class MElement
   }
 
   // IO routines
-  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, 
-                        int num=0, int elementary=1, int physical=1, int parentNum=0);
+  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
+                        int num=0, int elementary=1, int physical=1,
+                        int parentNum=0, std::vector<short> *ghosts=0);
   virtual void writePOS(FILE *fp, bool printElementary, bool printElementNumber, 
                         bool printGamma, bool printEta, bool printRho, 
                         bool printDisto,double scalingFactor=1.0, int elementary=1);
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index d121562d116138bcf7a3d3107246d8e94cfcb918..f82d2ba7f41fea6152075ce065d4002e821b388c 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -137,7 +137,8 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 }
 
 void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
-                           int elementary, int physical, int parentNum)
+                           int elementary, int physical, int parentNum,
+                           std::vector<short> *ghosts)
 {
   int type = getTypeForMSH();
 
@@ -156,14 +157,13 @@ void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
       fprintf(fp, " %d %d", abs(physical), elementary);
     else {
       if(parentNum)
-        fprintf(fp, " 4 %d %d %d %d", abs(physical), elementary, partE, parentNum);
+        fprintf(fp, " 5 %d %d 1 %d %d", abs(physical), elementary, partE, parentNum);
       else
-        fprintf(fp, " 3 %d %d %d", abs(physical), elementary, partE);
+        fprintf(fp, " 4 %d %d 1 %d", abs(physical), elementary, partE);
     }
   }
   else{
-    int tags[4] = {num ? num : numE, abs(physical), elementary, partE};
-    fwrite(tags, sizeof(int), 4, fp);
+    Msg::Error("Binary output not coded for polyhedra");
   }
 
   if(physical < 0) revert();
@@ -275,7 +275,8 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 }
 
 void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
-                        int elementary, int physical, int parentNum)
+                        int elementary, int physical, int parentNum,
+                        std::vector<short> *ghosts)
 {
   int type = getTypeForMSH();
 
@@ -294,14 +295,13 @@ void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
       fprintf(fp, " %d %d", abs(physical), elementary);
     else {
       if(parentNum)
-        fprintf(fp, " 4 %d %d %d %d", abs(physical), elementary, partE, parentNum);
+        fprintf(fp, " 5 %d %d 1 %d %d", abs(physical), elementary, partE, parentNum);
       else
-        fprintf(fp, " 3 %d %d %d", abs(physical), elementary, partE);
+        fprintf(fp, " 4 %d %d 1 %d", abs(physical), elementary, partE);
     }
   }
   else{
-    int tags[4] = {num ? num : numE, abs(physical), elementary, partE};
-    fwrite(tags, sizeof(int), 4, fp);
+    Msg::Error("Binary output not coded for polygons");
   }
 
   if(physical < 0) revert();
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index a9f355c15d4665256e53526cbd90b13e2430b16f..c4a59ceb093a3cd3337f82500a847a803b36b5f7 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -123,8 +123,9 @@ class MPolyhedron : public MElement {
   }
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
-  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, 
-                        int num=0, int elementary=1, int physical=1, int parentNum=0);
+  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
+                        int num=0, int elementary=1, int physical=1,
+                        int parentNum=0, std::vector<short> *ghosts=0);
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p) { _orig = p; }
   virtual int getNumChildren() const { return _parts.size(); }
@@ -225,8 +226,9 @@ class MPolygon : public MElement {
   }
   virtual bool isInside(double u, double v, double w);
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
-  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false, 
-                        int num=0, int elementary=1, int physical=1, int parentNum=0);
+  virtual void writeMSH(FILE *fp, double version=1.0, bool binary=false,
+                        int num=0, int elementary=1, int physical=1,
+                        int parentNum=0, std::vector<short> *ghosts=0);
   virtual MElement *getParent() const { return _orig; }
   virtual void setParent(MElement *p) { _orig = p; }
   virtual int getNumChildren() const { return _parts.size(); }
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 56aeef2bbe49e7405398512f8b06e2b17a0c8d60..3b9cf5e2385764cf721d3ef234d1baa709a22da9 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -730,16 +730,15 @@ struct MatchBoElemToGrVertex
   }
 };
 
-
 template <class ITERATOR>
 void fillit_(std::multimap<MFace, MElement*, Less_Face> &faceToElement,
 	     ITERATOR it_beg, ITERATOR it_end)
 {
   for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
     MElement *el = *IT;
-    for(int j=0;j < el->getNumFaces();j++) {
+    for(int j = 0; j < el->getNumFaces(); j++) {
       MFace e = el->getFace(j);
-      faceToElement.insert(std::make_pair(e,el));
+      faceToElement.insert(std::make_pair(e, el));
     }
   }
 }
@@ -748,11 +747,11 @@ template <class ITERATOR>
 void fillit_(std::multimap<MEdge, MElement*, Less_Edge> &edgeToElement,
 	     ITERATOR it_beg, ITERATOR it_end)
 {
-  for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
+  for (ITERATOR IT = it_beg; IT != it_end; ++IT){
     MElement *el = *IT;
-    for(int j=0;j < el->getNumEdges();j++) {
+    for(int j = 0; j < el->getNumEdges(); j++) {
       MEdge e = el->getEdge(j);
-      edgeToElement.insert(std::make_pair(e,el));
+      edgeToElement.insert(std::make_pair(e, el));
     }
   }
 }
@@ -763,15 +762,14 @@ void fillit_(std::multimap<MVertex*, MElement*> &vertexToElement,
 {
   for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
     MElement *el = *IT;
-    for(int j=0;j < el->getNumVertices();j++) {
+    for(int j = 0; j < el->getNumVertices(); j++) {
       MVertex* e = el->getVertex(j);
-      vertexToElement.insert(std::make_pair(e,el));
+      vertexToElement.insert(std::make_pair(e, el));
     }
   }
 }
 
-void assignPartitionBoundary(GModel *model,
-			     MFace &me,
+void assignPartitionBoundary(GModel *model, MFace &me,
 			     std::set<partitionFace*, Less_partitionFace> &pfaces,
 			     std::vector<MElement*> &v)
 {
@@ -1117,15 +1115,6 @@ int CreatePartitionBoundaries(GModel *model, bool createGhostCells)
       for(std::set<partitionFace*, Less_partitionFace>::iterator it = pfaces.begin();
           it != pfaces.end(); it++)
         addGhostCells(*it, vertexToElement, ghosts);
-#if 1
-    FILE *fp = fopen("ghosts.pos", "w");
-    fprintf(fp, "View \"ghosts\"{\n");
-    for(std::multimap<MElement*, short>::iterator it = ghosts.begin();
-        it != ghosts.end(); it++)
-      it->first->writePOS(fp, false, true, false, false, false, false);
-    fprintf(fp, "};\n");
-    fclose(fp);
-#endif
   }
 
   return 1;
diff --git a/doc/texinfo/gmsh.texi b/doc/texinfo/gmsh.texi
index 31503f03f605763c57cb2a3964f1284e785bb565..06e227e8c85fad858be92a56658a9140694a5495 100644
--- a/doc/texinfo/gmsh.texi
+++ b/doc/texinfo/gmsh.texi
@@ -3399,7 +3399,7 @@ for some examples.
 This chapter describes Gmsh's native ``MSH'' file format, used to store
 meshes and associated post-processing datasets. The MSH format exists in
 two flavors: ASCII and binary. The format has a version number
-(currently: 2.1) that is independent of Gmsh's main version number.
+(currently: 2.2) that is independent of Gmsh's main version number.
 
 (Remember that for small post-processing datasets you can also use
 human-readable ``parsed'' post-processing views, as described in
@@ -3505,7 +3505,7 @@ $ElementEndNodeData
 where
 @table @code
 @item @var{version-number}
-is a real number equal to 2.1
+is a real number equal to 2.2
 
 @item @var{file-type}
 is an integer equal to 0 in the ASCII file format.
@@ -3632,8 +3632,8 @@ gives the number of integer tags that follow for the @var{n}-th
 element. By default, the first @var{tag} is the number of the physical
 entity to which the element belongs; the second is the number of the
 elementary geometrical entity to which the element belongs; the third is
-the number of a mesh partition to which the element belongs. All tags
-must be postive integers, or zero. A zero tag is equivalent to no tag.
+the number of mesh partitions to which the element belongs, followed by
+the partition ids. A zero tag is equivalent to no tag.
 
 @item @var{node-number-list}
 is the list of the node numbers of the @var{n}-th element. The ordering of
@@ -3675,7 +3675,7 @@ file!):
 
 @smallexample
 $MeshFormat
-2.1 0 8
+2.2 0 8
 $EndMeshFormat
 $Nodes
 6                      @emph{six mesh nodes:}
@@ -3745,7 +3745,7 @@ $EndElements
 where
 @table @code
 @item @var{version-number}
-is a real number equal to 2.1.
+is a real number equal to 2.2.
 
 @item @var{file-type}
 is an integer equal to 1.
@@ -3797,7 +3797,7 @@ 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);
+fwrite(&header, sizeof(int), 2, file);
 @end example
 
 @item @var{elements-binary}
@@ -3809,13 +3809,12 @@ contain the tags, and the last (#@var{node-number-list} * 4) contain the
 node indices.
 
 Here is a pseudo C code to write @var{elements-binary} for triangles
-with the 3 standard tags (the physical and elementary regions, and the
-mesh partition):
+with the 2 standard tags (the physical and elementary regions):
 @example
 for(i = 0; i < number_of_triangles; i++)@{
-  int data[7] = @{num_i, physical, elementary, partition, 
+  int data[6] = @{num_i, physical, elementary, 
                  node_i_1, node_i_2, node_i_3@};
-  fwrite(data, sizeof(int), 7, file);
+  fwrite(data, sizeof(int), 6, file);
 @}
 @end example
 @end table