diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 5803463c09d6d63116d87a8dbf1dcdc68f9265a5..763cb786e9de6b756d61402f9949aefe77787580 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -146,7 +146,7 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
       break;
     }
     prevOne = &ges;
-    //    ges.print();
+    //        ges.print();
     loop.push_back(ges);
   }
 }
diff --git a/Geo/MElementOctree.cpp b/Geo/MElementOctree.cpp
index b9e344a76781f0ffc169b9941275483569522442..5c1c5b802edafde077700765ba1c08337025af87 100644
--- a/Geo/MElementOctree.cpp
+++ b/Geo/MElementOctree.cpp
@@ -122,6 +122,21 @@ MElementOctree::~MElementOctree()
   Octree_Delete(_octree);
 }
 
+
+std::vector<MElement *> MElementOctree::findAll(double x, double y, double z, int dim)
+{
+  double P[3] = {x, y, z};
+  std::list<void*> v;
+  std::vector<MElement*> e;
+  Octree_SearchAll(P, _octree,&v);
+  for (std::list<void*>::iterator it = v.begin();
+       it != v.end(); ++it){
+    MElement *el = (MElement*) *it;
+    if (dim == -1 || el->getDim() == dim)e.push_back(el);
+  }
+  return e;
+}
+
 MElement *MElementOctree::find(double x, double y, double z, int dim, bool strict)
 {
   double P[3] = {x, y, z};
diff --git a/Geo/MElementOctree.h b/Geo/MElementOctree.h
index 42e08aa0e2bc07d519189e29aaa4f7da1bd3ab37..29a44ad8ec4d1a6c3a5d47b01e95b497833c1f0d 100644
--- a/Geo/MElementOctree.h
+++ b/Geo/MElementOctree.h
@@ -22,6 +22,7 @@ class MElementOctree{
   ~MElementOctree();
   MElement *find(double x, double y, double z, int dim = -1, bool strict = false);
   Octree *getInternalOctree(){ return _octree; }
+  std::vector<MElement *> findAll(double x, double y, double z, int dim);
 };
 
 #endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 64b3e4d7292e5af6c9223c9d4c30fc001ebeac02..c2bc93cf202c1e5337fa14ed25a3b45fd595c44b 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -445,6 +445,37 @@ static void addOrRemove ( MVertex *v1, MVertex *v2, std::set<MEdge,Less_Edge> &
   else bedges.erase(it);
 }
 
+void filterOverlappingElements (int dim, std::vector<MElement*> &e, std::vector<MElement*> &eout,  std::vector<MElement*> &einter){
+  eout.clear();
+  MElementOctree octree (e);
+  for (int i=0;i<e.size();++i){
+    MElement *el = e[i];
+    bool intersection = false;
+    for (int j=0;j<el->getNumVertices();++j){
+      MVertex *v = el->getVertex(j);
+      std::vector<MElement *> inters = octree.findAll(v->x(),v->y(),v->z(),dim);
+      std::vector<MElement *> inters2;      
+      for (int k=0;k<inters.size();k++){
+	bool found = false;
+	for (int l=0;l<inters[k]->getNumVertices();l++){
+	  if (inters[k]->getVertex(l) == v)found = true;
+	}
+	if (!found)inters2.push_back(inters[k]);
+      }
+      if (inters2.size() >= 1 ){
+	intersection = true;
+      }
+    }
+    if (intersection){
+      printf("intersection found\n");
+      einter.push_back(el);       
+    }
+    else {
+      eout.push_back(el);
+    }
+  }
+}
+
 void modifyInitialMeshForTakingIntoAccountBoundaryLayers  (GFace *gf){
 
   BoundaryLayerColumns *_columns = buidAdditionalPoints2D (gf, M_PI/6.);
@@ -490,18 +521,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers  (GFace *gf){
 
 	  //avoid convergent errors
 	  if (dv2.length() < 0.5 * dv.length())break;
-	  //	printf("quadrangle generated\n");
 	  blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
-	  //blTris.push_back(new MTriangle(v11,v12,v22));
-	  //	  blTris.push_back(new MTriangle(v11,v22,v21));
-	  addOrRemove(v11,v12,bedges);
-	  addOrRemove(v12,v22,bedges);
-	  addOrRemove(v22,v21,bedges);
-	  addOrRemove(v21,v11,bedges);
-	  if(v11->onWhat() == gf)verts.insert(v11);
-	  if(v21->onWhat() == gf)verts.insert(v21);
-	  if(v12->onWhat() == gf)verts.insert(v12);
-	  if(v22->onWhat() == gf)verts.insert(v22);
 	  fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
 		  v11->x(),v11->y(),v11->z(),
 		  v12->x(),v12->y(),v12->z(),
@@ -510,95 +530,76 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers  (GFace *gf){
 	}
 	//	int M = std::max(c1._column.size(),c2._column.size());
 
-	/*
-	if (M>N) M = N+1;
-	// close with triangles
- 	for (int l=N-1;l < M-1 ;++l){
-	  MVertex *v11,*v12,*v21,*v22;
-	  v11 = c1._column[l>=c1._column.size() ? c1._column.size() -1 : l];
-	  v12 = c2._column[l>=c2._column.size() ? c2._column.size() -1 : l];
-	  v21 = c1._column[(l+1)>=c1._column.size() ? c1._column.size() -1 : l+1];
-	  v22 = c2._column[(l+1)>=c2._column.size() ? c2._column.size() -1 : l+1];
-
-	  MTriangle *nt = (v11 == v21) ? new MTriangle(v11,v12,v22) : new MTriangle(v11,v12,v21) ;
-	  blTris.push_back(nt);
-	  v11 = nt->getVertex(0);
-	  v12 = nt->getVertex(1);
-	  v21 = nt->getVertex(2);
-
-	  addOrRemove(v11,v12,bedges);
-	  addOrRemove(v12,v21,bedges);
-	  addOrRemove(v21,v11,bedges);
-	  if(v11->onWhat() == gf)verts.insert(v11);
-	  if(v21->onWhat() == gf)verts.insert(v21);
-	  if(v12->onWhat() == gf)verts.insert(v12);
-	  fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1};\n",
-		  v11->x(),v11->y(),v11->z(),
-		  v12->x(),v12->y(),v12->z(),
-		  v21->x(),v21->y(),v21->z());
-	}
-	*/
       }
    }
     ++ite;
   }
 
