diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index a495a7027741e7c98e6533b750857c0cbe09f81d..d04e042aa1ec19cc66fef19765c81189ceb68a7a 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -53,10 +53,19 @@ void GEdge::reverse()
 }
 
 unsigned int GEdge::getNumMeshElements()
-{ 
+{
   return lines.size();
 }
 
+unsigned int GEdge::getNumMeshParentElements()
+{
+  unsigned int n = 0;
+  for(unsigned int i = 0; i < lines.size(); i++)
+    if(lines[i]->ownsParent())
+      n++;
+  return n;
+}
+
 void GEdge::getNumMeshElements(unsigned *const c) const
 {
   c[0] += lines.size();
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index fd3f9a845e5984c5f1b1badc77e0525d6d57e491..7f98689c584e649ed9ad34de32b60f6e28012bce 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -141,6 +141,7 @@ class GEdge : public GEntity {
 
   // get total/by-type number of elements in the mesh
   unsigned int getNumMeshElements();
+  unsigned int getNumMeshParentElements();
   void getNumMeshElements(unsigned *const c) const;
 
   // get the start of the array of a type of element
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index 084d7b7b2fa28b9a25c97df74e49fb7ef5931ea8..ff0f3ee98159111e5f94ada8aeaa7dab8b8abe58 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -287,6 +287,7 @@ class GEntity {
 
   // get the number of mesh elements (total and by type) in the entity
   virtual unsigned int getNumMeshElements() { return 0; }
+  virtual unsigned int getNumMeshParentElements() { return 0; }
   virtual void getNumMeshElements(unsigned *const c) const { };
 
   // get the start of the array of a type of element
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 97645751eaa5c594a8772b6874fb59fafc66dee7..7321425206845b7cfbad23d64e1b4c94b7b212aa 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -93,6 +93,15 @@ unsigned int GFace::getNumMeshElements()
   return triangles.size() + quadrangles.size() + polygons.size();
 }
 
+unsigned int GFace::getNumMeshParentElements()
+{
+  unsigned int n = 0;
+  for(unsigned int i = 0; i < polygons.size(); i++)
+    if(polygons[i]->ownsParent())
+      n++;
+  return n;
+}
+
 void GFace::getNumMeshElements(unsigned *const c) const
 {
   c[0] += triangles.size();
diff --git a/Geo/GFace.h b/Geo/GFace.h
index bea0547df01240a268537d33d230fba17acca2b1..509eb6573ca1f76cc6d22753cd6588abfe70ff82 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -234,6 +234,7 @@ class GFace : public GEntity
 
   // get total/by-type number of elements in the mesh
   unsigned int getNumMeshElements();
+  unsigned int getNumMeshParentElements();
   void getNumMeshElements(unsigned *const c) const;
 
   // get the start of the array of a type of element
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index f19e13970cb22c4da5a830e277ccdfe35c352e8c..5d05bc505be428eae6c0aa71691ca7b80b3e4fd1 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -566,6 +566,16 @@ int GModel::getNumMeshElements()
   return n;
 }
 
+int GModel::getNumMeshParentElements()
+{
+  std::vector<GEntity*> entities;
+  getEntities(entities);
+  unsigned int n = 0;
+  for(unsigned int i = 0; i < entities.size(); i++)
+    n += entities[i]->getNumMeshParentElements();
+  return n;
+}
+
 int GModel::getNumMeshElements(unsigned c[5])
 {
   c[0] = 0; c[1] = 0; c[2] = 0; c[3] = 0; c[4] = 0;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 3a7936955086fde8db0a43e92e677eeb8a5f05a5..44a2b2ffaceceaf018e96961332c89df6a6d5603 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -288,6 +288,7 @@ class GModel
 
   // return the total number of elements in the mesh
   int getNumMeshElements();
+  int getNumMeshParentElements();
 
   // get the number of each type of element in the mesh at the largest
   // dimension and return the dimension
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index 57c7d0ba5098dafdf1803d1837e1ff7c890299ae..62849ad249526c583713ebfa0835fb6688acf426 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -56,6 +56,15 @@ unsigned int GRegion::getNumMeshElements()
     polyhedra.size();
 }
 
+unsigned int GRegion::getNumMeshParentElements()
+{
+  unsigned int n = 0;
+  for(unsigned int i = 0; i < polyhedra.size(); i++)
+    if(polyhedra[i]->ownsParent())
+      n++;
+  return n;
+}
+
 void GRegion::getNumMeshElements(unsigned *const c) const
 {
   c[0] += tetrahedra.size();
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index f869c133a9a6de69b6317fd011875a06b3b7d145..77acd62130f7ccb2988fb03368826259b34bad02 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -78,6 +78,7 @@ class GRegion : public GEntity {
 
   // get total/by-type number of elements in the mesh
   unsigned int getNumMeshElements();
+  unsigned int getNumMeshParentElements();
   void getNumMeshElements(unsigned *const c) const;
 
   // get the start of the array of a type of element
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 9e29a31d9ee1daec09dc6dc5775499fd737cd0a4..0d1da05eb7c5f533ca05adecb9b50ba41e65a90c 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1080,7 +1080,7 @@ int *MElement::getVerticesIdForMSH()
   return verts;
 }
 
-MElement *MElement::copy(int &num, std::map<int, MVertex*> &vertexMap,
+MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
                          std::map<MElement*, MElement*> &newParents,
                          std::map<MElement*, MElement*> &newDomains)
 {
@@ -1123,7 +1123,7 @@ MElement *MElement::copy(int &num, std::map<int, MVertex*> &vertexMap,
     std::map<MElement*, MElement*>::iterator it = newParents.find(eParent);
     MElement *newParent;
     if(it == newParents.end()) {
-      newParent = eParent->copy(++num, vertexMap, newParents, newDomains);
+      newParent = eParent->copy(vertexMap, newParents, newDomains);
       newParents[eParent] = newParent;
     }
     else
@@ -1140,7 +1140,7 @@ MElement *MElement::copy(int &num, std::map<int, MVertex*> &vertexMap,
     std::map<MElement*, MElement*>::iterator it = newDomains.find(dom);
     MElement *newDom;
     if(it == newDomains.end()) {
-      newDom = dom->copy(++num, vertexMap, newParents, newDomains);
+      newDom = dom->copy(vertexMap, newParents, newDomains);
       newDomains[dom] = newDom;
     }
     else
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 943b5ca5eb3fa85c72411f56ed25d926847e033d..2a901bb434be4d1d61a87e9ff383fbeee747e744 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -339,7 +339,7 @@ class MElement
   virtual int *getVerticesIdForMSH();
 
   // copy element and parent if any, vertexMap contains the new vertices
-  virtual MElement *copy(int &num, std::map<int, MVertex*> &vertexMap,
+  virtual MElement *copy(std::map<int, MVertex*> &vertexMap,
                          std::map<MElement*, MElement*> &newParents,
                          std::map<MElement*, MElement*> &newDomains);
 };
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 39f17442e718423d53815ec2f5f81ac5e017f2d5..f54b09c65b972c5618db6f81523e3b91c40cbff5 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -516,7 +516,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
   int eType = e->getTypeForMSH();
   std::vector<int> gePhysicals = ge->physicals;
 
-  MElement *copy = e->copy(numEle, vertexMap, newParents, newDomains);
+  MElement *copy = e->copy(vertexMap, newParents, newDomains);
 
   double lsMean = 0.;
   for(int k = 0; k < e->getNumVertices(); k++)
@@ -642,7 +642,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
   unsigned int nbTr = triangles.size();
   unsigned int nbTe = tetras.size();
 
-  MElement *copy = e->copy(numEle, vertexMap, newParents, newDomains);
+  MElement *copy = e->copy(vertexMap, newParents, newDomains);
   MElement *parent = eParent ? copy->getParent() : copy;
 
   double **nodeLs = new double*[e->getNumPrimaryVertices()];
@@ -1131,7 +1131,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
   int nbLs = (cutElem) ? ((primS > 1) ? primS + 1 : 1) : 1;
   fullMatrix<double> verticesLs(nbLs, numVert + 1);
 
- //Emi test compute all at once for POINTS (type = 11)
+ //compute all at once for ls POINTS (type = 11)
   std::vector<MVertex *> vert;
   for(unsigned int i = 0; i < gmEntities.size(); i++) {
     for(unsigned int j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
@@ -1139,7 +1139,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
     }
   }
   for(int k = 0; k < primS; k++){
-    if (primitives[k]->type() == 11){ //points
+    if (primitives[k]->type() == 11){
       ((gLevelsetPoints*)primitives[k])->computeLS(vert);
     }
   }
@@ -1167,7 +1167,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
     newPhysTags[d][0] = gm->getMaxPhysicalNumber(d); //max value at [dim][0]
   }
 
-  int numEle = gm->getNumMeshElements(); //element number increment
+  int numEle = gm->getNumMeshElements() + gm->getNumMeshParentElements(); //element number increment
   std::map<MElement*, MElement*> newParents; //map<oldParent, newParent>
   std::map<MElement*, MElement*> newDomains; //map<oldDomain, newDomain>