diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 0fa196375caa16551b859c2876d3c236b8748ff5..a188c8b6a3bc0f9d847f7c2d863ef243cf235486 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -617,6 +617,9 @@ static void writeElementMSH(FILE *fp, GModel *model, T *ele, bool saveAll,
   }
 
   model->setMeshElementIndex(ele, num); // should really be a multimap...
+
+  if(CTX::instance()->mesh.saveTri && ele->getNumChildren())
+    num += ele->getNumChildren() - 1;
 }
 
 template<class T>
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6f45c85c18548d9c02880ed744e595d8bba2f3ca..e5f2fed785e7590947436bb6106523e3c40471d9 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -630,13 +630,13 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
   if(CTX::instance()->mesh.saveTri && poly){
     for (int i = 0; i < getNumChildren() ; i++){
        MElement *t = getChild(i);
-       t->writeMSH(fp, version, binary, num, elementary, physical, 0, 0, 0, ghosts);
+       t->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
     }
     return;
   }
   else if(CTX::instance()->mesh.saveTri && polyl){
     MLine *l = new MLine(getVertex(0), getVertex(1));
-    l->writeMSH(fp, version, binary, num, elementary, physical, 0, 0, 0, ghosts);
+    l->writeMSH(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
     delete l;
     return;
   }
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index ee41a648accf17f6b2631f67926236b3d3c0d08f..be041390afd2bfed82bb03918dd4453f7b7179fc 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -732,7 +732,10 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
           elements[2][reg].push_back(tri);
         }
         else if(mf.getNumVertices() == 4){
-          MQuadrangle *quad = new MQuadrangle(mf.getVertex(0), mf.getVertex(1), mf.getVertex(2), mf.getVertex(3));
+          MQuadrangle *quad = new MQuadrangle(vertexMap[mf.getVertex(0)->getNum()],
+                                              vertexMap[mf.getVertex(1)->getNum()],
+                                              vertexMap[mf.getVertex(2)->getNum()],
+                                              vertexMap[mf.getVertex(3)->getNum()]);
           elements[3][reg].push_back(quad);
         }
         if(physTag)  assignLsPhysical(GM, reg, 2, physicals, physTag, gLsTag);
@@ -1375,13 +1378,10 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
         if(primS > 1)
           verticesLs(k, vi->getIndex()) = (*ls)(vi->x(), vi->y(), vi->z());
       }
-      else{
+      else
         verticesLs(0, vi->getIndex()) = (*ls)(vi->x(), vi->y(), vi->z());
-	//printf("xy = (%g,%g) val= %g(%g) ind=%d\n",vi->x(), vi->y(),  verticesLs(0, vi->getIndex()),-vi->x()+0.41, vi->getIndex() );
-      }
     }
   }
-  //exit(1);
 
   int numEle = gm->getNumMeshElements() + gm->getNumMeshParentElements(); //element number increment
   for(unsigned int i = 0; i < gmEntities.size(); i++) {