diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 86c554ab0849139e3bc1b27c97a23f43c2346abd..2618cce1c60ce231ef63a0f48e665fe45ddfda30 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -598,6 +598,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
 			     std::map<int, int> borderPhysTags[2],
 			     int gLsTag)
 {
+
   int elementary = ge->tag();
   int eType = e->getType();
   std::vector<int> gePhysicals = ge->physicals;
@@ -615,11 +616,9 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
   bool splitElem = false;
   double eps = 1.e-9;
   int lsTag =  (verticesLs(0, e->getVertex(0)->getIndex()) > eps) ? -1 : 1;
-  //printf("*** elementsplit %g(%d) \n", verticesLs(0, e->getVertex(0)->getIndex()) , e->getVertex(0)->getIndex());
-  for(int k = 1; k < e->getNumVertices(); k++){
+   for(int k = 1; k < e->getNumVertices(); k++){
     int lsTag2 = (verticesLs(0, e->getVertex(k)->getIndex()) > eps) ? -1: 1;
-    //printf("elementsplit %g(%d) \n", verticesLs(0, e->getVertex(k)->getIndex()) , e->getVertex(k)->getIndex());
-    if (lsTag*lsTag2 < 0.0) {
+     if (lsTag*lsTag2 < 0.0) {
       lsTag = -1;
       splitElem = true;
       break;
@@ -699,7 +698,6 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
       MVertex *v1N = vertexMap[v1->getNum()];
       double val0 = verticesLs(0, v0->getIndex()) - eps;
       double val1 = verticesLs(0, v1->getIndex()) - eps;
-      //printf("splitedge v0= (%g,%g) v1= (%g,%g) val0= %g(%d) val1= %g(%d) nums (%d %d)\n", v0->x(), v0->y(), v1->x(), v1->y(), val0, v0->getIndex(),  val1, v1->getIndex(), vertexMap[v0->getNum()]->getIndex(), vertexMap[v1->getNum()]->getIndex());
       if(val0*val1 > 0.0 && val0 < 0.0) {
         getPhysicalTag(-1, gePhysicals, newPhysTags[1]);
         int c = elements[1].count(gLsTag);
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 657d563b95d41daeb924fe32de6ffdc8ab9789e8..16f2e0192a3c26dde818c117878495b412350680 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -880,10 +880,9 @@ void gLevelsetMathEvalAll::hessian(double x, double y, double z,
     dfdzz = res[12];
   }
 }
-
-#if defined(HAVE_ANN)
-gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose)
-  :  _gm(gm), _nbClose(nbClose)
+#if defined(HAVE_ANN) 
+gLevelsetDistMesh::gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose, int tag)
+  : _gm(gm), _nbClose(nbClose), gLevelsetPrimitive(tag)
 {
   std::map<int, std::vector<GEntity*> > groups [4];
   gm->getPhysicalGroups(groups);
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 604fb0f8de75719cdd23eb72f3b7cf8879ba61f5..4f11c0e061e947498eae358e89eb2a9140c6eca4 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -353,7 +353,7 @@ class gLevelsetDistMesh: public gLevelsetPrimitive
   ANNidxArray _index;
   ANNdistArray _dist;
 public :
-  gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose = 5);
+  gLevelsetDistMesh(GModel *gm, std::string physical, int nbClose = 5, int tag=1);
   double operator () (double x, double y, double z) const;
   ~gLevelsetDistMesh();
   int type() const { return UNKNOWN; }
diff --git a/benchmarks/centerlines/aneurysm_centerlines.geo b/benchmarks/centerlines/aneurysm_centerlines.geo
index ba4987ebd668e066d7c5b549ffee3186ca71343e..b7986715fb26c54ab9e66b7776bf9c305b5cf3ed 100644
--- a/benchmarks/centerlines/aneurysm_centerlines.geo
+++ b/benchmarks/centerlines/aneurysm_centerlines.geo
@@ -17,7 +17,7 @@ Field[1].hLayer = 0.2;//percent of vessel radius
 
 Field[1].closeVolume =1;
 //Field[1].extrudeWall =1;
-Field[1].reMesh =1;
+//Field[1].reMesh =1;
 
 Field[1].run;