-  if (1){
-    for (BoundaryLayerColumns::iterf itf = _columns->beginf();
-	 itf != _columns->endf() ; ++itf){
-      MVertex *v = itf->first;
-      int nbCol = _columns->getNbColumns(v);
-      
-      for (int i=0;i<nbCol-1;i++){
-	//	printf("permut %d %d\n",permut[i],permut[i+1]);
-	const BoundaryLayerData & c1 = _columns->getColumn(v,i);
-	const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
-	int N = std::min(c1._column.size(),c2._column.size());
-	for (int l=0;l < N ;++l){
-	  MVertex *v11,*v12,*v21,*v22;
-	  v21 = c1._column[l];
-	  v22 = c2._column[l];	    
-	  if (l == 0){
-	    v11 = v;
-	    v12 = v;
-	  }
-	  else {
-	    v11 = c1._column[l-1];
-	    v12 = c2._column[l-1];	    
-	  }
-	  //	printf("quadrangle generated\n");
-	  if (v11 != v12){
-	    addOrRemove(v11,v12,bedges);
-	    addOrRemove(v12,v22,bedges);
-	    addOrRemove(v22,v21,bedges);
-	    addOrRemove(v21,v11,bedges);
-	    blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
-	    fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
-		    v11->x(),v11->y(),v11->z(),
+  for (BoundaryLayerColumns::iterf itf = _columns->beginf();
+       itf != _columns->endf() ; ++itf){
+    MVertex *v = itf->first;
+    int nbCol = _columns->getNbColumns(v);
+    
+    for (int i=0;i<nbCol-1;i++){
+      const BoundaryLayerData & c1 = _columns->getColumn(v,i);
+      const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
+      int N = std::min(c1._column.size(),c2._column.size());
+      for (int l=0;l < N ;++l){
+	MVertex *v11,*v12,*v21,*v22;
+	v21 = c1._column[l];
+	v22 = c2._column[l];	    
+	if (l == 0){
+	  v11 = v;
+	  v12 = v;
+	}
+	else {
+	  v11 = c1._column[l-1];
+	  v12 = c2._column[l-1];	    
+	}
+	if (v11 != v12){
+	  blQuads.push_back(new MQuadrangle(v11,v12,v22,v21));
+	  fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+		  v11->x(),v11->y(),v11->z(),
 		    v12->x(),v12->y(),v12->z(),
-		    v22->x(),v22->y(),v22->z(),
-		    v21->x(),v21->y(),v21->z());
-	  }
-	  else {
-	    addOrRemove(v,v22,bedges);
-	    addOrRemove(v22,v21,bedges);
-	    addOrRemove(v21,v,bedges);
-	    blTris.push_back(new MTriangle(v,v22,v21));
-	    fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
-		    v->x(),v->y(),v->z(),
-		    v22->x(),v22->y(),v22->z(),
-		    v21->x(),v21->y(),v21->z());
-	  }
-	  if(v11->onWhat() == gf)verts.insert(v11);
-	  if(v21->onWhat() == gf)verts.insert(v21);
-	  if(v12->onWhat() == gf)verts.insert(v12);
-	  if(v22->onWhat() == gf)verts.insert(v22);
+		  v22->x(),v22->y(),v22->z(),
+		  v21->x(),v21->y(),v21->z());
+	}
+	else {
+	  blTris.push_back(new MTriangle(v,v22,v21));
+	  fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
+		  v->x(),v->y(),v->z(),
+		  v22->x(),v22->y(),v22->z(),
+		  v21->x(),v21->y(),v21->z());
 	}
       }
     }
   }
