diff --git a/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp b/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp
index 6d5e7539294f2789fd7cba7872f1cddecfd0e1ae..6b1feed3c029691dacc2006764eb51ff53398b58 100644
--- a/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp
+++ b/include/Algo/ImplicitHierarchicalMesh/subdivision3.hpp
@@ -85,9 +85,7 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position
 		++degree;
 		p += position[it] ;
 		if(!map.edgeIsSubdivided(it))							// first cut the edges (if they are not already)
-		{
 			Algo::IHM::subdivideEdge<PFP>(map, it, position) ;	// and compute the degree of the face
-		}
 		it = map.phi1(it) ;
 	} while(it != old) ;
 	p /= typename PFP::REAL(degree) ;
@@ -96,38 +94,38 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position
 
 	Dart res;
 
-//	if(degree == 3 && sType ==  S_TRI)		//subdiviser une face triangulaire
-//	{
-//		Dart dd = map.phi1(old) ;
-//		Dart e = map.phi1(map.phi1(dd)) ;
-//		map.splitFace(dd, e) ;					// insert a new edge
-//		unsigned int id = map.getNewEdgeId() ;
-//		map.setEdgeId(map.phi_1(dd), id, EDGE) ;		// set the edge id of the inserted edge to the next available id
-//
-//		unsigned int idface = map.getFaceId(old);
-//		map.setFaceId(dd, idface, FACE) ;
-//		map.setFaceId(e, idface, FACE) ;
-//
-//		dd = e ;
-//		e = map.phi1(map.phi1(dd)) ;
-//		map.splitFace(dd, e) ;
-//		id = map.getNewEdgeId() ;
-//		map.setEdgeId(map.phi_1(dd), id, EDGE) ;
-//
-//		map.setFaceId(dd, idface, FACE) ;
-//		map.setFaceId(e, idface, FACE) ;
-//
-//		dd = e ;
-//		e = map.phi1(map.phi1(dd)) ;
-//		map.splitFace(dd, e) ;
-//		id = map.getNewEdgeId() ;
-//		map.setEdgeId(map.phi_1(dd), id, EDGE) ;
-//
-//		map.setFaceId(dd, idface, FACE) ;
-//		map.setFaceId(e, idface, FACE) ;
-//	}
-//	else
-//	{
+	if(degree == 3 && sType ==  S_TRI)		//subdiviser une face triangulaire
+	{
+		Dart dd = map.phi1(old) ;
+		Dart e = map.phi1(map.phi1(dd)) ;
+		map.splitFace(dd, e) ;					// insert a new edge
+		unsigned int id = map.getNewEdgeId() ;
+		map.setEdgeId(map.phi_1(dd), id, EDGE) ;		// set the edge id of the inserted edge to the next available id
+
+		unsigned int idface = map.getFaceId(old);
+		map.setFaceId(dd, idface, FACE) ;
+		map.setFaceId(e, idface, FACE) ;
+
+		dd = e ;
+		e = map.phi1(map.phi1(dd)) ;
+		map.splitFace(dd, e) ;
+		id = map.getNewEdgeId() ;
+		map.setEdgeId(map.phi_1(dd), id, EDGE) ;
+
+		map.setFaceId(dd, idface, FACE) ;
+		map.setFaceId(e, idface, FACE) ;
+
+		dd = e ;
+		e = map.phi1(map.phi1(dd)) ;
+		map.splitFace(dd, e) ;
+		id = map.getNewEdgeId() ;
+		map.setEdgeId(map.phi_1(dd), id, EDGE) ;
+
+		map.setFaceId(dd, idface, FACE) ;
+		map.setFaceId(e, idface, FACE) ;
+	}
+	else
+	{
 		Dart dd = map.phi1(old) ;
 		map.splitFace(dd, map.phi1(map.phi1(dd))) ;
 
@@ -164,7 +162,7 @@ void subdivideFace(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position
 			dd = map.phi2(map.phi1(dd));
 		}
 		while(dd != ne);
-//	}
+	}
 
 	map.setCurrentLevel(cur) ;
 }
@@ -182,8 +180,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 	map.setCurrentLevel(vLevel);
 
 
-	std::cout << "vLevel = " << vLevel << std::endl;
-
 	/*
 	 * au niveau du volume courant i
 	 * stockage d'un brin de chaque face de celui-ci
@@ -191,7 +187,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 	 */
 
 	DartMarkerStore mf(map);		// Lock a face marker to save one dart per face
-	//DartMarkerStore mv(map);		// Lock a vertex marker to compute volume center
 	CellMarker mv(map, VERTEX);
 
 	typename PFP::VEC3 volCenter;
