diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 3f16f02dc3e90208ca5ebaebb9ec865fbc73ee4e..da8627e9e8bab05b7d466e5c5d8c9cb409782df9 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1093,7 +1093,7 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
   if(getNumChildren() == 0) {
     for(int i = 0; i < getNumVertices(); i++) {
       MVertex *v = getVertex(i);
-      int numV = v->getNum();
+      int numV = v->getIndex();
       if(vertexMap.count(numV))
         vmv.push_back(vertexMap[numV]);
       else {
@@ -1107,7 +1107,7 @@ MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
     for(int i = 0; i < getNumChildren(); i++) {
       for(int j = 0; j < getChild(i)->getNumVertices(); j++) {
         MVertex *v = getChild(i)->getVertex(j);
-        int numV = v->getNum();
+        int numV = v->getIndex();
         if(vertexMap.count(numV))
           vmv.push_back(vertexMap[numV]);
         else {
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 6ef38440f58528c03a8041fd706cfa6cacb1fe7a..423e0157619f9baa6966cbc84d5258bc00007b89 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -612,9 +612,12 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
   //EMI : better for embedded dirichlet with smoothed properties
   //split according to values of vertices (keep +)
   bool splitElem = false;
-  int lsTag =  (verticesLs(0, e->getVertex(0)->getIndex()) > 0.0) ? -1 : 1;
+  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++){
-    int lsTag2 = (verticesLs(0, e->getVertex(k)->getIndex()) > 0.0) ? -1: 1;
+    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) {
       lsTag = -1;
       splitElem = true;
@@ -687,13 +690,17 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
 
   //create level set interface (pt in 1D, line in 2D or face in 3D)
   if (splitElem && copy->getDim()==2){
+    //printf("*** SPLITELEM \n");
     for(int k = 0; k < copy->getNumEdges(); k++){
       MEdge me  = copy->getEdge(k);
       MVertex *v0 = me.getVertex(0);
       MVertex *v1 = me.getVertex(1);
-      double val0 = verticesLs(0, v0->getIndex());
-      double val1 = verticesLs(0, v1->getIndex());
-      if (val0*val1 > 0.0 && val0 < 0.0) { 
+      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, val1, v0->getIndex(), v1->getIndex(), v0->getNum(), v1->getNum());
+      if ( (val0*val1 > 0.0 && val0 < 0.0 ) ) { 
+	//printf("split edge \n");
         MLine *lin = new MLine(v0, v1); 
         getPhysicalTag(-1, gePhysicals, newPhysTags[1]);
         int c = elements[1].count(gLsTag);
@@ -701,7 +708,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
         int physTag = (!gePhysicals.size()) ? 0 : getBorderTag(gLsTag, c, newPhysTags[1][0], borderPhysTags[0]);
         elements[1][reg].push_back(lin);
         if(physTag) assignLsPhysical(GM, reg, 1, physicals, physTag, gLsTag);
-      }   
+      }
     }
   }
   else if (splitElem && copy->getDim()==3){
@@ -1366,10 +1373,13 @@ 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++) {
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index ccd5058f1d90bedb08e7deecb6171bcacc6193eb..ddf98a0527d2c688aa70ca1c2b6a4ffc227c4fca 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -11,6 +11,7 @@
 #include <queue>
 #include <stack>
 #include "fullMatrix.h"
+#include "GModel.h"
 
 
 inline double det3(double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) {
@@ -494,6 +495,26 @@ double gLevelsetMathEval::operator() (const double x, const double y, const doub
     return 1.;
 }
 
+//EMI TO DO 
+gLevelsetDistGeom::gLevelsetDistGeom(std::string name, int tag) : gLevelsetPrimitive(tag) {
+  _model = new GModel();
+  _model->load(name);
+   for (GModel::fiter fit = _model->firstFace(); fit != _model->lastFace(); fit++){
+     if((*it)->geomType() == GEntity::DiscreteSurface){
+
+     }
+     else{
+       for (int i = 0; i < (*fit)->stl_triangles.size(); i += 3){
+       }
+     }
+   }
+}
+
+double gLevelsetDistGeom::operator() (const double x, const double y, const double z) const {
+  double dist = 1.0;
+  return dist;
+}
+
 
 #if defined (HAVE_POST)
 gLevelsetPostView::gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(tag), _viewIndex(index){
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 7cdd4b85706ebcfc0b4034f9cc654f43157d488e..de5974559e7fe7415677e67cb7c93bd0edce91d7 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -254,6 +254,17 @@ public:
   int type() const { return UNKNOWN; }
 };
 
+class GModel;
+class gLevelsetDistGeom: public gLevelsetPrimitive
+{
+  GModel *_model;
+public:
+  gLevelsetDistGeom(std::string f, int tag=1);
+  ~gLevelsetDistGeom(){ if (_model) delete _model; }
+  double operator () (const double x, const double y, const double z) const;
+  int type() const { return UNKNOWN; }
+};
+
 #if defined(HAVE_POST)
 class gLevelsetPostView : public gLevelsetPrimitive
 {