diff --git a/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp b/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp
new file mode 100644
index 00000000000..2d34f596cd8
--- /dev/null
+++ b/Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp
@@ -0,0 +1,191 @@
+// #define CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
+#define CGAL_PERIODIC_3_MESH_3_VERBOSE
+// #define CGAL_MESH_3_C3T3_HELPERS_VERBOSE
+// #define CGAL_MESH_3_EXUDER_VERBOSE
+#define CGAL_MESH_3_PROTECTION_DEBUG 31
+
+#include <CGAL/Periodic_3_mesh_3/config.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Surface_mesh.h>
+
+#include <CGAL/make_periodic_3_mesh_3.h>
+#include <CGAL/Periodic_3_mesh_3/IO/File_medit.h>
+#include <CGAL/Periodic_3_mesh_triangulation_3.h>
+
+#include <CGAL/boost/graph/Euler_operations.h>
+#include <CGAL/Labeled_mesh_domain_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+#include <CGAL/Mesh_domain_with_polyline_features_3.h>
+#include <CGAL/Polygon_mesh_processing/detect_features.h>
+#include <CGAL/Side_of_triangle_mesh.h>
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+
+namespace CGAL {
+
+// This is a wrapper to convert a polyhedral surface to a periodic polyhedral domain
+// over a user-provided canonical domain.
+//
+// It is the user's responsability to ensure that the polyhedral domain is actually periodic
+// over the canonical domain, i.e. there is periodic continuity at the boundaries
+// of the canonical domain.
+template<class TriangleMesh, class K>
+class Polyhedral_to_periodic_labeling_function_wrapper
+{
+public:
+  using return_type = int;
+  using Point_3 = typename K::Point_3;
+  using GT = typename details::Periodic_3_mesh_geom_traits_generator<K>::type;
+
+private:
+  const TriangleMesh& m_tmesh;
+  CGAL::Side_of_triangle_mesh<TriangleMesh, GT> m_sotm;
+  GT m_gt;
+
+public:
+  explicit Polyhedral_to_periodic_labeling_function_wrapper(const TriangleMesh& tmesh,
+                                                            const CGAL::Iso_cuboid_3<K>& domain)
+    : m_tmesh(tmesh), m_sotm(m_tmesh), m_gt(domain)
+  {
+    CGAL_precondition(CGAL::is_closed(tmesh));
+  }
+
+  Polyhedral_to_periodic_labeling_function_wrapper(const Polyhedral_to_periodic_labeling_function_wrapper& other)
+    : m_tmesh(other.m_tmesh), m_sotm(m_tmesh), m_gt(other.m_gt)
+  { }
+
+  return_type operator()(const Point_3& p) const
+  {
+    const Point_3 cp = P3T3::internal::robust_canonicalize_point(p, m_gt);
+    CGAL::Bounded_side res = m_sotm(cp);
+    return ((res == ON_BOUNDED_SIDE) ? 1 : 2); // set a region to '0' if it is not to be meshed
+  }
+};
+
+} // namespace CGAL
+
+// Kernel
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using FT = K::FT;
+using Point = K::Point_3;
+using Iso_cuboid = K::Iso_cuboid_3;
+
+// Domain
+using Mesh = CGAL::Surface_mesh<Point>;
+using Periodic_polyhedral_domain = CGAL::Polyhedral_to_periodic_labeling_function_wrapper<Mesh, K>;
+
+// Optional: add polyline features to be protected (i.e., preserved in the output)
+// using Periodic_mesh_domain = CGAL::Labeled_mesh_domain_3<K>; // no feature protection
+using Periodic_mesh_domain = CGAL::Mesh_domain_with_polyline_features_3<CGAL::Labeled_mesh_domain_3<K> >;
+using Polyline_3 = std::vector<Point>;
+using Polylines = std::list<Polyline_3>;
+
+// Triangulation
+using Tr = CGAL::Periodic_3_mesh_triangulation_3<Periodic_mesh_domain>::type;
+using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Tr>;
+
+// Criteria
+using Periodic_mesh_criteria = CGAL::Mesh_criteria_3<Tr>;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+// An arbitrary, simple polyhedral shape
+void generate_periodic_diamond(const Point& origin, // bottom, front, left point of the canonical domain
+                               const FT xs, const FT ys, const FT zs,
+                               Mesh& sm,
+                               Polylines& polylines)
+{
+  using vertex_descriptor = boost::graph_traits<Mesh>::vertex_descriptor;
+  using edge_descriptor = boost::graph_traits<Mesh>::edge_descriptor;
+
+  auto vpm = get(CGAL::vertex_point, sm);
+
+  // Three points on in the middle of the canonical domain
+  vertex_descriptor v0 = add_vertex(sm);
+  vertex_descriptor v1 = add_vertex(sm);
+  vertex_descriptor v2 = add_vertex(sm);
+  put(vpm, v0, Point(origin.x() + 0.25 * xs, origin.y() + 0.25 * ys, origin.z() + 0.5 * zs));
+  put(vpm, v1, Point(origin.x() + 0.75 * xs, origin.y() + 0.25 * ys, origin.z() + 0.5 * zs));
+  put(vpm, v2, Point(origin.x() + 0.50 * xs, origin.y() + 0.75 * ys, origin.z() + 0.5 * zs));
+
+  // two points out of the domain, but the intersection of this diamond
+  // with the canonical domain forms a periodic domain
+  vertex_descriptor v_up = add_vertex(sm);
+  vertex_descriptor v_do = add_vertex(sm);
+  put(vpm, v_up, Point(origin.x() + 0.5 * xs, origin.y() + 0.5 * ys, origin.z() - 0.5 * zs));
+  put(vpm, v_do, Point(origin.x() + 0.5 * xs, origin.y() + 0.5 * ys, origin.z() + 1.5 * zs));
+
+  CGAL::Euler::add_face(std::initializer_list<vertex_descriptor>{v1, v_up, v0}, sm);
+  CGAL::Euler::add_face(std::initializer_list<vertex_descriptor>{v2, v_up, v1}, sm);
+  CGAL::Euler::add_face(std::initializer_list<vertex_descriptor>{v0, v_up, v2}, sm);
+
+  CGAL::Euler::add_face(std::initializer_list<vertex_descriptor>{v0, v_do, v1}, sm);
+  CGAL::Euler::add_face(std::initializer_list<vertex_descriptor>{v1, v_do, v2}, sm);
+  CGAL::Euler::add_face(std::initializer_list<vertex_descriptor>{v2, v_do, v0}, sm);
+
+  CGAL::IO::write_OFF("periodic_spike.off", sm, CGAL::parameters::stream_precision(17));
+  assert(is_valid_polygon_mesh(sm));
+
+  // Get the sharp edges to protect, in a generic way
+  using EIFMap = boost::property_map<Mesh, CGAL::edge_is_feature_t>::type;
+  EIFMap eif = get(CGAL::edge_is_feature, sm);
+  CGAL::Polygon_mesh_processing::detect_sharp_edges(sm, 60 /*angle*/, eif);
+
+  std::cout << "Protect:" << std::endl;
+  int sharp_counter = 0;
+  for(edge_descriptor e : edges(sm))
+  {
+    if(get(eif, e))
+    {
+      ++sharp_counter;
+      polylines.emplace_back(std::initializer_list<Point>{get(vpm, source(e, sm)), get(vpm, target(e, sm))});
+      std::cout << get(vpm, source(e, sm)) << " " << get(vpm, target(e, sm)) << std::endl;
+    }
+  }
+
+  std::cout << sharp_counter << " sharp edges to protect" << std::endl;
+}
+
+int main(int argc, char** argv)
+{
+  const FT x_span = (argc > 1) ? atof(argv[1]) : 1;
+  const FT y_span = (argc > 2) ? atof(argv[2]) : x_span;
+  const FT z_span = (argc > 3) ? atof(argv[3]) : x_span;
+  const FT min_span = (std::min)({x_span, y_span, z_span});
+
+  const int number_of_copies_in_output = (argc > 4) ? atoi(argv[4]) : 4; // can be 1, 2, 4, or 8
+
+  Mesh sm;
+  Polylines polylines;
+  generate_periodic_diamond(CGAL::ORIGIN, x_span, y_span, z_span, sm, polylines);
+  Iso_cuboid canonical_cube(0, 0, 0, x_span, y_span, z_span);
+
+  Periodic_polyhedral_domain ppd(sm, canonical_cube);
+  Periodic_mesh_domain domain(ppd, canonical_cube);
+
+  // Insert the features in the domain
+// domain.add_features(polylines.begin(), polylines.end()); // @tmp
+
+  Periodic_mesh_criteria criteria(edge_size = 1 * min_span,
+                                  facet_angle = 30,
+                                  facet_size = 0.035 * min_span,
+                                  facet_distance = 0.025 * min_span,
+                                  cell_radius_edge_ratio = 2.,
+                                  cell_size = 0.05 * min_span);
+
+  // Mesh generation
+  C3t3 c3t3 = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria);
+
+  std::ofstream medit_file("output_periodic_polyhedral_shape.mesh");
+  CGAL::IO::output_periodic_mesh_to_medit(medit_file, c3t3, number_of_copies_in_output);
+
+  std::cout << "EXIT SUCCESS" << std::endl;
+  return 0;
+}