@@ -217,7 +212,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 			//compute volume centroid
 			if(!mv.isMarked(e))
 			{
-				//mv.markOrbit(VERTEX, e);
 				mv.mark(e);
 				volCenter += position[e];
 				++count;
@@ -246,7 +240,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 	std::vector<std::pair<Dart,Dart> > subdividedfaces;
 	subdividedfaces.reserve(25);
 
-	int i = 0;
 	//First step : subdivide edges and faces
 	//creates a i+1 edge level and i+1 face level
 	for (std::vector<Dart>::iterator face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
@@ -255,24 +248,12 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 
 		//if needed subdivide face
 		if(!map.faceIsSubdivided(d))
-		{
-			++i;
 			Algo::IHM::subdivideFace<PFP>(map, d, position);
-		}
-
-
-		std::cout << "CHECK FACES " << std::endl;
-		map.check();
 
 		//save a dart from the subdivided face
 		unsigned int cur = map.getCurrentLevel() ;
 
-		map.setCurrentLevel(map.getMaxLevel());
-		Dart old = map.faceOldestDart(d);
-		std::cout << "                  Oldest dart Face = " << old << std::endl;
-		unsigned int fLevel = map.faceLevel(old); //puisque dans tous les cas, la face est subdivisee
-		std::cout << "FACE LEVEL = " << fLevel << std::endl;
-
+		unsigned int fLevel = map.faceLevel(d) + 1; //puisque dans tous les cas, la face est subdivisee
 		map.setCurrentLevel(fLevel) ;
 
 
@@ -288,8 +269,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 		map.setCurrentLevel(cur);
 	}
 
-
-
 	map.setCurrentLevel(vLevel + 1) ; // go to the next level to perform volume subdivision
 
 	std::vector<Dart> newEdges;	//save darts from inner edges
@@ -301,33 +280,21 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 		Dart e = *edge;
 
 		Dart f1 = map.phi1(*edge);
-		int i  = 0;
+
 		do
 		{
-			std::cout << "corner #" << i << "  degree = " << map.faceDegree(e) << std::endl;
-			if(map.phi2(map.phi1(map.phi1(e))) != map.phi1(map.phi1(e)))
-				map.unsewFaces(map.phi1(map.phi1(e)));
+			map.unsewFaces(map.phi1(map.phi1(e)));
 
 			//TODO utile ?
-			if(map.phi2(map.phi1(e)) != map.phi1(e))
+			//if(map.phi2(map.phi1(e)) != map.phi1(e))
 				map.unsewFaces(map.phi1(e));
 
-			++i;
 			e = map.phi2(map.phi_1(e));
 		}
 		while(e != *edge);
 
 		map.closeHole(f1);
 
-		std::cout << "CHECK CLOSE HOLE " << std::endl;
-		map.check();
-
-		std::cout << "corner degree = " << i << std::endl;
-
-		std::cout << "face degree = " << map.faceDegree(map.phi2(f1)) << std::endl;
-//		if(map.faceDegree(map.phi2(f1)) > 6)
-//			return Dart::nil();
-
 		Dart old = map.phi2(map.phi1(e));
 		Dart dd = map.phi1(map.phi1(old)) ;
 		map.splitFace(old,dd) ;
@@ -395,8 +362,6 @@ Dart subdivideVolume(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& positi
 
 	map.setCurrentLevel(cur) ;
 
-	std::cout << std::endl;
-
 	return subdividedfaces.begin()->first;
 }
 
diff --git a/include/Geometry/bounding_box.h b/include/Geometry/bounding_box.h
index 06626778c33cc2968b9afc212b6b69b9eddaf116..bc1c639c8bbdfca117f81492c27c76b60b694d7d 100644
--- a/include/Geometry/bounding_box.h
+++ b/include/Geometry/bounding_box.h
@@ -86,6 +86,12 @@ class BoundingBox
 		// fusion with the given bounding box
 		void fusion(const BoundingBox<VEC>& bb) ;
 
+		//return true if the vector belongs strictly to a bounding box
+		bool contains(const VEC& p);
+
+		//return true if the bounding box belongs strictly to a bounding box
+		bool contains(const BoundingBox<VEC> & bb);
+
 		/**********************************************/
 		/*             STREAM OPERATORS               */
 		/**********************************************/
diff --git a/include/Geometry/bounding_box.hpp b/include/Geometry/bounding_box.hpp
index 5be1369dca3d1ecfa8572da72be7ad1ac542eca7..6db89968dd361c86f0c60e45b993d1c4082cceb1 100644
--- a/include/Geometry/bounding_box.hpp
+++ b/include/Geometry/bounding_box.hpp
@@ -178,6 +178,30 @@ void BoundingBox<VEC>::fusion(const BoundingBox<VEC>& bb)
 	}
 }
 
+
+template <typename VEC>
+bool BoundingBox<VEC>::contains(const VEC& p)
+{
+	for(unsigned int i = 0; i < m_pMin.dimension(); ++i)
+	{
+		if(m_pMin[i] > p[i])
+			return false ;
+		if(p[i] > m_pMax[i])
+			return false ;
+	}
+
+	return true;
+}
+
+template <typename VEC>
+bool BoundingBox<VEC>::contains(const BoundingBox<VEC>& bb)
+{
+	return this->contains(bb.min()) && this->contains(bb.max());
+}
+
+
+
+
 //template <typename VEC>
 //friend std::ostream& BoundingBox<VEC>::operator<<(std::ostream& out, const BoundingBox<VEC>& bb)
 //{
diff --git a/src/Topology/map/map3.cpp b/src/Topology/map/map3.cpp
index d87ba2c2e252ad7cfdc997b129883ebe2f395ffd..3a7d71bf369b744e6c48e691c85a833cdb79398f 100644
--- a/src/Topology/map/map3.cpp
+++ b/src/Topology/map/map3.cpp
@@ -1204,6 +1204,7 @@ bool Map3::check()
     {
         if (!m.isMarked(d)) // phi1 a au moins un antécédent ?
 		{
+        	std::cout << "dart = " << d << std::endl;
             CGoGNout << "Check: dart with no phi1 predecessor" << CGoGNendl;
             return false;
         }