diff --git a/Mesh_2/include/CGAL/Meshes/Double_map_container.h b/Mesh_2/include/CGAL/Meshes/Double_map_container.h
index 284121795e0..fa8ac113d76 100644
--- a/Mesh_2/include/CGAL/Meshes/Double_map_container.h
+++ b/Mesh_2/include/CGAL/Meshes/Double_map_container.h
@@ -74,7 +74,7 @@ namespace CGAL {
         m.pop_front();
       }
 
-      void remove_element(Element& e)
+      void remove_element(const Element& e)
       {
 #if CGAL_MESHES_DEBUG_DOUBLE_MAP
 	std::cerr << "remove_element(" << &*e << ")\n";
diff --git a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h
index 00ab4c4e478..aba5b56c898 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Dump_c3t3.h
@@ -58,19 +58,35 @@ template <typename C3t3>
 struct Dump_c3t3<C3t3, false> {
   void dump_c3t3(const C3t3&, std::string) {
     std::cerr << "Warning " << __FILE__ << ":" << __LINE__ << "\n"
-              << "  the c3t3 object cannot be dumped because some types are"
-              << " not streamable:\n";
-    if(!is_streamable<typename C3t3::Triangulation::Vertex>::value)
+              << "  the c3t3 object of following type:\n"
+              << typeid(C3t3).name() << std::endl
+              << "  cannot be dumped because some types are not streamable:\n";
+    if(!is_streamable<typename C3t3::Triangulation::Vertex>::value) {
       std::cerr << "     - C3t3::Triangulation::Vertex is not streamble\n";
+      std::cerr << "       "
+                << typeid(typename C3t3::Triangulation::Vertex).name()
+                << "\n";
+    }
 
-    if(!is_streamable<typename C3t3::Triangulation::Cell>::value)
+    if(!is_streamable<typename C3t3::Triangulation::Cell>::value) {
       std::cerr << "     - C3t3::Triangulation::Cell is not streamble\n";
+      std::cerr << "       "
+                << typeid(typename C3t3::Triangulation::Cell).name()
+                << "\n";
+    }
 
-    if(!is_streamable<typename C3t3::Surface_patch_index>::value)
+    if(!is_streamable<typename C3t3::Surface_patch_index>::value) {
       std::cerr << "     - C3t3::Surface_patch_index is not streamable\n";
-      
-    if(!is_streamable<typename C3t3::Subdomain_index>::value)
+      std::cerr << "       "
+                << typeid(typename C3t3::Surface_patch_index).name()
+                << "\n";
+    }
+    if(!is_streamable<typename C3t3::Subdomain_index>::value) {
       std::cerr << "     - C3t3::Subdomain_index is not streamable\n";      
+      std::cerr << "       "
+                << typeid(typename C3t3::Subdomain_index).name()
+                << "\n";
+    }
   }
 }; // end struct template specialization Dump_c3t3<C3t3, false>
 
diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h
index 0d0d9ed3a23..026db81e0bb 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_complex_3_in_triangulation_3_base.h
@@ -1,4 +1,5 @@
-// Copyright (c) 2009 INRIA Sophia-Antipolis (France).
+// Copyright (c) 2003-2009  INRIA Sophia-Antipolis (France).
+// Copyright (c) 2013       GeometryFactory Sarl (France).
 // All rights reserved.
 //
 // This file is part of CGAL (www.cgal.org).
@@ -16,7 +17,7 @@
 // $Id$
 //
 //
-// Author(s)     : Stéphane Tayeb
+// Author(s)     : Laurent Rineau, Stéphane Tayeb
 //
 //******************************************************************************
 // File Description : Implements class Mesh_complex_3_in_triangulation_3.
@@ -35,6 +36,7 @@
 #include <iostream>
 #include <fstream>
 #include <CGAL/Mesh_3/io_signature.h>
+#include <CGAL/Union_find.h>
 
 #ifdef CGAL_LINKED_WITH_TBB
   #include <tbb/atomic.h>
@@ -43,6 +45,46 @@
 namespace CGAL {
 namespace Mesh_3 {
 
+  namespace details {
+
+    template <typename Tr>
+    class C3t3_helper_class
+    {
+    protected:
+      typedef typename Tr::Vertex_handle Vertex_handle;
+      typedef typename Tr::Cell_handle   Cell_handle;
+      typedef typename Tr::Facet         Facet;
+      typedef typename Tr::Edge          Edge;
+
+      typedef std::pair<Vertex_handle, Vertex_handle> Pair_of_vertices;
+
+      // computes and return an ordered pair of Vertex
+      Pair_of_vertices
+      make_ordered_pair(const Vertex_handle vh1, const Vertex_handle vh2) const {
+	if (vh1 < vh2) {
+	  return std::make_pair(vh1, vh2);
+	}
+	else {
+	  return std::make_pair(vh2, vh1);
+	}
+      }
+
+      // same from an Edge
+      Pair_of_vertices
+      make_ordered_pair(const Edge e) const {
+        return make_ordered_pair(e.first->vertex(e.second),
+                                 e.first->vertex(e.third));
+      }
+
+      Facet canonical_facet(Cell_handle c, int i) const {
+        Cell_handle c2 = c->neighbor(i);
+        return (c2 < c) ? std::make_pair(c2,c2->index(c)) : std::make_pair(c,i);
+      }
+
+    }; // end class template C3t3_helper_class
+
+  } // end namespace Mesh_3::details
+
 /**
  * @class Mesh_complex_3_in_triangulation_3_base
  * @brief A data-structure to represent and maintain a 3D complex embedded
@@ -50,8 +92,10 @@ namespace Mesh_3 {
  */
 template<typename Tr, typename Concurrency_tag>
 class Mesh_complex_3_in_triangulation_3_base
+  : public details::C3t3_helper_class<Tr>
 {
   typedef Mesh_complex_3_in_triangulation_3_base<Tr, Concurrency_tag> Self;
+  typedef details::C3t3_helper_class<Tr> Base;
 
 public:
   // Triangulation types
@@ -67,6 +111,26 @@ public:
   typedef typename Tr::Cell::Surface_patch_index  Surface_patch_index;
   typedef typename Tr::Vertex::Index              Index;
 
+  enum Face_status{ NOT_IN_COMPLEX = 0,
+                    ISOLATED = 1, // - An ISOLATED edge is a marked edge,
+                                  //   without any incident facets.
+                    BOUNDARY,     // - An edge is on BOUNDARY if it has only
+                                  //   one incident facet.
+                                  // - A vertex is on BOUNDARY if all its
+                                  //   incident edges are REGULAR or on
+                                  //   BOUNDARY, at least one is on
+                                  //   BOUNDARY, and the incident facets
+                                  //   form only one connected component.
+                    REGULAR,      // - A facet that is in the complex is
+                                  //   REGULAR.
+                                  // - An edge is REGULAR if it has
+                                  //   exactly two incident facets.
+                                  // - A vertex is REGULAR if all it
+                                  //   incident edges are REGULAR, and the
+                                  //   incident facets form only one
+                                  //   connected component.
+                    SINGULAR};    // - SINGULAR is for all other cases.
+
   //-------------------------------------------------------
   // Constructors / Destructors
   //-------------------------------------------------------
@@ -75,7 +139,10 @@ public:
    * Builds an empty 3D complex.
    */
   Mesh_complex_3_in_triangulation_3_base()
-    : tr_()
+    : Base()
+    , tr_()
+    , edge_facet_counter_() //TODO: parallel!
+    , manifold_info_initialized_(false) //TODO: parallel!
   {
     // We don't put it in the initialization list because
     // tbb::atomic has no contructors
@@ -85,7 +152,10 @@ public:
 
   /// Copy constructor
   Mesh_complex_3_in_triangulation_3_base(const Self& rhs)
-    : tr_(rhs.tr_)
+    : Base()
+    , tr_(rhs.tr_)
+    , edge_facet_counter_(rhs.edge_facet_counter_)
+    , manifold_info_initialized_(rhs.manifold_info_initialized_)
   {
     number_of_facets_ = rhs.number_of_facets_;
     number_of_cells_ = rhs.number_of_cells_;
@@ -98,6 +168,8 @@ public:
     number_of_cells_ = 0;
     number_of_facets_ = 0;
     tr_.clear();
+    manifold_info_initialized_ = false;
+    edge_facet_counter_.clear();
   }
 
   /// Assignment operator
@@ -146,6 +218,87 @@ public:
     cell->set_surface_patch_index(i, index);
   }
 
+  /// Returns `NOT_IN_COMPLEX`, `BOUNDARY`, `REGULAR`, or `SINGULAR`,
+  /// depending on the number of incident facets in the complex, and the
+  /// number of connected components of its link
+  Face_status face_status(const Vertex_handle v) const
+  {
+    if(!manifold_info_initialized_) init_manifold_info();
+    const int n = v->cached_number_of_incident_facets();
+
+    if(n == 0) return NOT_IN_COMPLEX;
+
+    //test incident edges for REGULARITY and count BOUNDARY edges
+    typename std::vector<Edge> edges;
+    edges.reserve(64);
+    tr_.incident_edges(v, std::back_inserter(edges));
+    int number_of_boundary_incident_edges = 0; // could be a bool
+    for (typename std::vector<Edge>::iterator
+           eit=edges.begin(), end = edges.end();
+	 eit != end; eit++)
+    {
+      switch( face_status(*eit) )
+      {
+      case NOT_IN_COMPLEX: case REGULAR: break;
+      case BOUNDARY: ++number_of_boundary_incident_edges; break;
+      default :
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        std::cerr << "singular edge...\n";
+        std::cerr << v->point() << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        return SINGULAR;
+      }
+    }
+
+    // From here all incident edges (in complex) are REGULAR or BOUNDARY.
+    const int nb_components = union_find_of_incident_facets(v);
+    if(nb_components > 1) {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      std::cerr << "singular vertex: nb_components=" << nb_components << std::endl;
+      std::cerr << v->point() << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      return SINGULAR;
+    }
+    else { // REGULAR OR BOUNDARY
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      std::cerr << "regular or boundary: " << v->point() << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      if (number_of_boundary_incident_edges != 0)
+        return BOUNDARY;
+      else
+        return REGULAR;
+    }
+  }
+
+  /// This function should be called only when incident edges
+  /// are known to be REGULAR OR BOUNDARY
+  bool is_regular_or_boundary_for_vertices(Vertex_handle v) const {
+    return union_find_of_incident_facets(v) == 1;
+  }
+
+  /// Returns `NOT_IN_COMPLEX`, `BOUNDARY`, `REGULAR`, or `SINGULAR`,
+  /// depending on the number of incident facets in the complex
+  Face_status face_status(const Edge& edge) const
+  {
+    if(!manifold_info_initialized_) init_manifold_info();
+
+    switch(edge_facet_counter_[this->make_ordered_pair(edge)])
+    {
+    case 0: return NOT_IN_COMPLEX;
+    case 1: return BOUNDARY;
+    case 2: return REGULAR;
+    default: return SINGULAR;
+    }
+  }
+
+  /// Returns true if the vertex \c v has is incident to at least a facet
+  /// of the complex
+  bool has_incident_facets_in_complex(const Vertex_handle& v) const
+  {
+    if(!manifold_info_initialized_) init_manifold_info();
+    return v->cached_number_of_incident_facets() > 0;
+  }
+
   /// Returns true if facet \c facet is in complex
   bool is_in_complex(const Facet& facet) const
   {
@@ -314,6 +467,94 @@ public:
   /// Returns bbox
   Bbox_3 bbox() const;
 
+private:
+  void init_manifold_info() const {
+    for(typename Tr::All_vertices_iterator
+          vit = triangulation().finite_vertices_begin(),
+          end = triangulation().finite_vertices_end();
+        vit != end; ++vit)
+    {
+      vit->set_c2t3_cache(0, -1);
+    }
+    for(typename Tr::Finite_facets_iterator
+          fit = triangulation().finite_facets_begin(),
+          end = triangulation().finite_facets_end();
+        fit != end; ++fit)
+    {
+      if ( is_in_complex(*fit) ) {
+        const Cell_handle cell = fit->first;
+        const int i = fit->second;
+        for(int j = 0; j < 3; ++j)
+        {
+          const int edge_index_va = tr_.vertex_triple_index(i, j);
+          const int edge_index_vb = tr_.vertex_triple_index(i, (j == 2) ? 0 : (j+1));
+          const Vertex_handle edge_va = cell->vertex(edge_index_va);
+          const Vertex_handle edge_vb = cell->vertex(edge_index_vb);
+          ++edge_facet_counter_[this->make_ordered_pair(edge_va, edge_vb)];
+
+          const int n = edge_va->cached_number_of_incident_facets();
+          edge_va->set_c2t3_cache(n+1, -1);
+        }
+      }
+    }
+    manifold_info_initialized_ = true;
+  }
+
+  /// Extract the subset `F` of facets of the complex incident to `v` and
+  /// return the number of connected component of the adjacency graph of `F`.
+  int union_find_of_incident_facets(const Vertex_handle v) const
+  {
+    if( v->is_c2t3_cache_valid() )
+    {
+      const int n = v->cached_number_of_components();
+      if(n >= 0) return v->cached_number_of_components();
+    }
+
+    Union_find<Facet> facets;
+    { // fill the union find
+      std::vector<Facet> non_filtered_facets;
+      tr_.incident_facets(v, std::back_inserter(non_filtered_facets));
+
+      for(typename std::vector<Facet>::iterator
+            fit = non_filtered_facets.begin(),
+            end = non_filtered_facets.end();
+          fit != end; ++fit)
+      {
+        if(is_in_complex(*fit)) facets.push_back(*fit);
+      }
+    }
+
+    typedef std::map<Vertex_handle,
+                     typename Union_find<Facet>::handle> Vertex_set_map;
+    typedef typename Vertex_set_map::iterator Vertex_set_map_iterator;
+
+    Vertex_set_map vsmap;
+
+    for(typename Union_find<Facet>::iterator
+          it = facets.begin(), end = facets.end();
+        it != end; ++it)
+    {
+      const Cell_handle& ch = (*it).first;
+      const int& i = (*it).second;
+      for(int j=0; j < 3; ++j) {
+	const Vertex_handle w = ch->vertex(tr_.vertex_triple_index(i,j));
+	if(w != v){
+	  Vertex_set_map_iterator vsm_it = vsmap.find(w);
+	  if(vsm_it != vsmap.end()){
+	    facets.unify_sets(vsm_it->second, it);
+	  } else {
+	    vsmap.insert(std::make_pair(w, it));
+	  }
+	}
+      }
+    }
+    const int nb_components = static_cast<int>(facets.number_of_sets());
+
+    const int n = v->cached_number_of_incident_facets();
+    v->set_c2t3_cache(n, nb_components);
+    return nb_components;
+  }
+  
   //-------------------------------------------------------
   // Traversal
   //-------------------------------------------------------
@@ -525,8 +766,16 @@ private:
 
   // Private date members
   Triangulation tr_;
+
+  typedef typename Base::Pair_of_vertices Pair_of_vertices;
+  typedef std::map<Pair_of_vertices, int> Edge_facet_counter;
+
+  mutable Edge_facet_counter edge_facet_counter_;
+
   typename Number_of_elements<Concurrency_tag>::type number_of_facets_;
   typename Number_of_elements<Concurrency_tag>::type number_of_cells_;
+
+  mutable bool manifold_info_initialized_;
 };  // end class Mesh_complex_3_in_triangulation_3_base
 
 
@@ -545,6 +794,31 @@ Mesh_complex_3_in_triangulation_3_base<Tr,Ct>::add_to_complex(
     set_surface_patch_index(cell, i, index);
     set_surface_patch_index(mirror.first, mirror.second, index);
     ++number_of_facets_;
+    if(manifold_info_initialized_) {
+      for(int j = 0; j < 3; ++j)
+      {
+        int edge_index_va = tr_.vertex_triple_index(i, j);
+        int edge_index_vb = tr_.vertex_triple_index(i, (j == 2) ? 0 : (j+1));
+        Vertex_handle edge_va = cell->vertex(edge_index_va);
+        Vertex_handle edge_vb = cell->vertex(edge_index_vb);
+        ++edge_facet_counter_[this->make_ordered_pair(edge_va, edge_vb)];
+
+        const int n = edge_va->cached_number_of_incident_facets();
+        const int m = edge_va->cached_number_of_components();
+        edge_va->set_c2t3_cache(n+1, m);
+      }
+      const int dimension_plus_1 = tr_.dimension() + 1;
+      // update c2t3 for vertices of f
+      for (int j = 0; j < dimension_plus_1; j++) {
+        if (j != i) {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+          if(cell->vertex(j)->is_c2t3_cache_valid())
+            std::cerr << "(" << cell->vertex(j)->point() << ")->invalidate_c2t3_cache()\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+          cell->vertex(j)->invalidate_c2t3_cache();
+        }
+      }
+    }
   }
 }
 
@@ -559,6 +833,34 @@ Mesh_complex_3_in_triangulation_3_base<Tr,Ct>::remove_from_complex(const Facet&
     set_surface_patch_index(facet.first, facet.second, Surface_patch_index());
     set_surface_patch_index(mirror.first, mirror.second, Surface_patch_index());
     --number_of_facets_;
+    if(manifold_info_initialized_) {
+      const Cell_handle cell = facet.first;
+      const int i = facet.second;
+      for(int j = 0; j < 3; ++j)
+      {
+        const int edge_index_va = tr_.vertex_triple_index(i, j);
+        const int edge_index_vb = tr_.vertex_triple_index(i, (j == 2) ? 0 : (j+1));
+        const Vertex_handle edge_va = cell->vertex(edge_index_va);
+        const Vertex_handle edge_vb = cell->vertex(edge_index_vb);
+        --edge_facet_counter_[this->make_ordered_pair(edge_va, edge_vb)];
+
+        const int n = edge_va->cached_number_of_incident_facets();
+        CGAL_assertion(n>0);
+        const int m = edge_va->cached_number_of_components();
+        edge_va->set_c2t3_cache(n-1, m);
+      }
+      const int dimension_plus_1 = tr_.dimension() + 1;
+      // update c2t3 for vertices of f
+      for (int j = 0; j < dimension_plus_1; j++) {
+        if (j != facet.second) {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+          if(cell->vertex(j)->is_c2t3_cache_valid())
+            std::cerr << "(" << cell->vertex(j)->point() << ")->invalidate_c2t3_cache()\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+          cell->vertex(j)->invalidate_c2t3_cache();
+        }
+      }
+    }
   }
 }
 
diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h
index b9ef1c603c8..1a62495d315 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h
@@ -32,6 +32,7 @@
 #include <CGAL/Mesh_3/Dump_c3t3.h>
 
 #include<CGAL/Mesh_3/Refine_facets_3.h>
+#include<CGAL/Mesh_3/Refine_facets_manifold_base.h>
 #include<CGAL/Mesh_3/Refine_cells_3.h>
 #include <CGAL/Mesh_3/Refine_tets_visitor.h>
 #include <CGAL/Mesher_level_visitors.h>
diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
index cb2f91484b8..c181e4e92f1 100644
--- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
+++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h
@@ -337,10 +337,12 @@ template<class Tr,
                                        typename Criteria::Facet_quality>
 # endif
 #endif // CGAL_LINKED_WITH_TBB
+         , class Base_ = Refine_facets_3_base<typename MeshDomain::Index,
+                                              typename Tr::Facet,
+                                              Concurrency_tag>
 >
 class Refine_facets_3
-: public Refine_facets_3_base<typename MeshDomain::Index, typename Tr::Facet,
-                              Concurrency_tag>
+: public Base_
 , public Mesh_3::Mesher_level<Tr,
                       Refine_facets_3<Tr,
                                       Criteria,
@@ -348,7 +350,8 @@ class Refine_facets_3
                                       Complex3InTriangulation3,
                                       Previous_level_,
                                       Concurrency_tag,
-                                      Container_>,
+                                      Container_,
+                                      Base_>,
                       typename Tr::Facet,
                       Previous_level_,
                       Triangulation_mesher_level_traits_3<Tr>,
@@ -364,11 +367,10 @@ class Refine_facets_3
                           Complex3InTriangulation3,
                           Previous_level_,
                           Concurrency_tag,
-                          Container_>                  Self;
+                          Container_,
+                          Base_> Self;
 
-  typedef Refine_facets_3_base<typename MeshDomain::Index,
-                               typename Tr::Facet,
-                               Concurrency_tag>        Base;
+  typedef Base_ Base;
 
   typedef Mesher_level<Tr,
                       Refine_facets_3<Tr,
@@ -836,8 +838,8 @@ private:
 
 
 // For sequential
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 Refine_facets_3(Tr& triangulation,
                 const Cr& criteria,
                 const MD& oracle,
@@ -857,8 +859,8 @@ Refine_facets_3(Tr& triangulation,
 }
 
 // For parallel
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 Refine_facets_3(Tr& triangulation,
                 const Cr& criteria,
                 const MD& oracle,
@@ -881,9 +883,9 @@ Refine_facets_3(Tr& triangulation,
 }
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 scan_triangulation_impl()
 {
   typedef typename Tr::Finite_facets_iterator Finite_facet_iterator;
@@ -960,9 +962,9 @@ scan_triangulation_impl()
 
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 int
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 number_of_bad_elements_impl()
 {
   typedef typename Tr::Finite_facets_iterator Finite_facet_iterator;
@@ -1120,9 +1122,9 @@ number_of_bad_elements_impl()
 }
 
 // For sequential
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 Mesher_level_conflict_status
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 test_point_conflict_from_superior_impl(const Point& point, Zone& zone)
 {
   typedef typename Zone::Facets_iterator Facet_iterator;
@@ -1166,10 +1168,10 @@ test_point_conflict_from_superior_impl(const Point& point, Zone& zone)
 }
 
 // For parallel
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 template <typename Mesh_visitor>
 Mesher_level_conflict_status
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 test_point_conflict_from_superior_impl(const Point& point, Zone& zone,
                                        Mesh_visitor &visitor)
 {
@@ -1216,9 +1218,9 @@ test_point_conflict_from_superior_impl(const Point& point, Zone& zone,
 }
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
-typename Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::Zone
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
+typename Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::Zone
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 conflicts_zone_impl(const Point& point
                     , const Facet& facet
                     , bool &facet_is_in_its_cz)
@@ -1275,9 +1277,9 @@ conflicts_zone_impl(const Point& point
   return zone;
 }
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
-typename Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::Zone
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
+typename Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::Zone
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 conflicts_zone_impl(const Point& point
                     , const Facet& facet
                     , bool &facet_is_in_its_cz
@@ -1338,9 +1340,9 @@ conflicts_zone_impl(const Point& point
 
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 before_insertion_impl(const Facet& facet,
                       const Point& point,
                       Zone& zone)
@@ -1418,9 +1420,9 @@ before_insertion_impl(const Facet& facet,
 
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
-typename Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::Vertex_handle
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
+typename Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::Vertex_handle
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 insert_impl(const Point& point,
             const Zone& zone)
 {
@@ -1447,9 +1449,9 @@ insert_impl(const Point& point,
 
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 restore_restricted_Delaunay(const Vertex_handle& vertex)
 {
   typedef std::vector<Cell_handle> Cell_handle_vector;
@@ -1481,9 +1483,9 @@ restore_restricted_Delaunay(const Vertex_handle& vertex)
 //-------------------------------------------------------
 // Private methods
 //-------------------------------------------------------
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 treat_new_facet(Facet& facet)
 {
   // Treat facet
@@ -1525,9 +1527,9 @@ treat_new_facet(Facet& facet)
   set_facet_visited(mirror);
 }
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 dual_segment(const Facet & facet, Bare_point& p, Bare_point& q) const
 {
   Cell_handle c = facet.first;
@@ -1546,9 +1548,9 @@ dual_segment(const Facet & facet, Bare_point& p, Bare_point& q) const
       n->vertex(3)->point());
 }
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 dual_segment_exact(const Facet & facet, Bare_point& p, Bare_point& q) const
 {
   Cell_handle c = facet.first;
@@ -1569,9 +1571,9 @@ dual_segment_exact(const Facet & facet, Bare_point& p, Bare_point& q) const
       true);
 }
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 dual_ray(const Facet & facet, Ray_3& ray) const
 {
   Cell_handle c = facet.first;
@@ -1604,9 +1606,9 @@ dual_ray(const Facet & facet, Ray_3& ray) const
         n->vertex(3)->point()), l);
 }
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 dual_ray_exact(const Facet & facet, Ray_3& ray) const
 {
   Cell_handle c = facet.first;
@@ -1644,9 +1646,9 @@ dual_ray_exact(const Facet & facet, Ray_3& ray) const
 
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 compute_facet_properties(const Facet& facet,
                          Facet_properties& fp,
                          bool force_exact) const
@@ -1759,9 +1761,9 @@ compute_facet_properties(const Facet& facet,
 }
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 bool
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 is_facet_encroached(const Facet& facet,
                     const Point& point) const
 {
@@ -1783,9 +1785,9 @@ is_facet_encroached(const Facet& facet,
   return ( compare_distance(center, reference_point, point) != CGAL::SMALLER );
 }
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 bool
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 is_encroached_facet_refinable(Facet& facet) const
 {
   typedef typename Gt::Point_3 Point_3;
@@ -1872,9 +1874,9 @@ is_encroached_facet_refinable(Facet& facet) const
   * \c facet is an internal facet we are going to remove
   * \c source_facet is the facet we want to refine by inserting a new point
   */
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 bool
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 before_insertion_handle_facet_in_conflict_zone(Facet& facet,
                                                const Facet& source_facet)
 {
@@ -1899,9 +1901,9 @@ before_insertion_handle_facet_in_conflict_zone(Facet& facet,
 
 
 
-template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_>
+template<class Tr, class Cr, class MD, class C3T3_, class P_, class Ct, class C_, class B_>
 void
-Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
+Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,C_,B_>::
 after_insertion_handle_incident_facet(Facet& facet)
 {
   // If the facet is infinite or has been already visited,
diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h
new file mode 100644
index 00000000000..d47fc3b27d1
--- /dev/null
+++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_manifold_base.h
@@ -0,0 +1,627 @@
+// Copyright (c) 2003-2007  INRIA Sophia-Antipolis (France).
+// Copyright (c) 2008,2013  GeometryFactory, Sophia Antipolis (France)
+//
+// This file is part of CGAL (www.cgal.org).
+// You can redistribute it and/or modify it under the terms of the GNU
+// General Public License as published by the Free Software Foundation,
+// either version 3 of the License, or (at your option) any later version.
+//
+// Licensees holding a valid commercial license may use this file in
+// accordance with the commercial license agreement provided with the software.
+//
+// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+//
+// $URL$
+// $Id$
+//
+//
+// Author(s)     : Steve Oudot, David Rey, Mariette Yvinec, Laurent Rineau, Andreas Fabri
+
+#ifndef CGAL_MESH_3_REFINE_FACETS_MANIFOLD_BASE_H
+#define CGAL_MESH_3_REFINE_FACETS_MANIFOLD_BASE_H
+
+#include <CGAL/utility.h>
+#include <set>
+#include <vector>
+#include <boost/bimap.hpp>
+#include <boost/bimap/unordered_set_of.hpp>
+#include <boost/bimap/multiset_of.hpp>
+#include <boost/foreach.hpp>
+#include <boost/mpl/has_xxx.hpp>
+
+namespace CGAL {
+
+namespace Mesh_3 {
+
+BOOST_MPL_HAS_XXX_TRAIT_DEF(Has_manifold_criterion)
+
+template <typename Criteria,
+          bool has_Has_manifold_criterion =
+          has_Has_manifold_criterion<Criteria>::value>
+struct Has_manifold_criterion :
+    public CGAL::Boolean_tag<Criteria::Has_manifold_criterion::value>
+// when Criteria has the nested type Has_manifold_criterion
+{};
+
+template <typename Criteria>
+struct Has_manifold_criterion<Criteria, false> : public CGAL::Tag_false
+// when Criteria does not have the nested type Has_manifold_criterion
+{};
+
+
+
+template <typename Criteria>
+bool get_with_manifold(const Criteria& c, CGAL::Tag_false)
+{
+  return false;
+}
+
+template <typename Criteria>
+bool get_with_manifold(const Criteria& c, CGAL::Tag_true)
+{
+  return (c.topology() & MANIFOLD_WITH_BOUNDARY) != 0;
+}
+
+template <typename Criteria>
+bool get_with_manifold(const Criteria& c)
+{
+  return get_with_manifold(c, Has_manifold_criterion<Criteria>());
+}
+
+template <typename Criteria>
+bool get_with_boundary(const Criteria& c, CGAL::Tag_false)
+{
+  return false;
+}
+
+template <typename Criteria>
+bool get_with_boundary(const Criteria& c, CGAL::Tag_true)
+{
+  return (c.topology() & NO_BOUNDARY) == 0;
+}
+
+template <typename Criteria>
+bool get_with_boundary(const Criteria& c)
+{
+  return get_with_boundary(c, Has_manifold_criterion<Criteria>());
+}
+
+template <
+  typename Base_,
+  bool withBoundary = false
+  >
+class Refine_facets_manifold_base
+  : public Base_
+{
+public:
+
+  typedef Base_ Base ;
+
+  typedef typename Base::C3T3 C3t3;
+  typedef typename Base::Criteria Criteria;
+  typedef typename Base::Mesh_domain Mesh_domain;
+  typedef typename C3t3::Triangulation Tr;
+  typedef typename Tr::Point Point;
+  typedef typename Tr::Facet Facet;
+  typedef typename Tr::Vertex_handle Vertex_handle;
+
+  typedef typename Triangulation_mesher_level_traits_3<Tr>::Zone Zone;
+
+protected:
+  typedef typename Tr::Geom_traits GT;
+  typedef typename GT::FT FT;
+  typedef typename Tr::Edge Edge;
+  typedef typename Tr::Cell_handle Cell_handle;
+
+  typedef typename Tr::Facet_circulator Tr_facet_circulator;
+  typedef std::pair<Vertex_handle, Vertex_handle> EdgeVV;
+
+protected:
+  Tr& r_tr_;
+  C3t3& r_c3t3_;
+
+  typedef ::boost::bimap< EdgeVV,
+                          ::boost::bimaps::multiset_of<int> > Bad_edges;
+  typedef typename Bad_edges::value_type Bad_edge;
+
+  mutable Bad_edges m_bad_edges;
+  mutable std::set<Vertex_handle> m_bad_vertices;
+
+  mutable bool m_manifold_info_initialized;
+  mutable bool m_bad_vertices_initialized;
+  bool m_with_manifold_criterion;
+  bool m_with_boundary;
+
+private:
+  // computes and return an ordered pair of Vertex
+  EdgeVV make_edgevv(const Vertex_handle vh1,
+                     const Vertex_handle vh2) const {
+    if (vh1 < vh2) {
+      return std::make_pair(vh1, vh2);
+    }
+    else {
+      return std::make_pair(vh2, vh1);
+    }
+  }
+
+  // computes and returns the Edge type object from the EdgeVV object
+  // use tr.is_edge(...)
+  Edge edgevv_to_edge(const EdgeVV& arete) const {
+    Vertex_handle v1 = arete.first;
+    Vertex_handle v2 = arete.second;
+    Cell_handle c;
+    int index1, index2;
+
+    CGAL_assertion_code(bool is_edge =)
+    r_tr_.is_edge(v1, v2, c, index1, index2);
+    CGAL_assertion(is_edge);
+
+    return make_triple( c, index1, index2 );
+  }
+
+  // computes and returns the EdgeVV type object from the Edge object
+  EdgeVV edge_to_edgevv(const Edge& arete) const {
+    return make_edgevv(arete.first->vertex(arete.second),
+                       arete.first->vertex(arete.third) );
+  }
+
+  FT compute_distance_to_facet_center(const Facet& f,
+                                      const Vertex_handle v) const {
+    const Point& fcenter = f.first->get_facet_surface_center(f.second);
+    const Point& vpoint = v->point();
+
+    return
+      r_tr_.geom_traits().compute_squared_distance_3_object()(fcenter.point(),
+                                                              vpoint.point())
+      - vpoint.weight();
+  }
+
+  Facet
+  biggest_incident_facet_in_complex(const Vertex_handle v) const
+  {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "Bad vertex: " << v->point() << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::vector<Facet> facets;
+    facets.reserve(64);
+    r_tr_.incident_facets(v, std::back_inserter(facets));
+
+    typename std::vector<Facet>::iterator fit = facets.begin();
+    while(fit != facets.end() && !r_c3t3_.is_in_complex(*fit)) ++fit;
+    CGAL_assertion(fit!=facets.end());
+    CGAL_assertion_code(int facet_counter = 1);
+
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "  " << v->cached_number_of_incident_facets()
+              << " incident faces, with sizes:\n";
+    std::cerr << "    " << compute_distance_to_facet_center(*fit, v) << "\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    Facet biggest_facet = *fit++;
+    while(fit != facets.end() && !r_c3t3_.is_in_complex(*fit)) ++fit;
+
+    for (; fit != facets.end(); )
+    {
+      CGAL_assertion_code(++facet_counter);
+      Facet current_facet = *fit;
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      std::cerr << "    " << compute_distance_to_facet_center(*fit, v) << "\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      // is the current facet bigger than the current biggest one
+      if ( compute_distance_to_facet_center(current_facet, v) >
+           compute_distance_to_facet_center(biggest_facet, v) )
+      {
+        biggest_facet = current_facet;
+      }
+      ++fit;
+      while(fit != facets.end() && !r_c3t3_.is_in_complex(*fit)) ++fit;
+    }
+    CGAL_assertion(v->cached_number_of_incident_facets() ==
+                   facet_counter);
+    CGAL_assertion(r_c3t3_.is_in_complex(biggest_facet));
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "Biggest facet radius: "
+              << compute_distance_to_facet_center(biggest_facet, v)
+              << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+
+    return biggest_facet;
+  }
+
+  Facet biggest_incident_facet_in_complex(const Edge& arete) const {
+    // Find the first facet in the incident facets
+    // of the edge which is in the Complex
+    // use the list of incident facets in the complex
+    Vertex_handle fev = edge_to_edgevv(arete).first;
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "Bad edge: (" << fev->point()
+              << ", " << arete.first->vertex(arete.third)->point()
+              << ")\n  incident facets sizes:\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    Tr_facet_circulator fcirc = r_tr_.incident_facets(arete);
+    while(!r_c3t3_.is_in_complex(*fcirc)) ++fcirc;
+    Facet first_facet = *fcirc;
+    Facet biggest_facet = *fcirc;
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "    "
+              << compute_distance_to_facet_center(*fcirc, fev) << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+
+    for (++fcirc; *fcirc != first_facet; ++fcirc) 
+    {
+      while(!r_c3t3_.is_in_complex(*fcirc)) ++fcirc;
+      if(*fcirc == first_facet) break;
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      std::cerr << "    "
+                << compute_distance_to_facet_center(*fcirc, fev) << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+
+      // is the current facet bigger than the current biggest one
+      if ( compute_distance_to_facet_center(*fcirc, fev) >
+           compute_distance_to_facet_center(biggest_facet,
+                                            fev) ) {
+        biggest_facet = *fcirc;
+      }
+    }
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "Biggest facet radius: "
+              << compute_distance_to_facet_center(biggest_facet, fev)
+              << std::endl;
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+
+    return biggest_facet;
+  }
+
+  ///////////////////////
+  // For before_insertion
+
+  // Actions to perform on a facet inside the conflict zone
+  void
+  before_insertion_handle_facet_inside_conflict_zone (const Facet& f) 
+  {
+    if ( r_c3t3_.is_in_complex(f) ) {
+      // foreach edge of f
+      const Cell_handle cell = f.first;
+      const int i = f.second;
+      for(int j = 0; j < 3; ++j)
+      {
+        const int edge_index_va = r_tr_.vertex_triple_index(i, j);
+        const int edge_index_vb = r_tr_.vertex_triple_index(i, (j == 2) ? 0 : (j+1));
+        const Vertex_handle edge_va = cell->vertex(edge_index_va);
+        const Vertex_handle edge_vb = cell->vertex(edge_index_vb);
+        m_bad_edges.left.erase( make_edgevv(edge_va, edge_vb));
+      }
+    }
+  }
+
+  // Action to perform on a facet on the boundary of the conflict zone
+  void
+  before_insertion_handle_facet_on_boundary_of_conflict_zone(const Facet& f)
+  {
+    // perform the same operations as for an internal facet
+    before_insertion_handle_facet_inside_conflict_zone (f);
+
+    if(m_bad_vertices_initialized) {
+      const Cell_handle& c = f.first;
+      const int i = f.second;
+
+      // for each v of f
+      for (int j = 0; j < 4; j++)
+        if (i != j)
+          if(m_bad_vertices.erase(c->vertex(j)) > 0) {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+            std::cerr << "m_bad_vertices.erase("
+                      << c->vertex(j)->point() << ")\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+          }
+    }
+  }
+
+public:
+  Refine_facets_manifold_base(Tr& triangulation,
+                              const Criteria& criteria,
+                              const Mesh_domain& oracle,
+                              C3t3& c3t3,
+                              bool /*with_boundary*/ = withBoundary)
+    : Base(triangulation,
+           criteria,
+           oracle,
+           c3t3)
+    , r_tr_(triangulation)
+    , r_c3t3_(c3t3)
+    , m_manifold_info_initialized(false)
+    , m_bad_vertices_initialized(false)
+    , m_with_manifold_criterion(get_with_manifold(criteria))
+    , m_with_boundary(get_with_boundary(criteria))
+  {
+#ifdef CGAL_MESH_3_DEBUG_CONSTRUCTORS
+    std::cerr << "CONS: Refine_facets_manifold_base";
+    if(m_with_manifold_criterion) {
+      if(m_with_boundary)
+        std::cerr << " (with boundaries)\n";
+      else
+        std::cerr << " (without boundary)\n";
+    } else {
+      std::cerr << " (DEACTIVATED)\n";
+    }
+#endif
+  }
+
+private:
+  // Initialization function
+  void initialize_manifold_info() const {
+#ifdef CGAL_MESH_3_VERBOSE
+    std::cerr << "\nscanning edges ";
+    if(m_with_boundary)
+      std::cerr << "(boundaries allowed)";
+    std::cerr << "...\n";
+    int n = 0;
+#endif
+    for (typename Tr::Finite_edges_iterator
+           eit = r_tr_.finite_edges_begin(), end = r_tr_.finite_edges_end();
+         eit != end; ++eit)
+    {
+      if ( (r_c3t3_.face_status(*eit) == C3t3::SINGULAR) ||
+           ( (!m_with_boundary) &&
+             (r_c3t3_.face_status(*eit) == C3t3::BOUNDARY) ) )
+      {
+        m_bad_edges.insert(Bad_edge(edge_to_edgevv(*eit),
+                                    (r_c3t3_.face_status(*eit) ==
+                                     C3t3::SINGULAR ? 0 : 1)));
+#ifdef CGAL_MESH_3_VERBOSE
+        ++n;
+#endif
+      }
+    }
+    m_manifold_info_initialized = true;
+#ifdef CGAL_MESH_3_VERBOSE
+    std::cerr << "   -> found " << n << " bad edges\n";
+#endif
+  }
+
+  void initialize_bad_vertices() const
+  {
+    CGAL_assertion(m_bad_vertices_initialized == false);
+    CGAL_assertion(m_bad_vertices.empty());
+#ifdef CGAL_MESH_3_VERBOSE
+    std::cerr << "\nscanning vertices..." << std::endl;
+    int n = 0;
+#endif
+    for (typename Tr::Finite_vertices_iterator
+           vit = r_tr_.finite_vertices_begin(),
+           end = r_tr_.finite_vertices_end();
+         vit != end; ++vit)
+    {
+      if( r_c3t3_.face_status(vit) == C3t3::SINGULAR ) {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        std::cerr << "m_bad_edges.insert("
+                  << vit->point() << ")\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        m_bad_vertices.insert( vit );
+#ifdef CGAL_MESH_3_VERBOSE
+        ++n;
+#endif
+      }
+    }
+    m_bad_vertices_initialized = true;
+#ifdef CGAL_MESH_3_VERBOSE
+    std::cerr << "   -> found " << n << " bad vertices\n";
+#  ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+    std::cerr << "Bad vertices queue:\n";
+    BOOST_FOREACH(Vertex_handle v2, m_bad_vertices)
+    {
+      std::cerr << v2->point() << std::endl;
+    }
+    CGAL::dump_c3t3(r_c3t3_, "dump-at-scan-vertices");
+#  endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+#endif
+  }
+
+public:
+  void scan_triangulation_impl() {
+    Base::scan_triangulation_impl();
+#ifdef CGAL_MESH_3_VERBOSE
+    std::cerr << "scanning edges (lazy)" << std::endl;
+    std::cerr << "scanning vertices (lazy)" << std::endl;
+#endif
+  }
+
+  // Tells whether there remain elements to refine
+  bool no_longer_element_to_refine_impl() const {
+    if(Base::no_longer_element_to_refine_impl())
+    {
+      if(!m_with_manifold_criterion) return true;
+      // Disable the manifold criterion
+
+      if( ! m_manifold_info_initialized ) initialize_manifold_info();
+
+      if(m_bad_edges.left.empty())
+      {
+        if( ! m_bad_vertices_initialized ) initialize_bad_vertices();
+
+        return m_bad_vertices.empty();
+      }
+      else // m_bad_vertices is not empty
+        return false;
+    }
+    else // Base::no_longer_element_to_refine_impl() returned false
+      return false;
+  }
+
+  // Returns the next element to refine
+  Facet get_next_element_impl() {
+
+    if (!Base::no_longer_element_to_refine_impl()) {
+      return Base::get_next_element_impl();
+    }
+    else {
+      if(!m_bad_edges.left.empty()) {
+        Edge first_bad_edge = edgevv_to_edge(m_bad_edges.right.begin()->second);
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        const EdgeVV& edgevv = m_bad_edges.right.begin()->second;
+        std::cerr << "Bad edge "
+                  << edgevv.first->point()
+                  << " - "
+                  << edgevv.second->point()
+                  << "\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        return biggest_incident_facet_in_complex(first_bad_edge);
+      } else {
+        CGAL_assertion(!m_bad_vertices.empty());
+        const Vertex_handle& v = *m_bad_vertices.begin();
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        std::cerr << "Bad vertices queue:\n";
+        BOOST_FOREACH(Vertex_handle v2, m_bad_vertices)
+        {
+          std::cerr << v2->point() << std::endl;
+        }
+        std::cerr << "Bad vertex " << v->point() << "\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        CGAL_assertion(r_c3t3_.has_incident_facets_in_complex(v));
+        if(r_c3t3_.face_status(v) != C3t3::SINGULAR) {
+          dump_c3t3(r_c3t3_, "dump-crash");
+          CGAL_error_msg("r_c3t3_.face_status(v) != C3t3::SINGULAR");
+        }
+        return biggest_incident_facet_in_complex(v);
+      }
+    }
+  }
+
+  void before_insertion_impl(const Facet& f, const Point& s,
+                             Zone& zone) {
+    if( m_manifold_info_initialized )
+    {
+      for (typename Zone::Facets_iterator fit =
+             zone.internal_facets.begin();
+           fit != zone.internal_facets.end();
+           ++fit)
+        before_insertion_handle_facet_inside_conflict_zone (*fit);
+
+      for (typename Zone::Facets_iterator fit =
+             zone.boundary_facets.begin(); fit !=
+             zone.boundary_facets.end(); ++fit)
+        before_insertion_handle_facet_on_boundary_of_conflict_zone (*fit);
+    }
+    Base::before_insertion_impl(f, s, zone);
+  }
+
+  void after_insertion_impl(const Vertex_handle v) {
+    Base::after_insertion_impl(v);
+
+    if( ! m_manifold_info_initialized )
+      return;
+
+    // search for incident facets around v
+    typedef std::vector<Facet> Facets;
+    Facets facets;
+    facets.reserve(64);
+    r_tr_.incident_facets (v, std::back_inserter(facets));
+
+    // foreach f in star(v)
+    for (typename Facets::iterator fit = facets.begin();
+         fit != facets.end();
+         ++fit) 
+    {
+      // foreach edge of *fit
+      const Cell_handle cell = fit->first;
+      const int i = fit->second;
+      for(int j = 0; j < 3; ++j)
+      {
+        const int edge_index_va = r_tr_.vertex_triple_index(i, j);
+        const int edge_index_vb = r_tr_.vertex_triple_index(i, (j == 2) ? 0 : (j+1));
+        Edge edge(cell, edge_index_va, edge_index_vb);
+        // test if edge is in Complex
+        if ( r_c3t3_.face_status(edge) != C3t3::NOT_IN_COMPLEX ) {
+          // test if edge is not regular to store it as a "bad_edge"
+          // e.g. more than or equal to 3 incident facets (SINGULAR)
+          // or less than or equal to 1
+          // (BOUNDARY only, because ISOLATED is NA)
+          // This test is not efficient because
+          // edges are tried to be inserted several times
+          // TODO one day: test if the edge is still singular
+          if ( (r_c3t3_.face_status(edge) == C3t3::SINGULAR) ||
+               ( (!m_with_boundary) && 
+                 (r_c3t3_.face_status(edge) == C3t3::BOUNDARY) )
+               )
+          {
+            m_bad_edges.insert(Bad_edge(edge_to_edgevv(edge),
+                                        (r_c3t3_.face_status(edge) ==
+                                         C3t3::SINGULAR ? 0 : 1)));
+          }
+          else {
+            m_bad_edges.left.erase( edge_to_edgevv(edge) ); // @TODO: pourquoi?!
+          }
+        }
+      }
+    }
+
+    if(!m_bad_vertices_initialized) return;
+
+    // foreach v' in star of v
+    std::vector<Vertex_handle> vertices;
+    vertices.reserve(64);
+    r_tr_.incident_vertices(v, std::back_inserter(vertices));
+
+    // is_regular_or_boundary_for_vertices
+    // is used here also incident edges are not known to be
+    // REGULAR which may cause some singular vertices to be forgotten
+    // This causes no problem because
+    // those SINGULAR incident SINGULAR edges are going to be handled
+    for (typename std::vector<Vertex_handle>::iterator
+           vit = vertices.begin(), end = vertices.end();
+         vit != end; ++vit)
+    {
+      if ( r_c3t3_.has_incident_facets_in_complex(*vit)  &&
+           (r_c3t3_.face_status(*vit) == C3t3::SINGULAR)
+           // !r_c3t3_.is_regular_or_boundary_for_vertices(*vit)
+           )
+      {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        std::cerr << "m_bad_vertices.insert("
+                  << (*vit)->point() << ")\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+        m_bad_vertices.insert(*vit);
+      }
+    }
+
+    if ( r_c3t3_.has_incident_facets_in_complex(v) &&
+         (r_c3t3_.face_status(v) == C3t3::SINGULAR)
+         // !r_c3t3_.is_regular_or_boundary_for_vertices(v)
+         )
+    {
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      std::cerr << "m_bad_vertices.insert("
+                << v->point() << ")\n";
+#endif // CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+      m_bad_vertices.insert(v);
+    }
+  }
+
+  std::string debug_info() const
+  {
+    std::stringstream s;
+    s << Base::debug_info();
+    if(m_with_manifold_criterion) {
+      s << "," << m_bad_edges.left.size()
+        << "," << m_bad_vertices.size();
+    }
+    return s.str();
+  }
+
+  std::string debug_info_header() const
+  {
+    std::stringstream s;
+    s << Base::debug_info_header();
+    if(m_with_manifold_criterion) {
+      s << ",#bad edges,#bad vertices";
+    }
+    return s.str();
+  }
+};  // end Refine_facets_manifold_base
+
+}  // end namespace Mesh_3
+
+}  // end namespace CGAL
+
+
+#endif // CGAL_MESH_3_REFINE_FACETS_MANIFOLD_BASE_H
diff --git a/Mesh_3/include/CGAL/Mesh_facet_criteria_3.h b/Mesh_3/include/CGAL/Mesh_facet_criteria_3.h
index 69057089f10..7f8d2112dd8 100644
--- a/Mesh_3/include/CGAL/Mesh_facet_criteria_3.h
+++ b/Mesh_3/include/CGAL/Mesh_facet_criteria_3.h
@@ -50,6 +50,8 @@ private:
   typedef Mesh_facet_criteria_3<Tr> Self;
 
 public:
+  typedef CGAL::Tag_true Has_manifold_criterion;
+
   /**
    * @brief Constructor
    */
@@ -110,6 +112,10 @@ public:
     criteria_.add(criterion);
   }
 
+  Mesh_facet_topology topology() const {
+    return topology_;
+  }
+
 private:
   void init_aspect(const FT& angle_bound)
   {
@@ -138,7 +144,8 @@ private:
   
   void init_topo(const Mesh_facet_topology topology)
   {
-    switch ( topology )
+    topology_ = topology;
+    switch ( topology % MANIFOLD )
     {
       case FACET_VERTICES_ON_SURFACE:
       {
@@ -166,6 +173,7 @@ private:
   
 private:
   Criteria criteria_;
+  Mesh_facet_topology topology_;
 };  // end class Mesh_facet_criteria_3
 
 }  // end namespace CGAL
diff --git a/Mesh_3/include/CGAL/Mesh_facet_topology.h b/Mesh_3/include/CGAL/Mesh_facet_topology.h
index d15e9c54e8e..55718035c3c 100644
--- a/Mesh_3/include/CGAL/Mesh_facet_topology.h
+++ b/Mesh_3/include/CGAL/Mesh_facet_topology.h
@@ -29,9 +29,12 @@ namespace CGAL {
 
 enum Mesh_facet_topology {
   FACET_VERTICES_ON_SURFACE = 1,
-  FACET_VERTICES_ON_SAME_SURFACE_PATCH,
-  FACET_VERTICES_ON_SAME_SURFACE_PATCH_WITH_ADJACENCY_CHECK
-};  
+  FACET_VERTICES_ON_SAME_SURFACE_PATCH = 2,
+  FACET_VERTICES_ON_SAME_SURFACE_PATCH_WITH_ADJACENCY_CHECK = 3,
+  MANIFOLD_WITH_BOUNDARY = 8,
+  NO_BOUNDARY = 16,
+  MANIFOLD = 24
+};
   
 } // end namespace CGAL
 
diff --git a/Mesh_3/include/CGAL/Mesh_vertex_base_3.h b/Mesh_3/include/CGAL/Mesh_vertex_base_3.h
index 410fd65f4d0..0f0199d4835 100644
--- a/Mesh_3/include/CGAL/Mesh_vertex_base_3.h
+++ b/Mesh_3/include/CGAL/Mesh_vertex_base_3.h
@@ -212,12 +212,12 @@ public:
     cache_validity = true;
   }
 
-  int cached_number_of_incident_facets() const
+  std::size_t cached_number_of_incident_facets() const
   {
     return number_of_incident_facets_;
   }
     
-  int cached_number_of_components() const
+  std::size_t cached_number_of_components() const
   {
     return number_of_components_;
   }
@@ -232,8 +232,8 @@ public:
   }
 private:
 
-  int number_of_incident_facets_;
-  int number_of_components_; // number of components in the adjacency
+  std::size_t number_of_incident_facets_;
+  std::size_t number_of_components_; // number of components in the adjacency
   // graph of incident facets (in complex)
 
 
diff --git a/Surface_mesher/include/CGAL/Surface_mesher/Surface_mesher_manifold.h b/Surface_mesher/include/CGAL/Surface_mesher/Surface_mesher_manifold.h
index fbc4e390487..1d7b8b064cb 100644
--- a/Surface_mesher/include/CGAL/Surface_mesher/Surface_mesher_manifold.h
+++ b/Surface_mesher/include/CGAL/Surface_mesher/Surface_mesher_manifold.h
@@ -42,19 +42,12 @@ namespace CGAL {
       typedef C2T3 C2t3;
       typedef typename C2T3::Triangulation Tr;
       typedef typename Tr::Geom_traits GT;
-      typedef typename GT::FT FT;
       typedef typename Tr::Point Point;
       typedef typename Tr::Facet Facet;
-      typedef typename Tr::Cell_handle Cell_handle;
       typedef typename Tr::Vertex_handle Vertex_handle;
       typedef typename Triangulation_mesher_level_traits_3<Tr>::Zone Zone;
-      typedef std::list<Facet> Facets;
-      typedef std::list<Vertex_handle> Vertices;
-      typedef typename Facets::iterator Facets_iterator;
-      typedef typename Vertices::iterator Vertices_iterator;
-      typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator;
 
-  protected:
+  private:
     // because of the caching, these two members have to be mutable,
     // because they are actually updated in the const method
     // 'no_longer_element_to_refine_impl()'
@@ -62,7 +55,14 @@ namespace CGAL {
     mutable Bad_vertices bad_vertices;
     mutable bool bad_vertices_initialized;
 
-    private:
+  private:
+      typedef typename Tr::Cell_handle Cell_handle;
+      typedef std::list<Facet> Facets;
+      typedef std::list<Vertex_handle> Vertices;
+      typedef typename Facets::iterator Facets_iterator;
+      typedef typename Vertices::iterator Vertices_iterator;
+      typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator;
+
       Facet canonical_facet(const Facet& f) const {
 	const Cell_handle& c = f.first;
 	const Cell_handle& c2 = c->neighbor(f.second);
@@ -107,6 +107,27 @@ namespace CGAL {
 	return biggest_facet;
       }
 
+      void initialize_bad_vertices() const
+      {
+#ifdef CGAL_SURFACE_MESHER_VERBOSE
+	std::cerr << "scanning vertices" << std::endl;
+#endif
+	int n = 0;
+	for (Finite_vertices_iterator vit = SMMBB::tr.finite_vertices_begin();
+	     vit != SMMBB::tr.finite_vertices_end();
+	     ++vit) {
+	  if ( (SMMBB::c2t3.face_status(vit)  // @TODO: appeler is_regular
+		== C2t3::SINGULAR) ) {
+	    bad_vertices.insert( vit );
+	    ++n;
+	  }
+	}
+	bad_vertices_initialized = true;
+#ifdef CGAL_SURFACE_MESHER_VERBOSE
+	std::cerr << "   -> found " << n << " bad vertices\n";
+#endif
+      }
+
     public:
       Surface_mesher_manifold_base (C2T3& c2t3,
                                     const Surface& surface,
@@ -133,27 +154,6 @@ namespace CGAL {
 	return false;
       }
 
-      void initialize_bad_vertices() const
-      {
-#ifdef CGAL_SURFACE_MESHER_VERBOSE
-	std::cerr << "scanning vertices" << std::endl;
-#endif
-	int n = 0;
-	for (Finite_vertices_iterator vit = SMMBB::tr.finite_vertices_begin();
-	     vit != SMMBB::tr.finite_vertices_end();
-	     ++vit) {
-	  if ( (SMMBB::c2t3.face_status(vit)  // @TODO: appeler is_regular
-		== C2t3::SINGULAR) ) {
-	    bad_vertices.insert( vit );
-	    ++n;
-	  }
-	}
-	bad_vertices_initialized = true;
-#ifdef CGAL_SURFACE_MESHER_VERBOSE
-	std::cerr << "   -> found " << n << " bad vertices\n";
-#endif
-      }
-
       // Lazy initialization function
       void scan_triangulation_impl() {
 	SMMBB::scan_triangulation_impl();