+
   fprintf(ff2,"};\n");
   fclose(ff2);
+  
+  std::vector<MElement*> els,newels,oldels;
+  for (int i=0;i<blQuads.size();i++)els.push_back(blQuads[i]);
+  filterOverlappingElements (2,els,newels,oldels);
+  blQuads.clear();
+  for (int i=0;i<newels.size();i++)blQuads.push_back((MQuadrangle*)newels[i]);
+  for (int i=0;i<oldels.size();i++)delete oldels[i];
+  
+  for (int i=0;i<blQuads.size();i++){    
+    addOrRemove(blQuads[i]->getVertex(0),blQuads[i]->getVertex(1),bedges);
+    addOrRemove(blQuads[i]->getVertex(1),blQuads[i]->getVertex(2),bedges);
+    addOrRemove(blQuads[i]->getVertex(2),blQuads[i]->getVertex(3),bedges);
+    addOrRemove(blQuads[i]->getVertex(3),blQuads[i]->getVertex(0),bedges);
+    for (int j=0;j<4;j++) 
+      if(blQuads[i]->getVertex(j)->onWhat() == gf)verts.insert(blQuads[i]->getVertex(j));
+  }
+  for (int i=0;i<blTris.size();i++){    
+    addOrRemove(blTris[i]->getVertex(0),blTris[i]->getVertex(1),bedges);
+    addOrRemove(blTris[i]->getVertex(1),blTris[i]->getVertex(2),bedges);
+    addOrRemove(blTris[i]->getVertex(2),blTris[i]->getVertex(0),bedges);
+    for (int j=0;j<3;j++) 
+      if(blTris[i]->getVertex(j)->onWhat() == gf)verts.insert(blTris[i]->getVertex(j));
+  }
 
   discreteEdge ne (gf->model(), 444444,0,
 		   (*edges.begin())->getEndVertex());
@@ -1026,7 +1027,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
     outputScalarField(m->triangles, name, 1);
   }
   
