diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 0db6bd169210e9005e99a3c63a2df9351e75e9cc..4ed98e91e6e4d97f60a2de3fabfd815ca77a6c0e 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -29,6 +29,7 @@ set(SRC
   ExtrudeParams.cpp
   Geo.cpp
   GeoStringInterface.cpp GeoInterpolation.cpp
+  gmshLevelset.cpp
   findLinks.cpp
   SOrientedBoundingBox.cpp
   GeomMeshMatcher.cpp
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 5d05bc505be428eae6c0aa71691ca7b80b3e4fd1..6005e3357e71a7f17c230f3d09882eea3148b478 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1800,7 +1800,6 @@ GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
     cutGM->_storeElementsInEntities(elements[i]);
 
   cutGM->_associateEntityWithMeshVertices();
-
   cutGM->_storeVerticesInEntities(vertexMap);
 
   for(int i = 0; i < 4; i++)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index f54b09c65b972c5618db6f81523e3b91c40cbff5..3f2a3ca74463c25615973062421cccdc2f973873 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -8,9 +8,9 @@
 #include "GModel.h"
 #include "MElement.h"
 #include "MElementCut.h"
+#include "gmshLevelset.h"
 
 #if defined(HAVE_DINTEGRATION)
-#include "DILevelset.h"
 #include "Integration3D.h"
 #endif
 
diff --git a/contrib/DiscreteIntegration/DILevelset.cpp b/Geo/gmshLevelset.cpp
similarity index 94%
rename from contrib/DiscreteIntegration/DILevelset.cpp
rename to Geo/gmshLevelset.cpp
index 2ce725106798f1c5bfe2a8387534f78274a2877d..c784cf167be50e7af291692cbdcb9a0e14d3aa47 100644
--- a/contrib/DiscreteIntegration/DILevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -1,11 +1,8 @@
-#ifndef LEVELSET_CC
-#define LEVELSET_CC
-
-#include "DILevelset.h"
+#include "gmshLevelset.h"
 #include <queue>
 #include <stack>
 #include "fullMatrix.h"
-
+#include "mathEvaluator.h"
 
 inline double det3(double d11, double d12, double d13, double d21, double d22, double d23, double d31, double d32, double d33) {
   return d11 * (d22 * d33 - d23 * d32) - d21 * (d12 * d33 - d13 * d32) + d31 * (d12 * d23 - d13 * d22);
@@ -305,7 +302,7 @@ gLevelsetPoints::gLevelsetPoints(const gLevelsetPoints &lv) : gLevelsetPrimitive
   points = lv.points;
 }
 
-double gLevelsetPoints::operator()(const double &x, const double &y, const double &z) const{
+double gLevelsetPoints::operator()(const double x, const double y, const double z) const{
 
   if(mapP.empty()) printf("Levelset Points : call computeLS() before calling operator()\n");
 
@@ -465,11 +462,49 @@ void gLevelsetQuadric::init(){
   C = 0.;
 }
 
-double gLevelsetQuadric::operator()(const double &x, const double &y, const double &z) const{
+double gLevelsetQuadric::operator()(const double x, const double y, const double z) const{
   return(A[0][0] * x * x + 2. * A[0][1] * x * y + 2. * A[0][2] * x * z + A[1][1] * y * y 
         + 2. * A[1][2] * y * z + A[2][2] * z * z + B[0] * x + B[1] * y + B[2] * z + C);
 }
 
+gLevelsetMathEval::gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitive(tag) {
+    std::vector<std::string> expressions(1, f);
+    std::vector<std::string> variables(3);
+    variables[0] = "x";
+    variables[1] = "y";
+    variables[2] = "z";
+    _expr = new mathEvaluator(expressions, variables);
+}
+
+double gLevelsetMathEval::operator() (const double x, const double y, const double z) const {
+    std::vector<double> values(3), res(1);
+    values[0] = x;
+    values[1] = y;
+    values[2] = z;
+    if(_expr->eval(values, res)) return res[0];
+    return 1.;
+}
+
+#if defined (HAVE_POST)
+gLevelsetPostView::gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(tag), _viewIndex(index){
+  if(_viewIndex >= 0 && _viewIndex < (int)PView::list.size()){
+    PView *view = PView::list[_viewIndex];
+    _octree = new OctreePost(view);
+  }
+  else{
+    Msg::Error("Unknown View[%d] in PostView levelset", _viewIndex);
+    _octree = 0;
+  }
+}
+#endif
+
+double gLevelsetPostView::operator () (const double x, const double y, const double z) const  {
+  if(!_octree) return 1.;
+  double val = 1.;
+  _octree->searchScalar(x, y, z, &val, 0);
+  return val;
+}
+
 gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir, const double &R,
                                            int tag) : gLevelsetQuadric(tag) {
   A[0][0] = 1.;
@@ -634,4 +669,3 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const dou
 }
 
 gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv) : gLevelsetImproved(lv){}
-#endif
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index 75d39688e89efd83ea799f989c2cd5451c44adb8..9beef8585c70ec2d96d734c9a3b40636c0982878 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -7,74 +7,494 @@
 #define _GMSH_LEVELSET_H_
 
 #include <string>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h> // for abs()
 #include <vector>
+#include "fullMatrix.h"
+#include "MVertex.h"
 #include "GmshConfig.h"
+#include "mathEvaluator.h"
+
+#if defined(HAVE_POST)
+#include "PView.h"
+#include "OctreePost.h"
+#endif
 
-#if defined(HAVE_DINTEGRATION)
+// PRIMITIVE LEVELSET
+#define UNKNOWN      0
+#define SPHERE       1
+#define PLANE        2
+#define GENCYLINDER  3
+#define ELLIPS       4
+#define CONE         5
+#define QUADRIC      6
+#define BOX          7
+#define CYLINDER     8
+#define CONROD       9
+#define LSMESH      10
+#define LSPOINTS    11 // don't define 'POINTS' as it's reserved by win32
 
-#include "DILevelset.h"
-#include "mathEvaluator.h"
+// TOOLS
+#define CUT       12
+#define UNION     13
+#define INTER     14
+#define CRACK     15
+
+
+class gLevelset
+{
+protected:
+  static const short insideDomain = -1; // negative values of the levelset are inside the domain.
+  int tag_; // must be greater than 0
+public:
+  gLevelset() : tag_(-1) {}
+  gLevelset(const gLevelset &);
+  virtual ~gLevelset(){}
+  virtual gLevelset * clone() const{printf("Error virtual fct called gLevelset::clone()"); return 0;}
+  virtual double operator() (const double x, const double y, const double z) const = 0;
+  bool isInsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain > 0.;}
+  bool isOutsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain < 0.;}
+  bool isOnBorder      (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) == 0.;}
+  virtual std::vector<const gLevelset *> getChildren() const = 0;
+  virtual double choose (double d1, double d2) const = 0;
+  virtual int type() const = 0;
+  virtual bool isPrimitive() const = 0;
+  void setTag(int t) { tag_ = t; }
+  virtual int getTag() const { return tag_; }
+  void getPrimitives(std::vector<const gLevelset *> &primitives) const;
+  void getPrimitivesPO(std::vector<const gLevelset *> &primitives) const;
+  void getRPN(std::vector<const gLevelset *> &gLsRPN) const;
+  double H (const double &x, const double &y, const double &z) const {
+    if (isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
+    return 0.0;
+  }
+  void print() const {
+    printf("LS : ");
+    switch(type()) {
+    case SPHERE :      printf("SPHERE"); break;
+    case PLANE  :      printf("PLANE"); break;
+    case GENCYLINDER : printf("GENCYLINDER"); break;
+    case ELLIPS :      printf("ELLIPS"); break;
+    case CONE :        printf("CONE"); break;
+    case QUADRIC :     printf("QUADRIC"); break;
+    case BOX :         printf("BOX"); break;
+    case CYLINDER :    printf("CYLINDER"); break;
+    case CONROD :      printf("CONROD"); break;
+    case CUT :         printf("CUT"); break;
+    case UNION :       printf("UNION"); break;
+    case INTER :       printf("INTER"); break;
+    case LSMESH:       printf("LSMESH"); break;
+    case LSPOINTS:     printf("LSPOINTS"); break;
+    }
+    printf(" Tag=%d\n", getTag());
+  }
+};
+
+// ----------------------------------------------------------------------------------------------------------
+// PRIMITIVES
+class gLevelsetPrimitive : public gLevelset
+{
+public:
+  gLevelsetPrimitive() : gLevelset() {}
+  gLevelsetPrimitive(const gLevelsetPrimitive &lv) : gLevelset(lv) {}
+  gLevelsetPrimitive(int tag) {
+    if (tag < 1) {
+      printf("Tag of the levelset (%d) must be greater than 0.\n", tag);
+      tag = abs(tag);
+    }
+    tag_ = tag;
+  }
+  virtual double operator () (const double x, const double y, const double z) const = 0;
+  std::vector<const gLevelset *> getChildren() const { std::vector<const gLevelset *> p; return p; }
+  double choose (double d1, double d2) const {
+    printf("Cannot use function \"choose\" with a primitive!\n");
+    return d1;
+  }
+  virtual int type() const = 0;
+  bool isPrimitive() const {return true;}
+  
+};
+
+class gLevelsetSphere : public gLevelsetPrimitive
+{
+protected:
+  double xc, yc, zc, r;
+public:
+  gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
+  virtual double operator () (const double x, const double y, const double z) const
+    {return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;}
+  int type() const {return SPHERE;}
+};
+
+class gLevelsetPlane : public gLevelsetPrimitive
+{
+protected:
+  double a, b, c, d;
+public:
+  // define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
+  gLevelsetPlane (const double _a, const double _b, const double _c, const double _d, int tag) : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
+  // define the plane passing through the point pt and with outward normal norm
+  gLevelsetPlane (const double *pt, const double *norm, int tag);
+  // define the plane passing through the 3 points pt1,pt2,pt3 and with outward normal (pt1,pt2)x(pt1,pt3)
+  gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3, int tag);
+  // copy constructor
+  gLevelsetPlane(const gLevelsetPlane &lv);
+  virtual gLevelset * clone() const{return new gLevelsetPlane(*this);}
+  // return negative value inward and positive value outward
+  virtual double operator() (const double x, const double y, const double z) const
+    {return a * x + b * y + c * z + d;} 
+  int type() const {return PLANE;}
+};
+
+class gLevelsetPoints : public gLevelsetPrimitive
+{
+protected:
+  fullMatrix<double> points;
+  fullMatrix<double> surf;
+  fullMatrix<double> matAInv;
+  double delta;
+  std::map<SPoint3,double> mapP;
+  fullMatrix<double> generateRbfMat(int p, int index,
+				    const fullMatrix<double> &nodes1,
+				    const fullMatrix<double> &nodes2) const;
+  void RbfOp(int p, int index, 
+	     const fullMatrix<double> &cntrs,
+	     const fullMatrix<double> &nodes, 
+	     fullMatrix<double> &D, 
+	     bool isLocal = false) const;
+  void evalRbfDer(int p, int index,
+		  const fullMatrix<double> &cntrs,
+		  const fullMatrix<double> &nodes,
+		  const fullMatrix<double> &fValues, 
+		fullMatrix<double> &fApprox, bool isLocal = false) const;
+  void setup_level_set(const fullMatrix<double> &cntrs,
+		       fullMatrix<double> &level_set_nodes, 
+		       fullMatrix<double> &level_set_funvals);
+
+public:
+  // define the data points
+  gLevelsetPoints(fullMatrix<double> &_centers, int tag); 
+  // copy constructor
+  gLevelsetPoints(const gLevelsetPoints &lv);
+  virtual gLevelset * clone() const{return new gLevelsetPoints(*this);}
+  // return negative value inward and positive value outward
+  virtual double operator() (const double x, const double y, const double z) const;
+  void computeLS(std::vector<MVertex*> &vert);
+  int type() const {return LSPOINTS;}
+};
+
+class gLevelsetQuadric : public gLevelsetPrimitive
+{
+protected:
+  double A[3][3], B[3], C;
+  void translate (const double transl[3]);
+  void rotate (const double rotate[3][3]);  
+  void computeRotationMatrix (const double dir[3], double t[3][3]); 
+  void computeRotationMatrix (const double dir1[3], const double dir2[3] , double t[3][3]); 
+  void Ax (const double x[3], double res[3], double fact=1.0);
+  void xAx (const double x[3], double &res, double fact=1.0);
+  void init ();
+
+public:
+  gLevelsetQuadric() : gLevelsetPrimitive() {init(); }
+  gLevelsetQuadric(int tag) : gLevelsetPrimitive(tag) {init(); }
+  gLevelsetQuadric(const gLevelsetQuadric &);
+  virtual ~gLevelsetQuadric() {}
+  double operator () (const double x, const double y, const double z) const;
+  virtual int type() const = 0;
+};
+
+class gLevelsetGenCylinder : public gLevelsetQuadric
+{
+public:
+  gLevelsetGenCylinder (const double *pt, const double *dir, const double &R, int tag);
+  gLevelsetGenCylinder (const gLevelsetGenCylinder& );
+  virtual gLevelset * clone() const{return new gLevelsetGenCylinder(*this);}
+  int type() const {return GENCYLINDER;}
+};
+
+class gLevelsetEllipsoid : public gLevelsetQuadric
+{
+public:
+  gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int tag);
+  gLevelsetEllipsoid (const gLevelsetEllipsoid& );
+  virtual gLevelset * clone() const{return new gLevelsetEllipsoid(*this);}
+  int type() const {return ELLIPS;}
+};
+
+class gLevelsetCone : public gLevelsetQuadric
+{
+public:
+  gLevelsetCone (const double *pt, const double *dir, const double &angle, int tag);
+  gLevelsetCone (const gLevelsetCone& );
+  virtual gLevelset * clone() const{return new gLevelsetCone(*this);}
+  int type() const {return CONE;}
+};
+
+class gLevelsetGeneralQuadric : public gLevelsetQuadric
+{
+public:
+  gLevelsetGeneralQuadric (const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
+                           const double &z, const double &c, int tag);
+  gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& );
+  virtual gLevelset * clone() const{return new gLevelsetGeneralQuadric(*this);}
+  int type() const {return QUADRIC;}
+};
 
 class gLevelsetMathEval: public gLevelsetPrimitive
 {
   mathEvaluator *_expr;
 public:
-  gLevelsetMathEval(std::string f, int tag) : gLevelsetPrimitive(tag)
-  {
-    std::vector<std::string> expressions(1, f);
-    std::vector<std::string> variables(3);
-    variables[0] = "x";
-    variables[1] = "y";
-    variables[2] = "z";
-    _expr = new mathEvaluator(expressions, variables);
-  }
-  ~gLevelsetMathEval(){ delete _expr; }
-  virtual double operator () (const double &x, const double &y, const double &z) const
-  {
-    std::vector<double> values(3), res(1);
-    values[0] = x;
-    values[1] = y;
-    values[2] = z;
-    if(_expr->eval(values, res)) return res[0];
-    return 1.;
-  }
+  gLevelsetMathEval(std::string f, int tag);
+  ~gLevelsetMathEval(){ if (_expr) delete _expr; }
+  double operator () (const double x, const double y, const double z) const;
   int type() const { return UNKNOWN; }
 };
 
