diff --git a/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h b/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h
index 42b0de9a66f..13fd9f70083 100644
--- a/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h
+++ b/Mesh_3/include/CGAL/Compact_mesh_cell_base_3.h
@@ -98,12 +98,9 @@ namespace CGAL {
 // Base for Compact_mesh_cell_base_3, with specializations
 // for different values of Concurrency_tag
 // Sequential
-template <typename GT, typename Concurrency_tag>
+template <typename Point_3, typename Concurrency_tag>
 class Compact_mesh_cell_base_3_base
 {
-  typedef typename GT::Point_3                Point_3;
-  typedef typename GT::Weighted_point_3       Point;
-
 protected:
   Compact_mesh_cell_base_3_base()
     : bits_(0)
@@ -176,11 +173,9 @@ protected:
 #ifdef CGAL_LINKED_WITH_TBB
 // Class Compact_mesh_cell_base_3_base
 // Specialization for parallel
-template <typename GT>
-class Compact_mesh_cell_base_3_base<GT, Parallel_tag>
+template <typename Point_3>
+class Compact_mesh_cell_base_3_base<Point_3, Parallel_tag>
 {
-  typedef typename GT::Point_3  Point_3;
-
 protected:
   Compact_mesh_cell_base_3_base()
   {
@@ -260,14 +255,16 @@ protected:
 // Class Compact_mesh_cell_base_3
 // Cell base class used in 3D meshing process.
 // Adds information to Cb about the cell of the input complex containing it
-template< class GT,
-          class MD,
-          class TDS = void >
-class Compact_mesh_cell_base_3
-  : public Compact_mesh_cell_base_3_base<GT, typename TDS::Concurrency_tag>
+template< class Point_3,
+          class Weighted_point_3,
+          class Subdomain_index_,
+          class Surface_patch_index_,
+          class Index_,
+          class TDS>
+class Compact_mesh_cell_3
+  : public Compact_mesh_cell_base_3_base<Point_3, typename TDS::Concurrency_tag>
 {
-  typedef typename GT::FT FT;
-  typedef Compact_mesh_cell_base_3_base<GT,typename TDS::Concurrency_tag> Base;
+  typedef Compact_mesh_cell_base_3_base<Point_3,typename TDS::Concurrency_tag> Base;
   using Base::weighted_circumcenter_;
 
 public:
@@ -279,20 +276,12 @@ public:
   typedef typename TDS::Cell_data      TDS_data;
 
 
-  template <typename TDS2>
-  struct Rebind_TDS {
-    typedef Compact_mesh_cell_base_3<GT, MD, TDS2> Other;
-  };
-
-
   // Index Type
-  typedef typename MD::Subdomain_index      Subdomain_index;
-  typedef typename MD::Surface_patch_index  Surface_patch_index;
-  typedef typename MD::Index                Index;
+  typedef Subdomain_index_      Subdomain_index;
+  typedef Surface_patch_index_  Surface_patch_index;
+  typedef Index_                Index;
 
-  typedef GT                                Geom_traits;
-  typedef typename GT::Point_3              Point_3;
-  typedef typename GT::Weighted_point_3     Point;
+  typedef Weighted_point_3      Point;
 
 
   typedef Point*                            Point_container;
@@ -310,20 +299,10 @@ public:
 
 public:
   // Constructors
-  Compact_mesh_cell_base_3()
-    : surface_index_table_()
-    , surface_center_table_()
-#ifdef CGAL_INTRUSIVE_LIST
-    , next_intrusive_()
-    , previous_intrusive_()
-#endif
-    , surface_center_index_table_()
-    , sliver_value_(FT(0.))
-    , subdomain_index_()
-    , sliver_cache_validity_(false)
+  Compact_mesh_cell_3()
   {}
 
-  Compact_mesh_cell_base_3(const Compact_mesh_cell_base_3& rhs)
+  Compact_mesh_cell_3(const Compact_mesh_cell_3& rhs)
     : N(rhs.N)
     , V(rhs.V)
 #ifdef CGAL_INTRUSIVE_LIST
@@ -342,49 +321,29 @@ public:
     }
   }
 
-  Compact_mesh_cell_base_3 (Vertex_handle v0,
-                            Vertex_handle v1,
-                            Vertex_handle v2,
-                            Vertex_handle v3)
-    : surface_index_table_()
-    , surface_center_table_()
-    , V(CGAL::make_array(v0, v1, v2, v3))
-#ifdef CGAL_INTRUSIVE_LIST
-    , next_intrusive_()
-    , previous_intrusive_()
-#endif
-    , surface_center_index_table_()
-    , sliver_value_(FT(0.))
-    , subdomain_index_()
-    , sliver_cache_validity_(false)
+  Compact_mesh_cell_3 (Vertex_handle v0,
+                       Vertex_handle v1,
+                       Vertex_handle v2,
+                       Vertex_handle v3)
+    : V(CGAL::make_array(v0, v1, v2, v3))
   {
   }
 
 
-  Compact_mesh_cell_base_3 (Vertex_handle v0,
-                            Vertex_handle v1,
-                            Vertex_handle v2,
-                            Vertex_handle v3,
-                            Cell_handle n0,
-                            Cell_handle n1,
-                            Cell_handle n2,
-                            Cell_handle n3)
-    : surface_index_table_()
-    , surface_center_table_()
-    , N(CGAL::make_array(n0, n1, n2, n3))
+  Compact_mesh_cell_3 (Vertex_handle v0,
+                       Vertex_handle v1,
+                       Vertex_handle v2,
+                       Vertex_handle v3,
+                       Cell_handle n0,
+                       Cell_handle n1,
+                       Cell_handle n2,
+                       Cell_handle n3)
+    : N(CGAL::make_array(n0, n1, n2, n3))
     , V(CGAL::make_array(v0, v1, v2, v3))
-#ifdef CGAL_INTRUSIVE_LIST
-    , next_intrusive_()
-    , previous_intrusive_()
-#endif
-    , surface_center_index_table_()
-    , sliver_value_(FT(0.))
-    , subdomain_index_()
-    , sliver_cache_validity_(false)
   {
   }
 
-  ~Compact_mesh_cell_base_3()
+  ~Compact_mesh_cell_3()
   {
     if(!internal_tbb::is_null(weighted_circumcenter_)){
       internal_tbb::delete_circumcenter(weighted_circumcenter_);
@@ -550,11 +509,6 @@ public:
     return *weighted_circumcenter_;
   }
 
-  const Point_3& weighted_circumcenter() const
-  {
-    return weighted_circumcenter(Geom_traits());
-  }
-
   // Returns the index of the cell of the input complex that contains the cell
   Subdomain_index subdomain_index() const { return subdomain_index_; }
 
@@ -562,13 +516,13 @@ public:
   void set_subdomain_index(const Subdomain_index& index)
   { subdomain_index_ = index; }
 
-  void set_sliver_value(const FT& value)
+  void set_sliver_value(double value)
   {
     sliver_cache_validity_ = true;
     sliver_value_ = value;
   }
 
-  const FT& sliver_value() const
+  double sliver_value() const
   {
     CGAL_assertion(is_cache_valid());
     return sliver_value_;
@@ -646,6 +600,7 @@ public:
   static
   std::string io_signature()
   {
+    using Geom_traits = typename Kernel_traits<Point>::type;
     return
       Get_io_signature<Subdomain_index>()() + "+" +
       Get_io_signature<Regular_triangulation_cell_base_3<Geom_traits> >()()
@@ -683,81 +638,81 @@ private:
 
 
   /// Stores surface_index for each facet of the cell
-  std::array<Surface_patch_index, 4> surface_index_table_;
+  std::array<Surface_patch_index, 4> surface_index_table_ = {};
   /// Stores surface center of each facet of the cell
-  std::array<Point_3, 4> surface_center_table_;
+  std::array<Point_3, 4> surface_center_table_ = {};
   /// Stores surface center index of each facet of the cell
 
   std::array<Cell_handle, 4> N;
   std::array<Vertex_handle, 4> V;
 
 #ifdef CGAL_INTRUSIVE_LIST
-  Cell_handle next_intrusive_, previous_intrusive_;
+  Cell_handle next_intrusive_ = {}, previous_intrusive_ = {};
 #endif
   std::size_t time_stamp_;
 
-  std::array<Index, 4> surface_center_index_table_;
+  std::array<Index, 4> surface_center_index_table_ = {};
   /// Stores visited facets (4 first bits)
 
   //  Point_container _hidden;
 
-  FT sliver_value_;
+  double sliver_value_ = 0.;
 
   // The index of the cell of the input complex that contains me
-  Subdomain_index subdomain_index_;
+  Subdomain_index subdomain_index_ = {};
 
   TDS_data      _tds_data;
-  mutable bool sliver_cache_validity_;
+  mutable bool sliver_cache_validity_ = false;
 
+public:
 
-};  // end class Compact_mesh_cell_base_3
-
-template < class GT, class MT, class Cb >
-std::istream&
-operator>>(std::istream &is,
-           Compact_mesh_cell_base_3<GT, MT, Cb> &c)
-{
-  typename Compact_mesh_cell_base_3<GT, MT, Cb>::Subdomain_index index;
-  if(is_ascii(is))
-    is >> index;
-  else
-    read(is, index);
-  if(is) {
-    c.set_subdomain_index(index);
-    for(int i = 0; i < 4; ++i)
-    {
-      typename Compact_mesh_cell_base_3<GT, MT, Cb>::Surface_patch_index i2;
-      if(is_ascii(is))
-        is >> iformat(i2);
-      else
-      {
-        read(is, i2);
-      }
-      c.set_surface_patch_index(i, i2);
+  friend std::istream& operator>>(std::istream &is, Compact_mesh_cell_3 &c)
+  {
+    Subdomain_index index;
+    if(is_ascii(is))
+      is >> index;
+    else
+      read(is, index);
+    if(is) {
+      c.set_subdomain_index(index);
+      for(int i = 0; i < 4; ++i)
+        {
+          Surface_patch_index i2;
+          if(is_ascii(is))
+            is >> iformat(i2);
+          else
+            {
+              read(is, i2);
+            }
+          c.set_surface_patch_index(i, i2);
+        }
     }
+    return is;
   }
-  return is;
-}
 
-template < class GT, class MT, class Cb >
-std::ostream&
-operator<<(std::ostream &os,
-           const Compact_mesh_cell_base_3<GT, MT, Cb> &c)
-{
-  if(is_ascii(os))
-     os << c.subdomain_index();
-  else
-    write(os, c.subdomain_index());
-  for(int i = 0; i < 4; ++i)
+  friend
+  std::ostream& operator<<(std::ostream &os, const Compact_mesh_cell_3 &c)
   {
     if(is_ascii(os))
-      os << ' ' << oformat(c.surface_patch_index(i));
+      os << c.subdomain_index();
     else
-      write(os, c.surface_patch_index(i));
+      write(os, c.subdomain_index());
+    for(int i = 0; i < 4; ++i)
+      {
+        if(is_ascii(os))
+          os << ' ' << oformat(c.surface_patch_index(i));
+        else
+          write(os, c.surface_patch_index(i));
+      }
+    return os;
   }
-  return os;
-}
 
+};  // end class Compact_mesh_cell_3
+
+template< class GT,
+          class MD,
+          class TDS = void >
+class Compact_mesh_cell_base_3;
 
 // Specialization for void.
 template <typename GT, typename MD>
@@ -768,7 +723,35 @@ public:
   typedef Triangulation_data_structure::Vertex_handle   Vertex_handle;
   typedef Triangulation_data_structure::Cell_handle     Cell_handle;
   template <typename TDS2>
-  struct Rebind_TDS { typedef Compact_mesh_cell_base_3<GT, MD, TDS2> Other; };
+  struct Rebind_TDS {
+    typedef Compact_mesh_cell_3<typename GT::Point_3,
+                                typename GT::Weighted_point_3,
+                                typename MD::Subdomain_index,
+                                typename MD::Surface_patch_index,
+                                typename MD::Index,
+                                TDS2> Other;
+  };
+};
+
+template <typename GT,
+  typename Subdomain_index,
+  typename Surface_patch_index,
+  typename Index>
+class Compact_mesh_cell_generator_3
+{
+public:
+  typedef internal::Dummy_tds_3                         Triangulation_data_structure;
+  typedef Triangulation_data_structure::Vertex_handle   Vertex_handle;
+  typedef Triangulation_data_structure::Cell_handle     Cell_handle;
+  template <typename TDS2>
+  struct Rebind_TDS {
+    typedef Compact_mesh_cell_3<typename GT::Point_3,
+                                typename GT::Weighted_point_3,
+                                Subdomain_index,
+                                Surface_patch_index,
+                                Index,
+                                TDS2> Other;
+  };
 };
 
 }  // end namespace CGAL
diff --git a/Mesh_3/include/CGAL/Mesh_3/io_signature.h b/Mesh_3/include/CGAL/Mesh_3/io_signature.h
index d871b8ffe7d..2c0f02ab3b0 100644
--- a/Mesh_3/include/CGAL/Mesh_3/io_signature.h
+++ b/Mesh_3/include/CGAL/Mesh_3/io_signature.h
@@ -219,12 +219,10 @@ struct Get_io_signature<Weighted_point_3<K> >
 };
 
 #ifdef CGAL_TRIANGULATION_3_H
-template <class Gt, class Vb, class Cb, class C_tag>
+template <class Gt, class Tds>
 struct
-Get_io_signature<Triangulation_3<Gt, Triangulation_data_structure_3<Vb, Cb, C_tag> > >
+Get_io_signature<Triangulation_3<Gt, Tds > >
 {
-  typedef Triangulation_data_structure_3<Vb, Cb, C_tag> Tds;
-
   std::string operator()() {
     return std::string("Triangulation_3(") +
       Get_io_signature<typename Tds::Vertex::Point>()() +
diff --git a/Mesh_3/include/CGAL/Mesh_triangulation_3.h b/Mesh_3/include/CGAL/Mesh_triangulation_3.h
index 8d5373b1bd7..7bc643f4444 100644
--- a/Mesh_3/include/CGAL/Mesh_triangulation_3.h
+++ b/Mesh_3/include/CGAL/Mesh_triangulation_3.h
@@ -120,65 +120,40 @@ public:
 //
 template<class MD,
          class K_ = Default,
-         class Concurrency_tag = Sequential_tag,
+         class Concurrency_tag_ = Sequential_tag,
          class Vertex_base_ = Default,
          class Cell_base_   = Default>
-struct Mesh_triangulation_3;
-
-// Sequential version (default)
-template<class MD, class K_, class Concurrency_tag,
-         class Vertex_base_, class Cell_base_>
 struct Mesh_triangulation_3
 {
 private:
-  typedef typename Default::Lazy_get<K_, Kernel_traits<MD> >::type K;
+  using K = typename Default::Lazy_get<K_, Kernel_traits<MD> >::type;
 
-  typedef typename details::Mesh_geom_traits_generator<K>::type    Geom_traits;
+  using Geom_traits = typename details::Mesh_geom_traits_generator<K>::type;
 
-  typedef typename Default::Get<
+  using Indices_tuple = Mesh_3::internal::Indices_tuple_t<MD>;
+  using Vertex_base = typename Default::Get<
     Vertex_base_,
-    Mesh_vertex_base_3<Geom_traits, MD> >::type                    Vertex_base;
-  typedef typename Default::Get<
+    Mesh_vertex_generator_3<Geom_traits,
+                            Indices_tuple,
+                            typename MD::Index> >::type;
+  using Cell_base = typename Default::Get<
     Cell_base_,
-    Compact_mesh_cell_base_3<Geom_traits, MD> >::type              Cell_base;
-
-  typedef Triangulation_data_structure_3<Vertex_base,Cell_base>    Tds;
-  typedef Mesh_3_regular_triangulation_3_wrapper<Geom_traits, Tds> Triangulation;
+    Compact_mesh_cell_generator_3<Geom_traits,
+                                  typename MD::Subdomain_index,
+                                  typename MD::Surface_patch_index,
+                                  typename MD::Index> >::type;
+  using Concurrency_tag =
+      typename Default::Get<Concurrency_tag_, Sequential_tag>::type;
+  struct Tds : public Triangulation_data_structure_3<Vertex_base, Cell_base,
+                                                     Concurrency_tag> {};
+  using Triangulation =
+      Mesh_3_regular_triangulation_3_wrapper<Geom_traits, Tds>;
+  ;
 
 public:
-  typedef Triangulation type;
-  typedef type Type;
+  using type = Triangulation;
+  using Type = type;
 };  // end struct Mesh_triangulation_3
-
-#ifdef CGAL_LINKED_WITH_TBB
-// Parallel version (specialization)
-//
-template<class MD, class K_,
-         class Vertex_base_, class Cell_base_>
-struct Mesh_triangulation_3<MD, K_, Parallel_tag, Vertex_base_, Cell_base_>
-{
-private:
-  typedef typename Default::Get<K_, typename Kernel_traits<MD>::Kernel>::type K;
-
-  typedef typename details::Mesh_geom_traits_generator<K>::type       Geom_traits;
-
-  typedef typename Default::Get<
-    Vertex_base_,
-    Mesh_vertex_base_3<Geom_traits, MD> >::type                       Vertex_base;
-  typedef typename Default::Get<
-    Cell_base_,
-    Compact_mesh_cell_base_3<Geom_traits, MD> >::type                 Cell_base;
-
-  typedef Triangulation_data_structure_3<
-    Vertex_base, Cell_base, Parallel_tag>                             Tds;
-  typedef Mesh_3_regular_triangulation_3_wrapper<Geom_traits, Tds>    Triangulation;
-
-public:
-  typedef Triangulation type;
-  typedef type Type;
-};  // end struct Mesh_triangulation_3
-#endif // CGAL_LINKED_WITH_TBB
-
 }  // end namespace CGAL
 
 #include <CGAL/enable_warnings.h>
diff --git a/Mesh_3/include/CGAL/Mesh_vertex_base_3.h b/Mesh_3/include/CGAL/Mesh_vertex_base_3.h
index aeb3411a1e1..ecb1f9464c5 100644
--- a/Mesh_3/include/CGAL/Mesh_vertex_base_3.h
+++ b/Mesh_3/include/CGAL/Mesh_vertex_base_3.h
@@ -24,7 +24,7 @@
 
 
 #include <CGAL/Regular_triangulation_vertex_base_3.h>
-#include <CGAL/internal/Mesh_3/get_index.h>
+#include <CGAL/internal/Mesh_3/indices_management.h>
 #include <CGAL/Mesh_3/io_signature.h>
 #include <CGAL/Has_timestamp.h>
 #include <CGAL/tags.h>
@@ -92,9 +92,10 @@ protected:
 // Adds information to Vb about the localization of the vertex in regards
 // to the 3D input complex.
 template<class GT,
-         class MD,
-         class Vb = Regular_triangulation_vertex_base_3<GT> >
-class Mesh_vertex_base_3
+         class Indices_tuple,
+         class Index_,
+         class Vb>
+class Mesh_vertex_3
 : public Vb,
   public Mesh_vertex_base_3_base<
     typename Vb::Triangulation_data_structure::Concurrency_tag>
@@ -103,19 +104,12 @@ public:
   typedef Vb Cmvb3_base;
   typedef typename Vb::Vertex_handle  Vertex_handle;
 
-  // To get correct vertex type in TDS
-  template < class TDS3 >
-  struct Rebind_TDS {
-    typedef typename Vb::template Rebind_TDS<TDS3>::Other Vb3;
-    typedef Mesh_vertex_base_3 <GT, MD, Vb3> Other;
-  };
-
   // Types
-  typedef typename MD::Index                      Index;
-  typedef typename GT::FT                         FT;
+  typedef Index_                      Index;
+  typedef typename GT::FT             FT;
 
   // Constructor
-  Mesh_vertex_base_3()
+  Mesh_vertex_3()
     : Vb()
     , number_of_incident_facets_(0)
     , number_of_components_(0)
@@ -248,56 +242,77 @@ private:
   Vertex_handle previous_intrusive_;
 #endif
   std::size_t time_stamp_;
+public:
 
-};  // end class Mesh_vertex_base_3
+  friend std::istream& operator>>(std::istream &is, Mesh_vertex_3& v)
+  {
+    is >> static_cast<Cmvb3_base&>(v);
+    int dimension;
+    if(is_ascii(is)) {
+      is >> dimension;
+
+    } else {
+      CGAL::read(is, dimension);
+    }
+    v.set_dimension(dimension);
+    CGAL_assertion(v.in_dimension() >= -1);
+    CGAL_assertion(v.in_dimension() < 4);
+    Index index =
+      Mesh_3::internal::Read_write_index<Indices_tuple,
+                                         Index>()(is, v.in_dimension());
+    v.set_index(index);
+    return is;
+  }
+
+  friend std::ostream& operator<<(std::ostream &os, const Mesh_vertex_3& v)
+  {
+    os << static_cast<const Cmvb3_base&>(v);
+    if(is_ascii(os)) {
+      os << " " << v.in_dimension()
+         << " ";
+    } else {
+      CGAL::write(os, v.in_dimension());
+    }
+    Mesh_3::internal::Read_write_index<Indices_tuple,
+                                       Index>()(os,
+                                                v.in_dimension(),
+                                                v.index());
+    return os;
+  }
+};  // end class Mesh_vertex_3
 
 template<class GT,
          class MD,
-         class Vb>
-inline
-std::istream&
-operator>>(std::istream &is, Mesh_vertex_base_3<GT,MD,Vb>& v)
-{
-  typedef Mesh_vertex_base_3<GT,MD,Vb> Vertex;
-  typedef typename Vertex::Cmvb3_base Cmvb3_base;
-  is >> static_cast<Cmvb3_base&>(v);
-  int dimension;
-  if(is_ascii(is)) {
-    is >> dimension;
+         class Vb = Regular_triangulation_vertex_base_3<GT> >
+struct Mesh_vertex_base_3 {
+  using Triangulation_data_structure = internal::Dummy_tds_3;
+  using Vertex_handle = typename Triangulation_data_structure::Vertex_handle;
+  using Cell_handle = typename Triangulation_data_structure::Cell_handle;
 
-  } else {
-    CGAL::read(is, dimension);
-  }
-  v.set_dimension(dimension);
-  CGAL_assertion(v.in_dimension() >= -1);
-  CGAL_assertion(v.in_dimension() < 4);
-  typename Vertex::Index index =
-    Mesh_3::internal::Read_mesh_domain_index<MD>()(v.in_dimension(), is);
-  v.set_index(index);
-  return is;
-}
+  template < class TDS3 >
+  struct Rebind_TDS {
+    using Vb3 = typename Vb::template Rebind_TDS<TDS3>::Other;
+    using Other = Mesh_vertex_3 <GT,
+                                 Mesh_3::internal::Indices_tuple_t<MD>,
+                                 typename MD::Index, Vb3>;
+  };
+};
 
 template<class GT,
-         class MD,
-         class Vb>
-inline
-std::ostream&
-operator<<(std::ostream &os, const Mesh_vertex_base_3<GT,MD,Vb>& v)
-{
-  typedef Mesh_vertex_base_3<GT,MD,Vb> Vertex;
-  typedef typename Vertex::Cmvb3_base Cmvb3_base;
-  os << static_cast<const Cmvb3_base&>(v);
-  if(is_ascii(os)) {
-    os << " " << v.in_dimension()
-       << " ";
-  } else {
-    CGAL::write(os, v.in_dimension());
-  }
-  Mesh_3::internal::Write_mesh_domain_index<MD>()(os,
-                                                  v.in_dimension(),
-                                                  v.index());
-  return os;
-}
+         class Indices_tuple,
+         class Index,
+         class Vb = Regular_triangulation_vertex_base_3<GT> >
+struct Mesh_vertex_generator_3 {
+  using Triangulation_data_structure = internal::Dummy_tds_3;
+  using Vertex_handle = typename Triangulation_data_structure::Vertex_handle;
+  using Cell_handle = typename Triangulation_data_structure::Cell_handle;
+
+  template < class TDS3 >
+  struct Rebind_TDS {
+    typedef typename Vb::template Rebind_TDS<TDS3>::Other Vb3;
+    typedef Mesh_vertex_3 <GT, Indices_tuple, Index, Vb3> Other;
+  };
+};
 
 }  // end namespace CGAL
 
diff --git a/Mesh_3/include/CGAL/internal/Mesh_3/indices_management.h b/Mesh_3/include/CGAL/internal/Mesh_3/indices_management.h
index 62bde6e2135..fdb84a6c82b 100644
--- a/Mesh_3/include/CGAL/internal/Mesh_3/indices_management.h
+++ b/Mesh_3/include/CGAL/internal/Mesh_3/indices_management.h
@@ -23,10 +23,14 @@
 
 #include <boost/type_traits/is_same.hpp>
 #include <boost/utility/enable_if.hpp>
-#include <boost/variant.hpp>
+#include <boost/variant/variant.hpp>
+#include <boost/variant/get.hpp>
+#include <boost/variant/apply_visitor.hpp>
 #include <CGAL/Mesh_3/Has_features.h>
 #include <CGAL/IO/io.h>
 
+#include <tuple>
+
 namespace CGAL {
 namespace Mesh_3 {
 namespace internal {
@@ -49,6 +53,26 @@ struct Index_generator<T, T>
   typedef Index   type;
 };
 
+template <typename MD, bool has_feature = Has_features<MD>::value>
+struct Indices_tuple_generator
+{
+  using type = std::tuple<typename MD::Subdomain_index,
+                          typename MD::Surface_patch_index,
+                          typename MD::Curve_index,
+                          typename MD::Corner_index
+                          >;
+};
+
+template <typename MD>
+struct Indices_tuple_generator<MD, false>
+{
+  using type = std::tuple<typename MD::Subdomain_index,
+                          typename MD::Surface_patch_index>;
+};
+
+template <typename MD>
+using Indices_tuple_t = typename Indices_tuple_generator<MD>::type;
+
 // Nasty meta-programming to get a boost::variant of four types that
 // may not be all different.
 template <typename T0> struct seq1 {
@@ -245,6 +269,65 @@ struct Write_mesh_domain_index<Mesh_domain, false> {
 }; // end template partial specialization
    // Write_mesh_domain_index<Mesh_domain, false>
 
+template <typename, typename Index>
+struct Read_write_index {
+  void operator()(std::ostream& os, int, Index index) const {
+    if(is_ascii(os)) os << oformat(index);
+    else CGAL::write(os, index);
+  }
+  Index operator()(std::istream& is, int) const {
+    Index index;
+    if(is_ascii(is)) is >> iformat(index);
+    else CGAL::read(is, index);
+    return index;
+  }
+};
+
+struct Variant_write_visitor {
+  std::ostream& os;
+  template <typename T>
+  void operator()(T v) const {
+    if(is_ascii(os)) os << CGAL::oformat(v);
+    else CGAL::write(os, v);
+  }
+};
+
+template <typename Index>
+struct Variant_read_visitor {
+  std::istream& is;
+  Index& variant;
+  template <typename T>
+  void operator()(T) const {
+    T v;
+    if(is_ascii(is)) is >> CGAL::iformat(v);
+    else CGAL::read(is, v);
+    variant = v;
+  }
+};
+
+template <typename Indices_types, typename... Args>
+struct Read_write_index<Indices_types, boost::variant<Args...>> {
+  using Index = boost::variant<Args...>;
+  using index_seq = std::make_index_sequence<std::tuple_size<Indices_types>::value>;
+
+  template <std::size_t... Is>
+  Index get_index(int dimension, std::index_sequence<Is...>) const{
+    static const Index variants[] = { std::tuple_element_t<Is, Indices_types>{}... };
+    return variants[dimension < 0 ? 0 : 3-dimension];
+  }
+
+  void operator()(std::ostream& os, int, Index index) const {
+    Variant_write_visitor visitor{os};
+    apply_visitor(visitor, index);
+  }
+  Index operator()(std::istream& is, int dimension) const {
+    Index index = get_index(dimension, index_seq{});
+    Variant_read_visitor<Index> visitor{is, index};
+    apply_visitor(visitor, index);
+    return index;
+  }
+};
+
 } // end namespace internal
 } // end namespace Mesh_3
 } // end namespace CGAL
diff --git a/Mesh_3/test/Mesh_3/CMakeLists.txt b/Mesh_3/test/Mesh_3/CMakeLists.txt
index d731e12d4e4..98cfb1e2b23 100644
--- a/Mesh_3/test/Mesh_3/CMakeLists.txt
+++ b/Mesh_3/test/Mesh_3/CMakeLists.txt
@@ -118,6 +118,12 @@ if ( CGAL_FOUND )
     endforeach()
   endif()
 
+  if(BUILD_TESTING)
+    set_tests_properties(
+      execution___of__test_meshing_polyhedron_with_features
+      execution___of__test_meshing_verbose
+      PROPERTIES RESOURCE_LOCK Mesh_3_Tests_IO)
+  endif()
 else()
 
     message(STATUS "This program requires the CGAL library, and will not be compiled.")
diff --git a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/C3t3_rib_exporter_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/C3t3_rib_exporter_plugin.cpp
index f622cfa4862..795f6a247f9 100644
--- a/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/C3t3_rib_exporter_plugin.cpp
+++ b/Polyhedron/demo/Polyhedron/Plugins/Mesh_3/C3t3_rib_exporter_plugin.cpp
@@ -732,8 +732,8 @@ void
 C3t3_rib_exporter_plugin::
 write_surface_cells(const C3t3& c3t3, const Plane& /* plane */, std::ofstream& out)
 {
-  Geom_traits::Construct_point_3 wp2p
-    = c3t3.triangulation().geom_traits().construct_point_3_object();
+  const Geom_traits& gt = c3t3.triangulation().geom_traits();
+  Geom_traits::Construct_point_3 wp2p = gt.construct_point_3_object();
 
   for ( C3t3::Cells_in_complex_iterator it_cell = c3t3.cells_in_complex_begin(),
        end = c3t3.cells_in_complex_end() ; it_cell != end ; ++it_cell )
@@ -755,7 +755,7 @@ write_surface_cells(const C3t3& c3t3, const Plane& /* plane */, std::ofstream& o
 
     //const int TRANSPARENCY_ALPHA_VALUE = 100;
     CGAL::Bbox_3 bbox = c3t3.bbox();
-    float relPos = static_cast<float>((c->weighted_circumcenter().x() - bbox.xmin())
+    float relPos = static_cast<float>((c->weighted_circumcenter(gt).x() - bbox.xmin())
                                       / (bbox.xmax() - bbox.xmin()));
     float TRANSPARENCY_ALPHA_VALUE =
       1.f -