-  if(0){
+  if(1){
     std::list<BDS_Face*>::iterator itt = m->triangles.begin();
     while (itt != m->triangles.end()){
       BDS_Face *t = *itt;
@@ -1744,7 +1745,7 @@ void deMeshGFace::operator() (GFace *gf)
 }
 
 // for debugging, change value from -1 to -100;
-int debugSurface = -1; //-1; 
+int debugSurface = -100; //-1; 
 
 void meshGFace::operator() (GFace *gf)
 {
diff --git a/Mesh/meshGFaceBoundaryLayers.cpp b/Mesh/meshGFaceBoundaryLayers.cpp
index 09693f07692f0ad8760155c801ac8bc1066d31b0..d6dffcf635e1c625b58fb6f6ddf45ff5d3a30690 100644
--- a/Mesh/meshGFaceBoundaryLayers.cpp
+++ b/Mesh/meshGFaceBoundaryLayers.cpp
@@ -100,7 +100,6 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf, double _treshold) {
 
   FieldManager *fields = gf->model()->getFields();
   if(fields->boundaryLayer_field <= 0){
-    printf("no bl\n");
     return 0;
   }
   Field *bl_field = fields->get(fields->boundaryLayer_field);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 2cdaab0df3535bc0063c0e7674cd76cdebe352b9..db2b14c9d3a1c6c41d8a2068076c6d27020c2ef4 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -16,7 +16,6 @@
 #include "Numeric.h"
 #include "STensor3.h"
 #include "Context.h"
-#include "meshGFaceBoundaryLayers.h"
 #include "MLine.h"
 #include "MQuadrangle.h"
 
@@ -804,136 +803,6 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
   }
 }
 
-// initially adds some boundary layer points onto the mesh
-void addBoundaryLayerPoints  (GFace *gf){
-
-  std::set<MTri3*,compareTri3Ptr> AllTris;
-  std::vector<double> vSizes, vSizesBGM, Us, Vs;
-  std::vector<SMetric3> vMetricsBGM;
-
-  BoundaryLayerColumns *_columns = buidAdditionalPoints2D (gf, M_PI/6.);
-
-  buildMeshGenerationDataStructures
-    (gf, AllTris, vSizes, vSizesBGM, vMetricsBGM, Us, Vs);
-
-  for (BoundaryLayerColumns::iter it = _columns->begin(); it != _columns->end(); ++it){
-    for (int i=0;i<it->second._column.size();++i){
-      while(1){
-	MTri3 *worst = *AllTris.begin();
-	if (worst->isDeleted()){
-	  delete worst->tri();
-	  delete worst;
-	  AllTris.erase(AllTris.begin());
-	}
-	else {
-	  break;
-	}
-      }
-      MVertex *v = it->second._column[i];
-      SMetric3 mmm = it->second._metrics[i];
- 
-      double center[2],metric[3] = {1.,0,1.};
-      v->getParameter(0,center[0]);
-      v->getParameter(1,center[1]);
-      buildMeshMetric (gf,center,mmm,metric);
-      v->setIndex(Us.size());
-      vSizesBGM.push_back(vSizesBGM[((it->first))->getIndex()]);
-      vSizes.push_back(vSizes[((it->first))->getIndex()]);
-      Us.push_back(center[0]);
-      Vs.push_back(center[1]);
-      MTri3 *ptin =  search4Triangle (*AllTris.begin(), center, Us, Vs, AllTris);
-      if (ptin){
-	bool ok  = insertVertex
-	  (true, gf, v, center, ptin, AllTris,0, vSizes, vSizesBGM,vMetricsBGM, 
-	   Us, Vs, metric);
-	if (ok)gf->mesh_vertices.push_back(v);
-	else printf("failed !\n");
-      }
-      else {
-	printf("boundary layer point outside the domain ;-)\n");
-      }
-    }
-  } 
-
-  // Create all edges of the triangular mesh
-  std::set<MEdge,Less_Edge> allEdges;
-  {
-    std::set<MTri3*,compareTri3Ptr>::iterator it = AllTris.begin();
-    for (;it != AllTris.end();++it ){
-      MTri3 *worst = *it;
-      if (!worst->isDeleted()){
-	allEdges.insert(MEdge((*it)->tri()->getVertex(0),(*it)->tri()->getVertex(1)));
-	allEdges.insert(MEdge((*it)->tri()->getVertex(0),(*it)->tri()->getVertex(2)));
-	allEdges.insert(MEdge((*it)->tri()->getVertex(1),(*it)->tri()->getVertex(2)));
-      }
-    }
-  }
-  //  printf("%d edges generated\n",allEdges.size());
-
-  //  look for existing cavities and fill them up with quads
-  std::list<GEdge*> edges = gf->edges();
-  std::list<GEdge*>::iterator ite = edges.begin();
-  while(ite != edges.end()){
-    for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
-      MVertex *v1 = (*ite)->lines[i]->getVertex(0);
-      MVertex *v2 = (*ite)->lines[i]->getVertex(1);
-      int nbCol1 = _columns->getNbColumns(v1);
-      int nbCol2 = _columns->getNbColumns(v2);
-      if (nbCol1 >= 0 && nbCol2 >= 0){ 
-	const BoundaryLayerData & c1 = _columns->getColumn(v1,MEdge(v1,v2));
-	const BoundaryLayerData & c2 = _columns->getColumn(v2,MEdge(v1,v2));
-	int N = std::min(c1._column.size(),c2._column.size());
-	
-	// Look for the highest edge on the column
-	int highestHorizontal = -1;
-	int highestVertical = -1;
-	for (int l=0;l < N ;++l){
-	  MVertex *v11,*v12,*v21,*v22;
-	  v21 = c1._column[l];
-	  v22 = c2._column[l];	    
-	  if (l == 0){
-	    v11 = v1;
-	    v12 = v2;
-	  }
-	  else {
-	    v11 = c1._column[l-1];
-	    v12 = c2._column[l-1];	    
-	  }
-	  if (allEdges.find(MEdge(v11,v21)) != allEdges.end() && allEdges.find(MEdge(v12,v22)) != allEdges.end()){
-	    highestVertical = l;	    
-	  }
-	  else {
-	    break;
-	  }
-	  
-	  if (allEdges.find(MEdge(v21,v22))!= allEdges.end()){
-	    highestHorizontal = l;
-	  }
-	}
-	//	printf("%d vs %d %d\n",N,highestHorizontal,highestVertical);
-	N = std::min(highestHorizontal,highestVertical);
-		//N--;
-	for (int l=0;l <= N ;++l){
-	  MVertex *v11,*v12,*v21,*v22;
-	  v21 = c1._column[l];
-	  v22 = c2._column[l];	    
-	  if (l == 0){
-	    v11 = v1;
-	    v12 = v2;
-	  }
-	  else {
-	    v11 = c1._column[l-1];
-	    v12 = c2._column[l-1];	    
-	  }
-	  //	printf("quadrangle generated\n");
-	  gf->quadrangles.push_back(new MQuadrangle(v11,v12,v22,v21));
-	}
-      }
-    }
-    ++ite;
-  }
-  transferDataStructure(gf, AllTris, Us, Vs); 
-}
 
 void bowyerWatson(GFace *gf)
 {
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index ce4f189839760d02d98d1b418a7ab787ee69ee1c..6ac6654611af38ce2a4053121a3e6dabc6342a9c 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -98,7 +98,6 @@ void bowyerWatson(GFace *gf);
 void bowyerWatsonFrontal(GFace *gf);
 void bowyerWatsonFrontalLayers(GFace *gf, bool quad);
 void buildBackGroundMesh (GFace *gf);
-void addBoundaryLayerPoints  (GFace *gf);
 
 struct edgeXface
 {
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 73cd67ddaa74b1b25331334386bf148136946dee..296e607ffa9019253586369b45528ab9216661ec 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -512,7 +512,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> &regions)
     }
     else {
       sprintf(opts, "pe%c",  (Msg::GetVerbosity() < 3) ? 'Q': 
-	      (Msg::GetVerbosity() > 6) ? 'V': '\0');
+      	      (Msg::GetVerbosity() > 6) ? 'V': '\0');
     }
     try{
       tetrahedralize(opts, &in, &out);
@@ -936,8 +936,16 @@ bool buildFaceSearchStructure(GModel *model, fs_cont &search)
 {  
   search.clear();
 
-  GModel::fiter fit = model->firstFace();
-  while(fit != model->lastFace()){    
+  std::set<GFace*> faces_to_consider;
+  GModel::riter rit = model->firstRegion();
+  while(rit != model->lastRegion()){
+    std::list <GFace *> _faces = (*rit)->faces();
+    faces_to_consider.insert( _faces.begin(),_faces.end());
+    rit++;
+  }    
+
+  std::set<GFace*>::iterator fit = faces_to_consider.begin();
+  while(fit != faces_to_consider.end()){    
     for(unsigned int i = 0; i < (*fit)->triangles.size(); i++){
       MVertex *p1 = (*fit)->triangles[i]->getVertex(0);
       MVertex *p2 = (*fit)->triangles[i]->getVertex(1);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 8d52c8cb5f6e529aee1839403693c84231f0b6af..16dea86f4bb54e9dfcf294e667b43581fa639356 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -322,8 +322,19 @@ GRegion *getRegionFromBoundingFaces(GModel *model,
                                     std::set<GFace *> &faces_bound)
 {
   GModel::riter git = model->firstRegion();
+  //  for (std::set<GFace *>::iterator it = faces_bound.begin();
+  //       it != faces_bound.end(); ++it){
+  //    printf(" %d",(*it)->tag());
+  //  }
+  //  printf("\n");
   while (git != model->lastRegion()){
     std::list <GFace *> _faces = (*git)->faces();
+    //    printf("region %d %d faces\n",(*git)->tag(),_faces.size());
+    //    for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
+    //      printf(" %d",(*it)->tag());      
+    //    }
+    //  printf("\n");
+	   
     if(_faces.size() == faces_bound.size()){
       bool ok = true;
       for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
@@ -783,7 +794,7 @@ void insertVerticesInRegion (GRegion *gr)
       Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
                  theRegion.size(), faces_bound.size(), _t2 - _t1);
       GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
-      // Msg::Info("a region is found %p",myGRegion);
+       Msg::Info("a region is found %p",myGRegion);
       if(myGRegion) // a geometrical region associated to the list of faces has been found
         for(std::list<MTet4*>::iterator it2 = theRegion.begin(); 
             it2 != theRegion.end(); ++it2) (*it2)->setOnWhat(myGRegion);