-#if defined(HAVE_POST)
-
-#include "PView.h"
-#include "OctreePost.h"
 
+#if defined(HAVE_POST)
 class gLevelsetPostView : public gLevelsetPrimitive
 {
   int _viewIndex;
   OctreePost *_octree;
 public:
-  gLevelsetPostView(int index, int tag) : gLevelsetPrimitive(tag), _viewIndex(index)
-  {
-    if(_viewIndex >= 0 && _viewIndex < (int)PView::list.size()){
-      PView *view = PView::list[_viewIndex];
-      _octree = new OctreePost(view);
-    }
-    else{
-      Msg::Error("Unknown View[%d] in PostView levelset", _viewIndex);
-      _octree = 0;
+  gLevelsetPostView(int index, int tag) ;
+  ~gLevelsetPostView(){ if(_octree) delete _octree;}
+  double operator () (const double x, const double y, const double z) const;
+  int type() const { return UNKNOWN; }
+};
+#endif
+
+
+// --------------------------------------------------------------------------------------------------------------
+// TOOLS
+class gLevelsetTools : public gLevelset
+{
+protected:
+  std::vector<const gLevelset *> children;
+public:
+  gLevelsetTools () {}
+  gLevelsetTools (std::vector<const gLevelset *> &p) {children = p;}
+  gLevelsetTools (const gLevelsetTools &);
+  ~gLevelsetTools () {
+    for(int i = 0; i < (int)children.size(); i++)
+      delete children[i];
+  }
+  double operator () (const double x, const double y, const double z) const {
+    double d = (*children[0])(x, y, z);
+    for (int i = 1; i < (int)children.size(); i++){
+      double dt = (*children[i])(x, y, z);
+      d = choose(d, dt);
     }
+    return d;
   }
-  ~gLevelsetPostView(){ if(_octree) delete _octree; }
-  virtual double operator () (const double &x, const double &y, const double &z) const
-  {
-    if(!_octree) return 1.;
-    double val = 1.;
-    _octree->searchScalar(x, y, z, &val, 0);
-    return val;
+  std::vector<const gLevelset *> getChildren() const {
+    if(children.size() != 1) return children;
+    return children[0]->getChildren();
+  }
+  virtual double choose (double d1, double d2) const = 0;
+  virtual int type2() const = 0;
+  virtual int type() const {
+    if(children.size() != 1) return type2();
+    return children[0]->type();
+  }
+  bool isPrimitive() const {
+    if(children.size() != 1) return false;
+    return children[0]->isPrimitive();
+  }
+  int getTag() const {
+    if(children.size() != 1) return tag_;
+    return children[0]->getTag();
   }
-  int type() const { return UNKNOWN; }
 };
 
-#endif
+class gLevelsetReverse : public gLevelset
+{
+protected:
+  const gLevelset *ls;
+public:
+  gLevelsetReverse (const gLevelset *p) : ls(p){}
+  double operator () (const double x, const double y, const double z) const {
+    return -(*ls)(x, y, z);
+  }
+  std::vector<const gLevelset *> getChildren() const {return ls->getChildren();}
+  bool isPrimitive() const {return ls->isPrimitive();}
+  virtual double choose (double d1, double d2) const {return -ls->choose(d1,d2);}
+  virtual int type() const {return ls->type();}
+  int getTag() const { return ls->getTag(); }
+};
+
+
+// This levelset takes the first levelset in the list as the object and the others as tools that cut it
+class gLevelsetCut : public gLevelsetTools
+{
+public:
+  gLevelsetCut (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+  double choose (double d1, double d2) const {
+    return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
+  }
+  gLevelsetCut(const gLevelsetCut &lv):gLevelsetTools(lv){}
+  virtual gLevelset * clone() const{return new gLevelsetCut(*this);}
+  int type2() const {return CUT;}
+};
+
+// This levelset takes the minimum
+class gLevelsetUnion : public gLevelsetTools
+{
+public:
+  gLevelsetUnion (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+  gLevelsetUnion(const gLevelsetUnion &lv):gLevelsetTools(lv){}
+  virtual gLevelset * clone() const{return new gLevelsetUnion(*this);}
+  
+  double choose (double d1, double d2) const {
+    return (d1 < d2) ? d1 : d2; // lesser of d1 and d2
+  }
+  int type2() const {return UNION;}
+};
+
+// This levelset takes the maximum
+class gLevelsetIntersection : public gLevelsetTools
+{
+public:
+  gLevelsetIntersection (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+  gLevelsetIntersection(const gLevelsetIntersection &lv):gLevelsetTools(lv) { }
+  virtual gLevelset *clone() const { return new gLevelsetIntersection(*this); }
+
+  double choose (double d1, double d2) const {
+    return (d1 > d2) ? d1 : d2; // greater of d1 and d2
+  }
+  int type2() const {return INTER;}
+};
+
+// Crack defined by a normal and a tangent levelset
+class gLevelsetCrack : public gLevelsetTools
+{
+public:
+  gLevelsetCrack (std::vector<const gLevelset *> &p) {
+    if (p.size() != 2)
+      printf("Error : gLevelsetCrack needs 2 levelsets\n");
+    children.push_back(p[0]);
+    children.push_back(new gLevelsetReverse(p[0]));
+    if(p[1]) children.push_back(p[1]);
+  }
+  double choose (double d1, double d2) const {
+    return (d1 > d2) ? d1 : d2; // greater of d1 and d2
+  }
+  int type2() const {return CRACK;}
+};
+
+// --------------------------------------------------------------------------------------------------------------
+// IMPROVED LEVELSET
+class gLevelsetImproved : public gLevelset
+{
+protected:
+  gLevelset *Ls;
+public:
+  gLevelsetImproved(){}
+  gLevelsetImproved(const gLevelsetImproved &lv);
+  double operator() (const double x, const double y, const double z) const {return (*Ls)(x, y, z);}
+  std::vector<const gLevelset *> getChildren() const { return Ls->getChildren(); }
+  double choose (double d1, double d2) const { return Ls->choose(d1, d2); }
+  virtual int type() const = 0;
+  bool isPrimitive() const {return Ls->isPrimitive();}
+};
+
+class gLevelsetBox : public gLevelsetImproved
+{
+public:
+  // create a box with parallel faces : pt is a corner of the box,
+  //                                    dir1 is the direction of the first edge starting from pt,
+  //                                    dir2 is the direction of the second edge starting from pt,
+  //                                    dir3 is the direction of the third edge starting from pt,
+  //                                    a is the length of the first edge starting from pt,
+  //                                    b is the length of the second edge starting from pt,
+  //                                    c is the length of the third edge starting from pt.
+  // tags of the faces are : face normal to dir3 and not including pt : tag+0
+  //                         face normal to dir3 and     including pt : tag+1
+  //                         face normal to dir2 and     including pt : tag+2
+  //                         face normal to dir2 and not including pt : tag+3
+  //                         face normal to dir1 and not including pt : tag+4
+  //                         face normal to dir1 and     including pt : tag+5
+  gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3,
+               const double &a, const double &b, const double &c, int tag);
+  // create a box with the 8 vertices (pt1,...,pt8). 
+  // check if the faces are planar.
+  // tags of the faces are : face(pt5,pt6,pt7,pt8) : tag+0
+  //                         face(pt1,pt4,pt3,pt2) : tag+1
+  //                         face(pt1,pt2,pt6,pt5) : tag+2
+  //                         face(pt3,pt4,pt8,pt7) : tag+3
+  //                         face(pt2,pt3,pt7,pt6) : tag+4
+  //                         face(pt1,pt5,pt8,pt4) : tag+5
+  gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4,
+               const double *pt5, const double *pt6, const double *pt7, const double *pt8, int tag);
+  gLevelsetBox(const gLevelsetBox &);
+  virtual gLevelset * clone() const{return new gLevelsetBox(*this);}
+  int type() const {return BOX;}
+};
+
+class gLevelsetCylinder : public gLevelsetImproved
+{
+public:
+  // create a cylinder : pt is the point in the middle of the cylinder base,
+  //                     dir is the direction of the cylinder axis,
+  //                     R is the outer radius of the cylinder,
+  //                     H is the height of the cylinder.
+  // tags of the faces are : exterior face :             tag+0
+  //                         plane face including pt :   tag+1
+  //                         plane face opposite to pt : tag+2
+  gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &H, int tag);
+  // create a cylinder : pt is the point in the middle of the cylinder base,
+  //                     dir is the direction of the cylinder axis,
+  //                     R is the outer radius of the cylinder,
+  //                     r is the inner radius of the cylinder,
+  //                     H is the height of the cylinder.
+  // tags of the faces are : exterior face :             tag+0
+  //                         plane face including pt :   tag+1
+  //                         plane face opposite to pt : tag+2
+  //                         interior face :             tag+3
+  gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &r, const double &H, int tag);
+  gLevelsetCylinder(const gLevelsetCylinder &);
+  virtual gLevelset * clone() const{return new gLevelsetCylinder(*this);}
+  int type() const {return CYLINDER;}
+};
+
+class gLevelsetConrod : public gLevelsetImproved
+{
+public:
+  // create a connecting rod : pt is the point in the middle of the first bore,
+  //                           dir1 is the direction of the rod,
+  //                           dir2 is the direction of the axis of the bore,
+  //                           H1 is height of the first cylinder,
+  //                           H2 is the height of the second cylinder,
+  //                           H3 is the height of the rod,
+  //                           R1 is the outer radius of the first cylinder,
+  //                           r1 is the inner radius of the first cylinder,
+  //                           R2 is the outer radius of the second cylinder,
+  //                           r2 is the inner radius of the second cylinder,
+  //                           L1 is the width of the rod in the plane passing through the middle of the first bore,
+  //                           L2 is the width of the rod in the plane passing through the middle of the second bore,
+  //                           E is the distance between the axis of the cylinders.
+  // tags of the faces are : bottom face (+dir2) of the bore :      tag+2
+  //                         top    face (-dir2) of the bore :      tag+3
+  //                         rear   face (-dir1xdir2) of the bore : tag+4
+  //                         front  face (+dir1xdir2) of the bore : tag+5
+  //                         exterior face of the first cylinder :  tag+6
+  //                         bottom   face of the first cylinder :  tag+7
+  //                         top      face of the first cylinder :  tag+8
+  //                         exterior face of the second cylinder : tag+9
+  //                         bottom   face of the second cylinder : tag+10
+  //                         top      face of the second cylinder : tag+11
+  //                         interior face of the first  cylinder : tag+12
+  //                         interior face of the second cylinder : tag+13
+  gLevelsetConrod (const double *pt, const double *dir1, const double *dir2,
+                   const double &H1, const double &H2, const double &H3,
+                   const double &R1, const double &r1, const double &R2, const double &r2,
+                   const double &L1, const double &L2, const double &E, int tag);
+  gLevelsetConrod(const gLevelsetConrod &);
+  virtual gLevelset * clone() const{return new gLevelsetConrod(*this);}
+  int type() const {return CONROD;}
+};
 
-#endif
 
 #endif
diff --git a/Solver/FuncGradDisc.h b/Solver/FuncGradDisc.h
index 7021490625d2a166ce1c197c4d778198c30434f8..2b5b8c7e090402ab36aea485b60a7aa0e531ed05 100644
--- a/Solver/FuncGradDisc.h
+++ b/Solver/FuncGradDisc.h
@@ -12,7 +12,7 @@
 #define _FUNCGRADDISC_H_
 
 #include "simpleFunction.h"
-#include "DILevelset.h"
+#include "gmshLevelset.h"
 #include "MVertex.h"
 #include "GModel.h"
 //#include "gLevelSetMesh.cpp"
diff --git a/Solver/FuncHeaviside.h b/Solver/FuncHeaviside.h
index 6a24f67f2db2a18604c2b3a2e00b4629e354ab83..06e321e9c5a1cc8e706dde2bb73bbd248b8e9368 100644
--- a/Solver/FuncHeaviside.h
+++ b/Solver/FuncHeaviside.h
@@ -12,8 +12,7 @@
 #define _FUNCHEAVISIDE_H_
 
 #include "simpleFunction.h"
-#include "DILevelset.h"
-
+#include "gmshLevelset.h"
 
 class FuncHeaviside : public  simpleFunctionOnElement<double> {
  private :
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 5e610e4982a0e6ad38eef03ec1f8d9934063def7..4eda845f761184246efd502e6278206e6f94cb36 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -18,7 +18,7 @@
 #include "quadratureRules.h"
 #include "solverField.h"
 #include "MPoint.h"
-#include "DILevelset.h"
+#include "gmshLevelset.h"
 #if defined(HAVE_POST)
 #include "PView.h"
 #include "PViewData.h"
diff --git a/Solver/filters.h b/Solver/filters.h
index 5e05b1745dddacbde373b66a0f2b2de6381f50a5..038aaeae7ac410510af4fd5740f5e86bd46dfb42 100644
--- a/Solver/filters.h
+++ b/Solver/filters.h
@@ -15,7 +15,7 @@
 #include "dofManager.h"
 #include "GModel.h"
 #include "groupOfElements.h"
-#include "DILevelset.h"
+#include "gmshLevelset.h"
 
 class FilterNodeEnriched
 {
diff --git a/benchmarks/levelset/cube.geo b/benchmarks/levelset/cube.geo
index e93e0d127eb3565eb6cf6770e50bd592fd075d24..b69a57bbde2b4f9802c37e3123e087fffaf1e3d7 100644
--- a/benchmarks/levelset/cube.geo
+++ b/benchmarks/levelset/cube.geo
@@ -27,6 +27,6 @@ Levelset Cut (12) = {3,7};
 
 Levelset Point(13)={{0.5,0,0},{0.5,0.5,0.},{0.5,1,0},{0.5,0,0.5},{0.5,0,0.8},{0.5,0.5,0.5},{0.5,0.5,1},{0.5,1,0.5},{1.5,1,1}};
 
-Levelset CutMeshTri {13};
+Levelset CutMeshTri {5};
 
 
diff --git a/contrib/DiscreteIntegration/CMakeLists.txt b/contrib/DiscreteIntegration/CMakeLists.txt
index d1af4a09f53d90f4044f2dc58319c9aa6f63ab4c..b9cae1f5421fae44960e78c6e0ef012d1e0ecab4 100644
--- a/contrib/DiscreteIntegration/CMakeLists.txt
+++ b/contrib/DiscreteIntegration/CMakeLists.txt
@@ -4,7 +4,6 @@
 # bugs and problems to <gmsh@geuz.org>.
 
 set(SRC
-  DILevelset.cpp
   Integration3D.cpp
   recurCut.cpp
 )
diff --git a/contrib/DiscreteIntegration/DILevelset.h b/contrib/DiscreteIntegration/DILevelset.h
deleted file mode 100644
index b84d85760f7f5c6861d433feaecd70deba2be8ff..0000000000000000000000000000000000000000
--- a/contrib/DiscreteIntegration/DILevelset.h
+++ /dev/null
@@ -1,462 +0,0 @@
-#ifndef _gLEVELSET_H_
-#define _gLEVELSET_H_
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h> // for abs()
-#include <vector>
-#include "fullMatrix.h"
-#include "MVertex.h"
-
-// PRIMITIVE LEVELSET
-#define UNKNOWN      0
-#define SPHERE       1
-#define PLANE        2
-#define GENCYLINDER  3
-#define ELLIPS       4
-#define CONE         5
-#define QUADRIC      6
-#define BOX          7
-#define CYLINDER     8
-#define CONROD       9
-#define LSMESH      10
-#define LSPOINTS    11 // don't define 'POINTS' as it's reserved by win32
-
-// TOOLS
-#define CUT       12
-#define UNION     13
-#define INTER     14
-#define CRACK     15
-
-class gLevelset
-{
-protected:
-  static const short insideDomain = -1; // negative values of the levelset are inside the domain.
-  int tag_; // must be greater than 0
-public:
-  gLevelset() : tag_(-1) {}
-  gLevelset(const gLevelset &);
-  virtual ~gLevelset(){}
-  virtual gLevelset * clone() const{printf("Error virtual fct called gLevelset::clone()"); return 0;}
-  virtual double operator() (const double &x, const double &y, const double &z) const = 0;
-  // inline double operator () (const SPoint3 &p) const {return this->operator()(p.x(),p.y(),p.z());}
-  bool isInsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain > 0.;}
-  bool isOutsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain < 0.;}
-  bool isOnBorder      (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) == 0.;}
-  virtual std::vector<const gLevelset *> getChildren() const = 0;
-  virtual double choose (double d1, double d2) const = 0;
-  virtual int type() const = 0;
-  virtual bool isPrimitive() const = 0;
-  void setTag(int t) { tag_ = t; }
-  virtual int getTag() const { return tag_; }
-  void getPrimitives(std::vector<const gLevelset *> &primitives) const;
-  void getPrimitivesPO(std::vector<const gLevelset *> &primitives) const;
-  void getRPN(std::vector<const gLevelset *> &gLsRPN) const;
-  double H (const double &x, const double &y, const double &z) const {
-    if (isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
-    return 0.0;
-  }
-  void print() const {
-    printf("LS : ");
-    switch(type()) {
-    case SPHERE :      printf("SPHERE"); break;
-    case PLANE  :      printf("PLANE"); break;
-    case GENCYLINDER : printf("GENCYLINDER"); break;
-    case ELLIPS :      printf("ELLIPS"); break;
-    case CONE :        printf("CONE"); break;
-    case QUADRIC :     printf("QUADRIC"); break;
-    case BOX :         printf("BOX"); break;
-    case CYLINDER :    printf("CYLINDER"); break;
-    case CONROD :      printf("CONROD"); break;
-    case CUT :         printf("CUT"); break;
-    case UNION :       printf("UNION"); break;
-    case INTER :       printf("INTER"); break;
-    case LSMESH:       printf("LSMESH"); break;
-    case LSPOINTS:       printf("LSPOINTS"); break;
-    }
-    printf(" Tag=%d\n", getTag());
-  }
-};
-
-// ----------------------------------------------------------------------------------------------------------
-// PRIMITIVES
-class gLevelsetPrimitive : public gLevelset
-{
-public:
-  gLevelsetPrimitive() : gLevelset() {}
-  gLevelsetPrimitive(const gLevelsetPrimitive &lv) : gLevelset(lv) {}
-  gLevelsetPrimitive(int tag) {
-    if (tag < 1) {
-      printf("Tag of the levelset (%d) must be greater than 0.\n", tag);
-      tag = abs(tag);
-    }
-    tag_ = tag;
-  }
-  virtual double operator () (const double &x, const double &y, const double &z) const = 0;
-  std::vector<const gLevelset *> getChildren() const { std::vector<const gLevelset *> p; return p; }
-  double choose (double d1, double d2) const {
-    printf("Cannot use function \"choose\" with a primitive!\n");
-    return d1;
-  }
-  virtual int type() const = 0;
-  bool isPrimitive() const {return true;}
-  
-};
-
-class gLevelsetSphere : public gLevelsetPrimitive
-{
-protected:
-  double xc, yc, zc, r;
-public:
-  gLevelsetSphere (const double &x, const double &y, const double &z, const double &R, int tag) : gLevelsetPrimitive(tag), xc(x), yc(y), zc(z), r(R) {}
-  virtual double operator () (const double &x, const double &y, const double &z) const
-    {return sqrt((xc - x) * (xc - x) + (yc - y) * (yc - y) + (zc - z) * (zc - z)) - r;}
-  int type() const {return SPHERE;}
-};
-
-class gLevelsetPlane : public gLevelsetPrimitive
-{
-protected:
-  double a, b, c, d;
-public:
-  // define the plane _a*x+_b*y+_c*z+_d, with outward normal (a,b,c)
-  gLevelsetPlane (const double _a, const double _b, const double _c, const double _d, int tag) : gLevelsetPrimitive(tag), a(_a), b(_b), c(_c), d(_d) {}
-  // define the plane passing through the point pt and with outward normal norm
-  gLevelsetPlane (const double *pt, const double *norm, int tag);
-  // define the plane passing through the 3 points pt1,pt2,pt3 and with outward normal (pt1,pt2)x(pt1,pt3)
-  gLevelsetPlane (const double *pt1, const double *pt2, const double *pt3, int tag);
-  // copy constructor
-  gLevelsetPlane(const gLevelsetPlane &lv);
-  virtual gLevelset * clone() const{return new gLevelsetPlane(*this);}
-  // return negative value inward and positive value outward
-  virtual double operator() (const double &x, const double &y, const double &z) const
-    {return a * x + b * y + c * z + d;} 
-  int type() const {return PLANE;}
-};
-
-class gLevelsetPoints : public gLevelsetPrimitive
-{
-protected:
-  fullMatrix<double> points;
-  fullMatrix<double> surf;
-  fullMatrix<double> matAInv;
-  double delta;
-  std::map<SPoint3,double> mapP;
-  fullMatrix<double> generateRbfMat(int p, int index,
-				    const fullMatrix<double> &nodes1,
-				    const fullMatrix<double> &nodes2) const;
-  void RbfOp(int p, int index, 
-	     const fullMatrix<double> &cntrs,
-	     const fullMatrix<double> &nodes, 
-	     fullMatrix<double> &D, 
-	     bool isLocal = false) const;
-  void evalRbfDer(int p, int index,
-		  const fullMatrix<double> &cntrs,
-		  const fullMatrix<double> &nodes,
-		  const fullMatrix<double> &fValues, 
-		fullMatrix<double> &fApprox, bool isLocal = false) const;
-  void setup_level_set(const fullMatrix<double> &cntrs,
-		       fullMatrix<double> &level_set_nodes, 
-		       fullMatrix<double> &level_set_funvals);
-
-public:
-  // define the data points
-  gLevelsetPoints(fullMatrix<double> &_centers, int tag); 
-  // copy constructor
-  gLevelsetPoints(const gLevelsetPoints &lv);
-  virtual gLevelset * clone() const{return new gLevelsetPoints(*this);}
-  // return negative value inward and positive value outward
-  virtual double operator() (const double &x, const double &y, const double &z) const;
-  void computeLS(std::vector<MVertex*> &vert);
-  int type() const {return LSPOINTS;}
-};
-
-class gLevelsetQuadric : public gLevelsetPrimitive
-{
-protected:
-  double A[3][3], B[3], C;
-  void translate (const double transl[3]);
-  void rotate (const double rotate[3][3]);  
-  void computeRotationMatrix (const double dir[3], double t[3][3]); 
-  void computeRotationMatrix (const double dir1[3], const double dir2[3] , double t[3][3]); 
-  void Ax (const double x[3], double res[3], double fact=1.0);
-  void xAx (const double x[3], double &res, double fact=1.0);
-  void init ();
-
-public:
-  gLevelsetQuadric() : gLevelsetPrimitive() {init(); }
-  gLevelsetQuadric(int tag) : gLevelsetPrimitive(tag) {init(); }
-  gLevelsetQuadric(const gLevelsetQuadric &);
-  virtual ~gLevelsetQuadric() {}
-  double operator () (const double &x, const double &y, const double &z) const;
-  virtual int type() const = 0;
-};
-
-class gLevelsetGenCylinder : public gLevelsetQuadric
-{
-public:
-  gLevelsetGenCylinder (const double *pt, const double *dir, const double &R, int tag);
-  gLevelsetGenCylinder (const gLevelsetGenCylinder& );
-  virtual gLevelset * clone() const{return new gLevelsetGenCylinder(*this);}
-  int type() const {return GENCYLINDER;}
-};
-
-class gLevelsetEllipsoid : public gLevelsetQuadric
-{
-public:
-  gLevelsetEllipsoid (const double *pt, const double *dir, const double &a, const double &b, const double &c, int tag);
-  gLevelsetEllipsoid (const gLevelsetEllipsoid& );
-  virtual gLevelset * clone() const{return new gLevelsetEllipsoid(*this);}
-  int type() const {return ELLIPS;}
-};
-
-class gLevelsetCone : public gLevelsetQuadric
-{
-public:
-  gLevelsetCone (const double *pt, const double *dir, const double &angle, int tag);
-  gLevelsetCone (const gLevelsetCone& );
-  virtual gLevelset * clone() const{return new gLevelsetCone(*this);}
-  int type() const {return CONE;}
-};
-
-class gLevelsetGeneralQuadric : public gLevelsetQuadric
-{
-public:
-  gLevelsetGeneralQuadric (const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
-                           const double &z, const double &c, int tag);
-  gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& );
-  virtual gLevelset * clone() const{return new gLevelsetGeneralQuadric(*this);}
-  int type() const {return QUADRIC;}
-};
-
-// --------------------------------------------------------------------------------------------------------------
-// TOOLS
-class gLevelsetTools : public gLevelset
-{
-protected:
-  std::vector<const gLevelset *> children;
-public:
-  gLevelsetTools () {}
-  gLevelsetTools (std::vector<const gLevelset *> &p) {children = p;}
-  gLevelsetTools (const gLevelsetTools &);
-  ~gLevelsetTools () {
-    for(int i = 0; i < (int)children.size(); i++)
-      delete children[i];
-  }
-  double operator () (const double &x, const double &y, const double &z) const {
-    double d = (*children[0])(x, y, z);
-    for (int i = 1; i < (int)children.size(); i++){
-      double dt = (*children[i])(x, y, z);
-      d = choose(d, dt);
-    }
-    return d;
-  }
-  std::vector<const gLevelset *> getChildren() const {
-    if(children.size() != 1) return children;
-    return children[0]->getChildren();
-  }
-  virtual double choose (double d1, double d2) const = 0;
-  virtual int type2() const = 0;
-  virtual int type() const {
-    if(children.size() != 1) return type2();
-    return children[0]->type();
-  }
-  bool isPrimitive() const {
-    if(children.size() != 1) return false;
-    return children[0]->isPrimitive();
-  }
-  int getTag() const {
-    if(children.size() != 1) return tag_;
-    return children[0]->getTag();
-  }
-};
-
-class gLevelsetReverse : public gLevelset
-{
-protected:
-  const gLevelset *ls;
-public:
-  gLevelsetReverse (const gLevelset *p) : ls(p){}
-  double operator () (const double &x, const double &y, const double &z) const {
-    return -(*ls)(x, y, z);
-  }
-  std::vector<const gLevelset *> getChildren() const {return ls->getChildren();}
-  bool isPrimitive() const {return ls->isPrimitive();}
-  virtual double choose (double d1, double d2) const {return -ls->choose(d1,d2);}
-  virtual int type() const {return ls->type();}
-  int getTag() const { return ls->getTag(); }
-};
-
-
-// This levelset takes the first levelset in the list as the object and the others as tools that cut it
-class gLevelsetCut : public gLevelsetTools
-{
-public:
-  gLevelsetCut (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
-  double choose (double d1, double d2) const {
-    return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
-  }
-  gLevelsetCut(const gLevelsetCut &lv):gLevelsetTools(lv){}
-  virtual gLevelset * clone() const{return new gLevelsetCut(*this);}
-  int type2() const {return CUT;}
-};
-
-// This levelset takes the minimum
-class gLevelsetUnion : public gLevelsetTools
-{
-public:
-  gLevelsetUnion (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
-  gLevelsetUnion(const gLevelsetUnion &lv):gLevelsetTools(lv){}
-  virtual gLevelset * clone() const{return new gLevelsetUnion(*this);}
-  
-  double choose (double d1, double d2) const {
-    return (d1 < d2) ? d1 : d2; // lesser of d1 and d2
-  }
-  int type2() const {return UNION;}
-};
-
-// This levelset takes the maximum
-class gLevelsetIntersection : public gLevelsetTools
-{
-public:
-  gLevelsetIntersection (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
-  gLevelsetIntersection(const gLevelsetIntersection &lv):gLevelsetTools(lv) { }
-  virtual gLevelset *clone() const { return new gLevelsetIntersection(*this); }
-
-  double choose (double d1, double d2) const {
-    return (d1 > d2) ? d1 : d2; // greater of d1 and d2
-  }
-  int type2() const {return INTER;}
-};
-
-// Crack defined by a normal and a tangent levelset
-class gLevelsetCrack : public gLevelsetTools
-{
-public:
-  gLevelsetCrack (std::vector<const gLevelset *> &p) {
-    if (p.size() != 2)
-      printf("Error : gLevelsetCrack needs 2 levelsets\n");
-    children.push_back(p[0]);
-    children.push_back(new gLevelsetReverse(p[0]));
-    if(p[1]) children.push_back(p[1]);
-  }
-  double choose (double d1, double d2) const {
-    return (d1 > d2) ? d1 : d2; // greater of d1 and d2
-  }
-  int type2() const {return CRACK;}
-};
-
-// --------------------------------------------------------------------------------------------------------------
-// IMPROVED LEVELSET
-class gLevelsetImproved : public gLevelset
-{
-protected:
-  gLevelset *Ls;
-public:
-  gLevelsetImproved(){}
-  gLevelsetImproved(const gLevelsetImproved &lv);
-  double operator() (const double &x, const double &y, const double &z) const {return (*Ls)(x, y, z);}
-  std::vector<const gLevelset *> getChildren() const { return Ls->getChildren(); }
-  double choose (double d1, double d2) const { return Ls->choose(d1, d2); }
-  virtual int type() const = 0;
-  bool isPrimitive() const {return Ls->isPrimitive();}
-};
-
-class gLevelsetBox : public gLevelsetImproved
-{
-public:
-  // create a box with parallel faces : pt is a corner of the box,
-  //                                    dir1 is the direction of the first edge starting from pt,
-  //                                    dir2 is the direction of the second edge starting from pt,
-  //                                    dir3 is the direction of the third edge starting from pt,
-  //                                    a is the length of the first edge starting from pt,
-  //                                    b is the length of the second edge starting from pt,
-  //                                    c is the length of the third edge starting from pt.
-  // tags of the faces are : face normal to dir3 and not including pt : tag+0
-  //                         face normal to dir3 and     including pt : tag+1
-  //                         face normal to dir2 and     including pt : tag+2
-  //                         face normal to dir2 and not including pt : tag+3
-  //                         face normal to dir1 and not including pt : tag+4
-  //                         face normal to dir1 and     including pt : tag+5
-  gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3,
-               const double &a, const double &b, const double &c, int tag);
-  // create a box with the 8 vertices (pt1,...,pt8). 
-  // check if the faces are planar.
-  // tags of the faces are : face(pt5,pt6,pt7,pt8) : tag+0
-  //                         face(pt1,pt4,pt3,pt2) : tag+1
-  //                         face(pt1,pt2,pt6,pt5) : tag+2
-  //                         face(pt3,pt4,pt8,pt7) : tag+3
-  //                         face(pt2,pt3,pt7,pt6) : tag+4
-  //                         face(pt1,pt5,pt8,pt4) : tag+5
-  gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4,
-               const double *pt5, const double *pt6, const double *pt7, const double *pt8, int tag);
-  gLevelsetBox(const gLevelsetBox &);
-  virtual gLevelset * clone() const{return new gLevelsetBox(*this);}
-  int type() const {return BOX;}
-};
-
-class gLevelsetCylinder : public gLevelsetImproved
-{
-public:
-  // create a cylinder : pt is the point in the middle of the cylinder base,
-  //                     dir is the direction of the cylinder axis,
-  //                     R is the outer radius of the cylinder,
-  //                     H is the height of the cylinder.
-  // tags of the faces are : exterior face :             tag+0
-  //                         plane face including pt :   tag+1
-  //                         plane face opposite to pt : tag+2
-  gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &H, int tag);
-  // create a cylinder : pt is the point in the middle of the cylinder base,
-  //                     dir is the direction of the cylinder axis,
-  //                     R is the outer radius of the cylinder,
-  //                     r is the inner radius of the cylinder,
-  //                     H is the height of the cylinder.
-  // tags of the faces are : exterior face :             tag+0
-  //                         plane face including pt :   tag+1
-  //                         plane face opposite to pt : tag+2
-  //                         interior face :             tag+3
-  gLevelsetCylinder (const double *pt, const double *dir, const double &R, const double &r, const double &H, int tag);
-  gLevelsetCylinder(const gLevelsetCylinder &);
-  virtual gLevelset * clone() const{return new gLevelsetCylinder(*this);}
-  int type() const {return CYLINDER;}
-};
-
-class gLevelsetConrod : public gLevelsetImproved
-{
-public:
-  // create a connecting rod : pt is the point in the middle of the first bore,
-  //                           dir1 is the direction of the rod,
-  //                           dir2 is the direction of the axis of the bore,
-  //                           H1 is height of the first cylinder,
-  //                           H2 is the height of the second cylinder,
-  //                           H3 is the height of the rod,
-  //                           R1 is the outer radius of the first cylinder,
-  //                           r1 is the inner radius of the first cylinder,
-  //                           R2 is the outer radius of the second cylinder,
-  //                           r2 is the inner radius of the second cylinder,
-  //                           L1 is the width of the rod in the plane passing through the middle of the first bore,
-  //                           L2 is the width of the rod in the plane passing through the middle of the second bore,
-  //                           E is the distance between the axis of the cylinders.
-  // tags of the faces are : bottom face (+dir2) of the bore :      tag+2
-  //                         top    face (-dir2) of the bore :      tag+3
-  //                         rear   face (-dir1xdir2) of the bore : tag+4
-  //                         front  face (+dir1xdir2) of the bore : tag+5
-  //                         exterior face of the first cylinder :  tag+6
-  //                         bottom   face of the first cylinder :  tag+7
-  //                         top      face of the first cylinder :  tag+8
-  //                         exterior face of the second cylinder : tag+9
-  //                         bottom   face of the second cylinder : tag+10
-  //                         top      face of the second cylinder : tag+11
-  //                         interior face of the first  cylinder : tag+12
-  //                         interior face of the second cylinder : tag+13
-  gLevelsetConrod (const double *pt, const double *dir1, const double *dir2,
-                   const double &H1, const double &H2, const double &H3,
-                   const double &R1, const double &r1, const double &R2, const double &r2,
-                   const double &L1, const double &L2, const double &E, int tag);
-  gLevelsetConrod(const gLevelsetConrod &);
-  virtual gLevelset * clone() const{return new gLevelsetConrod(*this);}
-  int type() const {return CONROD;}
-};
-
-
-#endif
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index f8a8afa14d862afcf4cb9964405e5f445d94fc50..c2c338fc077050a45928ae62cdd7e19cb87373d7 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -6,7 +6,7 @@
 #include <set>
 #include <map>
 #include <cmath>
-#include "DILevelset.h"
+#include "gmshLevelset.h"
 
 // Element type
 #define DI_LIN  1
diff --git a/contrib/DiscreteIntegration/recurCut.h b/contrib/DiscreteIntegration/recurCut.h
index 132d333695bc09b42b079e43a9052c4bb8dd4207..aca3b61c4562cdbef859161db71315d208136a5b 100644
--- a/contrib/DiscreteIntegration/recurCut.h
+++ b/contrib/DiscreteIntegration/recurCut.h
@@ -1,7 +1,7 @@
 #ifndef _RECUR_H_
 #define _RECUR_H_
 
-#include "DILevelset.h"
+#include "gmshLevelset.h"
 #include "Integration3D.h"
 
 class RecurElement
diff --git a/gmshpy/gmshGeo.i b/gmshpy/gmshGeo.i
index a5eb3008a11bc50efe20793132158c699775441a..8d64807c13d22d3bb23769e55f579c2e284e5746 100644
--- a/gmshpy/gmshGeo.i
+++ b/gmshpy/gmshGeo.i
@@ -20,10 +20,7 @@
   #include "discreteFace.h"
   #include "discreteEdge.h"
   #include "discreteVertex.h"
-#ifdef HAVE_DINTEGRATION
   #include "gmshLevelset.h"
-  #include "DILevelset.h"
-#endif
   #include "MElement.h"
   #include "MVertex.h"
   #include "MTriangle.h"
@@ -81,8 +78,5 @@ namespace std {
 %include "SPoint2.h"
 %include "SBoundingBox3d.h"
 %include "Curvature.h"
-#ifdef HAVE_DINTEGRATION
 %include "gmshLevelset.h"
-%include "DILevelset.h"
-#endif