diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index ec300ed0ff7646e0a3317430d643a75958ba91d3..676cee9c348a402c0153487e79b8e5584167d219 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -950,30 +950,31 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
           else
             poly[1].push_back(mt);
         }
-	// //if quads
-        // for (unsigned int i = nbQ; i < quads.size(); i++){
-        //   MVertex *mv[4] = {NULL, NULL, NULL, NULL};
-        //   for(int j = 0; j < 4; j++){
-        //   int numV = getElementVertexNum(triangles[i]->pt(j), e);
-        //     if(numV == -1) {
-        //       MVertex *newv = new MVertex(quads[i]->x(j), quads[i]->y(j), quads[i]->z(j));
-        //       std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
-        //       mv[j] = *(it.first);
-        //       if (!it.second) newv->deleteLast();
-        //     }
-        //     else {
-	//       std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
-	//       if(it == vertexMap.end()) {
-        //         mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
-        //                             quads[i]->z(j), 0, numV);
-        //         vertexMap[numV] = mv[j];
-        //       }
-        //       else mv[j] = it->second;
-        //     }
-	//   }
-        //   MQuadrangle *mq = new MQuadrangle(mv[0], mv[1], mv[2], mv[3], 0, 0);
-	//   elements[3][reg].push_back(mq);
-        // }
+	//if quads
+        for (unsigned int i = nbQ; i < quads.size(); i++){
+          MVertex *mv[4] = {NULL, NULL, NULL, NULL};
+          for(int j = 0; j < 4; j++){
+          int numV = getElementVertexNum(quads[i]->pt(j), e);
+            if(numV == -1) {
+              MVertex *newv = new MVertex(quads[i]->x(j), quads[i]->y(j), quads[i]->z(j));
+              std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
+              mv[j] = *(it.first);
+              if (!it.second) newv->deleteLast();
+            }
+            else {
+	      std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
+	      if(it == vertexMap.end()) {
+                mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
+                                    quads[i]->z(j), 0, numV);
+                vertexMap[numV] = mv[j];
+              }
+              else mv[j] = it->second;
+            }
+	  }
+          MQuadrangle *mq = new MQuadrangle(mv[0], mv[1], mv[2], mv[3], 0, 0);
+	  int reg = getElementaryTag(quads[i]->lsTag(), elementary, newElemTags[2]);
+	  elements[3][reg].push_back(mq);
+        }
 
         bool own = (eParent && !e->ownsParent()) ? false : true;
         if(poly[0].size()) {
@@ -1337,7 +1338,9 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
     }
 
     for(DI_Point::Container::iterator it = cp.begin(); it != cp.end(); it++) delete *it;
+    printf("coucou delet dim =%d\n", gmEntities[i]->dim());
     for(unsigned int k = 0; k < lines.size(); k++) delete lines[k];
+    printf("coucou delet \n");
     for(unsigned int k = 0; k < triangles.size(); k++) delete triangles[k];
     for(unsigned int k = 0; k < quads.size(); k++) delete quads[k];
     for(unsigned int k = 0; k < tetras.size(); k++) delete tetras[k];
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index b8d68b5b6443243af97877982ec447193ea4a895..76fb9977c67e1675b7cdb8523d9e631ab77ac995 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -5,7 +5,7 @@
 #include "recurCut.h"
 #include "Gauss.h"
 
-#define ZERO_LS_TOL  1.e-7
+#define ZERO_LS_TOL  1.e-3
 #define EQUALITY_TOL 1.e-15
 
 // cross product
@@ -729,7 +729,7 @@ bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
   double sideLength = distance(pt(s1), pt(s2));
   if (dist / sideLength < 1.e-5) {
     if(order0 == 1) setPolynomialOrder(1);
-    printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quad edge\n", dist, sideLength, dist/sideLength);
+    printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quadratic edge\n", dist, sideLength, dist/sideLength);
     return true; // do not add the quadratic edge if xm is very close from the middle of the edge
   }
   DI_Point *p0 = &mid_[edge];
@@ -1447,7 +1447,8 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<gLevelset *>
                              std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
                              std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
 {
-  bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // if true, use quadratic sub triangles
+  //bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // if true, use quadratic sub triangles
+  bool quadT = false;
   int LStag = RPNi.back()->getTag();
 
   int nbZe = 0, ze[3];
@@ -1565,7 +1566,8 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<gLevelset *> &R
                        std::vector<DI_Tetra *> &subTetras, std::vector<DI_Triangle *> &surfTriangles,
                        std::vector<DI_CuttingPoint *> &cuttingPoints, std::vector<DI_QualError *> &QError) const
 {
-  bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // use of quadratic surf triangles and quadratic sub tetrahedra
+  //bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // use of quadratic surf triangles and quadratic sub tetrahedra
+  bool quadT = false;
   int tag = RPNi.back()->getTag();
 
   int nbZe = 0, ze[4];