Merge remote-tracking branch 'cgal/master' into Mesh_3-doc_fixes-GF

BGL/include/CGAL/boost/graph/Euler_operations.h

@@ -1432,6 +1432,106 @@ add_face_to_border(typename boost::graph_traits<Graph>::halfedge_descriptor h1,
   return newh;
 }
 
+/**
+ *  \returns `true` if `e` satisfies the *link condition* \cgalCite{degn-tpec-98}, which guarantees that the surface is also 2-manifold after the edge collapse.
+ */
+template<typename Graph>
+bool
+does_satisfy_link_condition(typename boost::graph_traits<Graph>::edge_descriptor e,
+                            const Graph& g)
+{
+  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
+  typedef typename boost::graph_traits<Graph>::halfedge_descriptor halfedge_descriptor;
+  typedef CGAL::Halfedge_around_source_iterator<Graph> out_edge_iterator;
+
+  halfedge_descriptor v0_v1 = halfedge(e,g);
+  halfedge_descriptor v1_v0 = opposite(v0_v1,g);
+
+  vertex_descriptor v0 = target(v1_v0,g), v1 = target(v0_v1,g);
+
+  vertex_descriptor vL = target(next(v0_v1,g),g);
+  vertex_descriptor vR = target(next(v1_v0,g),g);
+
+  out_edge_iterator eb1, ee1 ;
+  out_edge_iterator eb2, ee2 ;
+
+  // The following loop checks the link condition for v0_v1.
+  // Specifically, that for every vertex 'k' adjacent to both 'p and 'q', 'pkq' is a face of the mesh.
+  //
+  for ( boost::tie(eb1,ee1) = halfedges_around_source(v0,g) ;  eb1 != ee1 ; ++ eb1 )
+  {
+    halfedge_descriptor v0_k = *eb1;
+
+    if ( v0_k != v0_v1 )
+    {
+      vertex_descriptor k = target(v0_k,g);
+
+      for ( boost::tie(eb2,ee2) =  halfedges_around_source(k,g) ; eb2 != ee2 ; ++ eb2 )
+      {
+        halfedge_descriptor k_v1 = *eb2;
+
+        if ( target(k_v1,g) == v1 )
+        {
+          // At this point we know p-q-k are connected and we need to determine if this triangle is a face of the mesh.
+          //
+          // Since the mesh is known to be triangular there are at most two faces sharing the edge p-q.
+          //
+          // If p->q is NOT a border edge, the top face is p->q->t where t is target(next(p->q))
+          // If q->p is NOT a border edge, the bottom face is q->p->b where b is target(next(q->p))
+          //
+          // If k is either t or b then p-q-k *might* be a face of the mesh. It won't be if k==t but p->q is border
+          // or k==b but q->b is a border (because in that case even though there exists triangles p->q->t (or q->p->b)
+          // they are holes, not faces)
+          //
+
+          bool lIsFace =   ( vL == k && (! is_border(v0_v1,g)) )
+            || ( vR == k && (! is_border(v1_v0,g)) ) ;
+
+          if ( !lIsFace )
+          {
+            // CGAL_ECMS_TRACE(3,"  k=V" << get(Vertex_index_map,k) << " IS NOT in a face with p-q. NON-COLLAPSABLE edge." ) ;
+            return false ;
+          }
+          else
+          {
+            //CGAL_ECMS_TRACE(4,"  k=V" << get(Vertex_index_map,k) << " is in a face with p-q") ;
+          }
+        }
+      }
+    }
+  }
+
+  // detect isolated triangle (or triangle attached to a mesh with non-manifold vertices)
+  if (!is_border(v0_v1,g) && is_border(opposite(next(v0_v1,g), g), g)
+                          && is_border(opposite(prev(v0_v1,g), g), g) ) return false;
+  if (!is_border(v1_v0,g) && is_border(opposite(next(v1_v0,g), g), g)
+                          && is_border(opposite(prev(v1_v0,g), g), g) ) return false;
+
+  if ( !is_border(v0_v1,g) && !is_border(v1_v0,g) )
+  {
+    if ( is_border(v0,g) && is_border(v1,g) )
+    {
+      //CGAL_ECMS_TRACE(3,"  both p and q are boundary vertices but p-q is not. NON-COLLAPSABLE edge." ) ;
+      return false ;
+    }
+    else
+    {
+      if ( is_tetrahedron(v0_v1,g) )
+      {
+        //CGAL_ECMS_TRACE(3,"  p-q belongs to a tetrahedron. NON-COLLAPSABLE edge." ) ;
+        return false ;
+      }
+      if ( next(v0_v1, g) == opposite(prev(v1_v0, g), g) &&
+           prev(v0_v1, g) == opposite(next(v1_v0, g), g) )
+      {
+        //CGAL_ECMS_TRACE(3,"  degenerate volume." ) ;
+        return false ;
+      }
+    }
+  }
+
+  return true ;
+}
 
 /**
  * collapses an edge in a graph.
@@ -1464,6 +1564,7 @@ collapse_edge(typename boost::graph_traits<Graph>::edge_descriptor e,
   typedef typename Traits::halfedge_descriptor            halfedge_descriptor;
 
   CGAL_precondition(is_valid_edge_descriptor(e, g));
+  CGAL_precondition(does_satisfy_link_condition(e,g));
 
   halfedge_descriptor pq = halfedge(e,g);
   halfedge_descriptor qp = opposite(pq, g);
@@ -1584,6 +1685,7 @@ collapse_edge(typename boost::graph_traits<Graph>::edge_descriptor v0v1,
   typedef typename Traits::halfedge_descriptor            halfedge_descriptor;
 
   CGAL_precondition(is_valid_edge_descriptor(v0v1, g));
+  CGAL_precondition(does_satisfy_link_condition(v0v1,g));
   CGAL_precondition(!get(Edge_is_constrained_map, v0v1));
 
   halfedge_descriptor pq = halfedge(v0v1,g);
@@ -1754,109 +1856,6 @@ flip_edge(typename boost::graph_traits<Graph>::halfedge_descriptor h,
   set_halfedge(foh,oh,g);
 }
 
-/**
- *  \returns `true` if `e` satisfies the *link condition* \cgalCite{degn-tpec-98}, which guarantees that the surface is also 2-manifold after the edge collapse.
- */
-template<typename Graph>
-bool
-does_satisfy_link_condition(typename boost::graph_traits<Graph>::edge_descriptor e,
-                            const Graph& g)
-{
-  typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
-  typedef typename boost::graph_traits<Graph>::halfedge_descriptor halfedge_descriptor;
-  typedef CGAL::Halfedge_around_source_iterator<Graph> out_edge_iterator;
-
-  CGAL_precondition(is_valid_edge_descriptor(e, g));
-
-  halfedge_descriptor v0_v1 = halfedge(e,g);
-  halfedge_descriptor v1_v0 = opposite(v0_v1,g);
-
-  vertex_descriptor v0 = target(v1_v0,g), v1 = target(v0_v1,g);
-
-  vertex_descriptor vL = target(next(v0_v1,g),g);
-  vertex_descriptor vR = target(next(v1_v0,g),g);
-
-  out_edge_iterator eb1, ee1 ;
-  out_edge_iterator eb2, ee2 ;
-
-  // The following loop checks the link condition for v0_v1.
-  // Specifically, that for every vertex 'k' adjacent to both 'p and 'q', 'pkq' is a face of the mesh.
-  //
-  for ( boost::tie(eb1,ee1) = halfedges_around_source(v0,g) ;  eb1 != ee1 ; ++ eb1 )
-  {
-    halfedge_descriptor v0_k = *eb1;
-
-    if ( v0_k != v0_v1 )
-    {
-      vertex_descriptor k = target(v0_k,g);
-
-      for ( boost::tie(eb2,ee2) =  halfedges_around_source(k,g) ; eb2 != ee2 ; ++ eb2 )
-      {
-        halfedge_descriptor k_v1 = *eb2;
-
-        if ( target(k_v1,g) == v1 )
-        {
-          // At this point we know p-q-k are connected and we need to determine if this triangle is a face of the mesh.
-          //
-          // Since the mesh is known to be triangular there are at most two faces sharing the edge p-q.
-          //
-          // If p->q is NOT a border edge, the top face is p->q->t where t is target(next(p->q))
-          // If q->p is NOT a border edge, the bottom face is q->p->b where b is target(next(q->p))
-          //
-          // If k is either t or b then p-q-k *might* be a face of the mesh. It won't be if k==t but p->q is border
-          // or k==b but q->b is a border (because in that case even though there exists triangles p->q->t (or q->p->b)
-          // they are holes, not faces)
-          //
-
-          bool lIsFace =   ( vL == k && (! is_border(v0_v1,g)) )
-            || ( vR == k && (! is_border(v1_v0,g)) ) ;
-
-          if ( !lIsFace )
-          {
-            // CGAL_ECMS_TRACE(3,"  k=V" << get(Vertex_index_map,k) << " IS NOT in a face with p-q. NON-COLLAPSABLE edge." ) ;
-            return false ;
-          }
-          else
-          {
-            //CGAL_ECMS_TRACE(4,"  k=V" << get(Vertex_index_map,k) << " is in a face with p-q") ;
-          }
-        }
-      }
-    }
-  }
-
-  // detect isolated triangle (or triangle attached to a mesh with non-manifold vertices)
-  if (!is_border(v0_v1,g) && is_border(opposite(next(v0_v1,g), g), g)
-                          && is_border(opposite(prev(v0_v1,g), g), g) ) return false;
-  if (!is_border(v1_v0,g) && is_border(opposite(next(v1_v0,g), g), g)
-                          && is_border(opposite(prev(v1_v0,g), g), g) ) return false;
-
-  if ( !is_border(v0_v1,g) && !is_border(v1_v0,g) )
-  {
-    if ( is_border(v0,g) && is_border(v1,g) )
-    {
-      //CGAL_ECMS_TRACE(3,"  both p and q are boundary vertices but p-q is not. NON-COLLAPSABLE edge." ) ;
-      return false ;
-    }
-    else
-    {
-      if ( is_tetrahedron(v0_v1,g) )
-      {
-        //CGAL_ECMS_TRACE(3,"  p-q belongs to a tetrahedron. NON-COLLAPSABLE edge." ) ;
-        return false ;
-      }
-      if ( next(v0_v1, g) == opposite(prev(v1_v0, g), g) &&
-           prev(v0_v1, g) == opposite(next(v1_v0, g), g) )
-      {
-        //CGAL_ECMS_TRACE(3,"  degenerate volume." ) ;
-        return false ;
-      }
-    }
-  }
-
-  return true ;
-}
-
 #ifndef CGAL_NO_DEPRECATED_CODE
 /// \cond SKIP_IN_MANUAL
 template<typename Graph>

Documentation/doc/resources/1.8.13/BaseDoxyfile.in

@@ -461,18 +461,6 @@ HTML_HEADER            = ${CGAL_DOC_HEADER_PACKAGE}
 
 HTML_FOOTER            = ${CGAL_DOC_RESOURCE_DIR}/footer.html
 
-# The HTML_STYLESHEET tag can be used to specify a user-defined cascading style
-# sheet that is used by each HTML page. It can be used to fine-tune the look of
-# the HTML output. If left blank doxygen will generate a default style sheet.
-# See also section "Doxygen usage" for information on how to generate the style
-# sheet that doxygen normally uses.
-# Note: It is recommended to use HTML_EXTRA_STYLESHEET instead of this tag, as
-# it is more robust and this tag (HTML_STYLESHEET) will in the future become
-# obsolete.
-# This tag requires that the tag GENERATE_HTML is set to YES.
-
-HTML_STYLESHEET        = ${CGAL_DOC_RESOURCE_DIR}/stylesheet.css
-
 # If the HTML_TIMESTAMP tag is set to YES then the footer of each generated HTML
 # page will contain the date and time when the page was generated. Setting this
 # to YES can help to show when doxygen was last run and thus if the

GraphicsView/include/CGAL/Qt/vec.h

@@ -13,6 +13,8 @@
 #ifndef QGLVIEWER_VEC_H
 #define QGLVIEWER_VEC_H
 
+#include <QtGlobal>
+
 #include <iostream>
 #include <math.h>
 

Installation/CHANGES.md

@@ -112,6 +112,7 @@ extracted from labeled images.
 -   Added new meshing criteria `facet_min_size` and `cell_min_size` to prevent Delaunay refinement from creating simplices smaller than the prescribed bound.
 
 ### [3D Periodic Mesh Generation](https://doc.cgal.org/5.6/Manual/packages.html#PkgPeriodic3Mesh3)
+-   Periodic Mesh Generation now supports non-cubic domains.
 -   Deprecated usage of boost parameters in favor of function named parameters.
 
 ### [2D Hyperbolic Triangulations](https://doc.cgal.org/5.6/Manual/packages.html#PkgHyperbolicTriangulation2)

Installation/CMakeLists.txt

@@ -930,6 +930,9 @@ if(CGAL_BRANCH_BUILD)
     find_package(Eigen3 REQUIRED)
     find_package(Qt5 COMPONENTS Core Widgets OpenGL Gui REQUIRED)
     find_package(VTK COMPONENTS vtkImagingGeneral vtkIOImage NO_MODULE)
+    if(VTK_FOUND)
+      get_target_property(VTK_INCLUDE_DIRS VTK::IOImage INTERFACE_INCLUDE_DIRECTORIES)
+    endif()
     find_package(IPE)
     find_package(RS3)
     find_package(LEDA)
@@ -1197,6 +1200,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}/../${package}/include/${header}"
       endif()
     endforeach() # loop on packages
     #Now check that a cpp file including all documented headers compiles
+    file(WRITE ${CGAL_BINARY_DIR}/test_headers.cpp "#include <CGAL/Periodic_3_mesh_3/config.h>\n")
     foreach(header ${list_of_headers_to_test})
       file(APPEND ${CGAL_BINARY_DIR}/test_headers.cpp "#include <${header}>\n")
     endforeach() #loop on headers to include in test file

Kernel_23/doc/Kernel_23/Concepts/FunctionObjectConcepts.h

@@ -7254,6 +7254,20 @@ public:
                              const Kernel::Weighted_point_3& q,
                              const Kernel::Weighted_point_3& r,
                              const Kernel::Weighted_point_3& s);
+
+  /*!
+  constructs the point which is the center of the smallest orthogonal sphere to the input weighted points.
+  */
+  Kernel::Point_3 operator()(const Kernel::Weighted_point_3& p,
+                             const Kernel::Weighted_point_3& q,
+                             const Kernel::Weighted_point_3& r);
+
+  /*!
+  constructs the point which is the center of the smallest orthogonal sphere to the input weighted points.
+  */
+  Kernel::Point_3 operator()(const Kernel::Weighted_point_3& p,
+                             const Kernel::Weighted_point_3& q);
+
   /// @}
 }; /* end Kernel::ConstructWeightedCircumcenter_3 */
 

Kernel_d/doc/Kernel_d/Concepts/Kernel--Contained_in_affine_hull_d.h

@@ -19,7 +19,7 @@ affine hull of the points in `A = tuple [first,last)`.
 
 \tparam ForwardIterator has `Kernel_d::Point_d` as value type.
 */
-template <class ForwardIterator> Bounded_side
+template <class ForwardIterator> bool
 operator()( ForwardIterator first, ForwardIterator last, const
 Kernel_d::Point_d& p);
 

Kernel_d/doc/Kernel_d/Concepts/Kernel--Contained_in_linear_hull_d.h

@@ -18,7 +18,7 @@ linear hull of the vectors in `A = tuple [first,last)`.
 \pre The objects are of the same dimension.
 \tparam ForwardIterator has `Kernel_d::Vector_d` as value type.
 */
-template <class ForwardIterator> Bounded_side
+template <class ForwardIterator> bool
 operator()( ForwardIterator first, ForwardIterator last, const
 Kernel_d::Vector_d& v);
 

Kernel_d/doc/Kernel_d/Concepts/Kernel--Contained_in_simplex_d.h

@@ -19,7 +19,7 @@ simplex of the points in `A = tuple [first,last)`.
 \pre The objects in \f$ A\f$ are of the same dimension and affinely independent.
 \tparam ForwardIterator has `Kernel_d::Point_d` as value type.
 */
-template <class ForwardIterator> Bounded_side
+template <class ForwardIterator> bool
 operator()( ForwardIterator first, ForwardIterator last, const
 Kernel_d::Point_d& p);
 

Mesh_2/include/CGAL/Mesh_2/Refine_edges.h

@@ -569,7 +569,7 @@ public:
     std::cerr << "E";
 #endif
     // @todo Perhaps we should remove destroyed edges too.
-    // @warning This code has been rewroten!
+    // @warning This code has been rewritten!
 
     Face_circulator fc = tr.incident_faces(v), fcbegin(fc);
     if( fc == 0 ) return;
@@ -586,7 +586,7 @@ public:
     Face_handle fh;
     int index = 0; // Avoids a warning.
                    // We know that is_edge must return true, and is_edge will assign something to index
-                   // but the compiler does not so it will issue a maybe uninitialized warning
+                   // but the compiler does not so it will issue a maybe-uninitialized warning
 
     CGAL_assume_code(bool is_edge = )
     tr.is_edge(va, v, fh, index);

Mesh_3/doc/Mesh_3/Concepts/MeshTriangulationTraits_3.h

@@ -78,6 +78,15 @@ public:
   */
   typedef unspecified_type Is_degenerate_3;
 
+  /*!
+  A predicate object that must provide the function operator:
+
+  `bool operator()(Point_3 p, Point_3 q, Point_3 r)`
+
+  which returns `true` iff `p`, `q`, and `r` are collinear.
+  */
+  typedef unspecified_type Collinear_3;
+
   /*!
   A constructor object that must provide the function operators:
 
@@ -228,6 +237,19 @@ public:
   /*!
   A constructor object that must provide the function operators:
 
+  `Point_3 operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s)`,
+
+  `Point_3 operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s)`,
+
+  `Point_3 operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s)`,
+
+  which return the center of the smallest orthogonal sphere to the input weighted points.
+  */
+  typedef unspecified_type Construct_weighted_circumcenter_3;
+
+  /*!
+  A constructor object that must provide the function operators:
+
   `Point_3 operator()(Point_3 p, Point_3 q, Point_3 r)`,
 
   `Point_3 operator()(Tetrahedron_3 t)`,

Mesh_3/include/CGAL/Labeled_mesh_domain_3.h

@@ -1153,7 +1153,7 @@ public:
     Intersection operator()(const Segment_3& s) const
     {
 #ifndef CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
-      CGAL_precondition(r_domain_.do_intersect_surface_object()(s));
+      CGAL_precondition(r_domain_.do_intersect_surface_object()(s) != boost::none);
 #endif // NOT CGAL_MESH_3_NO_LONGER_CALLS_DO_INTERSECT_3
       return this->operator()(s.source(),s.target());
     }

Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h

@@ -1776,15 +1776,15 @@ private:
 
   template <typename OutputIterator>
   OutputIterator
-  get_conflict_zone_topo_change(const Vertex_handle& vertex,
-                                const Weighted_point& conflict_point,
-                                OutputIterator conflict_cells) const;
+  get_conflict_zone_after_move_topo_change(const Vertex_handle& new_vertex,
+                                           const Weighted_point& old_position,
+                                           OutputIterator conflict_cells) const;
 
   template <typename CellsOutputIterator,
             typename FacetsOutputIterator>
   void
-  get_conflict_zone_topo_change(const Vertex_handle& v,
-                                const Weighted_point& conflict_point,
+  get_conflict_zone_topo_change(const Vertex_handle& old_vertex,
+                                const Weighted_point& new_position,
                                 CellsOutputIterator insertion_conflict_cells,
                                 FacetsOutputIterator insertion_conflict_boundary,
                                 CellsOutputIterator removal_conflict_cells,
@@ -2617,12 +2617,11 @@ update_mesh_topo_change(const Vertex_handle& old_vertex,
   }
   else
   {
-    // Removing from c3t3 cells which will be destroyed by revert_move
+    // Removing from c3t3 the cells which will be destroyed by revert_move
     // is done by move_point_topo_change_conflict_zone_known, called by revert_move
 
 #ifdef CGAL_MESH_3_C3T3_HELPERS_VERBOSE
-     std::cerr << "update_mesh_topo_change: revert move to "
-               << old_position << "\n";
+     std::cerr << "update_mesh_topo_change: revert move to " << old_position << "\n";
 #endif
 
     //reset caches in case cells are re-used by the compact container
@@ -2632,7 +2631,7 @@ update_mesh_topo_change(const Vertex_handle& old_vertex,
 
     // Revert move
     Vertex_handle revert_vertex = revert_move(new_vertex, old_position,
-                          std::inserter(outdated_cells, outdated_cells.end()));
+                                              std::inserter(outdated_cells, outdated_cells.end()));
 
     //restore meta-data (cells should have same connectivity as before move)
     //cells should be the same (connectivity-wise) as before initial move
@@ -2662,8 +2661,8 @@ update_mesh(const Vertex_handle& old_vertex,
   Vertex_handle new_vertex = move_point(old_vertex, move,
                                         std::back_inserter(outdated_cells),
                                         CGAL::Emptyset_iterator());
-  // move_point has invalidated caches
 
+  // move_point has invalidated caches
   restore_mesh(outdated_cells.begin(), outdated_cells.end());
 
   // Fill modified vertices
@@ -3150,7 +3149,8 @@ move_point_topo_change(const Vertex_handle& old_vertex,
                                 could_lock_zone);
   if (insertion_conflict_cells.empty())
     return old_vertex;//new_position coincides with an existing vertex (not old_vertex)
-                      //and old_vertex should not be removed of the nb_vertices will change
+                      //and old_vertex should not be removed if the nb_vertices will change
+
   reset_circumcenter_cache(removal_conflict_cells);
   reset_sliver_cache(removal_conflict_cells);
   reset_circumcenter_cache(insertion_conflict_cells);
@@ -3260,7 +3260,9 @@ move_point_topo_change_conflict_zone_known(
   // Remove conflict zone cells from c3t3 (they will be deleted by insert/remove)
   remove_cells_and_facets_from_c3t3(conflict_zone.begin(), conflict_zone.end());
 
-// Start Move point // Insert new_vertex, remove old_vertex
+// Start Move point
+
+  // Insert new_vertex, remove old_vertex
   int dimension = c3t3_.in_dimension(old_vertex);
   Index vertex_index = c3t3_.index(old_vertex);
   FT meshing_info = old_vertex->meshing_info();
@@ -3285,14 +3287,16 @@ move_point_topo_change_conflict_zone_known(
   c3t3_.set_dimension(new_vertex,dimension);
   c3t3_.set_index(new_vertex,vertex_index);
   new_vertex->set_meshing_info(meshing_info);
-  // End Move point
+
+// End Move point
 
   //// Fill outdated_cells
-  // Get conflict zone in new triangulation and set cells outdated
+  // Get the union of the cells impacted by the insertion and the removal
   Cell_vector new_conflict_cells;
   new_conflict_cells.reserve(64);
-  get_conflict_zone_topo_change(new_vertex, old_position,
-                                std::back_inserter(new_conflict_cells));
+  get_conflict_zone_after_move_topo_change(new_vertex, old_position,
+                                           std::back_inserter(new_conflict_cells));
+
   std::copy(new_conflict_cells.begin(),new_conflict_cells.end(),outdated_cells);
 
   // Fill deleted_cells
@@ -3432,6 +3436,7 @@ get_least_square_surface_plane(const Vertex_handle& v,
                                Surface_patch_index patch_index) const
 {
   typedef typename C3T3::Triangulation::Triangle Triangle;
+
   typename GT::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
 
   // Get incident facets
@@ -3461,17 +3466,10 @@ get_least_square_surface_plane(const Vertex_handle& v,
          (patch_index == Surface_patch_index() ||
           c3t3_.surface_patch_index(f) == patch_index) )
     {
-      ref_facet = f;
+      if(ref_facet.first == Cell_handle())
+        ref_facet = f;
 
-      // In the case of a periodic triangulation, the incident facets of a point
-      // do not necessarily have the same offsets. Worse, the surface centers
-      // might not have the same offset as their facet. Thus, no solution except
-      // calling a function 'get_closest_triangle(p, t)' that simply returns t
-      // for a non-periodic triangulation, and checks all possible offsets for
-      // periodic triangulations
-
-      Triangle t = c3t3_.triangulation().triangle(f);
-      Triangle ct = tr_.get_closest_triangle(cp(position), t);
+      const Triangle ct = tr_.get_incident_triangle(f, v);
       triangles.push_back(ct);
     }
   }
@@ -3483,7 +3481,6 @@ get_least_square_surface_plane(const Vertex_handle& v,
   // Compute least square fitting plane
   Plane_3 plane;
   Bare_point point;
-
   CGAL::linear_least_squares_fitting_3(triangles.begin(),
                                        triangles.end(),
                                        plane,
@@ -3492,12 +3489,12 @@ get_least_square_surface_plane(const Vertex_handle& v,
                                        tr_.geom_traits(),
                                        Default_diagonalize_traits<double, 3>());
 
-   return std::make_pair(plane,
-     ref_facet.first->get_facet_surface_center(ref_facet.second));
+  // The surface center of a facet might have an offset in periodic triangulations
+  const Bare_point& ref_facet_scp = ref_facet.first->get_facet_surface_center(ref_facet.second);
+  const Bare_point& ref_point = tr_.get_closest_point(cp(position), ref_facet_scp);
+  return std::make_pair(plane, ref_point);
 }
 
-
-
 template <typename C3T3, typename MD>
 typename C3T3_helpers<C3T3,MD>::Bare_point
 C3T3_helpers<C3T3,MD>::
@@ -3519,6 +3516,7 @@ project_on_surface_if_possible(const Vertex_handle& v,
                                const Bare_point& p,
                                Surface_patch_index index) const
 {
+  // @todo should call below if it's available...
   // return domain_.project_on_surface(p);
 
   typename GT::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
@@ -3848,8 +3846,8 @@ template <typename CellsOutputIterator,
           typename FacetsOutputIterator>
 void
 C3T3_helpers<C3T3,MD>::
-get_conflict_zone_topo_change(const Vertex_handle& v,
-                              const Weighted_point& conflict_point,
+get_conflict_zone_topo_change(const Vertex_handle& old_vertex,
+                              const Weighted_point& new_position,
                               CellsOutputIterator insertion_conflict_cells,
                               FacetsOutputIterator insertion_conflict_boundary,
                               CellsOutputIterator removal_conflict_cells,
@@ -3863,21 +3861,20 @@ get_conflict_zone_topo_change(const Vertex_handle& v,
 // Parallel
   if (boost::is_convertible<Concurrency_tag, Parallel_tag>::value)
   {
-    tr_.incident_cells_threadsafe(v, removal_conflict_cells);
+    tr_.incident_cells_threadsafe(old_vertex, removal_conflict_cells);
   }
   // Sequential
   else
 # endif // CGAL_LINKED_WITH_TBB
   {
-    tr_.incident_cells(v, removal_conflict_cells);
+    tr_.incident_cells(old_vertex, removal_conflict_cells);
   }
 
   // Get conflict_point conflict zone
   int li=0;
   int lj=0;
   typename Tr::Locate_type lt;
-  Cell_handle cell = tr_.locate(
-    conflict_point, lt, li, lj, v->cell(), could_lock_zone);
+  Cell_handle cell = tr_.locate(new_position, lt, li, lj, old_vertex->cell(), could_lock_zone);
 
   if (could_lock_zone && *could_lock_zone == false)
     return;
@@ -3886,7 +3883,7 @@ get_conflict_zone_topo_change(const Vertex_handle& v,
     return;
 
   // Find conflict zone
-  tr_.find_conflicts(conflict_point,
+  tr_.find_conflicts(new_position,
                      cell,
                      insertion_conflict_boundary,
                      insertion_conflict_cells,
@@ -3897,34 +3894,45 @@ template <typename C3T3, typename MD>
 template <typename OutputIterator>
 OutputIterator
 C3T3_helpers<C3T3,MD>::
-get_conflict_zone_topo_change(const Vertex_handle& vertex,
-                              const Weighted_point& conflict_point,
-                              OutputIterator conflict_cells) const
+get_conflict_zone_after_move_topo_change(const Vertex_handle& new_vertex,
+                                         const Weighted_point& old_position,
+                                         OutputIterator conflict_cells) const
 {
-  // Get triangulation_vertex incident cells
+  // Gather the impacted cells: the union of `old_point` conflict zone and `new_vertex` incident cells
+
+  // Get the incident cells of `new_vertex`
   Cell_vector incident_cells_;
   incident_cells_.reserve(64);
-  tr_.incident_cells(vertex, std::back_inserter(incident_cells_));
+  tr_.incident_cells(new_vertex, std::back_inserter(incident_cells_));
 
-  // Get conflict_point conflict zone
+  // Get the conflict zone of `old_point`
   Cell_vector deleted_cells;
   deleted_cells.reserve(64);
 
-  // Vertex removal is forbidden
   int li=0;
   int lj=0;
-  typename Tr::Locate_type locate_type;
-  Cell_handle cell = tr_.locate(conflict_point, locate_type, li, lj, vertex->cell());
+  typename Tr::Locate_type lt;
+  Cell_handle cell = tr_.locate(old_position, lt, li, lj, new_vertex->cell());
 
-  if ( Tr::VERTEX == locate_type )
-    return conflict_cells;
+  // `Periodic_mesh_triangulation::remove()` can refuse to remove a point if this removal
+  // would compromise the 1-cover property (i.e., no too-long edges).
+  // The cells incident to `old_vertex` have not been modified at the TDS level if removal
+  // was rejected, but they still must be gathered here because of the call to
+  // `remove_cells_and_facets_from_c3t3()`
+  if (lt == Tr::VERTEX)
+  {
+    CGAL_assertion((std::is_same<typename Tr::Periodic_tag, CGAL::Tag_true>::value));
+    tr_.incident_cells(cell->vertex(li), std::back_inserter(deleted_cells));
+  }
+  else
+  {
+    tr_.find_conflicts(old_position,
+                       cell,
+                       CGAL::Emptyset_iterator(),
+                       std::back_inserter(deleted_cells),
+                       CGAL::Emptyset_iterator());
+  }
 
-  // Find conflict zone
-  tr_.find_conflicts(conflict_point,
-                     cell,
-                     CGAL::Emptyset_iterator(),
-                     std::back_inserter(deleted_cells),
-                     CGAL::Emptyset_iterator());
 
   // Compute union of conflict_point conflict zone and triangulation_vertex
   // incident cells
@@ -3938,7 +3946,6 @@ get_conflict_zone_topo_change(const Vertex_handle& vertex,
   return conflict_cells;
 }
 
-
 template <typename C3T3, typename MD>
 typename C3T3_helpers<C3T3,MD>::Facet_boundary
 C3T3_helpers<C3T3,MD>::

Mesh_3/include/CGAL/Mesh_3/Protect_edges_sizing_field.h

@@ -561,6 +561,12 @@ insert_corners()
   Initial_corners corners;
   domain_.get_corners(std::back_inserter(corners));
 
+#if CGAL_MESH_3_PROTECTION_DEBUG & 1
+  std::cout << corners.size() << " corners to treat" << std::endl;
+  for(const auto& e : corners)
+    std::cout << "Corner #" << CGAL::IO::oformat(e.first) << ", " <<  e.second << std::endl;
+#endif
+
   Dt dt;
   for ( typename Initial_corners::iterator it = corners.begin(),
        end = corners.end() ; it != end ; ++it )

Mesh_3/include/CGAL/Mesh_3/Refine_cells_3.h

@@ -469,6 +469,7 @@ public:
   std::string debug_info_element_impl(const Cell_handle &ch) const
   {
     std::stringstream sstr;
+    sstr.precision(17);
     sstr << "Cell " << (void*)(ch.operator->()) << " { " << std::endl
     << "  " << *ch->vertex(0) << std::endl
     << "  " << *ch->vertex(1) << std::endl
@@ -749,8 +750,9 @@ int
 Refine_cells_3<Tr,Cr,MD,C3T3_,P_,Ct,C_>::
 number_of_bad_elements_impl()
 {
-  typedef typename MD::Subdomain Subdomain;
-  typedef typename Tr::Finite_cells_iterator Finite_cell_iterator;
+  typedef typename MD::Subdomain_index        Subdomain_index;
+  typedef boost::optional<Subdomain_index>    Subdomain;
+  typedef typename Tr::Finite_cells_iterator  Finite_cell_iterator;
 
   int count = 0;
 #if defined(CGAL_MESH_3_VERBOSE) || defined(CGAL_MESH_3_PROFILING)

Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h

@@ -348,6 +348,7 @@ public:
   std::string debug_info_element_impl(const Facet &facet) const
   {
     std::stringstream sstr;
+    sstr.precision(17);
     sstr << "Facet { " << std::endl
     << "  " << *facet.first->vertex((facet.second+1)%4) << std::endl
     << "  " << *facet.first->vertex((facet.second+2)%4) << std::endl
@@ -477,7 +478,7 @@ protected:
   /// Insert facet into refinement queue
   void insert_bad_facet(Facet facet, const Quality& quality)
   {
-#if CGAL_MESH_3_VERY_VERBOSE
+#ifdef CGAL_MESH_3_VERY_VERBOSE
     std::stringstream s;
     s << "insert_bad_facet(" << debug_info_element_impl(facet) << ", ...) by thread "
       << std::this_thread::get_id() << '\n';
@@ -1051,7 +1052,8 @@ int
 Refine_facets_3<Tr,Cr,MD,C3T3_,P_,Ct,B_,C_>::
 number_of_bad_elements_impl()
 {
-  typedef typename MD::Subdomain Subdomain;
+  typedef typename MD::Subdomain_index        Subdomain_index;
+  typedef boost::optional<Subdomain_index>    Subdomain;
   typedef typename Tr::Finite_facets_iterator Finite_facet_iterator;
 
   int count = 0, count_num_bad_surface_facets = 0;
@@ -1732,6 +1734,9 @@ Refine_facets_3_base<Tr,Cr,MD,C3T3_,Ct,C_>::
 is_facet_encroached(const Facet& facet,
                     const Weighted_point& point) const
 {
+  typedef typename MD::Subdomain_index        Subdomain_index;
+  typedef boost::optional<Subdomain_index>    Subdomain;
+
   if ( r_tr_.is_infinite(facet) || ! this->is_facet_on_surface(facet) )
     return false;
 
@@ -1741,9 +1746,40 @@ is_facet_encroached(const Facet& facet,
   const Bare_point& center = get_facet_surface_center(facet);
   const Weighted_point& reference_point = r_tr_.point(cell, (facet_index+1)&3);
 
+#ifdef CGAL_MESHES_DEBUG_REFINEMENT_POINTS
+  std::cout << "---------------------------------------------------" << std::endl;
+  std::cout << "Facet " << r_tr_.point(cell, (facet_index+1)%4) << " "
+                        << r_tr_.point(cell, (facet_index+2)%4) << " "
+                        << r_tr_.point(cell, (facet_index+3)%4) << std::endl;
+  std::cout << "center: " << center << std::endl;
+  std::cout << "cell point: " << reference_point << std::endl;
+  std::cout << "refinement point: " << point << std::endl;
+  std::cout << "greater or equal? " << r_tr_.greater_or_equal_power_distance(center, reference_point, point) << std::endl;
+  std::cout << "greater or equal (other way)? " << r_tr_.greater_or_equal_power_distance(center, point, reference_point) << std::endl;
+  std::cout << "index of cell " << r_c3t3_.subdomain_index(cell) << std::endl;
+#endif
+
   // the facet is encroached if the new point is closer to the center than
   // any vertex of the facet
-  return r_tr_.greater_or_equal_power_distance(center, reference_point, point);
+  if(r_tr_.greater_or_equal_power_distance(center, reference_point, point))
+    return true;
+
+  // In an ideal (exact) world, when the predicate above returns true then the insertion
+  // of the refinement point will shorten the dual of the facet but that dual will
+  // still intersects the surface (domain), otherwise the power distance to the surface
+  // center would be shorter.
+  //
+  // In the real world, we can make an error both when we switch back to the inexact kernel
+  // and when we evaluate the domain (e.g. trigonometry-based implicit functions).
+  //
+  // An issue can then arise when we update the restricted Delaunay due to the insertion
+  // of another point, and we do not notice that a facet should in fact have been encroached
+  // by a previous insertion.
+  Bare_point cc;
+  r_tr_.dual_exact(facet, point, cc);
+
+  const Subdomain subdomain = r_oracle_.is_in_domain_object()(cc);
+  return (!subdomain || *subdomain != r_c3t3_.subdomain_index(cell));
 }
 
 template<class Tr, class Cr, class MD, class C3T3_, class Ct, class C_>

Mesh_3/include/CGAL/Mesh_3/Slivers_exuder.h

@@ -964,9 +964,9 @@ pump_vertices(FT sliver_criterion_limit,
       bool vertex_pumped = false;
       for( int i = 0; i < 4; ++i )
       {
-        // pump_vertices_on_surfaces is a boolean template parameter.  The
-        // following condition is pruned at compiled time, if
-        // pump_vertices_on_surfaces==false.
+        // pump_vertices_on_surfaces is a Boolean template parameter.
+        // The following condition is pruned at compile time,
+        // if pump_vertices_on_surfaces is `false`.
         if( pump_vertices_on_surfaces || c3t3_.in_dimension(c->vertex(i)) > 2 )
         {
           if( pump_vertex<pump_vertices_on_surfaces>(c->vertex(i)) )
@@ -976,7 +976,9 @@ pump_vertices(FT sliver_criterion_limit,
             break;
           }
           else
+          {
             ++num_of_ignored_vertices_;
+          }
 
           ++num_of_treated_vertices_;
         }
@@ -987,14 +989,14 @@ pump_vertices(FT sliver_criterion_limit,
         this->cells_queue_pop_front();
 
       visitor.after_cell_pumped(this->cells_queue_size());
-  #ifdef CGAL_MESH_3_EXUDER_VERBOSE
+#ifdef CGAL_MESH_3_EXUDER_VERBOSE
       std::cerr << boost::format("\r             \r"
                                  "(%1%,%2%,%3%) (%|4$.1f| vertices/s)")
         % this->cells_queue_size()
         % num_of_pumped_vertices_
         % num_of_ignored_vertices_
         % (num_of_treated_vertices_ / running_time_.time());
-  #endif // CGAL_MESH_3_EXUDER_VERBOSE
+#endif // CGAL_MESH_3_EXUDER_VERBOSE
     }
   }
 
@@ -1054,12 +1056,12 @@ pump_vertex(const Vertex_handle& pumped_vertex,
   {
     typename GT::Construct_point_3 cp = tr_.geom_traits().construct_point_3_object();
 
-    const Weighted_point& pwp = tr_.point(pumped_vertex);
-    Weighted_point wp(cp(pwp), best_weight);
+    const Weighted_point& old_position = tr_.point(pumped_vertex);
+    Weighted_point new_point(cp(old_position), best_weight);
 
     // Insert weighted point into mesh
     // note it can fail if the mesh is non-manifold at pumped_vertex
-    return update_mesh<pump_vertices_on_surfaces>(wp,
+    return update_mesh<pump_vertices_on_surfaces>(new_point,
                                                   pumped_vertex,
                                                   could_lock_zone);
   }

Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h

@@ -414,7 +414,7 @@ get_sq_distance_to_closest_vertex(const Tr& tr,
   typedef std::vector<Vertex_handle>              Vertex_container;
 
   // There is no need to use tr.min_squared_distance() here because we are computing
-  // distances between 'v' and a neighbor within their common cell, which means
+  // distances between 'v' and a neighboring vertex within a common cell, which means
   // that even if we are using a periodic triangulation, the distance is correctly computed.
   typename GT::Compute_squared_distance_3 csqd = tr.geom_traits().compute_squared_distance_3_object();
   typename GT::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
@@ -475,7 +475,7 @@ get_sq_distance_to_closest_vertex(const Tr& tr,
   typedef typename Vertex_container::iterator                VC_it;
 
   // There is no need to use tr.min_squared_distance() here because we are computing
-  // distances between 'v' and a neighbor within their common cell, which means
+  // distances between 'v' and a neighboring vertex within a common cell, which means
   // that even if we are using a periodic triangulation, the distance is correctly computed.
   typename GT::Compute_squared_distance_3 csqd = tr.geom_traits().compute_squared_distance_3_object();
   typename GT::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();

Mesh_3/include/CGAL/Mesh_3/radius_ratio.h

@@ -29,8 +29,7 @@ radius_ratio(const typename K::Point_3& p0,
              K k = K())
 {
   typedef typename K::FT FT;
-  typename K::Compute_squared_distance_3 sq_distance =
-    k.compute_squared_distance_3_object();
+
   typename K::Compute_squared_radius_3 comp_sq_circumradius =
     k.compute_squared_radius_3_object();
   typename K::Compute_volume_3 volume =

Mesh_3/include/CGAL/Mesh_3/vertex_perturbation.h

@@ -119,10 +119,13 @@ edge_sq_length(const typename Tr::Edge& e,
   typedef typename Tr::Bare_point      Bare_point;
   typedef typename Tr::Weighted_point  Weighted_point;
 
-  typename GT::Construct_point_3 cp =
-    tr.geom_traits().construct_point_3_object();
+  // There is no need to use tr.min_squared_distance() here because we are computing
+  // distances between vertices within a common cell, which means that even
+  // if we are using a periodic triangulation, the distance is correctly computed.
   typename GT::Compute_squared_distance_3 sq_distance =
     tr.geom_traits().compute_squared_distance_3_object();
+  typename GT::Construct_point_3 cp =
+    tr.geom_traits().construct_point_3_object();
 
   const Weighted_point& wp = tr.point(e.first, e.second);
   const Weighted_point& wq = tr.point(e.first, e.third);
@@ -979,10 +982,13 @@ private:
     CGAL_assertion(cell->has_vertex(v));
     const typename C3T3::Triangulation& tr = c3t3.triangulation();
 
-    typename GT::Construct_point_3 cp =
-      tr.geom_traits().construct_point_3_object();
+    // There is no need to use tr.min_squared_distance() here because we are computing
+    // distances between vertices within a common cell, which means that even
+    // if we are using a periodic triangulation, the distance is correctly computed.
     typename GT::Compute_squared_distance_3 sq_distance =
       tr.geom_traits().compute_squared_distance_3_object();
+    typename GT::Construct_point_3 cp =
+      tr.geom_traits().construct_point_3_object();
 
     const int i = cell->index(v);
     const Weighted_point& wp0 = tr.point(cell, i);

Mesh_3/include/CGAL/Mesh_domain_with_polyline_features_3.h

@@ -883,7 +883,7 @@ public:
   }
 
   void build_curves_aabb_tree() const {
-#if CGAL_MESH_3_VERBOSE
+#ifdef CGAL_MESH_3_VERBOSE
     std::cerr << "Building curves AABB tree...";
     CGAL::Real_timer timer;
     timer.start();
@@ -909,7 +909,7 @@ public:
     }
     curves_aabb_tree_ptr_->build();
     curves_aabb_tree_is_built = true;
-#if CGAL_MESH_3_VERBOSE
+#ifdef CGAL_MESH_3_VERBOSE
     timer.stop();
     std::cerr << " done (" << timer.time() * 1000 << " ms)" << std::endl;
 #endif

Mesh_3/include/CGAL/Mesh_triangulation_3.h

@@ -62,12 +62,14 @@ public:
   typedef typename Base::Triangle                             Triangle;
 
   typedef typename Base::Vertex_handle                        Vertex_handle;
+  typedef typename Base::Facet                                Facet;
   typedef typename Base::Cell_handle                          Cell_handle;
 
   typedef typename Geom_traits::Vector_3                      Vector;
 
   using Base::geom_traits;
   using Base::point;
+  using Base::triangle;
 
   static std::string io_signature() { return Get_io_signature<Base>()(); }
 
@@ -83,9 +85,9 @@ public:
     return q;
   }
 
-  const Triangle& get_closest_triangle(const Bare_point& /*p*/, const Triangle& t) const
+  Triangle get_incident_triangle(const Facet& f, const Vertex_handle) const
   {
-    return t;
+    return triangle(f);
   }
 
   void set_point(const Vertex_handle v,

Mesh_3/test/Mesh_3/test_meshing_polyhedron_with_features.cpp

@@ -26,7 +26,7 @@
 #include <sstream>
 
 static constexpr bool verbose =
-#if CGAL_MESH_3_VERBOSE
+#ifdef CGAL_MESH_3_VERBOSE
   true;
 #else
   false;

Periodic_3_mesh_3/doc/Periodic_3_mesh_3/CGAL/Implicit_to_labeled_subdomains_function_wrapper.h

@@ -13,7 +13,7 @@ the subdomain index `2`.
 Note that for the 3D mesh generator [`f(p)=0`] means that p is outside the domain.
 Since this wrapper has values into `{1, 2}`, both the interior and the exterior of
 the periodic domain described by the input implicit function are meshed,
-thus yielding a periodic mesh of the entire canonical cube.
+thus yielding a periodic mesh of the entire canonical cuboid.
 
 \tparam Function provides the definition of the function.
         This parameter stands for a model of the concept `ImplicitFunction`

Periodic_3_mesh_3/doc/Periodic_3_mesh_3/CGAL/Periodic_3_function_wrapper.h

@@ -5,50 +5,50 @@ namespace CGAL {
 
 The class `Periodic_3_function_wrapper` is a helper class designed to wrap
 an (a priori non-periodic) implicit function describing a domain through the relationship
-[`p` is inside if `f(p)<0`] and defined over the canonical cube to a function
-defined over the whole Euclidean space and periodic, with the same period as the canonical cube.
+[`p` is inside if `f(p)<0`] and defined over the canonical cuboid to a function
+defined over the whole Euclidean space and periodic, with the same period as the canonical cuboid.
 
 More precisely, if `f` is the real function defined either over \f$ \mathbb R^3\f$
-or over the canonical cube, we construct the periodic real function \f$ f^{\ast} \f$ defined over
+or over the canonical cuboid, we construct the periodic real function \f$ f^{\ast} \f$ defined over
 \f$ \mathbb R^3\f$ as follows:
-- For any point \f$(x,y,z)\f$ in the canonical cube, \f$ f^{\ast}(x,y,z) = f(x,y,z)\f$
-- For any point \f$(x,y,z)\f$ outside the canonical cube, there is a unique canonical representative
-\f$(x_0,y_0,z_0)\f$ of \f$(x,y,z)\f$ in the canonical cube, i.e.,
+- For any point \f$(x,y,z)\f$ in the canonical cuboid, \f$ f^{\ast}(x,y,z) = f(x,y,z)\f$
+- For any point \f$(x,y,z)\f$ outside the canonical cuboid, there is a unique canonical representative
+\f$(x_0,y_0,z_0)\f$ of \f$(x,y,z)\f$ in the canonical cuboid, i.e.,
 \f$(x,y,z)=(x_0 + k_x c, y_0 + k_y c, z_0 + k_z c)\f$ with \f$(k_x,k_y,k_z)\f$ in \f$ \mathbb Z^3\f$,
 and \f$ f^{\ast}(x,y,z) = f(x_0,y_0,z_0) \f$.
 
-For example, if considering the unit cube as canonical cube, an oracle answering a
+For example, if considering the unit cube as canonical cuboid, an oracle answering a
 query such as <I>"what is the value of the implicit function at this point?"</I>
 at the point `(2.5, 2.5, 2.5)` will be in fact evaluated at the canonical representative, that is
 `(0.5, 0.5, 0.5)`.
 Consequently, it is then not required to provide an input domain that is defined over the whole
-space or periodic, but only defined over the canonical cube.
+space or periodic, but only defined over the canonical cuboid.
 
 \cgalFigureBegin{Periodic_3_mesh_3FromCanonicalToWhole, periodicity_base.svg}
 Illustration in 2D (cut view) of a domain defined by an implicit function that is transformed
 into a periodic implicit function.
-Only the values of the implicit function that are in the canonical cube are used:
+Only the values of the implicit function that are in the canonical cuboid are used:
 the values of the implicit function at \f$ P \f$ and \f$ Q \f$ are obtained by evaluating
 instead at \f$ P' \f$ and \f$ Q' \f$, as shown on the right.
 \cgalFigureEnd
 
 In practice, the implicit function provided by the user is likely defined
-over a larger domain than the canonical cube (in general, it is \f$ \mathbb R^3\f$).
+over a larger domain than the canonical cuboid (in general, it is \f$ \mathbb R^3\f$).
 Note that -- when constructing artificially periodic functions -- all the values of the implicit function
-for points outside this canonical cube are unused since queries are always answered by looking at the canonical representative.
+for points outside this canonical cuboid are unused since queries are always answered by looking at the canonical representative.
 \cgalFigureRef{Periodic_3_mesh_3FromCanonicalToWholeDiscard} gives an example of such domain where some information is discarded.
 
 \cgalFigureBegin{Periodic_3_mesh_3FromCanonicalToWholeDiscard, periodicity.svg}
 Illustration in 2D (cut view) of a domain defined by an implicit function artificially made periodic.
-Any value of the function outside of the canonical cube is ignored.
+Any value of the function outside of the canonical cuboid is ignored.
 \cgalFigureEnd
 
 Note also that when constructing artificially periodic functions, it is the responsibility of the user
-to provide an input function that is compatible with the canonical cube (that is, whose isovalues
+to provide an input function that is compatible with the canonical cuboid (that is, whose isovalues
 are <em>periodically</em> continuous and without intersections).
 \cgalFigureRef{Periodic_3_mesh_3ContinuityIssue} is an example of a bad choice
-of input function and canonical cube: there is no continuity of the isovalues
-at the border of the canonical cube. In such configuration, the mesher might
+of input function and canonical cuboid: there is no continuity of the isovalues
+at the border of the canonical cuboid. In such configuration, the mesher might
 or might not finish and the result is likely to be non-manifold and to contain self-intersections.
 
 \cgalFigureBegin{Periodic_3_mesh_3ContinuityIssue, periodicity_issue.svg}
@@ -82,7 +82,7 @@ public:
   /// \name Creation
   /// @{
   /*!
-   * \brief Construction from an implicit function and the canonical cube.
+   * \brief Construction from an implicit function and the canonical cuboid.
    */
   Periodic_3_function_wrapper(Function f, const Iso_cuboid_3& domain);
   /// @}

Periodic_3_mesh_3/doc/Periodic_3_mesh_3/Concepts/Periodic_3MeshDomain_3.h

@@ -12,7 +12,7 @@ domain is defined over the three-dimensional flat torus.
 From a syntactic point of view, it defines almost the same requirements
 as the concept `MeshDomain_3` and thus `Periodic_3MeshDomain_3` refines `MeshDomain_3`:
 the concept `Periodic_3MeshDomain_3` additionally requires an access to the user-defined
-canonical cube via the function `bounding_box`.
+canonical cuboid via the function `bounding_box`.
 However, the oracle must take into account the periodicity of the domain (see Section
 \ref Periodic_3_mesh_3InputDomain).
 
@@ -32,12 +32,12 @@ class Periodic_3MeshDomain_3
 {
 public:
   /*!
-  The canonical cube type.
+  The canonical cuboid type.
   */
   typedef unspecified_type Iso_cuboid_3;
 
   /*!
-  returns the user-chosen cube that is the canonical instance of the flat torus.
+  returns the user-chosen cuboid that is the canonical instance of the flat torus.
   */
   const Iso_cuboid_3& bounding_box();
 

Periodic_3_mesh_3/doc/Periodic_3_mesh_3/PackageDescription.txt

@@ -50,9 +50,9 @@ at the beginning of the meshing process, thus ensuring that the underlying perio
 triangulation forms at all times a simplicial complex within a single copy of the periodic space \f$ \mathbb T_c^3\f$
 (see Sections \ref P3Triangulation3secspace and \ref P3Triangulation3secintro
 of the package \ref PkgPeriodic3Triangulation3).
-By identifying a single copy of the flat torus \f$ \mathbb T_c^3\f$ (where `c`
-denotes the period) with a cube of side `c` in \f$ \mathbb R^3\f$, the meshing process
-can be exclusively conducted within a cube of side `c`.
+By identifying a single copy of the flat torus \f$ \mathbb T_{c_x,c_y,c_y}^3\f$ (where `c_x`,`c_y`,
+and `c_z` denote the period along the canonical axes) with a cuboid of side lengths `c_x`,`c_y`,
+and`c_z` in \f$ \mathbb R^3\f$, the meshing process can be exclusively conducted within a cuboid.
 The mesh within a single copy is created using \cgal's \ref PkgMesh3 package, but
 because \cgal's \ref PkgMesh3 package aims to mesh traditional (non-periodic)
 domains, an interface is necessary between \cgal's \ref PkgMesh3 package
@@ -80,7 +80,7 @@ and \cgal's three-dimensional mesh generator:
 The following class allows to construct a periodic implicit function from an implicit function that is not periodic:
 - `CGAL::Periodic_3_function_wrapper<Function,BGT>`
 
-The following class allows to split the canonical cube in two subdomains,
+The following class allows to split the canonical cuboid in two subdomains,
 separated by the zero-level of an implicit function:
 - `CGAL::Implicit_to_labeled_subdomains_function_wrapper<Function,BGT>`
 

Periodic_3_mesh_3/doc/Periodic_3_mesh_3/Periodic_3_mesh_3.txt

@@ -24,7 +24,7 @@ The domain to be meshed is a subset of the three-dimensional flat torus
 The domain may be connected or composed of multiple components
 and/or subdivided in several subdomains.
 The current implementation provides classes to represent
-domains bounded by isosurfaces of implicit functions defined over a cube.
+domains bounded by isosurfaces of implicit functions defined over a cuboid.
 
 Boundary and subdivision surfaces are either
 smooth or piecewise-smooth surfaces, formed with planar or curved surface patches.
@@ -66,16 +66,20 @@ meshes discretizing (non-periodic) 3D domains and to the \ref PkgPeriodic3Triang
 which are used as underlying triangulation structures of the mesh.
 
 A periodic mesh extends, by definition, infinitely in space. We consider the flat torus \f$ \mathbb T_c^3\f$,
-whose <I>canonical cube</I> has side length `c` (this canonical cube is named <I>original domain</I>
-in Chapter \ref PkgPeriodic3Triangulation3; we rename it here to avoid the confusion
-with the domain defined in Chapter \ref PkgMesh3). Well-chosen "dummy" points
+whose <I>canonical cuboid</I> has side lengths `(c_x, c_y, c_z)`. This canonical cuboid is named
+<I>original domain</I> in Chapter \ref PkgPeriodic3Triangulation3; we rename it here to avoid the confusion
+with the domain defined in Chapter \ref PkgMesh3. Well-chosen "dummy" points
 are inserted at the beginning of the meshing process, ensuring that the projection
 of the periodic triangulation into the flat torus \f$ \mathbb T_c^3\f$ forms at all times a simplicial complex
 (see Sections \ref P3Triangulation3secspace and \ref P3Triangulation3secintro
 of the manual of 3D periodic triangulations).
-Thanks to this construction, the meshing process can be exclusively conducted
-within the canonical cube.
-The mesh can then be created using the \ref PkgMesh3 package of \cgal.
+
+\warning These dummy points are spread on a grid whose density depends on the shortest side
+of the canonical cuboid. If the anisotropy of the canonical cuboid is too important, the grid may be
+very dense, and the meshing process may be slow and memory consuming.
+
+Using dummy points, one can constrain the periodic mesh generation to be exclusively conducted
+within the canonical cuboid. The mesh can then be created using the \ref PkgMesh3 package of \cgal.
 As this package originally aims to mesh non-periodic domains of \f$ \mathbb R^3\f$, an interface
 is necessary between the packages \ref PkgMesh3 and \ref PkgPeriodic3Triangulation3.
 This package provides this interface.
@@ -149,23 +153,23 @@ the function has negative values.
 \subsubsection Periodic_3_mesh_3InputDomainPeriodicity Periodicity of the Input Domain
 
 As described in Section \ref Periodic_3_mesh_3Mesh_3, the periodic mesh is in fact
-constructed over a single cube of side `c` in \f$ \mathbb R^3\f$, the <I>canonical cube</I>
+constructed over a single cube of sides `(c_x, c_y, c_z)` in \f$ \mathbb R^3\f$, the <I>canonical cuboid</I>
 of the flat torus \f$ \mathbb T_c^3\f$. The origin (given by three coordinates \f$ \alpha\f$, \f$ \beta\f$,
-and \f$ \gamma\f$) of this cube and the period `c` are input parameters chosen
-by the user. The cube \f$ [\alpha,\alpha+c)\times[\beta,\beta+c)\times[\gamma,\gamma+c)\f$
+and \f$ \gamma\f$) of this cube and the periods `c_x`, `c_y` and `c_z` are input parameters chosen
+by the user. The cuboid \f$ [\alpha,\alpha+c_x)\times[\beta,\beta+c_y)\times[\gamma,\gamma+c_z)\f$
 contains exactly one representative of each element in \f$ \mathbb T_c^3\f$.
-Although the mesh is only constructed over the canonical cube, some of the oracles
-used during the generation of the mesh must sometimes be evaluated outside of the canonical cube.
+Although the mesh is only constructed over the canonical cuboid, some of the oracles
+used during the generation of the mesh must sometimes be evaluated outside of the canonical cuboid.
 The implicit function describing the domain to be meshed must thus be defined
-over the whole Euclidean space and be periodic, with a period compatible with the canonical cube.
+over the whole Euclidean space and be periodic, with a period compatible with the canonical cuboid.
 
 \subsubsection Periodic_3_mesh_3ArtificialPeriodicity Enforcing Domain Periodicity
 
 The specifications of the input implicit function described in the previous section are quite restrictive.
 To relax these requirements, this package also offers a wrapper class,
 `CGAL::Periodic_3_function_wrapper`, to artificially construct
-periodic functions compatible with the user-defined canonical cube,
-from the values of an implicit function over the canonical cube.
+periodic functions compatible with the user-defined canonical cuboid,
+from the values of an implicit function over the canonical cuboid.
 It is thus possible to construct periodic domains described by implicit functions
 that are not intrinsically periodic, for example a sphere (see \cgalFigureRef{Periodic_3_mesh_3Periodic_implicit_sphere})
 or a cone (see Section \ref Periodic_3_mesh_3ConeWithSharpFeatures).
@@ -187,7 +191,7 @@ whose dual Voronoi edges intersect the surface patch. Such mesh facets are calle
 are approximated by sequences of mesh edges and the 0-dimensional exposed features
 are represented by mesh vertices.
 
-It is possible to extract the facets of the complex (restricted to the canonical cube)
+It is possible to extract the facets of the complex (restricted to the canonical cuboid)
 as a `FaceGraph`, using the function `facets_in_complex_3_to_triangle_mesh()`.
 
 \subsection Periodic_3_mesh_3DelaunayRefinement Delaunay Refinement
@@ -456,7 +460,7 @@ Cut view (middle). Another cut is shown, using 8 copies (right).
 
 While the implicit function used in the previous example is defined and periodic
 over the complete space, it is also possible to consider non-periodic implicit functions
-defined entirely within the canonical cube (or over the whole space) by using the wrapper class
+defined entirely within the canonical cuboid (or over the whole space) by using the wrapper class
 `CGAL::Periodic_3_function_wrapper`. Values will then be periodically duplicated, creating
 a periodic function. For example, replacing the previous domain with the following
 non-periodic implicit function (a sphere):
@@ -482,7 +486,7 @@ will yield the mesh shown on \cgalFigureRef{Periodic_3_mesh_3Periodic_implicit_s
 <img src="periodic_implicit_spheres.png" style="max-width:70%;"/>
 </center>
 \cgalFigureCaptionBegin{Periodic_3_mesh_3Periodic_implicit_sphere}
-Periodic mesh of an implicit sphere that is entirely contained in the input cube,
+Periodic mesh of an implicit sphere that is entirely contained in the input cuboid,
 shown with 8 copies (left). A cut view along one of the axes (right).
 \cgalFigureCaptionEnd
 
@@ -502,7 +506,7 @@ and wrapping the function as follows:
 \code{.cpp}
 ...
 Function_wrapper wrapper(schwarz_p);
-Periodic_mesh_domain domain(wrapper, canonical_cube);
+Periodic_mesh_domain domain(wrapper, canonical_cuboid);
 ...
 \endcode
 
@@ -587,9 +591,9 @@ and corners that are specified by the user:
 - Two curves cannot intersect, except at a common endpoint.
 
 \warning For conveniency, curves and corners do not need to be restricted
-to the canonical cube, but users should be mindful that curves and corners will exist
+to the canonical cuboid, but users should be mindful that curves and corners will exist
 in all periodic copies and the requirements described above must be satisfied.
-For example, if considering the unit cube as canonical cube, it is not valid
+For example, if considering the unit cube as canonical cuboid, it is not valid
 to add the segment `(2,2,2)--(3,3,3)` as feature and the point `(1.5, 1.5, 1.5)`
 as corner: once the periodicity is taken in account, the point is actually the middle
 of the segment.
@@ -625,11 +629,11 @@ see Section \ref Periodic_3_mesh_3SubMultipleCopies for more details.)
 It is possible to prescribe features that will, due to periodicity, form a
 continuous polyline that extends infinitely in space.
 A simple example of such occurrence is the segment `(0,0,0) -- (1,1,1)` when
-considering the unit cube as canonical cube.
+considering the unit cube as canonical cuboid.
 \cgalFigureRef{Periodic_3_mesh_3Periodic_protection} shows an implicit function
 describing a square-based prism such that the axis of the prism is the (Oz) axis.
 (The prism is thus 'cut' into 4 pieces when considered within a single copy of the periodic space.)
-The canonical cube is the unit cube and the following polylines have been specified to protect edges:
+The canonical cuboid is the unit cube and the following polylines have been specified to protect edges:
 \code
   // These are four vertical edges (orthogonal to xOy), strictly in the domain.
   // They correspond to the four sharp edges of the prism.

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/CMakeLists.txt

@@ -20,11 +20,15 @@ create_single_source_cgal_program("mesh_implicit_multi_domain.cpp")
 create_single_source_cgal_program("mesh_implicit_shape_with_subdomains.cpp")
 create_single_source_cgal_program("mesh_implicit_shape_with_optimizers.cpp")
 create_single_source_cgal_program("mesh_implicit_shape_with_features.cpp")
+create_single_source_cgal_program("mesh_periodic_polyhedral_domain.cpp")
 
 foreach(
   target
-  mesh_implicit_shape mesh_implicit_multi_domain
-  mesh_implicit_shape_with_subdomains mesh_implicit_shape_with_optimizers
-  mesh_implicit_shape_with_features)
+  mesh_implicit_shape
+  mesh_implicit_multi_domain
+  mesh_implicit_shape_with_subdomains
+  mesh_implicit_shape_with_optimizers
+  mesh_implicit_shape_with_features
+  mesh_periodic_polyhedral_domain)
   target_link_libraries(${target} PUBLIC CGAL::Eigen3_support)
 endforeach()

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_multi_domain.cpp

@@ -78,11 +78,11 @@ int main(int argc, char** argv)
   Multi_domain_wrapper multi_domain_function(funcs, vps);
   Periodic_mesh_domain domain(multi_domain_function, canonical_cube);
 
-  Periodic_mesh_criteria criteria(params::facet_angle(30).
-                                          facet_size(0.04).
-                                          facet_distance(0.025).
-                                          cell_radius_edge_ratio(2.).
-                                          cell_size(0.04));
+  Periodic_mesh_criteria criteria(params::facet_angle(30)
+                                         .facet_size(0.04)
+                                         .facet_distance(0.025)
+                                         .cell_radius_edge_ratio(2.)
+                                         .cell_size(0.04));
 
   // Mesh generation
   C3t3 c3t3 = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria);

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape.cpp

@@ -47,25 +47,32 @@ FT schwarz_p(const Point& p)
 
 int main(int argc, char** argv)
 {
-  // 'int' because the 'schwarz_p' function is periodic over the domain only if
-  // the length of the side of the domain is an integer.
-  int domain_size = (argc > 1) ? atoi(argv[1]) : 1;
-  int number_of_copies_in_output = (argc > 2) ? atoi(argv[2]) : 4; // can be 1, 2, 4, or 8
+  // 'atoi' because the 'schwarz_p' function is periodic over the domain
+  // only if the length of the side of the domain is an integer.
+  const int x_span = (argc > 1) ? atoi(argv[1]) : 1;
+  const int y_span = (argc > 2) ? atoi(argv[2]) : x_span;
+  const int z_span = (argc > 3) ? atoi(argv[3]) : x_span;
+  const int min_span = (std::min)({x_span, y_span, z_span});
 
-  Iso_cuboid canonical_cube(0, 0, 0, domain_size, domain_size, domain_size);
+  const int number_of_copies_in_output = (argc > 4) ? atoi(argv[4]) : 4; // can be 1, 2, 4, or 8
 
-  Periodic_mesh_domain domain =
+  const Iso_cuboid canonical_cube(0, 0, 0, x_span, y_span, z_span);
+
+  const Periodic_mesh_domain domain =
     Periodic_mesh_domain::create_implicit_mesh_domain(schwarz_p, canonical_cube);
 
-  Periodic_mesh_criteria criteria(params::facet_angle(30).
-                                          facet_size(0.035 * domain_size).
-                                          facet_distance(0.025 * domain_size).
-                                          cell_radius_edge_ratio(2.).
-                                          cell_size(0.05));
+  Periodic_mesh_criteria criteria(params::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::cout << "Created mesh with " << c3t3.number_of_vertices_in_complex() << " vertices"
+                        << " and " <<  c3t3.number_of_cells_in_complex() << " cells" << std::endl;
+
   std::ofstream medit_file("output_implicit_shape.mesh");
   CGAL::IO::output_periodic_mesh_to_medit(medit_file, c3t3, number_of_copies_in_output);
 

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape_with_features.cpp

@@ -91,11 +91,11 @@ int main(int argc, char** argv)
       Periodic_function(cone_function, canonical_cube), canonical_cube);
 
   // Mesh criteria
-  Periodic_mesh_criteria criteria(params::edge_size(0.02 * domain_size).
-                                          facet_angle(0.05 * domain_size).
-                                          facet_size(0.02 * domain_size).
-                                          cell_radius_edge_ratio(2).
-                                          cell_size(0.5));
+  Periodic_mesh_criteria criteria(params::edge_size(0.02 * domain_size)
+                                         .facet_angle(30)
+                                         .facet_size(0.02 * domain_size)
+                                         .cell_radius_edge_ratio(2)
+                                         .cell_size(0.5 * domain_size));
 
   // Create the features that we want to preserve
   Polylines polylines;

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape_with_optimizers.cpp

@@ -62,11 +62,11 @@ int main(int argc, char** argv)
   Periodic_mesh_domain domain =
     Periodic_mesh_domain::create_implicit_mesh_domain(double_p, canonical_cube);
 
-  Periodic_mesh_criteria criteria(params::facet_angle(30).
-                                          facet_size(0.05 * domain_size).
-                                          facet_distance(0.025 * domain_size).
-                                          cell_radius_edge_ratio(2.).
-                                          cell_size(0.05));
+  Periodic_mesh_criteria criteria(params::facet_angle(30)
+                                         .facet_size(0.05 * domain_size)
+                                         .facet_distance(0.025 * domain_size)
+                                         .cell_radius_edge_ratio(2.)
+                                         .cell_size(0.05 * domain_size));
 
   // Mesh generation with optimizers
   C3t3 c3t3 = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria,

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_implicit_shape_with_subdomains.cpp

@@ -67,11 +67,11 @@ int main(int argc, char** argv)
   size.set_size(0.1, volume_dimension, domain.index_from_subdomain_index(2)); // exterior
   size.set_size(0.03, volume_dimension, domain.index_from_subdomain_index(1)); // interior
 
-  Periodic_mesh_criteria criteria(params::facet_angle(30.).
-                                          facet_size(0.05).
-                                          facet_distance(0.025).
-                                          cell_radius_edge_ratio(2.).
-                                          cell_size(size));
+  Periodic_mesh_criteria criteria(params::facet_angle(30.)
+                                         .facet_size(0.05)
+                                         .facet_distance(0.025)
+                                         .cell_radius_edge_ratio(2.)
+                                         .cell_size(size));
 
   // Mesh generation
   C3t3 c3t3 = CGAL::make_periodic_3_mesh_3<C3t3>(domain, criteria);

Periodic_3_mesh_3/examples/Periodic_3_mesh_3/mesh_periodic_polyhedral_domain.cpp

@@ -0,0 +1,159 @@
+#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/Polyhedral_mesh_domain_3.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::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
+namespace params = 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)
+{
+  using vertex_descriptor = boost::graph_traits<Mesh>::vertex_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));
+}
+
+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;
+  generate_periodic_diamond(CGAL::ORIGIN, x_span, y_span, z_span, sm);
+  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);
+
+  Periodic_mesh_criteria criteria(params::edge_size(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;
+}

Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/IO/File_medit.h

@@ -89,9 +89,14 @@ void output_to_medit(std::ostream& os,
                      const Vertex_index_property_map& vertex_pmap,
                      const Facet_index_property_map& facet_pmap,
                      const Cell_index_property_map& cell_pmap,
-                     const Facet_index_property_map_twice& facet_twice_pmap = Facet_index_property_map_twice(),
-                     const bool print_each_facet_twice = false)
+                     const Facet_index_property_map_twice& = Facet_index_property_map_twice())
 {
+#ifdef CGAL_MESH_3_IO_VERBOSE
+  std::cout << "Output to medit;\n"
+            << "\toccurrences = " << occurrence_count << "\n"
+            << "\tdistinguish_copies = " << distinguish_copies << std::endl;
+#endif
+
   // if occurrence_count equals:
   // "1" --> only draws a single instance of the domain.
   // "2" --> draws 2 occurrences of the domain, displaying an additional domain
@@ -110,11 +115,10 @@ void output_to_medit(std::ostream& os,
   typedef typename Tr::Weighted_point            Weighted_point;
 
   typedef typename C3T3::Vertex_handle           Vertex_handle;
+  typedef typename C3T3::Facet                   Facet;
   typedef typename C3T3::Cell_handle             Cell_handle;
 
   typedef typename Tr::Vertex_iterator           Vertex_iterator;
-  typedef typename C3T3::Facet_iterator          Facet_iterator;
-  typedef typename C3T3::Cell_iterator           Cell_iterator;
 
   typedef typename Tr::Offset                    Offset;
 
@@ -126,21 +130,21 @@ void output_to_medit(std::ostream& os,
   int Oz_rn = 1 + (((occurrence_count - 1) >> 2) & 1);
 
   int number_of_vertices = static_cast<int>(tr.number_of_vertices());
-  int number_of_facets = static_cast<int>(c3t3.number_of_facets());
-  int number_of_cells = static_cast<int>(c3t3.number_of_cells());
+  int number_of_facets = static_cast<int>(c3t3.number_of_facets_in_complex());
+  int number_of_cells = static_cast<int>(c3t3.number_of_cells_in_complex());
 
   // Hardcoded values can be passed here to force more copies
   // Ox_rn = 20; Oy_rn = 20; Oz_rn = 1;
 
   int occ_mult = Ox_rn * Oy_rn * Oz_rn;
 
-#ifdef CGAL_PERIODIC_3_MESH_3_VERBOSE
-  std::cerr << "Outputting mesh to medit... " << std::endl;
-  std::cerr << "occurrences over each axis: "
+#ifdef CGAL_MESH_3_IO_VERBOSE
+  std::cout << "Outputting mesh to medit... " << std::endl;
+  std::cout << "occurrences over each axis: "
             << Ox_rn << " " << Oy_rn << " " << Oz_rn << std::endl;
-  std::cerr << number_of_vertices << " vertices" << std::endl;
-  std::cerr << number_of_facets << " facets" << std::endl;
-  std::cerr << number_of_cells << " cells" << std::endl;
+  std::cout << number_of_vertices << " vertices" << std::endl;
+  std::cout << number_of_facets << " facets" << std::endl;
+  std::cout << number_of_cells << " cells" << std::endl;
 #endif
 
   os << std::setprecision(17);
@@ -155,7 +159,7 @@ void output_to_medit(std::ostream& os,
 
   // On each axis, we repeat n+1 times the point, where 'n' is the number of
   // instances of the mesh that will be printed over that axis. This is because
-  // a cell 'c' might have point(c,i) that is equal to v with an offset 2
+  // a cell 'c' might have point(c,i) that lives in the +1 (in x, y, or z) offset
 
   boost::unordered_map<Vertex_handle, int> V;
   int inum = 1; // '1' because medit ids start at 1
@@ -194,65 +198,45 @@ void output_to_medit(std::ostream& os,
     }
   }
 
-  int medit_number_of_triangles = occ_mult * number_of_facets;
-  if(print_each_facet_twice)
-    medit_number_of_triangles *= 2;
+  int medit_number_of_triangles = 2 * occ_mult * number_of_facets;
 
   os << "Triangles\n" << medit_number_of_triangles << std::endl;
-  for(int i=0; i<Oz_rn; ++i)
+  for(auto cit = c3t3.triangulation().cells_begin();
+           cit != c3t3.triangulation().cells_end(); ++cit)
   {
-    for(int j=0; j<Oy_rn; ++j)
+    Cell_handle c = cit;
+    for(int s=0; s<4; ++s)
     {
-      for(int k=0; k<Ox_rn; ++k)
+      if(!c3t3.is_in_complex(c, s))
+        continue;
+
+      // Offsets in x/y/z
+      for(int i=0; i<Oz_rn; ++i)
       {
-        const Offset off(k, j, i);
-        for(Facet_iterator fit = c3t3.facets_begin(); fit != c3t3.facets_end(); ++fit)
+        for(int j=0; j<Oy_rn; ++j)
         {
-          for(int l=0; l<4; ++l)
+          for(int k=0; k<Ox_rn; ++k)
           {
-            if(l == fit->second)
-              continue;
+            const Offset off(k, j, i);
 
-            Cell_handle c = fit->first;
-            Vertex_handle v = c->vertex(l);
-            const Offset combined_off = tr.combine_offsets(
-                                          off, tr.int_to_off(c->offset(l)));
-            const int vector_offset = combined_off.x() +
-                                      combined_off.y() * (Ox_rn + 1) +
-                                      combined_off.z() * (Ox_rn + 1) * (Oy_rn + 1);
-            const int id = vector_offset * number_of_vertices + V[v];
-            CGAL_assertion(1 <= id && id <= medit_number_of_vertices);
-            os << id << " ";
-          }
-
-          // For multiple copies, color to distinguish copies rather than to distinguish subdomains
-          if(!distinguish_copies || occ_mult == 1)
-            os << get(facet_pmap, *fit) << '\n';
-          else
-            os << 1 + k + 3*j + 9*i << '\n';
-
-          // Print triangle again if needed
-          if(print_each_facet_twice)
-          {
-            for(int l=0; l<4; ++l)
+            for(int vi=1; vi<4; ++vi) // vertices of the facet in complex
             {
-              if(l == fit->second)
-                continue;
-
-              Cell_handle c = fit->first;
-              Vertex_handle v = c->vertex(l);
-              const Offset combined_off = tr.combine_offsets(
-                                            off, tr.int_to_off(c->offset(l)));
+              const int pos = (s+vi) % 4;
+              const Offset combined_off = tr.combine_offsets(off, tr.int_to_off(c->offset(pos)));
               const int vector_offset = combined_off.x() +
                                         combined_off.y() * (Ox_rn + 1) +
                                         combined_off.z() * (Ox_rn + 1) * (Oy_rn + 1);
+
+              const Vertex_handle v = c->vertex(pos);
               const int id = vector_offset * number_of_vertices + V[v];
+
               CGAL_assertion(1 <= id && id <= medit_number_of_vertices);
               os << id << " ";
             }
 
+            // For multiple copies, color to distinguish copies rather than to distinguish subdomains
             if(!distinguish_copies || occ_mult == 1)
-              os << get(facet_twice_pmap, *fit) << '\n';
+              os << get(facet_pmap, Facet(c, s)) << '\n';
             else
               os << 1 + k + 3*j + 9*i << '\n';
           }
@@ -269,12 +253,12 @@ void output_to_medit(std::ostream& os,
       for(int k=0; k<Ox_rn; ++k)
       {
         const Offset off(k, j, i);
-        for(Cell_iterator cit = c3t3.cells_begin(); cit !=c3t3.cells_end(); ++cit)
+        for(auto cit = c3t3.cells_in_complex_begin();
+                 cit !=c3t3.cells_in_complex_end(); ++cit)
         {
           for(int l=0; l<4; ++l)
           {
-            const Offset combined_off = tr.combine_offsets(
-                                          off, tr.int_to_off(cit->offset(l)));
+            const Offset combined_off = tr.combine_offsets(off, tr.int_to_off(cit->offset(l)));
             const int vector_offset = combined_off.x() +
                                       combined_off.y() * (Ox_rn + 1) +
                                       combined_off.z() * (Ox_rn + 1) * (Oy_rn + 1);
@@ -297,18 +281,21 @@ void output_to_medit(std::ostream& os,
   os << "End" << std::endl;
 }
 
-template <class C3T3, bool rebind, bool no_patch>
+template <class C3T3, bool rebind>
 void output_to_medit(std::ostream& os,
                      const C3T3& c3t3,
                      const int occurrence_count,
                      const bool distinguish_copies)
 {
 #ifdef CGAL_MESH_3_IO_VERBOSE
-  std::cerr << "Output to medit:\n";
+  std::cout << "Output to medit:\n";
 #endif
+
   CGAL_precondition(c3t3.triangulation().is_1_cover());
 
-  typedef CGAL::SMDS_3::Medit_pmap_generator<C3T3, rebind, no_patch>  Generator;
+  // periodic meshes always print facets twice because the facet in complex
+  // might be on the boundary of the domain
+  typedef CGAL::SMDS_3::Medit_pmap_generator<C3T3, rebind, true>      Generator;
   typedef typename Generator::Cell_pmap                               Cell_pmap;
   typedef typename Generator::Facet_pmap                              Facet_pmap;
   typedef typename Generator::Facet_pmap_twice                        Facet_pmap_twice;
@@ -320,11 +307,10 @@ void output_to_medit(std::ostream& os,
   Vertex_pmap vertex_pmap(c3t3, cell_pmap, facet_pmap);
 
   Periodic_3_mesh_3::output_to_medit(os, c3t3, occurrence_count, distinguish_copies,
-                                     vertex_pmap, facet_pmap, cell_pmap, facet_pmap_twice,
-                                     Generator().print_twice());
+                                     vertex_pmap, facet_pmap, cell_pmap, facet_pmap_twice);
 
 #ifdef CGAL_MESH_3_IO_VERBOSE
-  std::cerr << "done.\n";
+  std::cout << "done.\n";
 #endif
 }
 
@@ -341,32 +327,19 @@ namespace IO {
  * \param distinguish_copies if set to `true`, each copy is assigned a unique color.
  *                           Otherwise, all domains are drawn with subdomain index-based colors.
  * \param rebind if set to `true`, labels of cells are rebinded into [1..nb_of_labels]
- * \param show_patches if set to `true`, patches are labeled with different labels than
- *                     cells. If `false`, each surface facet is written twice, using label of
- *                     each adjacent cell.
  */
 template <class C3T3>
 void output_periodic_mesh_to_medit(std::ostream& os,
                                    const C3T3& c3t3,
                                    const int occurrence_count = 8,
                                    const bool distinguish_copies = true,
-                                   bool rebind = false,
-                                   bool show_patches = false)
+                                   const bool rebind = false,
+                                   const bool /*show_patches*/ = false) // for backward compatibility
 {
   if(rebind)
-  {
-    if(show_patches)
-      Periodic_3_mesh_3::output_to_medit<C3T3, true, false>(os, c3t3, occurrence_count, distinguish_copies);
-    else
-      Periodic_3_mesh_3::output_to_medit<C3T3, true, true>(os, c3t3, occurrence_count, distinguish_copies);
-  }
+    Periodic_3_mesh_3::output_to_medit<C3T3, true>(os, c3t3, occurrence_count, distinguish_copies);
   else
-  {
-    if(show_patches)
-      Periodic_3_mesh_3::output_to_medit<C3T3, false, false>(os, c3t3, occurrence_count, distinguish_copies);
-    else
-      Periodic_3_mesh_3::output_to_medit<C3T3, false, true>(os, c3t3, occurrence_count, distinguish_copies);
-  }
+    Periodic_3_mesh_3::output_to_medit<C3T3, false>(os, c3t3, occurrence_count, distinguish_copies);
 }
 
 } // namespace IO

Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/Protect_edges_sizing_field.h

@@ -10,10 +10,6 @@
 //
 // Author(s)     : Stephane Tayeb, Laurent Rineau, Mael Rouxel-Labbé
 //
-//******************************************************************************
-// File Description :
-//******************************************************************************
-
 #ifndef CGAL_PERIODIC_3_MESH_3_PROTECT_EDGES_SIZING_FIELD_H
 #define CGAL_PERIODIC_3_MESH_3_PROTECT_EDGES_SIZING_FIELD_H
 
@@ -38,6 +34,9 @@
 #include <CGAL/Mesh_3/Protect_edges_sizing_field.h>
 #include <CGAL/SMDS_3/utilities.h>
 #include <CGAL/Mesh_3/Triangulation_helpers.h>
+#if CGAL_MESH_3_PROTECTION_DEBUG
+#  include <CGAL/Mesh_3/Dump_c3t3.h>
+#endif
 
 #include <CGAL/enum.h>
 #include <CGAL/STL_Extension/internal/Has_member_visited.h>
@@ -463,7 +462,16 @@ private:
       if(dim == 0) msg << "corner (";
       else msg << "point (";
       msg << p << ")";
+#if CGAL_MESH_3_PROTECTION_DEBUG & 4
+      CGAL_error_msg(([this, str = msg.str()]()
+                      {
+                        CGAL_USE(this);
+                        dump_c3t3(this->c3t3_, "dump-bug");
+                        return str.c_str();
+                      }()));
+#else
       CGAL_error_msg(msg.str().c_str());
+#endif
     }
     return s;
   }
@@ -1029,9 +1037,8 @@ get_maximum_weight(const Vertex_handle protection_vertex, const FT intended_weig
   if(max_possible_weight < minimal_weight_)
     max_possible_weight = minimal_weight_;
 
-  CGAL_assertion_code(const Weighted_point& pvwp = c3t3_.triangulation().point(protection_vertex);)
-  CGAL_assertion_code(const Bare_point& pvp =
-    c3t3_.triangulation().geom_traits().construct_point_3_object()(pvwp);)
+  CGAL_assertion_code(const Weighted_point pvwp = c3t3_.triangulation().point(protection_vertex);)
+  CGAL_assertion_code(const Bare_point pvp = c3t3_.triangulation().geom_traits().construct_point_3_object()(pvwp);)
   CGAL_assertion_code(const int dim = get_dimension(protection_vertex);)
   CGAL_assertion_code(const Index index = c3t3_.index(protection_vertex);)
   CGAL_assertion_code(const FT w_max = CGAL::square(query_size(pvp, dim, index));)
@@ -1671,7 +1678,7 @@ smart_insert_point(const Bare_point& p, Weight w, int dim, const Index& index,
   CGAL_assertion(nearest_vh != Vertex_handle());
   CGAL_assertion(tr.point(nearest_vh) != cwp(tr.canonicalize_point(p)));
 
-#if CGAL_MESH_3_PROTECTION_DEBUG & 2
+#if CGAL_MESH_3_PROTECTION_DEBUG & 16
   std::cerr << "Nearest power vertex of (" << p << ") is "
             << &*nearest_vh << " (" << c3t3_.triangulation().point(nearest_vh) << ") "
             << "at distance: " << sq_d << std::endl;
@@ -2091,6 +2098,7 @@ insert_balls(const Vertex_handle& vp,
   //   n = 2(d-sq) / (sp+sq)
   // =======================
 
+  const FT d_signF = static_cast<FT>(d_sign);
   int n = static_cast<int>(std::floor(FT(2)*(d-sq) / (sp+sq))+.5);
   // if(minimal_weight_ != 0 && n == 0) return;
 
@@ -2113,10 +2121,8 @@ insert_balls(const Vertex_handle& vp,
                                                 curve_index, d_sign)
                 << ")\n";
 #endif
-      const FT sgn = (d_sign == CGAL::POSITIVE) ? 1.
-                   : (d_sign == CGAL::NEGATIVE ? -1. : 0.);
       const Bare_point new_point =
-        domain_.construct_point_on_curve(vpp, curve_index, sgn * d / 2);
+        domain_.construct_point_on_curve(vpp, curve_index, d_signF * d / 2);
       const int dim = 1; // new_point is on edge
       const Index index = domain_.index_from_curve_index(curve_index);
       const FT point_weight = CGAL::square(size_(new_point, dim, index));
@@ -2162,7 +2168,6 @@ insert_balls(const Vertex_handle& vp,
   FT norm_step_size = dleft_frac * step_size;
 
   // Initial distance
-  FT d_signF = static_cast<FT>(d_sign);
   FT pt_dist = d_signF * norm_step_size;
   Vertex_handle prev = vp;
 
@@ -2189,8 +2194,18 @@ insert_balls(const Vertex_handle& vp,
   else
   {
     CGAL_assertion_code(using boost::math::float_prior);
-    CGAL_assertion(n==0 ||
-                   dleft_frac >= float_prior(float_prior(1.)));
+#if CGAL_MESH_3_PROTECTION_DEBUG & 4
+    CGAL_assertion_msg(n==0 ||
+                       dleft_frac >= float_prior(float_prior(1.)),
+                       ([this]()
+                        {
+                          CGAL_USE(this);
+                          dump_c3t3(this->c3t3_, "dump-bug");
+                          return "the sampling of protecting balls is not possible";
+                        }()));
+#else
+    CGAL_assertion(n==0 || dleft_frac >= float_prior(float_prior(1.)));
+#endif
   }
 
   // Launch balls
@@ -2248,7 +2263,7 @@ void
 Protect_edges_sizing_field<C3T3, MD, Sf>::
 refine_balls()
 {
-#if CGAL_MESH_3_PROTECTION_DEBUG & 4
+#if CGAL_MESH_3_PROTECTION_DEBUG & 8
   dump_c3t3(c3t3_, "dump-before-refine_balls");
   dump_c3t3_edges(c3t3_, "dump-before-refine_balls");
 #endif
@@ -2384,7 +2399,7 @@ refine_balls()
       }
     }
 
-#if CGAL_MESH_3_PROTECTION_DEBUG & 4
+#if CGAL_MESH_3_PROTECTION_DEBUG & 8
     dump_c3t3(c3t3_, "dump-before-check_and_repopulate_edges");
     dump_c3t3_edges(c3t3_, "dump-before-check_and_repopulate_edges");
 #endif

Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_3/config.h

@@ -20,6 +20,10 @@
   #define CGAL_MESH_3_VERBOSE
 #endif
 
+#ifdef CGAL_MESH_3_CONFIG_H
+  #error "The header '<CGAL/Periodic_3_mesh_3/config.h>' must be included before '<CGAL/Mesh_3/config.h>'"
+#endif
+
 #include <CGAL/Mesh_3/config.h>
 
 // Whether to remove dummy points or not during the protection of sharp features

Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h

@@ -30,6 +30,7 @@
 #include <CGAL/Periodic_3_regular_triangulation_3.h>
 
 // vertex and cell bases
+#include <CGAL/Triangulation_vertex_base_with_info_3.h> // to mark dummy vertices
 #include <CGAL/Mesh_vertex_base_3.h>
 #include <CGAL/Mesh_cell_base_3.h>
 
@@ -44,6 +45,7 @@
 #include <iostream>
 #include <iterator>
 #include <limits>
+#include <unordered_set>
 #include <utility>
 #include <vector>
 
@@ -51,7 +53,7 @@ namespace CGAL {
 
 /// This class currently provides an interface between the class
 /// `CGAL::Periodic_3_regular_triangulation_3` and the mesher `Mesh_3`.
-/// As periodic triangulations are parallelized, a lot of these functions will
+/// If periodic triangulations are parallelized, a lot of these functions will
 /// become obsolete.
 template<class Gt_, class Tds_>
 class Periodic_3_regular_triangulation_3_wrapper
@@ -140,12 +142,18 @@ public:
 
   /// Concurrency related
   template <typename Cell_handle>
-  bool try_lock_cell(const Cell_handle &, int = 0) const { return true; }
+  bool try_lock_cell(const Cell_handle &, int = 0) const
+  {
+    std::cerr << "ERROR: P3M3 does not yet support parallel execution" << std::endl;
+    CGAL_assertion(false);
+    return true;
+  }
 
   bool try_lock_and_get_incident_cells(Vertex_handle /* v */,
                                        std::vector<Cell_handle>& /* cells */) const
   {
-    std::cerr << "ERROR: implement try_lock_and_get_incident_cells()"<< std::endl;
+    std::cerr << "ERROR: P3M3 does not yet support parallel execution" << std::endl;
+    CGAL_assertion(false);
     return true;
   }
 
@@ -156,7 +164,7 @@ public:
   bool is_infinite(const Facet&) const { return false; }
   bool is_infinite(const Cell_handle) const { return false; }
   bool is_infinite(const Cell_handle, int) const { return false; }
-  bool is_infinite(const Cell_handle c, int i, int j) const;
+  bool is_infinite(const Cell_handle, int, int) const { return false; }
 
   Cell_handle infinite_cell() const
   {
@@ -180,9 +188,6 @@ public:
 
   void set_domain(const Iso_cuboid& domain)
   {
-    CGAL_precondition_msg(domain.xmax() - domain.xmin() == domain.ymax() - domain.ymin() &&
-                          domain.xmax() - domain.xmin() == domain.zmax() - domain.zmin(),
-                          "The fundamental domain must be a cube.");
     Base::set_domain(domain);
   }
 
@@ -191,13 +196,19 @@ public:
     return P3T3::internal::robust_canonicalize_point(p, geom_traits());
   }
 
-  // @todo it might be dangerous to call robust_canonicalize without also changing
+  // @fixme it might be dangerous to call robust_canonicalize() without also changing
   // <p, offset> = construct_periodic_point(p) (lack of consistency in the result)
   Weighted_point canonicalize_point(const Weighted_point& p) const
   {
     return P3T3::internal::robust_canonicalize_point(p, geom_traits());
   }
 
+  // 1-cover, so we can take a const&
+  const Weighted_point& point(const Vertex_handle v) const
+  {
+    return v->point();
+  }
+
   Triangle triangle(const Facet& f) const
   {
     Periodic_triangle ptri = periodic_triangle(f);
@@ -297,51 +308,52 @@ public:
     this->v_offsets.clear();
   }
 
-  FT compute_power_distance_to_power_sphere(const Cell_handle& c, const int i) const
+  FT compute_power_distance_to_power_sphere(const Cell_handle c, const int i) const
   {
     typename Geom_traits::Compute_power_distance_to_power_sphere_3 cr =
       geom_traits().compute_power_distance_to_power_sphere_3_object();
 
     Offset o_nb = this->neighbor_offset(c, i, c->neighbor(i));
-    Offset o_vt = this->get_offset(c->neighbor(i), c->neighbor(i)->index(c));
+    Offset o_vt = get_offset(c->neighbor(i), c->neighbor(i)->index(c));
 
-    const Weighted_point& wp0 = this->point(c->vertex(0)); // need the canonical point
-    const Weighted_point& wp1 = this->point(c->vertex(1));
-    const Weighted_point& wp2 = this->point(c->vertex(2));
-    const Weighted_point& wp3 = this->point(c->vertex(3));
-    const Weighted_point& wq = this->point(c->neighbor(i)->vertex(c->neighbor(i)->index(c)));
-    const Offset& op0 = this->get_offset(c, 0);
-    const Offset& op1 = this->get_offset(c, 1);
-    const Offset& op2 = this->get_offset(c, 2);
-    const Offset& op3 = this->get_offset(c, 3);
-    const Offset& oq = o_vt - o_nb;
+    const Weighted_point& wp0 = point(c->vertex(0)); // need the canonical point
+    const Weighted_point& wp1 = point(c->vertex(1));
+    const Weighted_point& wp2 = point(c->vertex(2));
+    const Weighted_point& wp3 = point(c->vertex(3));
+    const Weighted_point& wq = point(c->neighbor(i)->vertex(c->neighbor(i)->index(c)));
+    const Offset& op0 = get_offset(c, 0);
+    const Offset& op1 = get_offset(c, 1);
+    const Offset& op2 = get_offset(c, 2);
+    const Offset& op3 = get_offset(c, 3);
+    const Offset oq = o_vt - o_nb;
 
     return cr(wp0, wp1, wp2, wp3, wq, op0, op1, op2, op3, oq);
   }
 
-  // The functions below are needed by Mesh_3 but need a specific implementation
-  // for the periodic case because we need to try with different offsets
-  // to get a result
-  FT compute_power_distance_to_power_sphere(const Cell_handle& c,
+  // The functions below are used in Mesh_3 and need a specific implementation
+  // for the periodic case because we need to try with different offsets to get the result
+  FT compute_power_distance_to_power_sphere(const Cell_handle c,
                                             const Vertex_handle v) const
   {
+    // @fixme need to introduce Compare_power_distances_to_power_sphere_3(4 points, query)
     typename Geom_traits::Compute_power_distance_to_power_sphere_3 cr =
       geom_traits().compute_power_distance_to_power_sphere_3_object();
 
     FT min_power_dist = std::numeric_limits<FT>::infinity();
 
+    const Weighted_point& wp0 = point(c->vertex(0)); // need the canonical point
+    const Weighted_point& wp1 = point(c->vertex(1));
+    const Weighted_point& wp2 = point(c->vertex(2));
+    const Weighted_point& wp3 = point(c->vertex(3));
+    const Weighted_point& wq = point(v);
+    const Offset& op0 = get_offset(c, 0);
+    const Offset& op1 = get_offset(c, 1);
+    const Offset& op2 = get_offset(c, 2);
+    const Offset& op3 = get_offset(c, 3);
+
     for(int i = 0; i < 3; ++i) {
       for(int j = 0; j < 3; ++j) {
         for(int k = 0; k < 3; ++k) {
-          const Weighted_point& wp0 = this->point(c->vertex(0)); // need the canonical point
-          const Weighted_point& wp1 = this->point(c->vertex(1));
-          const Weighted_point& wp2 = this->point(c->vertex(2));
-          const Weighted_point& wp3 = this->point(c->vertex(3));
-          const Weighted_point& wq = this->point(v);
-          const Offset& op0 = this->get_offset(c, 0);
-          const Offset& op1 = this->get_offset(c, 1);
-          const Offset& op2 = this->get_offset(c, 2);
-          const Offset& op3 = this->get_offset(c, 3);
           const Offset oq(i-1, j-1, k-1);
 
           FT power_dist = cr(wp0, wp1, wp2, wp3, wq, op0, op1, op2, op3, oq);
@@ -371,11 +383,14 @@ public:
           const Offset off(i-1, j-1, k-1);
           if(tester(c, off))
           {
-            return construct_tetrahedron(
-                       canonic_wp, this->point(c->vertex((index+1)&3)),
-                       this->point(c->vertex((index+2)&3)), this->point(c->vertex((index+3)&3)),
-                       off, this->get_offset(c, (index+1)&3),
-                       this->get_offset(c, (index+2)&3), this->get_offset(c, (index+3)&3));
+            return construct_tetrahedron(canonic_wp,
+                                         point(c->vertex((index+1)&3)),
+                                         point(c->vertex((index+2)&3)),
+                                         point(c->vertex((index+3)&3)),
+                                         off,
+                                         get_offset(c, (index+1)&3),
+                                         get_offset(c, (index+2)&3),
+                                         get_offset(c, (index+3)&3));
           }
         }
       }
@@ -385,7 +400,7 @@ public:
     return Tetrahedron();
   }
 
-  Bounded_side side_of_power_sphere(const Cell_handle& c, const Weighted_point& p,
+  Bounded_side side_of_power_sphere(const Cell_handle c, const Weighted_point& p,
                                     bool perturb = false) const
   {
     Weighted_point canonical_p = canonicalize_point(p);
@@ -404,62 +419,85 @@ public:
     }
 
     return bs;
-
-    return Base::side_of_power_sphere(c, canonical_p, Offset(), perturb);
   }
 
-  // Warning : This is a periodic version that computes the smallest possible
+  // Warning: This is a periodic version that computes the smallest possible distance
   // between 'p' and 'q', for all possible combinations of offsets
   FT min_squared_distance(const Bare_point& p, const Bare_point& q) const
   {
-    typename Geom_traits::Compute_squared_distance_3 csd =
+    typename Geom_traits::Compare_squared_distance_3 compare_sd =
+      geom_traits().compare_squared_distance_3_object();
+    typename Geom_traits::Compute_squared_distance_3 compute_sd =
       geom_traits().compute_squared_distance_3_object();
 
-    const Bare_point cp = canonicalize_point(p);
-    const Bare_point cq = canonicalize_point(q);
+    bool used_exact = false;
+    std::pair<Bare_point, Offset> pp_p = P3T3::internal::construct_periodic_point(p, used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(q, used_exact, geom_traits());
 
-    FT min_sq_dist = std::numeric_limits<FT>::infinity();
+    Offset min_off;
 
     for(int i = 0; i < 3; ++i) {
       for(int j = 0; j < 3; ++j) {
-        for(int k = 0; k < 3; ++k) {
-          FT sq_dist = csd(cq, construct_point(std::make_pair(cp, Offset(i-1, j-1, k-1))));
+        for(int k = 0; k < 3; ++k)
+        {
+          const Offset o(i-1, j-1, k-1);
 
-          if(sq_dist < min_sq_dist)
-            min_sq_dist = sq_dist;
-        }
-      }
-    }
-
-    return min_sq_dist;
-  }
-
-  // Warning : This function finds which offset 'Oq' should be applied to 'q' so
-  // that the distance between 'p' and '(q, Oq)' is minimal.
-  //
-  // \pre 'p' lives in the canonical instance.
-  Bare_point get_closest_point(const Bare_point& p, const Bare_point& q) const
-  {
-    Bare_point rq;
-    const Bare_point cq = canonicalize_point(q);
-    FT min_sq_dist = std::numeric_limits<FT>::infinity();
-
-    for(int i = -1; i < 2; ++i) {
-      for(int j = -1; j < 2; ++j) {
-        for(int k = -1; k < 2; ++k) {
-          const Bare_point tcq = construct_point(std::make_pair(cq, Offset(i, j, k)));
-          FT sq_dist = geom_traits().compute_squared_distance_3_object()(p, tcq);
-
-          if(sq_dist < min_sq_dist)
+          if((i == 0 && j == 0 && k == 0) ||
+              compare_sd(q, p, q, p,
+                         pp_q.second, pp_p.second + o,
+                         pp_q.second, pp_p.second + min_off) == SMALLER)
           {
-            rq = tcq;
-            min_sq_dist = sq_dist;
+            min_off = o;
           }
         }
       }
     }
 
-    return rq;
+    return compute_sd(q, p, pp_q.second, pp_p.second + min_off);
+  }
+
+  // Warning: This function finds which offset 'Oq' should be applied to 'q' such
+  // that the distance between 'p' and '(q, Oq)' is minimal.
+  //
+  // \pre 'p' lives in the canonical instance.
+  Bare_point get_closest_point(const Bare_point& p, const Bare_point& q) const
+  {
+    CGAL_precondition(p.x() < domain().xmax());
+    CGAL_precondition(p.y() < domain().ymax());
+    CGAL_precondition(p.z() < domain().zmax());
+    CGAL_precondition(p.x() >= domain().xmin());
+    CGAL_precondition(p.y() >= domain().ymin());
+    CGAL_precondition(p.z() >= domain().zmin());
+
+    typename Geom_traits::Compare_squared_distance_3 compare_sd =
+      geom_traits().compare_squared_distance_3_object();
+    typename Geom_traits::Construct_point_3 cp =
+      geom_traits().construct_point_3_object();
+
+    bool used_exact = false;
+    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(q, used_exact, geom_traits());
+
+    Offset min_off;
+    Offset null_offset(0,0,0);
+
+    for(int i = 0; i < 3; ++i) {
+      for(int j = 0; j < 3; ++j) {
+        for(int k = 0; k < 3; ++k)
+        {
+          const Offset o(i-1, j-1, k-1);
+
+          if((i == 0 && j == 0 && k == 0) ||
+             compare_sd(p, q, p, q,
+                        null_offset, pp_q.second + o,
+                        null_offset, pp_q.second + min_off) == SMALLER)
+          {
+            min_off = o;
+          }
+        }
+      }
+    }
+
+    return cp(q, pp_q.second + min_off);
   }
 
   Weighted_point get_closest_point(const Weighted_point& wp, const Weighted_point& wq) const
@@ -471,46 +509,28 @@ public:
     return cwp(get_closest_point(cp(wp), cp(wq)), cw(wq));
   }
 
-  Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const
+  // returns the triangle corresponding to f, with a geometric shift
+  // so that it is incident to ref_v's canonical position
+  Triangle get_incident_triangle(const Facet& f, const Vertex_handle ref_v) const
   {
-    typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_3_object();
+    typename Geom_traits::Construct_point_3 cp = geom_traits().construct_point_3_object();
     typename Geom_traits::Construct_translated_point_3 tr = geom_traits().construct_translated_point_3_object();
-    typename Geom_traits::Compute_squared_distance_3 csd = geom_traits().compute_squared_distance_3_object();
+    typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_3_object();
+    typename Geom_traits::Construct_triangle_3 ct = geom_traits().construct_triangle_3_object();
 
-    // It doesn't matter which point we use to canonicalize the triangle as P3M3 is necessarily
-    // in one cover and we have to look at all the neighboring copies anyway since we do not
-    // have control of 'p'.
-    Bare_point canon_p0 = canonicalize_point(t[0]);
-    Vector_3 move_to_canonical = cv(t[0], canon_p0);
-    const std::array<Bare_point, 3> ct = { canon_p0,
-                                           tr(t[1], move_to_canonical),
-                                           tr(t[2], move_to_canonical) };
+    CGAL_precondition(f.first != Cell_handle() && f.first->has_vertex(ref_v));
+    const int ref_v_pos = f.first->index(ref_v);
+    const Bare_point& ref_p = cp(point(ref_v));
+    const Bare_point ref_p_in_f = cp(point(f.first, ref_v_pos));
+    Vector_3 move_to_canonical = cv(ref_p_in_f, ref_p);
 
-    FT min_sq_dist = std::numeric_limits<FT>::infinity();
-    Triangle rt;
-    for(int i = -1; i < 2; ++i) {
-      for(int j = -1; j < 2; ++j) {
-        for(int k = -1; k < 2; ++k) {
+    const int s = f.second;
+    const Bare_point mp0 = tr(cp(point(f.first, (s+1)%4)), move_to_canonical);
+    const Bare_point mp1 = tr(cp(point(f.first, (s+2)%4)), move_to_canonical);
+    const Bare_point mp2 = tr(cp(point(f.first, (s+3)%4)), move_to_canonical);
+    const Triangle t = ct(mp0, mp1, mp2);
 
-          const Triangle tt(
-            construct_point(std::make_pair(ct[0], Offset(i, j, k))),
-            construct_point(std::make_pair(ct[1], Offset(i, j, k))),
-            construct_point(std::make_pair(ct[2], Offset(i, j, k))));
-
-          const FT sq_dist = csd(p, tt);
-
-          if(sq_dist == FT(0))
-            return rt;
-
-          if(sq_dist < min_sq_dist) {
-            rt = tt;
-            min_sq_dist = sq_dist;
-          }
-        }
-      }
-    }
-
-    return rt;
+    return t;
   }
 
   // Warning: This is a periodic version that computes the smallest possible
@@ -522,37 +542,51 @@ public:
                                        const Weighted_point& q,
                                        const Weighted_point& r) const
   {
-    typename Geom_traits::Compute_power_product_3 power_distance =
-      geom_traits().compute_power_product_3_object();
+    CGAL_precondition(this->is_1_cover());
 
-    // canonicalize the points
-    const Weighted_point cp =
-      geom_traits().construct_weighted_point_3_object()(canonicalize_point(p));
-    const Weighted_point cq = canonicalize_point(q);
-    const Weighted_point cr = canonicalize_point(r);
+    typename Geom_traits::Construct_point_3 cp =
+      geom_traits().construct_point_3_object();
+    typename Geom_traits::Compare_power_distance_3 compare_power_distance =
+      geom_traits().compare_power_distance_3_object();
 
-    FT min_power_distance_to_q = std::numeric_limits<FT>::infinity();
-    FT min_power_distance_to_r = std::numeric_limits<FT>::infinity();
+    // Compute the offsets that would bring p, q, and r into the canonical domain
+    bool used_exact = false;
+    std::pair<Bare_point, Offset> pp_p = P3T3::internal::construct_periodic_point(p, used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_q = P3T3::internal::construct_periodic_point(cp(q), used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_r = P3T3::internal::construct_periodic_point(cp(r), used_exact, geom_traits());
 
-    for(int i = 0; i < 3; ++i) {
-      for(int j = 0; j < 3; ++j) {
-        for(int k = 0; k < 3; ++k) {
-          const Weighted_point cp_copy =
-            construct_weighted_point(std::make_pair(cp, Offset(i-1, j-1, k-1)));
-
-          const FT power_distance_to_q = power_distance(cp_copy, cq);
-          if(power_distance_to_q < min_power_distance_to_q)
-            min_power_distance_to_q = power_distance_to_q;
-
-          const FT power_distance_to_r = power_distance(cp_copy, cr);
-          if (power_distance_to_r < min_power_distance_to_r)
-            min_power_distance_to_r = power_distance_to_r;
+    // To compare pp(p, q) to pp(p, r), we first need to know the best offsets that minimize these distances
+    auto get_offset_minimizing_power_product = [&](const Weighted_point& wp,
+                                                   const Offset& base_wp_offset) -> Offset
+    {
+      Offset min_wp_offset;
+      for(int i = 0; i < 3; ++i) {
+        for(int j = 0; j < 3; ++j) {
+          for(int k = 0; k < 3; ++k)
+          {
+            const Offset off(i-1, j-1, k-1);
+            if((i == 0 && j == 0 && k == 0) ||
+               compare_power_distance(p, wp, wp,
+                                      pp_p.second,
+                                      base_wp_offset + off,
+                                      base_wp_offset + min_wp_offset) == SMALLER)
+            {
+              min_wp_offset = off;
+            }
+          }
         }
       }
-    }
 
-    CGAL_postcondition(min_power_distance_to_r < 0.5 && min_power_distance_to_q < 0.5);
-    return min_power_distance_to_q >= min_power_distance_to_r;
+      return min_wp_offset;
+    };
+
+    Offset min_q_off = get_offset_minimizing_power_product(q, pp_q.second);
+    Offset min_r_off = get_offset_minimizing_power_product(r, pp_r.second);
+
+    return !(compare_power_distance(p, q, r,
+                                    pp_p.second,
+                                    pp_q.second + min_q_off,
+                                    pp_r.second + min_r_off) == SMALLER);
   }
 
   /// \name Locate functions
@@ -574,8 +608,8 @@ public:
     return Base::nearest_power_vertex(canonicalize_point(p), start);
   }
 
-  /// Return the squared distance (note: _NOT_ the power distance) between the
-  /// 'p' and the closest vertex for the power distance.
+  /// Return the squared distance (note: _NOT_ the power distance)
+  /// between 'p' and the closest vertex for the power distance.
   std::pair<Vertex_handle, FT>
   nearest_power_vertex_with_sq_distance(const Bare_point& p, Cell_handle start) const
   {
@@ -607,14 +641,14 @@ public:
     query_offset = this->combine_offsets(Offset(), -query_offset);
 
     Vertex_handle nearest = Base::nearest_vertex_in_cell(c, canonical_p, query_offset);
-    const Weighted_point& nearest_wp = this->point(nearest);
+    const Weighted_point& nearest_wp = point(nearest);
     Offset offset_of_nearest = Base::get_min_dist_offset(canonical_p, query_offset, nearest);
     FT min_sq_dist = csd(canonical_p, cp(nearest_wp), query_offset, offset_of_nearest);
 
     std::vector<Vertex_handle> vs;
     vs.reserve(32);
 
-    while(true)
+    for(;;)
     {
       Vertex_handle tmp = nearest;
 
@@ -622,7 +656,7 @@ public:
       for(typename std::vector<Vertex_handle>::const_iterator vsit = vs.begin();
                                                               vsit != vs.end(); ++vsit)
       {
-        // Can happen in 27-sheeted triangulations composed of few points
+        // Can happen in periodic triangulations composed of few points
         if(point(*vsit) == point(nearest))
           continue;
 
@@ -632,7 +666,7 @@ public:
         {
           tmp = *vsit;
           offset_of_nearest = min_dist_offset;
-          const Weighted_point& vswp = this->point(tmp);
+          const Weighted_point& vswp = point(tmp);
           min_sq_dist = csd(canonical_p, cp(vswp), query_offset, min_dist_offset);
         }
       }
@@ -647,6 +681,41 @@ public:
     return std::make_pair(nearest, min_sq_dist);
   }
 
+  std::pair<Vertex_handle, FT>
+  nearest_power_vertex_with_sq_distance(const Vertex_handle v) const
+  {
+    typename Geom_traits::Construct_point_3 cp = geom_traits().construct_point_3_object();
+    typename Geom_traits::Compute_squared_distance_3 csd = geom_traits().compute_squared_distance_3_object();
+
+    Vertex_handle min_v {};
+    FT min_sq_dist = -1;
+
+    std::vector<Cell_handle> inc_cells;
+    std::unordered_set<Vertex_handle> visited_vertices;
+    incident_cells(v, std::back_inserter(inc_cells));
+    for(Cell_handle ch : inc_cells)
+    {
+      CGAL_assertion(ch->has_vertex(v));
+      int v_pos = ch->index(v);
+      for(int i=1; i<4; ++i)
+      {
+        int vi_pos = (v_pos + i) % 4;
+        Vertex_handle vi = ch->vertex(vi_pos);
+        if(!visited_vertices.insert(vi).second) // already visited
+          continue;
+
+        FT sq_dist_i = csd(cp(point(ch, v_pos)), cp(point(ch, vi_pos)));
+        if(min_v == Vertex_handle() || sq_dist_i < min_sq_dist)
+        {
+          min_v = vi;
+          min_sq_dist = sq_dist_i;
+        }
+      }
+    }
+
+    return { min_v, min_sq_dist };
+  }
+
   Cell_handle locate(const Weighted_point& p,
                      Cell_handle start = Cell_handle(),
                      bool* CGAL_assertion_code(could_lock_zone) = nullptr) const
@@ -847,6 +916,76 @@ public:
     return make_object(s);
   }
 
+  void dual_exact(const Facet& f, const Weighted_point& ws,
+                  Bare_point& cc) const
+  {
+    // first find the offset minimizing the distance between the facet and the fourth point
+    typename Geom_traits::Construct_point_3 construct_point =
+      geom_traits().construct_point_3_object();
+
+    // @fixme need to introduce Compare_power_distances_to_power_sphere_3(3 points, query)
+    typename Geom_traits::Construct_weighted_circumcenter_3 wcc =
+      geom_traits().construct_weighted_circumcenter_3_object();
+    typename Geom_traits::Compare_squared_distance_3 compare_sd =
+      geom_traits().compare_squared_distance_3_object();
+
+    const Cell_handle c = f.first;
+    const int i = f.second;
+
+    const Bare_point fcc = wcc(point(c, (i+1)%4), point(c, (i+2)%4), point(c, (i+3)%4));
+    const Bare_point& s = construct_point(ws);
+
+    bool used_exact = false;
+    std::pair<Bare_point, Offset> pp_fcc = P3T3::internal::construct_periodic_point(fcc, used_exact, geom_traits());
+    std::pair<Bare_point, Offset> pp_s = P3T3::internal::construct_periodic_point(s, used_exact, geom_traits());
+
+    Offset min_off;
+
+    for(int i = 0; i < 3; ++i) {
+      for(int j = 0; j < 3; ++j) {
+        for(int k = 0; k < 3; ++k)
+        {
+          const Offset o(i-1, j-1, k-1);
+
+          if((i == 0 && j == 0 && k == 0) ||
+             compare_sd(fcc, s, fcc, s,
+                        pp_fcc.second, pp_s.second + o,
+                        pp_fcc.second, pp_s.second + min_off) == SMALLER)
+          {
+            min_off = o;
+          }
+        }
+      }
+    }
+
+    typedef typename Kernel_traits<Bare_point>::Kernel                Kernel;
+    typedef Exact_predicates_exact_constructions_kernel               EKernel;
+
+    typedef Cartesian_converter<Kernel, EKernel>                      To_exact;
+    typedef Cartesian_converter<EKernel, Kernel>                      Back_from_exact;
+
+    typedef CGAL::Periodic_3_regular_triangulation_traits_3<EKernel>  Exact_Rt;
+    typedef typename Exact_Rt::Weighted_point_3                       EWeighted_point_3;
+
+    To_exact to_exact;
+    Back_from_exact back_from_exact;
+
+    Exact_Rt etraits(to_exact(domain()));
+    Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter =
+      etraits.construct_weighted_circumcenter_3_object();
+
+    const EWeighted_point_3& cp = to_exact(c->vertex((i+1)%4)->point());
+    const EWeighted_point_3& cq = to_exact(c->vertex((i+2)%4)->point());
+    const EWeighted_point_3& cr = to_exact(c->vertex((i+3)%4)->point());
+    const EWeighted_point_3& cs = to_exact(ws);
+
+    cc = back_from_exact(exact_weighted_circumcenter(cp, cq, cr, cs,
+                                                     pp_fcc.second + get_offset(c, (i+1)%4),
+                                                     pp_fcc.second + get_offset(c, (i+2)%4),
+                                                     pp_fcc.second + get_offset(c, (i+3)%4),
+                                                     pp_s.second + min_off));
+  }
+
   void dual_segment(const Facet& facet, Bare_point& p, Bare_point& q) const
   {
     typename Base::Periodic_segment_3 ps = Base::dual(facet);
@@ -890,12 +1029,14 @@ public:
                            to_exact(point(n->vertex(2))), to_exact(point(n->vertex(3))),
                            get_offset(n, 0), get_offset(n, 1),
                            get_offset(n, 2), get_offset(n, 3));
+
     typename EKernel::Point_3 dp;
 
     // get the offset of the first weighted circumcenter
     Offset transl_wc1;
-    while(true) /* while not in */
+    for(;;) /* while not in */
     {
+      // can safely perform a construction here because the kernel has exact constructions
       dp = etraits.construct_point_3_object()(exact_wc1, transl_wc1);
 
       if(dp.x() < dom.xmin())
@@ -916,7 +1057,7 @@ public:
 
     // get the offset of the second weighted circumcenter
     Offset transl_wc2;
-    while(true) /* while not in */
+    for(;;) /* while not in */
     {
       dp = etraits.construct_point_3_object()(exact_wc2, transl_wc2);
 
@@ -942,6 +1083,7 @@ public:
     Offset cumm_off((std::min)(o1.x(), o2.x()),
                     (std::min)(o1.y(), o2.y()),
                     (std::min)(o1.z(), o2.z()));
+
     EPoint_3 ewc1 = exact_construct_point(exact_wc1, transl_wc1);
     EPoint_3 ewc2 = exact_construct_point(exact_wc2, transl_wc2);
     p = back_from_exact(exact_construct_point(ewc1, o1 - cumm_off));
@@ -984,6 +1126,13 @@ template<class MD,
          class Cell_base_ = Default>
 class Periodic_3_mesh_triangulation_3
 {
+  // Triangulation_vertex_base_with_info_3 only does default initialization
+  // and not value initialization, but we cannot initialize info() during Mesh_3's refinement
+  struct Boolean_with_def_value
+  {
+    bool is_dummy_vertex = false;
+  };
+
   // default K
   typedef typename Default::Get<K_, typename Kernel_traits<MD>::Kernel>::type K;
 
@@ -993,12 +1142,13 @@ class Periodic_3_mesh_triangulation_3
   // Periodic vertex and cell bases
   typedef Periodic_3_triangulation_ds_vertex_base_3<> VbDS;
   typedef Regular_triangulation_vertex_base_3<Geom_traits, VbDS> PVb;
+  typedef Triangulation_vertex_base_with_info_3<Boolean_with_def_value, Geom_traits, PVb> Vb;
 
   typedef Periodic_3_triangulation_ds_cell_base_3<> CbDS;
   typedef Regular_triangulation_cell_base_3<Geom_traits, CbDS> RCb;
   typedef Regular_triangulation_cell_base_with_weighted_circumcenter_3<Geom_traits, RCb> PCb;
 
-  typedef Mesh_vertex_base_3<Geom_traits, MD, PVb> Default_Vb;
+  typedef Mesh_vertex_base_3<Geom_traits, MD, Vb> Default_Vb;
   typedef Mesh_cell_base_3<Geom_traits, MD, PCb> Default_Cb;
 
   // default Vb/Cb

Periodic_3_mesh_3/include/CGAL/make_periodic_3_mesh_3.h

@@ -30,22 +30,61 @@
 #include <CGAL/Mesh_3/C3T3_helpers.h>
 #include <CGAL/Named_function_parameters.h>
 
+#include <boost/random/random_number_generator.hpp>
+#include <boost/random/linear_congruential.hpp>
+
 namespace CGAL {
 namespace Periodic_3_mesh_3 {
 namespace internal {
 
-template<typename C3T3>
-void mark_dummy_points(C3T3& c3t3)
+template<typename C3T3, typename MeshDomain>
+bool insert_dummy_points(C3T3& c3t3,
+                         const MeshDomain& domain)
 {
-  CGAL_precondition(c3t3.triangulation().is_1_cover());
+  typedef typename C3T3::Triangulation::Vertex_handle                         Vertex_handle;
 
-  typedef typename C3T3::Triangulation::Vertex_iterator       Vertex_iterator;
+  typename C3T3::Triangulation::Geom_traits::Construct_point_3 cp =
+    c3t3.triangulation().geom_traits().construct_point_3_object();
 
-  for(Vertex_iterator vit = c3t3.triangulation().vertices_begin();
-                      vit != c3t3.triangulation().vertices_end(); ++vit)
+  std::vector<Vertex_handle> dummy_vertices = c3t3.triangulation().insert_generic_dummy_points();
+  CGAL_postcondition(c3t3.triangulation().is_1_cover());
+
+  // Abuse the multi cover function to check if cells have a small-enough orthoradius
+  c3t3.triangulation().update_cover_data_after_converting_to_27_sheeted_covering();
+  if(!c3t3.triangulation().can_be_converted_to_1_sheet())
   {
-    c3t3.set_index(vit, 0);
+    std::cerr << "Error: dummy points do not create a 1-cover" << std::endl;
+    CGAL_postcondition(false);
+    return false;
   }
+
+  for(Vertex_handle dvh : dummy_vertices)
+  {
+    dvh->info().is_dummy_vertex = true;
+
+    // @fixme
+    // The real in-dimension of a dummy vertex is impossible to evaluate: it can fall
+    // on a surface or even a curve, but there is no way to know this because the concept (and models)
+    // of MeshDomain_3 do not have a `Patch_index Do_intersect()(Point_3)`.
+    //
+    // It seems that setting the dimension wrongly, to a too-low dimension does not cause
+    // problems, whereas on the other way is problematic: if setting it to '3' but the dummy vertex
+    // is part of the restricted Delaunay, then the refinement will refine until the vertex
+    // is no longer part of the restricted Delaunay, creating ugly spots in the mesh...
+    c3t3.set_dimension(dvh, 0);
+
+    // @fixme
+    // This should be setting patch indices if the dummy vertex is not in dim 3.
+    // Could maybe do it later, re-assigning indices of dummy vertices once the c3t3
+    // has been scanned, but still before refinement starts.
+    auto opt_si = domain.is_in_domain_object()(cp(c3t3.triangulation().point(dvh)));
+    if(opt_si.has_value())
+      c3t3.set_index(dvh, domain.index_from_subdomain_index(*opt_si));
+    else
+      c3t3.set_index(dvh, 0);
+  }
+
+  return true;
 }
 
 template <typename C3T3, typename MeshDomain, typename MeshCriteria>
@@ -60,7 +99,6 @@ void init_c3t3_with_features(C3T3& c3t3,
   CGAL::Periodic_3_mesh_3::Protect_edges_sizing_field<C3T3, MeshDomain, Sizing_field>
     protect_edges(c3t3, domain, Sizing_field(criteria.edge_criteria_object()));
   protect_edges.set_nonlinear_growth_of_balls(nonlinear);
-
   protect_edges(true);
 }
 
@@ -85,8 +123,10 @@ struct C3t3_initializer_base
                   const parameters::internal::Mesh_3_options& mesh_options)
   {
     c3t3.triangulation().set_domain(domain.bounding_box());
-    c3t3.triangulation().insert_dummy_points();
-    mark_dummy_points(c3t3);
+
+    // Enforce 1-cover by adding dummy points
+    if(!insert_dummy_points(c3t3, domain))
+      return;
 
     // Call the basic initialization from c3t3, which handles features and
     // adds a bunch of points on the surface

Periodic_3_mesh_3/include/CGAL/refine_periodic_3_mesh_3.h

@@ -17,13 +17,13 @@
 
 #include <CGAL/license/Periodic_3_mesh_3.h>
 
-#include <CGAL/Mesh_3/config.h>
 #include <CGAL/Periodic_3_mesh_3/config.h>
+#include <CGAL/optimize_periodic_3_mesh_3.h>
 
 #include <CGAL/Mesh_3/C3T3_helpers.h>
-#include <CGAL/SMDS_3/Dump_c3t3.h>
 #include <CGAL/Mesh_3/Triangulation_helpers.h>
 #include <CGAL/refine_mesh_3.h>
+#include <CGAL/SMDS_3/Dump_c3t3.h>
 #include <CGAL/Time_stamper.h>
 
 #include <CGAL/Named_function_parameters.h>
@@ -34,22 +34,8 @@
 #include <iterator>
 
 namespace CGAL {
-
 namespace internal {
 
-template<class C3T3, class MeshDomain>
-void project_dummy_points_of_surface(C3T3& c3t3, const MeshDomain& domain)
-{
-  typedef typename C3T3::Vertex_handle                     Vertex_handle;
-  typedef CGAL::Hash_handles_with_or_without_timestamps    Hash_fct;
-  typedef boost::unordered_set<Vertex_handle, Hash_fct>    Vertex_container;
-
-  Vertex_container vertex_container;
-  find_points_to_project(c3t3, std::insert_iterator<Vertex_container>(vertex_container, vertex_container.begin()));
-
-  project_points(c3t3, domain, vertex_container.begin(), vertex_container.end());
-}
-
 template<class C3T3, class OutputIterator>
 void find_points_to_project(C3T3& c3t3, OutputIterator vertices)
 {
@@ -65,21 +51,22 @@ void find_points_to_project(C3T3& c3t3, OutputIterator vertices)
     int ind = face_it->second;
     Cell_handle c = face_it->first;
 
-    for(int i = 1; i < 4; i++) {
+    for(int i = 1; i < 4; i++)
+    {
       Vertex_handle v = c->vertex((ind+i)&3);
-
-      typename C3T3::Index index = c3t3.index(v);
-      if(const int* i = boost::get<int>(&index))
+      if(v->info().is_dummy_vertex)
       {
-        if(*i == 0) // '0' is the index of dummies
-          *vertices++ = v;
+#ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
+        std::cout << c3t3.triangulation().point(v) << " must be projected" << std::endl;
+#endif
+        *vertices++ = v;
       }
     }
   }
 }
 
 template<class C3T3, class MeshDomain, class InputIterator>
-void project_points(C3T3& c3t3,
+bool project_points(C3T3& c3t3,
                     const MeshDomain& domain,
                     InputIterator vertex_begin,
                     InputIterator vertex_end)
@@ -95,40 +82,115 @@ void project_points(C3T3& c3t3,
 
   typename C3T3::Triangulation::Geom_traits::Construct_point_3 cp =
     c3t3.triangulation().geom_traits().construct_point_3_object();
+  typename C3T3::Triangulation::Geom_traits::Compute_squared_distance_3 csd =
+    c3t3.triangulation().geom_traits().compute_squared_distance_3_object();
 
   CGAL::Mesh_3::C3T3_helpers<C3T3, MeshDomain> helper(c3t3, domain);
   CGAL::Mesh_3::Triangulation_helpers<Tr> tr_helpers;
 
+  bool did_something = false;
+
   for(InputIterator it = vertex_begin; it != vertex_end; ++it)
   {
-    Vertex_handle vh = *it;
+    Vertex_handle old_vertex = *it;
 
-    const Weighted_point& vh_wp = c3t3.triangulation().point(vh);
-    const Bare_point& vh_p = cp(vh_wp);
-    const Bare_point new_point = helper.project_on_surface(vh, vh_p);
+    const Weighted_point& weighted_old_position = c3t3.triangulation().point(old_vertex);
+    CGAL_assertion(weighted_old_position.weight() == FT(0)); // point projection happens before optimizers
 
-    const FT sq_d = CGAL::squared_distance(new_point, vh_p);
+    const Bare_point& old_position = cp(weighted_old_position);
+    const Bare_point new_position = helper.project_on_surface(old_vertex, old_position);
+    const FT sq_d = csd(new_position, old_position);
 
 #ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
-    std::cerr << "vh: " << &*vh << std::endl;
-    std::cerr << "vhp: " << vh_p << std::endl;
-    std::cerr << "projected: " << new_point << std::endl;
+    std::cerr << "\n\nMove dummy vertex" << std::endl;
+    std::cerr << "old_vertex: " << &*old_vertex << std::endl;
+    std::cerr << "old_position: " << old_position << std::endl;
+    std::cerr << "new_position: " << new_position << std::endl;
     std::cerr << "squared distance from dummy to surface: " << sq_d << std::endl;
 #endif
 
     // Skip tiny moves for efficiency
-    if(sq_d < 1e-10) // arbitrary value, maybe compare it to the surface distance criterium ?
+    auto min_v_and_sqd = c3t3.triangulation().nearest_power_vertex_with_sq_distance(old_vertex);
+    CGAL_postcondition(min_v_and_sqd.first != Vertex_handle() && min_v_and_sqd.second != FT(-1));
+
+    if(sq_d < 0.01 * min_v_and_sqd.second)
+    {
+#ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
+      std::cout << "REJECTED because dummy point is close enough to the surface" << std::endl;
+#endif
       continue;
+    }
 
     // Do not project if the projected point is in a protection ball
-    if(tr_helpers.inside_protecting_balls(c3t3.triangulation(), vh, new_point))
+    if(tr_helpers.inside_protecting_balls(c3t3.triangulation(), old_vertex, new_position))
+    {
+#ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
+      std::cout << "REJECTED because new pos is within protection ball" << std::endl;
+#endif
       continue;
+    }
 
-    const Vector_3 move(vh_p, new_point);
-    Vertex_handle new_vertex = helper.update_mesh(vh, move);
-    if(new_vertex != vh) // if the move has successfully been performed
-      c3t3.set_dimension(new_vertex, 2);
+    // For periodic triangulations, the move is always performed using insert+remove,
+    // so new_vertex cannot be old_vertex if the move has succeeded
+    const Vector_3 move(old_position, new_position);
+    Vertex_handle new_vertex = helper.update_mesh(old_vertex, move);
+
+    // if the move has successfully been performed
+    if(new_vertex != old_vertex && new_vertex != Vertex_handle())
+    {
+      new_vertex->info().is_dummy_vertex = false;
+      c3t3.set_dimension(new_vertex, 2); // on the surface
+
+      // @fixme
+      // This actually should be the index from the surface patch index...
+      // It can be obtained either by modifying project_on_surface to return the surface_patch index
+      auto opt_si = domain.is_in_domain_object()(cp(c3t3.triangulation().point(new_vertex)));
+      if(opt_si.has_value())
+        c3t3.set_index(new_vertex, domain.index_from_subdomain_index(*opt_si));
+      else
+        c3t3.set_index(new_vertex, 0);
+    }
+    else
+    {
+#ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
+      std::cerr << "Warning: failed to create projection" << std::endl;
+#endif
+    }
+
+    // The vertex `old_vertex` can still exist in the P3RT3:
+    // - if the target already existed
+    // - if its removal would have compromised the 1-cover property of the periodic triangulation
+    // It's (almost) pointless to try and move it again, so fix it
+    if(c3t3.triangulation().tds().is_vertex(old_vertex))
+    {
+#ifdef CGAL_PERIODIC_3_MESH_3_DEBUG_DUMMY_PROJECTION
+      std::cerr << "Warning: failed to remove pre-projection: " << c3t3.triangulation().point(old_vertex) << std::endl;
+#endif
+      old_vertex->info().is_dummy_vertex = false;
+    }
+
+    did_something = true;
   }
+
+  return did_something;
+}
+
+template<class C3T3, class MeshDomain>
+void project_dummy_points_of_surface(C3T3& c3t3,
+                                     const MeshDomain& domain)
+{
+  typedef typename C3T3::Vertex_handle                     Vertex_handle;
+  typedef CGAL::Hash_handles_with_or_without_timestamps    Hash_fct;
+  typedef boost::unordered_set<Vertex_handle, Hash_fct>    Vertex_container;
+
+  bool did_something = false;
+  do
+  {
+    Vertex_container vertex_container;
+    find_points_to_project(c3t3, std::insert_iterator<Vertex_container>(vertex_container, vertex_container.begin()));
+    did_something = project_points(c3t3, domain, vertex_container.begin(), vertex_container.end());
+  }
+  while(did_something);
 }
 
 } // namespace internal
@@ -410,21 +472,21 @@ void refine_periodic_3_mesh_3_impl(C3T3& c3t3,
   // Odt
   if(odt)
   {
-    odt_optimize_mesh_3(c3t3, domain,
-                        parameters::time_limit = odt.time_limit(),
-                        parameters::max_iteration_number = odt.max_iteration_number(),
-                        parameters::convergence = odt.convergence(),
-                        parameters::freeze_bound = odt.bound());
+    odt_optimize_periodic_3_mesh_3(c3t3, domain,
+                                   parameters::time_limit = odt.time_limit(),
+                                   parameters::max_iteration_number = odt.max_iteration_number(),
+                                   parameters::convergence = odt.convergence(),
+                                   parameters::freeze_bound = odt.bound());
   }
 
   // Lloyd
   if(lloyd)
   {
-    lloyd_optimize_mesh_3(c3t3, domain,
-                          parameters::time_limit = lloyd.time_limit(),
-                          parameters::max_iteration_number = lloyd.max_iteration_number(),
-                          parameters::convergence = lloyd.convergence(),
-                          parameters::freeze_bound = lloyd.bound());
+    lloyd_optimize_periodic_3_mesh_3(c3t3, domain,
+                                     parameters::time_limit = lloyd.time_limit(),
+                                     parameters::max_iteration_number = lloyd.max_iteration_number(),
+                                     parameters::convergence = lloyd.convergence(),
+                                     parameters::freeze_bound = lloyd.bound());
   }
 
   if(odt || lloyd)
@@ -440,9 +502,9 @@ void refine_periodic_3_mesh_3_impl(C3T3& c3t3,
     if(perturb.is_time_limit_set())
       perturb_time_limit = perturb.time_limit();
 
-    perturb_mesh_3(c3t3, domain,
-                   parameters::time_limit = perturb_time_limit,
-                   parameters::sliver_bound = perturb.bound());
+    perturb_periodic_3_mesh_3(c3t3, domain,
+                              parameters::time_limit = perturb_time_limit,
+                              parameters::sliver_bound = perturb.bound());
 
     dump_c3t3(c3t3, mesh_options.dump_after_perturb_prefix);
   }
@@ -455,9 +517,9 @@ void refine_periodic_3_mesh_3_impl(C3T3& c3t3,
     if(exude.is_time_limit_set())
       exude_time_limit = exude.time_limit();
 
-    exude_mesh_3(c3t3,
-                 parameters::time_limit = exude_time_limit,
-                 parameters::sliver_bound = exude.bound());
+    exude_periodic_3_mesh_3(c3t3,
+                            parameters::time_limit = exude_time_limit,
+                            parameters::sliver_bound = exude.bound());
 
     dump_c3t3(c3t3, mesh_options.dump_after_perturb_prefix);
   }

Periodic_3_mesh_3/test/Periodic_3_mesh_3/test_implicit_shapes_bunch.cpp

@@ -202,7 +202,7 @@ int main()
   Iso_cuboid canonical_cube(0, 0, 0, domain_size, domain_size, domain_size);
 
   // Array of the functions
-  const int functions_count = 11;
+  const unsigned int functions_count = 11;
 
   std::array<Periodic_function, functions_count> implicit_functions =
   {{
@@ -219,7 +219,7 @@ int main()
     Periodic_function(split_p, canonical_cube)
   }};
 
-  for(int i=0; i<functions_count; ++i)
+  for(unsigned int i=0; i<functions_count; ++i)
   {
     // Periodic mesh domain
     Periodic_mesh_domain domain =

Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/CGAL/Periodic_3_triangulation_3.h

@@ -313,8 +313,7 @@ The elements of the enum have the following meaning:
 /// @{
 
 /*!
-Introduces an empty triangulation `t` with `domain` as
-original domain.
+Introduces an empty triangulation `t` with `domain` as original domain.
 \pre `domain` is a cube.
 */
 Periodic_3_triangulation_3(const Iso_cuboid & domain = Iso_cuboid(0,0,0,1,1,1),

Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3RegularTriangulationTraits_3.h

@@ -81,10 +81,33 @@ typedef unspecified_type Power_side_of_oriented_power_sphere_3;
 /*!
 A predicate object that must provide the function operators:
 
-`Orientation operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s, FT w)`,
+` operator()(Point_3 p, Point_3 q, Point_3 r, Point_3 s,
+             Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q, Periodic_3_offset_3 o_r, Periodic_3_offset_3 o_s)`,
+
+which compares the squared distance between `(p, o_p)` and `(q, o_q)` and the squared distance
+between `(r, o_r)` and `(s, o_s)`, and returns `SMALLER`, `EQUAL`, or `LARGER`.
+
+\pre `p`, `q`, `r`, and `s` lie inside the domain.
+
+*/
+typedef unspecified_type Compare_squared_distance_3;
+
+/// @}
+
+/// \name
+/// @{
+
+/*!
+A predicate object that must provide the function operators:
+
+`Comparison_result operator()(Weighted_point_3 p, Weighted_point_3 q, Weighted_point_3 r, Weighted_point_3 s, FT w)`,
+
+and
+
+`Comparison_result operator()(Weighted_point_3 p, FT w)`,
 
 which compares the weight of the smallest sphere orthogonal to the input weighted
-points with the input weight `w` and returns a `SMALLER`, `EQUAL`, or `LARGER`.
+points with the input weight `w` and returns `SMALLER`, `EQUAL`, or `LARGER`.
 
 \pre `p`, `q`, `r`, and `s` lie inside the domain.
 

Periodic_3_triangulation_3/doc/Periodic_3_triangulation_3/Concepts/Periodic_3TriangulationTraits_3.h

@@ -84,6 +84,41 @@ A predicate object that must provide the function operator
 
 `Comparison_result operator()(Point_3 p, Point_3 q, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`,
 
+which returns `SMALLER` (`EQUAL`, `LARGER`) if the x-coordinate of `p` is smaller (equal, larger) than
+the x-coordinate of `q`.
+
+\pre `p`, `q` lie inside the domain.
+*/
+typedef unspecified_type Compare_x_3;
+
+/*!
+A predicate object that must provide the function operator
+
+`Comparison_result operator()(Point_3 p, Point_3 q, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`,
+
+which returns `SMALLER` (`EQUAL`, `LARGER`) if the y-coordinate of `p` is smaller (equal, larger) than
+the y-coordinate of `q`.
+\pre `p`, `q` lie inside the domain.
+*/
+typedef unspecified_type Compare_y_3;
+
+/*!
+A predicate object that must provide the function operator
+
+`Comparison_result operator()(Point_3 p, Point_3 q, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`,
+
+which returns `SMALLER` (`EQUAL`, `LARGER`) if the z-coordinate of `p` is smaller (equal, larger) than
+the z-coordinate of `q`.
+
+\pre `p`, `q` lie inside the domain.
+*/
+typedef unspecified_type Compare_z_3;
+
+/*!
+A predicate object that must provide the function operator
+
+`Comparison_result operator()(Point_3 p, Point_3 q, Periodic_3_offset_3 o_p, Periodic_3_offset_3 o_q)`,
+
 which returns `EQUAL` if the two point-offset pairs are equal.
 Otherwise it must return a consistent order for any two points chosen
 in a same line.

Periodic_3_triangulation_3/include/CGAL/Periodic_3_Delaunay_triangulation_3.h

@@ -665,29 +665,6 @@ private:
   }
 
 public:
-  /** @name Geometric access functions */
-
-  Point point(const Periodic_point& pp) const
-  {
-    return point(pp, geom_traits().construct_point_3_object());
-  }
-
-  // The following functions return the "real" position in space (unrestrained
-  // to the fundamental domain) of the vertices v and c->vertex(idx),
-  // respectively
-
-  Point point(Vertex_handle v) const
-  {
-    return point(v, geom_traits().construct_point_3_object());
-  }
-
-  Point point(Cell_handle c, int idx) const
-  {
-    return point(c, idx, geom_traits().construct_point_3_object());
-  }
-
-  // end of geometric functions
-
   Periodic_point periodic_circumcenter(Cell_handle c) const {
     return Base::periodic_circumcenter(c, geom_traits().construct_circumcenter_3_object());
   }
@@ -926,17 +903,8 @@ Periodic_3_Delaunay_triangulation_3<Gt,Tds>::
 move_point(Vertex_handle v, const Point& p)
 {
   CGAL_expensive_precondition(is_vertex(v));
-  // Remember an incident vertex to restart
-  // the point location after the removal.
-  // Cell_handle c = v->cell();
-  //Vertex_handle old_neighbor = c->vertex(c->index(v) == 0 ? 1 : 0);
-  //  CGAL_assertion(old_neighbor != v);
-
   remove(v);
-
-  if(number_of_vertices() == 0)
-    return insert(p);
-  return insert(p);//, old_neighbor->cell());
+  return insert(p);
 }
 
 template < class Gt, class Tds >

Periodic_3_triangulation_3/include/CGAL/Periodic_3_offset_3.h

@@ -20,6 +20,9 @@
 #include <CGAL/assertions.h>
 #include <CGAL/Cartesian.h>
 
+#include <iostream>
+#include <type_traits>
+
 namespace CGAL {
 
 class Periodic_3_offset_3
@@ -43,18 +46,27 @@ public:
   int& z() { return _offz; }
   int z() const { return _offz; }
 
-  int &operator[](int i) {
+  // Use sfinae on the operator[] to accept only integral types as argument
+  template <typename T,
+            typename std::enable_if<std::is_integral<T>::value>::type* = nullptr>
+  int& operator[](T i)
+  {
     if (i==0) return _offx;
     if (i==1) return _offy;
     CGAL_assertion(i==2);
     return _offz;
   }
-  int operator[](int i) const {
+
+  template <typename T,
+            typename std::enable_if<std::is_integral<T>::value>::type* = nullptr>
+  int operator[](T i) const
+  {
     if (i==0) return _offx;
     if (i==1) return _offy;
     CGAL_assertion(i==2);
     return _offz;
   }
+
   void operator+=(const Periodic_3_offset_3 &other) {
     _offx += other._offx;
     _offy += other._offy;

Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_3.h

@@ -220,13 +220,20 @@ public:
   };
 
 public:
+  FT compute_cover_threshold() const
+  {
+    FT min_span = (std::min)({ domain().xmax() - domain().xmin(),
+                               domain().ymax() - domain().ymin(),
+                               domain().zmax() - domain().zmin() });
+    return FT(0.015625) * CGAL::square(min_span);
+  }
+
   /** @name Creation */
   Periodic_3_regular_triangulation_3(const Iso_cuboid& domain = Iso_cuboid(0, 0, 0, 1, 1, 1),
                                      const Geometric_traits& gt = Geometric_traits())
     : Tr_Base(domain, gt)
   {
-    orthosphere_radius_threshold = FT(0.015625) * (domain.xmax() - domain.xmin())
-                                                * (domain.xmax() - domain.xmin());
+    orthosphere_radius_threshold = compute_cover_threshold();
   }
 
   template < typename InputIterator >
@@ -236,8 +243,7 @@ public:
                                      bool is_large_point_set = false)
     : Tr_Base(domain, gt)
   {
-    orthosphere_radius_threshold = FT(0.015625) * (domain.xmax() - domain.xmin())
-                                                * (domain.xmax() - domain.xmin());
+    orthosphere_radius_threshold = compute_cover_threshold();
 
     insert(first, last, is_large_point_set);
   }
@@ -383,8 +389,7 @@ public:
 
   virtual void update_cover_data_after_setting_domain ()
   {
-    orthosphere_radius_threshold = FT(0.015625) * (domain().xmax() - domain().xmin())
-                                                * (domain().xmax() - domain().xmin());
+    orthosphere_radius_threshold = compute_cover_threshold();
   }
 
   // the function below is used in `convert_to_1_sheeted_covering()` of P3T3
@@ -609,6 +614,8 @@ public:
     return geom_traits().compare_power_distance_3_object()(p, q, r, o1, o2, o3) == SMALLER;
   }
 
+  // @fixme the overloads with offset might run into an issue if the intermediate construction
+  // does not preserve the orientation... See Robust_periodic_weighted_circumcenter_traits_3.h
   Bare_point construct_weighted_circumcenter(const Weighted_point &p, const Weighted_point &q,
                                              const Weighted_point &r) const
   {
@@ -1746,12 +1753,10 @@ operator>> (std::istream& is, Periodic_3_regular_triangulation_3<GT, TDS>& tr)
 {
   typedef Periodic_3_regular_triangulation_3<GT,TDS>   P3RT3;
   typedef typename P3RT3::Tr_Base                      Tr_Base;
-  typedef typename GT::FT                              FT;
 
   is >> static_cast<Tr_Base&>(tr);
 
-  tr.orthosphere_radius_threshold = FT(0.015625) * (tr.domain().xmax() - tr.domain().xmin())
-                                                 * (tr.domain().xmax() - tr.domain().xmin());
+  tr.orthosphere_radius_threshold = tr.compute_cover_threshold();
 
   tr.insert_cells_with_too_big_orthoball(tr.cells_begin(), tr.cells_end());
 

Periodic_3_triangulation_3/include/CGAL/Periodic_3_regular_triangulation_traits_3.h

@@ -85,6 +85,8 @@ public:
       Compute_power_distance_to_power_sphere_3;
   typedef Functor_with_offset_weighted_points_adaptor_3<Self, typename Kernel::Compute_squared_distance_3>
       Compute_squared_distance_3;
+  typedef Functor_with_offset_weighted_points_adaptor_3<Self, typename Kernel::Compare_squared_distance_3>
+      Compare_squared_distance_3;
 
   // Operations
   Construct_weighted_point_3 construct_weighted_point_3_object() const {
@@ -111,6 +113,12 @@ public:
       this->construct_point_3_object(), construct_weighted_point_3_object());
   }
 
+  Compare_squared_distance_3 compare_squared_distance_3_object() const {
+    return Compare_squared_distance_3(
+      this->Base::compare_squared_distance_3_object(),
+      this->construct_point_3_object(), construct_weighted_point_3_object());
+  }
+
   // predicates
   Power_side_of_bounded_power_sphere_3 power_side_of_bounded_power_sphere_3_object() const {
     return Power_side_of_bounded_power_sphere_3(

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3.h

@@ -22,16 +22,19 @@
 #include <CGAL/basic.h>
 
 #include <CGAL/Periodic_3_triangulation_3/internal/Periodic_3_triangulation_iterators_3.h>
+#include <CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h>
+
 #include <CGAL/Periodic_3_triangulation_ds_cell_base_3.h>
 #include <CGAL/Periodic_3_triangulation_ds_vertex_base_3.h>
 #include <CGAL/Periodic_3_triangulation_traits_3.h>
 #include <CGAL/Triangulation_data_structure_3.h>
 #include <CGAL/Triangulation_cell_base_3.h>
 #include <CGAL/Triangulation_vertex_base_3.h>
-#include <CGAL/assertions.h>
-#include <CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h>
+#include <CGAL/Triangulation_vertex_base_with_info_3.h>
+#include <CGAL/Delaunay_triangulation_3.h>
 
 #include <CGAL/array.h>
+#include <CGAL/assertions.h>
 #include <CGAL/Number_types/internal/Exact_type_selector.h>
 #include <CGAL/NT_converter.h>
 #include <CGAL/Unique_hash_map.h>
@@ -123,7 +126,7 @@ public:
 
   typedef typename GT::Periodic_3_offset_3     Offset;
   typedef typename GT::Iso_cuboid_3            Iso_cuboid;
-  typedef std::array<int, 3>           Covering_sheets;
+  typedef std::array<int, 3>                   Covering_sheets;
 
   // point types
   typedef typename TDS::Vertex::Point          Point;
@@ -259,17 +262,6 @@ public:
                              const Geometric_traits& gt = Geometric_traits())
     : _gt(gt), _tds()
   {
-    typedef typename internal::Exact_field_selector<FT>::Type EFT;
-    typedef NT_converter<FT,EFT> NTC;
-    CGAL_USE_TYPE(NTC);
-    CGAL_precondition_code( NTC ntc; )
-    CGAL_precondition(ntc(domain.xmax())-ntc(domain.xmin())
-                                    == ntc(domain.ymax())-ntc(domain.ymin()));
-    CGAL_precondition(ntc(domain.ymax())-ntc(domain.ymin())
-                                    == ntc(domain.zmax())-ntc(domain.zmin()));
-    CGAL_precondition(ntc(domain.zmax())-ntc(domain.zmin())
-                                    == ntc(domain.xmax())-ntc(domain.xmin()));
-
     _gt.set_domain(domain);
     _cover = CGAL::make_array(3,3,3);
     init_tds();
@@ -487,6 +479,8 @@ public:
     return i;
   }
 
+  // The first 3 bits are a 0 or 1 offsets in the xyz directions
+  // For example, 6 = 4 + 2 + 0 <=> 110 in binary <=> Offset(1, 1, 0)
   Offset int_to_off(int i) const
   {
     return Offset((i>>2)&1,(i>>1)&1,i&1);
@@ -694,7 +688,8 @@ public:
   // to the fundamental domain) of the vertices v and c->vertex(idx),
   // respectively
   template <class ConstructPoint>
-  Point point(Vertex_handle v, ConstructPoint cp) const {
+  Point point(Vertex_handle v, ConstructPoint cp) const
+  {
     return point(periodic_point(v), cp);
   }
 
@@ -755,6 +750,25 @@ public:
     return Point();
   }
 
+  Point point(const Periodic_point& pp) const
+  {
+    return point(pp, geom_traits().construct_point_3_object());
+  }
+
+  // The following functions return the "real" position in space (unrestrained
+  // to the fundamental domain) of the vertices v and c->vertex(idx),
+  // respectively
+
+  Point point(Vertex_handle v) const
+  {
+    return point(v, geom_traits().construct_point_3_object());
+  }
+
+  Point point(Cell_handle c, int idx) const
+  {
+    return point(c, idx, geom_traits().construct_point_3_object());
+  }
+
   Periodic_point periodic_point(const Vertex_handle v) const
   {
     if(is_1_cover())
@@ -1288,6 +1302,8 @@ private:
 public:
   std::vector<Vertex_handle> insert_dummy_points();
 
+  std::vector<Vertex_handle> insert_generic_dummy_points();
+
 protected:
   // this is needed for compatibility reasons
   template <class Conflict_test, class OutputIteratorBoundaryFacets,
@@ -1762,6 +1778,8 @@ protected:
   {
     CGAL_precondition(c != Cell_handle());
 
+    // @fixme this might run into an issue if the intermediate construction does not preserve
+    // the orientation... See Robust_periodic_weighted_circumcenter_traits_3.h
     Point_3 p = construct_circumcenter(c->vertex(0)->point(), c->vertex(1)->point(),
                                        c->vertex(2)->point(), c->vertex(3)->point(),
                                        get_offset(c, 0), get_offset(c, 1),
@@ -2290,14 +2308,10 @@ inexact_periodic_locate(const Point& p, const Offset& o_p,
   }
 
   CGAL_postcondition(start!=Cell_handle());
-  CGAL_assertion(start->neighbor(0)->neighbor(
-      start->neighbor(0)->index(start))==start);
-  CGAL_assertion(start->neighbor(1)->neighbor(
-      start->neighbor(1)->index(start))==start);
-  CGAL_assertion(start->neighbor(2)->neighbor(
-      start->neighbor(2)->index(start))==start);
-  CGAL_assertion(start->neighbor(3)->neighbor(
-      start->neighbor(3)->index(start))==start);
+  CGAL_assertion(start->neighbor(0)->neighbor(start->neighbor(0)->index(start)) == start);
+  CGAL_assertion(start->neighbor(1)->neighbor(start->neighbor(1)->index(start)) == start);
+  CGAL_assertion(start->neighbor(2)->neighbor(start->neighbor(2)->index(start)) == start);
+  CGAL_assertion(start->neighbor(3)->neighbor(start->neighbor(3)->index(start)) == start);
 
   // We implement the remembering visibility/stochastic walk.
 
@@ -2523,7 +2537,6 @@ Periodic_3_triangulation_3<GT,TDS>::periodic_insert(
     Locate_type /*lt*/, Cell_handle c, const Conflict_tester& tester,
     Point_hider& hider, CoverManager& cover_manager, Vertex_handle vh)
 {
-  Vertex_handle v;
   CGAL_precondition(number_of_vertices() != 0);
   CGAL_assertion_code(
       Locate_type lt_assert; int i_assert; int j_assert;);
@@ -2568,7 +2581,7 @@ Periodic_3_triangulation_3<GT,TDS>::periodic_insert(
   // Insertion. Warning: facets[0].first MUST be in conflict!
   // Compute the star and put it into the data structure.
   // Store the new cells from the star in nbs.
-  v = _tds._insert_in_hole(cells.begin(), cells.end(), facet.first, facet.second);
+  Vertex_handle v = _tds._insert_in_hole(cells.begin(), cells.end(), facet.first, facet.second);
   v->set_point(p);
 
   //TODO: this could be done within the _insert_in_hole without losing any
@@ -2577,6 +2590,7 @@ Periodic_3_triangulation_3<GT,TDS>::periodic_insert(
   //- Find the modified _insert_in_hole in the branch svn history of TDS
   std::vector<Cell_handle> nbs;
   incident_cells(v, std::back_inserter(nbs));
+
   // For all neighbors of the newly added vertex v: fetch their offsets from
   // the tester and reset them in the triangulation data structure.
   for(typename std::vector<Cell_handle>::iterator cit = nbs.begin();
@@ -2751,8 +2765,10 @@ Periodic_3_triangulation_3<GT,TDS>::create_initial_triangulation(const Point& p)
 
   return vir_vertices[0][0][0];
 }
+
 #define CGAL_INCLUDE_FROM_PERIODIC_3_TRIANGULATION_3_H
 #include <CGAL/Periodic_3_triangulation_3/internal/Periodic_3_triangulation_dummy_36.h>
+#include <CGAL/Periodic_3_triangulation_3/internal/Periodic_3_triangulation_dummy_generator.h>
 #undef CGAL_INCLUDE_FROM_PERIODIC_3_TRIANGULATION_3_H
 
 /** finds all cells that are in conflict with the currently added point
@@ -2848,13 +2864,11 @@ template < class Conflict_tester, class Point_hider, class CoverManager >
 inline typename Periodic_3_triangulation_3<GT,TDS>::Vertex_handle
 Periodic_3_triangulation_3<GT,TDS>::insert_in_conflict(const Point& p,
     Locate_type lt, Cell_handle c, int li, int lj,
-    const Conflict_tester& tester, Point_hider& hider, CoverManager& cover_manager) {
-  CGAL_assertion((domain().xmin() <= p.x())
-                               && (p.x() < domain().xmax()));
-  CGAL_assertion((domain().ymin() <= p.y())
-                               && (p.y() < domain().ymax()));
-  CGAL_assertion((domain().zmin() <= p.z())
-                               && (p.z() < domain().zmax()));
+    const Conflict_tester& tester, Point_hider& hider, CoverManager& cover_manager)
+{
+  CGAL_assertion((domain().xmin() <= p.x()) && (p.x() < domain().xmax()));
+  CGAL_assertion((domain().ymin() <= p.y()) && (p.y() < domain().ymax()));
+  CGAL_assertion((domain().zmin() <= p.z()) && (p.z() < domain().zmax()));
 
   if(number_of_vertices() == 0) {
     Vertex_handle vh = create_initial_triangulation(p);
@@ -2874,10 +2888,8 @@ Periodic_3_triangulation_3<GT,TDS>::insert_in_conflict(const Point& p,
       vstart = c->vertex(0);
     else
       vstart = vvmit->second.first;
-    CGAL_assertion(virtual_vertices.find(vstart)
-                                 == virtual_vertices.end());
-    CGAL_assertion(virtual_vertices_reverse.find(vstart)
-                                 != virtual_vertices_reverse.end());
+    CGAL_assertion(virtual_vertices.find(vstart) == virtual_vertices.end());
+    CGAL_assertion(virtual_vertices_reverse.find(vstart) != virtual_vertices_reverse.end());
   }
 
   CGAL_assertion( number_of_vertices() != 0 );
@@ -2971,8 +2983,8 @@ is_valid(bool verbose, int level) const
   }
 
   bool error = false;
-  for(Cell_iterator cit = cells_begin();
-       cit != cells_end(); ++cit) {
+  for(Cell_iterator cit = cells_begin(); cit != cells_end(); ++cit) {
+    CGAL_assertion(cit != Cell_handle());
     for(int i=0; i<4; i++) {
       CGAL_assertion(cit != cit->neighbor(i));
       for(int j=i+1; j<4; j++) {
@@ -2986,14 +2998,15 @@ is_valid(bool verbose, int level) const
       p[i] = &cit->vertex(i)->point();
       off[i] = get_offset(cit,i);
     }
+
     if(orientation(*p[0], *p[1], *p[2], *p[3],
                    off[0], off[1], off[2], off[3]) != POSITIVE) {
       if(verbose) {
-        std::cerr<<"Periodic_3_triangulation_3: wrong orientation:"<<std::endl;
-        std::cerr<<off[0]<<'\t'<<*p[0]<<'\n'
-                         <<off[1]<<'\t'<<*p[1]<<'\n'
-                         <<off[2]<<'\t'<<*p[2]<<'\n'
-                         <<off[3]<<'\t'<<*p[3]<<std::endl;
+        std::cerr << "Periodic_3_triangulation_3: wrong orientation:"<<std::endl;
+        std::cerr << off[0] << '\t' << *p[0] << '\n'
+                  << off[1] << '\t' << *p[1] << '\n'
+                  << off[2] << '\t' << *p[2] << '\n'
+                  << off[3] << '\t' << *p[3] << std::endl;
       }
       error = true;
     }
@@ -3145,7 +3158,9 @@ template < class GT, class TDS >
 template < class PointRemover, class CoverManager >
 inline bool
 Periodic_3_triangulation_3<GT,TDS>::
-periodic_remove(Vertex_handle v, PointRemover& remover, CoverManager& cover_manager,
+periodic_remove(Vertex_handle v,
+                PointRemover& remover,
+                CoverManager& cover_manager,
                 const bool abort_if_cover_change)
 {
   // Construct the set of vertex triples on the boundary
@@ -3278,6 +3293,8 @@ periodic_remove(Vertex_handle v, PointRemover& remover, CoverManager& cover_mana
                         vh_off_map[vmap[i_ch->vertex(2)]],
                         vh_off_map[vmap[i_ch->vertex(3)]]);
 
+    // cells created above are deleted in update_cover_data_during_management() even if we abort
+
     // Update the edge length management
     if(cover_manager.update_cover_data_during_management(new_ch, new_cells,
                                                          abort_if_cover_change))

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Periodic_3_regular_triangulation_dummy_288.h

@@ -37,7 +37,7 @@ std::vector<Weighted_point> dummy_points()
 
         FT z = (FT(k) * domain_z / FT(8)) + domain().zmin();
 
-        dummy_points.push_back(Weighted_point(Bare_point(x, y, z), 0));
+        dummy_points.emplace_back(Bare_point(x, y, z), 0);
       }
     }
   }

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Periodic_3_triangulation_dummy_generator.h

@@ -0,0 +1,636 @@
+// Copyright (c) 2009   INRIA Sophia-Antipolis (France).
+// All rights reserved.
+//
+// This file is part of CGAL (www.cgal.org).
+//
+// $URL$
+// $Id$
+// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
+//
+//
+// Author(s)     : Mael Rouxel-Labbé
+
+#ifdef CGAL_INCLUDE_FROM_PERIODIC_3_TRIANGULATION_3_H
+
+template < class GT, class TDS >
+inline std::vector<typename Periodic_3_triangulation_3<GT,TDS>::Vertex_handle >
+Periodic_3_triangulation_3<GT, TDS>::
+insert_generic_dummy_points()
+{
+  // the "info" field is the index of the neighbor offset in `neigh_offsets`
+  using DT3_Vb = CGAL::Triangulation_vertex_base_with_info_3<std::size_t, GT>;
+  using DT3_Cb = CGAL::Delaunay_triangulation_cell_base_3<GT>;
+  using DT3_TDS = CGAL::Triangulation_data_structure_3<DT3_Vb, DT3_Cb>;
+  using DT3 = CGAL::Delaunay_triangulation_3<GT, DT3_TDS>;
+  using DT3_VH = typename DT3::Vertex_handle;
+  using DT3_CH = typename DT3::Cell_handle;
+
+  // Compute the number of subdivisions in all directions ------------------------------------------
+
+  const std::array<FT, 3> spans = { domain().xmax() - domain().xmin(),
+                                    domain().ymax() - domain().ymin(),
+                                    domain().zmax() - domain().zmin() };
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 1)
+  std::cout << "Domain spans:\n"
+            << domain().xmin() << " " << domain().xmax() << "\n"
+            << domain().ymin() << " " << domain().ymax() << "\n"
+            << domain().zmin() << " " << domain().zmax() << std::endl;
+#endif
+
+  std::size_t min_pos = std::min_element(std::cbegin(spans), std::cend(spans)) - std::cbegin(spans);
+  std::size_t max_pos = std::max_element(std::cbegin(spans), std::cend(spans)) - std::cbegin(spans);
+  std::size_t mid_pos;
+
+  if(min_pos == max_pos) // cubic
+  {
+    min_pos = 0; // x
+    mid_pos = 1; // y
+    max_pos = 2; // z
+  }
+  else
+  {
+    mid_pos = (min_pos + 1) % 3;
+    if(mid_pos == max_pos)
+      mid_pos = (max_pos + 1) % 3;
+
+    if(min_pos > mid_pos && spans[min_pos] == spans[mid_pos])
+      std::swap(min_pos, mid_pos); // just for convenience
+  }
+
+  CGAL_assertion(min_pos < 3 && mid_pos < 3 && max_pos < 3);
+  CGAL_assertion(min_pos != max_pos && mid_pos != min_pos && mid_pos != max_pos);
+
+  std::array<int, 3> nums_steps;
+  std::array<FT, 3> steps;
+
+  // Min:
+  // this doesn't work for P3M3 + sharp features due to weights creating a harder constraint of 1/64th domain_size²
+// #define CGAL_P3T3_USE_EXPERIMENTAL_LARGE_STEP_IN_DUMMY_GENERATION
+#ifdef CGAL_P3T3_USE_EXPERIMENTAL_LARGE_STEP_IN_DUMMY_GENERATION
+  nums_steps[min_pos] = 3;
+#else
+  nums_steps[min_pos] = 6;
+#endif
+  steps[min_pos] = spans[min_pos] / nums_steps[min_pos];
+
+  // Mid: do not use the min step, but redistribute the error between nums_steps[mid_pos] * min_step and the actual span
+  nums_steps[mid_pos] = int(to_interval(spans[mid_pos] / steps[min_pos]).first); // flooring, "min" is not a typo
+  steps[mid_pos] = spans[mid_pos] / nums_steps[mid_pos];
+
+  // Max: smaller step in the max span direction as to avoid cospherical configurations
+#ifdef CGAL_P3T3_USE_EXPERIMENTAL_LARGE_STEP_IN_DUMMY_GENERATION
+  const FT minor_step = spans[min_pos] / FT(4); // a ratio of 3:4 makes for nicely shaped tetrahedra
+#else
+  const FT minor_step = spans[min_pos] / FT(8); // a ratio of 6:8 makes for nicely shaped tetrahedra
+#endif
+  nums_steps[max_pos] = int(to_interval(spans[max_pos] / minor_step).first); // flooring
+  CGAL_assertion(nums_steps[max_pos] >= steps[min_pos]);
+
+  // Important! Consecutive levels in the max length have a shift (that's the `k % 2 != 0` part).
+  // Hence, to get proper periodicity, the number of steps needs to be even, otherwise the level
+  // `max - 1` and `0` would have the same shift and the Delaunay stars would not be identical everywhere.
+  if(nums_steps[max_pos] % 2 != 0)
+    ++nums_steps[max_pos];
+
+  // now that the number of steps is known, re-adjust the step length so that it matches the height
+  steps[max_pos] = spans[max_pos] / nums_steps[max_pos];
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 1)
+  std::cout << "min|mid|max: " << min_pos << " " << mid_pos << " " << max_pos << std::endl;
+  std::cout << "nums_steps[min_pos] = " << nums_steps[min_pos] << std::endl;
+  std::cout << "steps[min_pos] = " << steps[min_pos] << std::endl;
+  std::cout << "nums_steps[mid_pos] = " << nums_steps[mid_pos] << std::endl;
+  std::cout << "steps[mid_pos] = " << steps[mid_pos] << std::endl;
+  std::cout << "minor step = " << minor_step << std::endl;
+  std::cout << "nums_steps[max_pos] = " << nums_steps[max_pos] << std::endl;
+  std::cout << "steps[max_pos] = " << steps[max_pos] << std::endl;
+#endif
+
+  auto grid_offset_to_lattice_offset = [&](const Offset& loff) -> Offset
+  {
+    const int hloz = int(std::ceil(double(loff[max_pos]) / 2.));
+
+    Offset goff;
+    goff[min_pos] = loff[min_pos] + hloz;
+    goff[mid_pos] = loff[mid_pos] + hloz;
+    goff[max_pos] = loff[max_pos];
+    return goff;
+  };
+
+  auto lattice_offset_to_grid_offset = [&](const Offset& loff) -> Offset
+  {
+    const int hloz = int(std::ceil(double(loff[max_pos]) / 2.));
+
+    Offset goff;
+    goff[min_pos] = loff[min_pos] - hloz;
+    goff[mid_pos] = loff[mid_pos] - hloz;
+    goff[max_pos] = loff[max_pos];
+    return goff;
+  };
+
+  // This should be prefered from building points from the lattice offset
+  // because the grid offset (by construction) aligns with the canonical domain,
+  // and there is thus fewer numerical errors.
+  auto construct_point_from_grid_offset = [&](const Offset& goff) -> Point_3
+  {
+    std::array<FT, 3> coords;
+    coords[min_pos] = goff[min_pos] * steps[min_pos];
+    coords[mid_pos] = goff[mid_pos] * steps[mid_pos];
+    coords[max_pos] = goff[max_pos] * steps[max_pos];
+
+    if(goff[max_pos] % 2 != 0)
+    {
+      coords[min_pos] += 0.5 * steps[min_pos];
+      coords[mid_pos] += 0.5 * steps[mid_pos];
+    }
+
+    return { coords[0], coords[1], coords[2] };
+  };
+
+  auto construct_point_from_lattice_offset = [&](const Offset& loff) -> Point_3
+  {
+    // @fixme? verify that the lattice is indeed already reduced
+
+    // For `min_pos = 0`, `mid_pos = 1`, `max_pos = 2`, this returns:
+    //   CGAL::Origin + loff[0]*lv0 + loff[1]*lv1 + loff[2]*lv2
+    // with
+    //   Vector lv0 {        x_step,             0,      0 };
+    //   Vector lv1 {             0,        y_step,      0 };
+    //   Vector lv2 { -0.5 * x_step, -0.5 * y_step, z_step };
+
+    std::array<FT, 3> coords;
+    coords[min_pos] = loff[min_pos] * steps[min_pos] - loff[max_pos] * steps[min_pos] / FT(2);
+    coords[mid_pos] = loff[mid_pos] * steps[mid_pos] - loff[max_pos] * steps[mid_pos] / FT(2);
+    coords[max_pos] = loff[max_pos] * steps[max_pos];
+
+    return { coords[0], coords[1], coords[2] };
+  };
+
+  CGAL_USE(construct_point_from_lattice_offset);
+
+  // Create the dummy points -----------------------------------------------------------------------
+
+  std::vector<Point_3> dummy_points;
+  dummy_points.reserve(nums_steps[0] * nums_steps[1] * nums_steps[2]);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 32)
+  DT3 dt3;
+#endif
+
+  std::array<FT, 3> coords;
+  for(int i=0; i<nums_steps[min_pos]; ++i)
+  {
+    for(int j=0; j<nums_steps[mid_pos]; ++j)
+    {
+      for (int k=0; k<nums_steps[max_pos]; ++k)
+      {
+        coords[min_pos] = domain().min_coord(static_cast<int>(min_pos)) + i * steps[min_pos];
+        coords[mid_pos] = domain().min_coord(static_cast<int>(mid_pos)) + j * steps[mid_pos];
+
+        if(k % 2 != 0)
+        {
+          coords[min_pos] += steps[min_pos] / FT(2);
+          coords[mid_pos] += steps[mid_pos] / FT(2);
+        }
+
+        coords[max_pos] = domain().min_coord(static_cast<int>(max_pos)) + k * steps[max_pos];
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 2)
+        std::cout << "Add dummy: " << coords[0] << " " << coords[1] << " " << coords[2] << std::endl;
+#endif
+        dummy_points.emplace_back(coords[0], coords[1], coords[2]);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 32)
+        DT3_VH dt3_v = dt3.insert(dummy_points.back());
+#endif
+      }
+    }
+  }
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 1)
+  std::cout << dummy_points.size() << " dummy points" << std::endl;
+#endif
+
+  // Build the periodic triangulation --------------------------------------------------------------
+
+  clear();
+  tds().set_dimension(3);
+  set_cover(CGAL::make_array(1,1,1));
+
+  const std::size_t nv = static_cast<std::size_t>(nums_steps[0] * nums_steps[1] * nums_steps[2]);
+
+  // Create vertices
+  std::vector<Vertex_handle> vertices;
+  vertices.reserve(nv);
+
+  std::size_t id = 0;
+  for(int i=0; i<nums_steps[min_pos]; ++i)
+  {
+    for(int j=0; j<nums_steps[mid_pos]; ++j)
+    {
+      for(int k=0; k<nums_steps[max_pos]; ++k)
+      {
+        Vertex_handle vh = tds().create_vertex();
+        vertices.push_back(vh);
+        vh->set_point(Point(dummy_points[id++])); // wrap with Point for regular triangulations
+      }
+    }
+  }
+
+  // Create cells
+  std::vector<Offset> star_lnoffs;
+  std::vector<std::array<std::size_t, 4> > star_cells;
+
+  auto compute_periodic_star = [&]() -> void
+  {
+    // Lattice offsets sufficient to get the complete periodic star
+    //
+    // Considering all offsets is likely an overkill because it's a special lattice with bounded angles
+    // and ratio, but it's a tiny triangulation so it does not matter
+    const std::array<Offset, 75> neigh_offsets =
+    {{
+      {0, 0, 0},
+
+      // 14 offsets that are entirely contained within the scaled domain
+      {-1, -1, -1}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}, {1, 1, 0}, {1, 0, 1}, {0, -1, -1},
+      {0, 1, 1}, {-1, 0, -1}, {-1, -1, 0}, {1, 1, 1}, {0, 0, -1}, {0, -1, 0}, {-1, 0, 0},
+
+      // 36 offsets that have a guaranteed intersection with the scaled domain
+      {1, 2, 0}, {1, 0, 2}, {-1, -2, -2}, {2, 1, 0}, {2, 0, 1}, {1, -1, -1}, {0, 1, 2},
+      {-2, -1, -2}, {0, 2, 1}, {-1, 1, -1}, {-2, -2, -1}, {-1, -1, 1}, {1, 1, 2},
+      {-1, -1, -2}, {1, 2, 1}, {0, 1, -1}, {-1, -2, -1}, {0, -1, 1}, {2, 1, 1},
+      {1, 0, -1}, {1, -1, 0}, {-2, -1, -1}, {-1, 0, 1}, {-1, 1, 0}, {1, 2, 2},
+      {-1, 0, -2}, {-1, -2, 0}, {2, 1, 2}, {0, -1, -2}, {2, 2, 1}, {1, 1, -1},
+      {0, -2, -1}, {1, -1, 1}, {-2, -1, 0}, {-2, 0, -1}, {-1, 1, 1},
+
+      // 24 offsets that might have an intersection with the scaled domain (6 of them will)
+      {3, 2, 1}, {2, 1, -1}, {3, 1, 2}, {2, -1, 1}, {1, -1, -2}, {1, -2, -1},
+      {2, 3, 1}, {1, 2, -1}, {1, 3, 2}, {-1, 2, 1}, {-1, 1, -2}, {-2, 1, -1},
+      {2, 1, 3}, {1, -1, 2}, {1, 2, 3}, {-1, 1, 2}, {-1, -2, 1}, {-2, -1, 1},
+      {-1, -2, -3}, {-1, -3, -2}, {-2, -1, -3}, {-3, -1, -2}, {-2, -3, -1}, {-3, -2, -1}
+    }};
+
+    star_lnoffs.reserve(neigh_offsets.size());
+    star_lnoffs.push_back(neigh_offsets[0]);
+
+    DT3 dt3;
+    DT3_VH v0 = dt3.insert(CGAL::ORIGIN);
+    CGAL_assertion(v0 != DT3_VH());
+    v0->info() = 0;
+
+    for(std::size_t oi=1; oi<neigh_offsets.size(); ++oi)
+    {
+      const Offset& lnoff = neigh_offsets[oi];
+
+      Offset lidx;
+      lidx[min_pos] = lnoff[0]; // not a typo, this aligns the neigh offsets with the min/mid/max
+      lidx[mid_pos] = lnoff[1];
+      lidx[max_pos] = lnoff[2];
+
+      // Using gp because it's less likely to produce numerical errors
+      const Offset gidx = lattice_offset_to_grid_offset(lidx);
+      const Point_3 gp = construct_point_from_grid_offset(gidx);
+
+      CGAL_warning_code(const Point_3 lp = construct_point_from_lattice_offset(lidx);)
+      CGAL_warning_code(if(CGAL::squared_distance(gp, lp) > 1e-10 * spans[min_pos])) // might fail with inexact constructions
+      CGAL_warning_code({ std::cout << "WARNING: gp/lp " << gp << " ||| " << lp << std::endl; })
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 8)
+      std::cout << "Lattice Offset: " << lidx[min_pos] << " " << lidx[mid_pos] << " " << lidx[max_pos]
+                << " yields p: " << gp << std::endl;
+#endif
+
+      DT3_VH vh = dt3.insert(gp);
+      CGAL_assertion(vh != DT3_VH() && vh != v0);
+      vh->info() = oi;
+
+      star_lnoffs.push_back(lidx);
+    }
+
+    std::list<DT3_CH> cells;
+    dt3.incident_cells(v0, std::back_inserter(cells));
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 8)
+    std::set<Offset> neighbors_in_use;
+#endif
+
+    star_cells.reserve(cells.size());
+    for(DT3_CH ch : cells)
+    {
+      // To avoid keeping tracks of cells being already created or not, we only insert cells
+      // when they appear in the incident cells of the vertex with the smallest global ID within the call.
+      //
+      // Put the center vertex at pos 0 in the cell so that it is simple to find it
+      constexpr std::array<std::array<int, 4>, 4> permutations =
+      {{
+        { 0, 1, 2, 3 },
+        { 1, 2, 0, 3 },
+        { 2, 0, 1, 3 },
+        { 3, 1, 0, 2 },
+      }};
+
+      const int v0_pos = ch->index(v0);
+      CGAL_postcondition(v0_pos >= 0 && v0_pos <= 3);
+      CGAL_assertion(ch->vertex(permutations[v0_pos][0]) == v0);
+      CGAL_assertion(ch->vertex(permutations[v0_pos][0])->info() == 0);
+
+      star_cells.emplace_back(std::array<std::size_t, 4>{ ch->vertex(permutations[v0_pos][0])->info(),
+                                                          ch->vertex(permutations[v0_pos][1])->info(),
+                                                          ch->vertex(permutations[v0_pos][2])->info(),
+                                                          ch->vertex(permutations[v0_pos][3])->info() });
+
+      CGAL_assertion_code
+      (
+        for(int pos1=0; pos1<4; ++pos1)
+          for(int pos2=0; pos2<4; ++pos2)
+            if(pos1 != pos2)
+              CGAL_assertion(star_cells.back()[pos1] != star_cells.back()[pos2]);
+      )
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 8)
+      for(int i=0; i<4; ++i)
+        neighbors_in_use.insert(neigh_offsets[ch->vertex(i)->info()]);
+#endif
+    }
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 8)
+    std::cout << neighbors_in_use.size() << " lattice offsets in the canonical star:\n";
+    for(const Offset& loff : neighbors_in_use)
+      std::cout << loff[min_pos] << " " << loff[mid_pos] << " " << loff[max_pos] << std::endl;
+
+    std::cout << "as grid offsets:\n";
+    for(const Offset& loff : neighbors_in_use)
+    {
+      const Offset goff = lattice_offset_to_grid_offset(loff);
+      std::cout << goff[min_pos] << " " << goff[mid_pos] << " " << goff[max_pos] << std::endl;
+    }
+#endif
+  };
+
+  compute_periodic_star();
+
+  // @todo might be able to avoid a map by numbering uniquely cells
+  // with a clear vertex_id<->cell_id formula
+  std::map<std::set<std::size_t>, std::pair<Facet, Facet> > neighboring_facets;
+
+  for(int i=0; i<nums_steps[min_pos]; ++i)
+  {
+    for(int j=0; j<nums_steps[mid_pos]; ++j)
+    {
+      for(int k=0; k<nums_steps[max_pos]; ++k)
+      {
+        auto idx_to_global = [&](const Offset& v_idx) -> std::size_t
+        {
+          CGAL_assertion(v_idx[min_pos] >= 0 && v_idx[min_pos] < nums_steps[min_pos]);
+          CGAL_assertion(v_idx[mid_pos] >= 0 && v_idx[mid_pos] < nums_steps[mid_pos]);
+          CGAL_assertion(v_idx[max_pos] >= 0 && v_idx[max_pos] < nums_steps[max_pos]);
+
+          return v_idx[min_pos] * nums_steps[mid_pos] * nums_steps[max_pos]
+               + v_idx[mid_pos] * nums_steps[max_pos]
+               + v_idx[max_pos];
+        };
+
+        auto add_facet = [&](const std::set<std::size_t>& face_vertices_ids, const Facet& f) -> void
+        {
+          CGAL_assertion(f.first != Cell_handle() && f.second >= 0 && f.second <= 3);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+          std::cout << std::endl;
+          std::cout << "Add facet";
+          for(const std::size_t vi : face_vertices_ids)
+            std::cout << " " << vi;
+          std::cout << std::endl;
+
+          std::cout << "Positions\n";
+          for(const std::size_t vi : face_vertices_ids)
+            std::cout << vertices[vi]->point() << std::endl;
+          std::cout << std::endl;
+#endif
+
+          auto insertion_result = neighboring_facets.emplace(face_vertices_ids, std::make_pair(f, Facet()));
+          if(!insertion_result.second) // already exists in the facet map
+          {
+            std::pair<Facet, Facet>& fp = insertion_result.first->second;
+            CGAL_assertion(fp.first.first != Cell_handle()); // one existing neighbor since insertion failed
+            CGAL_assertion(fp.second.first == Cell_handle()); // two neighbors at most
+            CGAL_assertion(fp.second.first != f.first); // neighbors are different
+            fp.second = f;
+
+            CGAL_assertion(fp.second.first != Cell_handle()); // properly set up
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+            std::cout << "facets" << std::endl;
+            std::cout << &*(fp.first.first->vertex((fp.first.second+1)%4)) << std::endl;
+            std::cout << &*(fp.first.first->vertex((fp.first.second+2)%4)) << std::endl;
+            std::cout << &*(fp.first.first->vertex((fp.first.second+3)%4)) << std::endl;
+            std::cout << &*(fp.second.first->vertex((fp.second.second+1)%4)) << std::endl;
+            std::cout << &*(fp.second.first->vertex((fp.second.second+2)%4)) << std::endl;
+            std::cout << &*(fp.second.first->vertex((fp.second.second+3)%4)) << std::endl;
+            std::cout << fp.first.first->vertex((fp.first.second+1)%4)->point() << std::endl;
+            std::cout << fp.first.first->vertex((fp.first.second+2)%4)->point() << std::endl;
+            std::cout << fp.first.first->vertex((fp.first.second+3)%4)->point() << std::endl;
+            std::cout << fp.second.first->vertex((fp.second.second+1)%4)->point() << std::endl;
+            std::cout << fp.second.first->vertex((fp.second.second+2)%4)->point() << std::endl;
+            std::cout << fp.second.first->vertex((fp.second.second+3)%4)->point() << std::endl;
+
+            std::set<Vertex_handle> f1_vs = { fp.first.first->vertex((fp.first.second+1)%4),
+                                              fp.first.first->vertex((fp.first.second+2)%4),
+                                              fp.first.first->vertex((fp.first.second+3)%4) };
+            std::set<Vertex_handle> f2_vs = { fp.second.first->vertex((fp.second.second+1)%4),
+                                              fp.second.first->vertex((fp.second.second+2)%4),
+                                              fp.second.first->vertex((fp.second.second+3)%4) };
+            CGAL_assertion(f1_vs == f2_vs);
+#endif
+
+            // neighbors could be set here instead of storing pairs and doing it later,
+            // but it doesn't cost much and it enables checking that every face has exactly
+            // two incident cells
+          }
+
+          const std::pair<Facet, Facet>& fp = insertion_result.first->second;
+          CGAL_USE(fp);
+          CGAL_assertion(fp.first.first != Cell_handle()); // properly set up
+        };
+
+        auto add_cell = [&](const std::array<std::size_t, 4>& star_cell) -> Cell_handle
+        {
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 8)
+          std::cout << std::endl;
+          std::cout << "-- Considering cell " << star_cell[0] << " " << star_cell[1] << " " << star_cell[2] << " " << star_cell[3] << std::endl;
+#endif
+
+          // convert to lattice offset, add (lattice) neighbor offset, convert back to grid
+
+          Offset idx;
+          idx[min_pos] = i;
+          idx[mid_pos] = j;
+          idx[max_pos] = k;
+
+          const Offset lidx = grid_offset_to_lattice_offset(idx);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+          std::cout << "(Recall that min/mid/max are: " << min_pos << " " << mid_pos << " " << max_pos << ")" << std::endl;
+          std::cout << "(Recall that steps are: " << steps[min_pos] << " " << steps[mid_pos] << " " << steps[max_pos] << ")" << std::endl;
+          std::cout << "(Recall that nums steps are: " << nums_steps[min_pos] << " " << nums_steps[mid_pos] << " " << nums_steps[max_pos] << ")" << std::endl;
+          std::cout << "grid idx: " << idx[min_pos] << " " << idx[mid_pos] << " " << idx[max_pos] << std::endl;
+          std::cout << "lattice idx: " << lidx[min_pos] << " " << lidx[mid_pos] << " " << lidx[max_pos] << std::endl;
+#endif
+
+          std::array<std::size_t, 4> glob_idx;
+          std::array<CGAL::Periodic_3_offset_3, 4> offsets { };
+
+          for(int pos=0; pos<4; ++pos)
+          {
+            const Offset& lnoff = star_lnoffs[star_cell[pos]];
+
+            // lnoff has already been shuffled to match min/mid/max
+            Offset lnidx = lidx + lnoff;
+
+            // canonical (to be) ijk of the neighbor
+            Offset cidx = lattice_offset_to_grid_offset(lnidx);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+            std::cout << std::endl;
+            std::cout << "lattice neighbor offset[" << pos << "] = " << lnoff[min_pos] << " " << lnoff[mid_pos] << " " << lnoff[max_pos] << std::endl;
+            std::cout << "lattice neighbor idx[" << pos << "] = " << lnidx[min_pos] << " " << lnidx[mid_pos] << " " << lnidx[max_pos] << std::endl;
+            std::cout << "canonical idx[" << pos << "] = " << cidx[min_pos] << " " << cidx[mid_pos] << " " << cidx[max_pos] << " (before)" << std::endl;
+#endif
+
+            CGAL_warning_code(const Point_3 gp = construct_point_from_grid_offset(cidx);)
+            CGAL_warning_code(const Point_3 lp = construct_point_from_lattice_offset(lnidx);)
+            CGAL_warning_code(if(CGAL::squared_distance(gp, lp) > 1e-10 * spans[min_pos])) // might fail with inexact constructions
+            CGAL_warning_code({ std::cout << "WARNING: gp/lp " << gp << " ||| " << lp << std::endl; })
+
+            // canonicalize using indices to avoid constructions
+            for(int l=0; l<3; ++l)
+            {
+              offsets[pos][l] = 0;
+
+              while(cidx[l] >= nums_steps[l])
+              {
+                cidx[l] -= nums_steps[l];
+                offsets[pos][l] += 1;
+              }
+
+              while(cidx[l] < 0)
+              {
+                cidx[l] += nums_steps[l];
+                offsets[pos][l] -= 1;
+              }
+            }
+
+            glob_idx[pos] = idx_to_global(cidx);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+            std::cout << "final canon idx[" << pos << "] = " << cidx[min_pos] << " " << cidx[mid_pos] << " " << cidx[max_pos] << std::endl;
+            std::cout << "periodic offset[" << pos << "] = " << offsets[pos] << std::endl;
+            std::cout << "Glob ID[" << pos << "] = " << glob_idx[pos]<< std::endl;
+            std::cout << "Position[" << pos << "] = " << vertices[glob_idx[pos]]->point() << std::endl;
+#endif
+
+            CGAL_assertion(glob_idx[pos] < nv);
+          }
+
+          CGAL_assertion_code
+          (
+            for(int pos1=0; pos1<4; ++pos1)
+              for(int pos2=0; pos2<4; ++pos2)
+                if(pos1 != pos2)
+                  CGAL_assertion(glob_idx[pos1] != glob_idx[pos2]);
+          )
+
+          // trick to add facets only once per cell; can't really do that for all operations
+          // due to cells on the border (or, it could be done, but then
+          // is_not_on_top_right_up_corner must be reworked)
+          //
+          // The trick works because glob_idx[0] is always the center vertex
+          if(glob_idx[0] != (std::min)({glob_idx[0], glob_idx[1], glob_idx[2], glob_idx[3]}))
+            return Cell_handle();
+
+          // here and below, it's a new cell
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+          std::cout << "\nCell position around star center" << std::endl;
+          std::cout << construct_point(vertices[glob_idx[0]]->point(), offsets[0]) << std::endl;
+          std::cout << construct_point(vertices[glob_idx[1]]->point(), offsets[1]) << std::endl;
+          std::cout << construct_point(vertices[glob_idx[2]]->point(), offsets[2]) << std::endl;
+          std::cout << construct_point(vertices[glob_idx[3]]->point(), offsets[3]) << std::endl;
+#endif
+
+          CGAL_assertion(orientation(construct_point(vertices[glob_idx[0]]->point(), offsets[0]),
+                                     construct_point(vertices[glob_idx[1]]->point(), offsets[1]),
+                                     construct_point(vertices[glob_idx[2]]->point(), offsets[2]),
+                                     construct_point(vertices[glob_idx[3]]->point(), offsets[3])) == CGAL::POSITIVE);
+
+          Cell_handle ch = tds().create_cell();
+          CGAL_assertion(ch != Cell_handle());
+
+          CGAL_assertion_code(for(int pos=0; pos<4; ++pos) {)
+          CGAL_assertion(vertices.at(glob_idx[pos]) != Vertex_handle());
+          CGAL_assertion_code(})
+
+          ch->set_vertices(vertices[glob_idx[0]], vertices[glob_idx[1]], vertices[glob_idx[2]], vertices[glob_idx[3]]);
+          set_offsets(ch, offsets[0], offsets[1], offsets[2], offsets[3]);
+
+          for(int pos=0; pos<4; ++pos)
+            vertices[glob_idx[pos]]->set_cell(ch);
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 16)
+          std::cout << "\nCell position in P3T3" << std::endl;
+          std::cout << point(ch, 0) << std::endl;
+          std::cout << point(ch, 1) << std::endl;
+          std::cout << point(ch, 2) << std::endl;
+          std::cout << point(ch, 3) << std::endl;
+#endif
+
+          // add faces to the neighboring map
+          add_facet(std::set<std::size_t>{glob_idx[0], glob_idx[1], glob_idx[2]}, Facet(ch, 3));
+          add_facet(std::set<std::size_t>{glob_idx[0], glob_idx[1], glob_idx[3]}, Facet(ch, 2));
+          add_facet(std::set<std::size_t>{glob_idx[0], glob_idx[2], glob_idx[3]}, Facet(ch, 1));
+          add_facet(std::set<std::size_t>{glob_idx[1], glob_idx[2], glob_idx[3]}, Facet(ch, 0));
+
+          return ch;
+        };
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 4)
+        std::cout << std::endl;
+        std::cout << " ====================== VERTEX " << i << " " << j << " " << k << " ==============" << std::endl;
+#endif
+
+        for(const auto& star_cell : star_cells)
+          add_cell(star_cell);
+      }
+    }
+  }
+
+  // build neighboring info
+  for(const auto& e : neighboring_facets)
+  {
+    const std::pair<Facet, Facet>& fp = e.second;
+    const Cell_handle ch_a = fp.first.first, ch_b = fp.second.first;
+    const int s_a = fp.first.second, s_b = fp.second.second;
+
+    CGAL_assertion(ch_a != Cell_handle() && ch_b != Cell_handle() && ch_a != ch_b);
+    CGAL_assertion(ch_a->neighbor(s_a) == Cell_handle() && ch_b->neighbor(s_b) == Cell_handle());
+
+    ch_a->set_neighbor(s_a, ch_b);
+    ch_b->set_neighbor(s_b, ch_a);
+  }
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 8)
+  std::cout << "facet map size: " << neighboring_facets.size() << std::endl;
+  std::cout << number_of_vertices() << " nv" << std::endl;
+  std::cout << number_of_cells() << " nc" << std::endl;
+#endif
+
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 2)
+  CGAL_postcondition(is_valid(true));
+#else
+  CGAL_postcondition(is_valid());
+#endif
+
+  return vertices;
+}
+
+#endif // CGAL_INCLUDE_FROM_PERIODIC_3_TRIANGULATION_3_H

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/Robust_periodic_weighted_circumcenter_traits_3.h

@@ -67,16 +67,45 @@ public:
                      const Offset& o_p, const Offset& o_q,
                      const Offset& o_r, const Offset& o_s) const
   {
-    typename Base_traits::Construct_weighted_point_3 p2wp =
+    typename Base_traits::Compute_weight_3 weight =
+               base_traits.compute_weight_3_object();
+    typename Base_traits::Construct_point_3 point =
+               base_traits.construct_point_3_object();
+    typename Base_traits::Construct_weighted_point_3 weighted_point =
                base_traits.construct_weighted_point_3_object();
     typename Base_traits::Construct_weighted_circumcenter_3 cwc =
                base_traits.construct_weighted_circumcenter_3_object();
+    typename Base_traits::Orientation_3 orientation =
+               base_traits.orientation_3_object();
     typename Base_traits::Power_side_of_oriented_power_sphere_3 ps =
                base_traits.power_side_of_oriented_power_sphere_3_object();
 
-    Point_3 c = cwc(p, q, r, s, o_p, o_q, o_r, o_s);
+    // Calling Construct_weighted_circumcenter_3(p,q,r,s,o_p,o_q,o_r,o_s) will construct
+    // the 3D weighted points (in Functor_with_offset_adaptor) but that is bad because
+    // if we're unlucky pqrs is almost flat and the orientation might become wrong
+    // (i.e. != POSITIVE) when we construct the Euclidean, non-periodic 3D representations
+    // of the points
+    const Point_3 euc_p = point(p, o_p);
+    const Point_3 euc_q = point(q, o_q);
+    const Point_3 euc_r = point(r, o_r);
+    const Point_3 euc_s = point(s, o_s);
 
-    if(ps(p, q, r, s, p2wp(c), o_p, o_q, o_r, o_s, Offset()) != ON_POSITIVE_SIDE)
+    bool needs_exact = false;
+    Point_3 c;
+    if(orientation(euc_p, euc_q, euc_r, euc_s) != POSITIVE) // see comment above
+    {
+      needs_exact = true;
+    }
+    else
+    {
+      c = cwc(weighted_point(euc_p, weight(p)),
+              weighted_point(euc_q, weight(q)),
+              weighted_point(euc_r, weight(r)),
+              weighted_point(euc_s, weight(s)));
+    }
+
+    if(needs_exact ||
+       ps(p, q, r, s, weighted_point(c), o_p, o_q, o_r, o_s, Offset()) != ON_POSITIVE_SIDE)
     {
       // switch to exact
       typedef typename Base_traits::Kernel                         K;
@@ -95,28 +124,55 @@ public:
               to_exact(p), to_exact(q), to_exact(r), to_exact(s),
               o_p, o_q, o_r, o_s));
 
-      CGAL_assertion(ps(p, q, r, s, p2wp(c), o_p, o_q, o_r, o_s, Offset()) == ON_POSITIVE_SIDE);
+      CGAL_assertion(ps(p, q, r, s, weighted_point(c), o_p, o_q, o_r, o_s, Offset()) == ON_POSITIVE_SIDE);
     }
 
     return c;
   }
 
-  Point_3 operator()(const Weighted_point_3 & p,
-                     const Weighted_point_3 & q,
-                     const Weighted_point_3 & r,
-                     const Offset& o_p, const Offset& o_q,
+  Point_3 operator()(const Weighted_point_3& p,
+                     const Weighted_point_3& q,
+                     const Weighted_point_3& r,
+                     const Offset& o_p,
+                     const Offset& o_q,
                      const Offset& o_r) const
   {
-    typename Base_traits::Construct_weighted_point_3 p2wp =
+    typename Base_traits::Compute_weight_3 weight =
+               base_traits.compute_weight_3_object();
+    typename Base_traits::Construct_point_3 point =
+               base_traits.construct_point_3_object();
+    typename Base_traits::Construct_weighted_point_3 weighted_point =
                base_traits.construct_weighted_point_3_object();
     typename Base_traits::Construct_weighted_circumcenter_3 cwc =
                base_traits.construct_weighted_circumcenter_3_object();
+    typename Base_traits::Collinear_3 collinear =
+               base_traits.collinear_3_object();
     typename Base_traits::Power_side_of_bounded_power_sphere_3 ps =
                base_traits.power_side_of_bounded_power_sphere_3_object();
 
-    Point_3 c = cwc(p, q, r, o_p, o_q, o_r);
+    // Calling Construct_weighted_circumcenter_3(p,q,r,o_p,o_q,o_r) will construct
+    // the 3D weighted points (in Functor_with_offset_adaptor) but that is bad because
+    // if we're unlucky pqr is almost collinear and it becomes so when we construct the Euclidean,
+    // non-periodic 3D representations of the points
+    const Point_3 euc_p = point(p, o_p);
+    const Point_3 euc_q = point(q, o_q);
+    const Point_3 euc_r = point(r, o_r);
 
-    if(ps(p, q, r, p2wp(c), o_p, o_q, o_r, Offset()) != ON_BOUNDED_SIDE)
+    bool needs_exact = false;
+    Point_3 c;
+    if(collinear(euc_p, euc_q, euc_r)) // see comment above
+    {
+      needs_exact = true;
+    }
+    else
+    {
+      c = cwc(weighted_point(euc_p, weight(p)),
+              weighted_point(euc_q, weight(q)),
+              weighted_point(euc_r, weight(r)));
+    }
+
+    if(needs_exact ||
+       ps(p, q, r, weighted_point(c), o_p, o_q, o_r, Offset()) != ON_BOUNDED_SIDE)
     {
       // switch to exact
       typedef typename Base_traits::Kernel                         K;
@@ -134,26 +190,51 @@ public:
             etraits.construct_weighted_circumcenter_3_object()(
               to_exact(p), to_exact(q), to_exact(r), o_p, o_q, o_r));
 
-      CGAL_assertion(ps(p, q, r, p2wp(c), o_p, o_q, o_r, Offset()) == ON_BOUNDED_SIDE);
+      CGAL_assertion(ps(p, q, r, weighted_point(c), o_p, o_q, o_r, Offset()) == ON_BOUNDED_SIDE);
     }
 
     return c;
   }
 
-  Point_3 operator()(const Weighted_point_3 & p,
-                     const Weighted_point_3 & q,
-                     const Offset& o_p, const Offset& o_q) const
+  Point_3 operator()(const Weighted_point_3& p,
+                     const Weighted_point_3& q,
+                     const Offset& o_p,
+                     const Offset& o_q) const
   {
-    typename Base_traits::Construct_weighted_point_3 p2wp =
+    typename Base_traits::Equal_3 equal =
+               base_traits.equal_3_object();
+    typename Base_traits::Compute_weight_3 weight =
+               base_traits.compute_weight_3_object();
+    typename Base_traits::Construct_point_3 point =
+               base_traits.construct_point_3_object();
+    typename Base_traits::Construct_weighted_point_3 weighted_point =
                base_traits.construct_weighted_point_3_object();
     typename Base_traits::Construct_weighted_circumcenter_3 cwc =
                base_traits.construct_weighted_circumcenter_3_object();
     typename Base_traits::Power_side_of_bounded_power_sphere_3 ps =
                base_traits.power_side_of_bounded_power_sphere_3_object();
 
-    Point_3 c = cwc(p, q, o_p, o_q);
+    // Calling Construct_weighted_circumcenter_3(p,q,o_p,o_q) will construct
+    // the 3D weighted points (in Functor_with_offset_adaptor) but that is bad because
+    // if we're unlucky p and q become equal when we construct the Euclidean,
+    // non-periodic 3D representations of the points
+    const Point_3 euc_p = point(p, o_p);
+    const Point_3 euc_q = point(q, o_q);
 
-    if(ps(p, q, p2wp(c), o_p, o_q, Offset()) != ON_BOUNDED_SIDE)
+    bool needs_exact = false;
+    Point_3 c;
+    if(equal(euc_p, euc_q)) // see comment above
+    {
+      needs_exact = true;
+    }
+    else
+    {
+      c = cwc(weighted_point(euc_p, weight(p)),
+              weighted_point(euc_q, weight(q)));
+    }
+
+    if(needs_exact ||
+       ps(p, q, weighted_point(c), o_p, o_q, Offset()) != ON_BOUNDED_SIDE)
     {
       // switch to exact
       typedef typename Base_traits::Kernel                         K;
@@ -171,7 +252,7 @@ public:
             etraits.construct_weighted_circumcenter_3_object()(
               to_exact(p), to_exact(q), o_p, o_q));
 
-      CGAL_assertion(ps(p, q, p2wp(c), o_p, o_q, Offset()) == ON_BOUNDED_SIDE);
+      CGAL_assertion(ps(p, q, weighted_point(c), o_p, o_q, Offset()) == ON_BOUNDED_SIDE);
     }
 
     return c;

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_3/internal/canonicalize_helper.h

@@ -39,60 +39,13 @@ namespace CGAL {
 namespace P3T3 { // can't name it Periodic_3_triangulation_3 because it's already a class...
 namespace internal {
 
+// Given a point `p` in space, compute its offset `o` with respect to the canonical
+// domain (i.e., p + o * d is in the canonical domain) and returns `(p, o)`
 template <typename Gt_>
 std::pair<typename Gt_::Point_3, typename Gt_::Periodic_3_offset_3>
-construct_periodic_point_exact(const typename Gt_::Point_3& p,
-                               const Gt_& gt)
-{
-  typedef Gt_                                                  Geom_traits;
-  typedef typename Geom_traits::Periodic_3_offset_3            Offset;
-  typedef typename Geom_traits::Iso_cuboid_3                   Iso_cuboid;
-
-  const Iso_cuboid& domain = gt.get_domain();
-
-  typedef typename Geom_traits::Kernel                         K;
-  typedef typename Exact_kernel_selector<K>::Exact_kernel      EK;
-  typedef typename Exact_kernel_selector<K>::C2E               C2E;
-
-  C2E to_exact;
-
-  typedef Periodic_3_triangulation_traits_3<EK> Exact_traits;
-  Exact_traits etraits(to_exact(domain));
-
-  Offset transl(0, 0, 0);
-  typename EK::Point_3 ep = to_exact(p);
-  typename EK::Point_3 dp;
-
-  const typename EK::Iso_cuboid_3& exact_domain = etraits.get_domain();
-
-  while(true) /* while not in */
-  {
-    dp = etraits.construct_point_3_object()(ep, transl);
-
-    if(dp.x() < exact_domain.xmin())
-      transl.x() += 1;
-    else if(dp.y() < exact_domain.ymin())
-      transl.y() += 1;
-    else if(dp.z() < exact_domain.zmin())
-      transl.z() += 1;
-    else if(!(dp.x() < exact_domain.xmax()))
-      transl.x() -= 1;
-    else if(!(dp.y() < exact_domain.ymax()))
-      transl.y() -= 1;
-    else if(!(dp.z() < exact_domain.zmax()))
-      transl.z() -= 1;
-    else
-      break;
-  }
-
-  return std::make_pair(p, transl);
-}
-
-// Given a point `p` in space, compute its offset `o` with respect
-// to the canonical instance and returns (p, o)
-template <typename Gt_>
-std::pair<typename Gt_::Point_3, typename Gt_::Periodic_3_offset_3>
-construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue, const Gt_& gt)
+construct_periodic_point(const typename Gt_::Point_3& p,
+                         bool& encountered_issue,
+                         const Gt_& gt)
 {
   typedef Gt_                                                  Geom_traits;
   typedef typename Geom_traits::Point_3                        Point;
@@ -101,13 +54,18 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
 
   const Iso_cuboid& domain = gt.get_domain();
 
+  // Use these rather than Construct_point_3 to avoid construction inaccuracies
+  typename Geom_traits::Compare_x_3 cmp_x3 = gt.compare_x_3_object();
+  typename Geom_traits::Compare_y_3 cmp_y3 = gt.compare_y_3_object();
+  typename Geom_traits::Compare_z_3 cmp_z3 = gt.compare_z_3_object();
+
   // Check if p lies within the domain. If not, translate.
   if(!(p.x() < domain.xmin()) && p.x() < domain.xmax() &&
      !(p.y() < domain.ymin()) && p.y() < domain.ymax() &&
      !(p.z() < domain.zmin()) && p.z() < domain.zmax())
+  {
     return std::make_pair(p, Offset());
-
-  typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
+  }
 
   // Numerical approximations might create inconsistencies between the constructions
   // and the comparisons. For example in a cubic domain of size 2:
@@ -120,9 +78,10 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
   //
   // If this is happening the 'Last_change' enum will break this infinite
   // loop and return the wrong point and the 'encountered_issue' bool will be
-  // set to 'true'. An exact version of this function should then be called.
+  // set to 'true'. An exact version of this function is then be called.
 
-  enum Last_change {
+  enum Last_change
+  {
     NO_LAST_CHANGE,
     INCREASED_X, DECREASED_X, INCREASED_Y, DECREASED_Y, INCREASED_Z, DECREASED_Z
   };
@@ -131,13 +90,14 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
   bool in = false;
 
   Offset transl(0, 0, 0);
-  Point dp;
+  const Offset null_off(0, 0, 0);
+
+  Point domain_m(domain.xmin(), domain.ymin(), domain.zmin());
+  Point domain_M(domain.xmax(), domain.ymax(), domain.zmax());
 
   while(!in)
   {
-    dp = cp(p, transl);
-
-    if(dp.x() < domain.xmin())
+    if(cmp_x3(p, domain_m, transl, null_off) == SMALLER)
     {
       if(lc == DECREASED_X) // stuck in a loop
         break;
@@ -145,7 +105,7 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
       lc = INCREASED_X;
       transl.x() += 1;
     }
-    else if(dp.y() < domain.ymin())
+    else if(cmp_y3(p, domain_m, transl, null_off) == SMALLER)
     {
       if(lc == DECREASED_Y) // stuck in a loop
         break;
@@ -153,7 +113,7 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
       lc = INCREASED_Y;
       transl.y() += 1;
     }
-    else if(dp.z() < domain.zmin())
+    else if(cmp_z3(p, domain_m, transl, null_off) == SMALLER)
     {
       if(lc == DECREASED_Z) // stuck in a loop
         break;
@@ -161,7 +121,7 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
       lc = INCREASED_Z;
       transl.z() += 1;
     }
-    else if(!(dp.x() < domain.xmax()))
+    else if(!(cmp_x3(p, domain_M, transl, null_off) == SMALLER))
     {
       if(lc == INCREASED_X) // stuck in a loop
         break;
@@ -169,7 +129,7 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
       lc = DECREASED_X;
       transl.x() -= 1;
     }
-    else if(!(dp.y() < domain.ymax()))
+    else if(!(cmp_y3(p, domain_M, transl, null_off) == SMALLER))
     {
       if(lc == INCREASED_Y) // stuck in a loop
         break;
@@ -177,7 +137,7 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
       lc = DECREASED_Y;
       transl.y() -= 1;
     }
-    else if(!(dp.z() < domain.zmax()))
+    else if(!(cmp_z3(p, domain_M, transl, null_off) == SMALLER))
     {
       if(lc == INCREASED_Z) // stuck in a loop
         break;
@@ -193,157 +153,79 @@ construct_periodic_point(const typename Gt_::Point_3& p, bool& encountered_issue
 
   std::pair<Point, Offset> pp(p, transl);
 
-  if(dp.x() < domain.xmin() || !(dp.x() < domain.xmax()) ||
-     dp.y() < domain.ymin() || !(dp.y() < domain.ymax()) ||
-     dp.z() < domain.zmin() || !(dp.z() < domain.zmax()))
+  if(cmp_x3(p, domain_m, transl, null_off) == SMALLER || // < min
+     cmp_y3(p, domain_m, transl, null_off) == SMALLER ||
+     cmp_z3(p, domain_m, transl, null_off) == SMALLER ||
+     !(cmp_x3(p, domain_M, transl, null_off) == SMALLER) || // >= max
+     !(cmp_y3(p, domain_M, transl, null_off) == SMALLER) ||
+     !(cmp_z3(p, domain_M, transl, null_off) == SMALLER))
   {
     encountered_issue = true;
-    pp = construct_periodic_point_exact(p, gt);
   }
 
   return pp;
 }
 
 template <typename Gt_>
-bool
-is_point_too_close_to_border(const std::pair<typename Gt_::Point_3,
-                                             typename Gt_::Periodic_3_offset_3>& pbp,
-                             const Gt_& gt)
+typename Gt_::Point_3
+constrain_to_canonical_domain(const typename Gt_::Point_3& p,
+                              const Gt_& gt)
 {
   typedef Gt_                                         Geom_traits;
   typedef typename Geom_traits::FT                    FT;
-  typedef typename Geom_traits::Point_3               Bare_point;
   typedef typename Geom_traits::Iso_cuboid_3          Iso_cuboid;
 
   typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
 
-  const Bare_point p = cp(pbp.first /*point*/, pbp.second /*offset*/);
-  const FT px = p.x();
-  const FT py = p.y();
-  const FT pz = p.z();
-
   const Iso_cuboid& domain = gt.get_domain();
-  const FT dxm = domain.xmin();
-  const FT dym = domain.ymin();
-  const FT dzm = domain.zmin();
-  const FT dxM = domain.xmax();
-  const FT dyM = domain.ymax();
-  const FT dzM = domain.zmax();
+  FT x = p.x();
+  FT y = p.y();
+  FT z = p.z();
 
-  // simply comparing to FT::epsilon() is probably not completely satisfactory
-  const FT eps = std::numeric_limits<FT>::epsilon();
+  if(p.x() < domain.xmin() || p.x() >= domain.xmax())
+    x = domain.xmin();
+  if(p.y() < domain.ymin() || p.y() >= domain.ymax())
+    y = domain.ymin();
+  if(p.z() < domain.zmin() || p.z() >= domain.zmax())
+    z = domain.zmin();
 
-  FT diff = CGAL::abs(px - dxm);
-  if(diff < eps && diff > 0) return true;
-  diff = CGAL::abs(px - dxM);
-  if(diff < eps && diff > 0) return true;
-  diff = CGAL::abs(py - dym);
-  if(diff < eps && diff > 0) return true;
-  diff = CGAL::abs(py - dyM);
-  if(diff < eps && diff > 0) return true;
-  diff = CGAL::abs(pz - dzm);
-  if(diff < eps && diff > 0) return true;
-  diff = CGAL::abs(pz - dzM);
-  if(diff < eps && diff > 0) return true;
-  return false;
-}
-
-template <typename Gt_>
-typename Gt_::Point_3
-snap_to_domain_border(const typename Gt_::Point_3& p, const Gt_& gt)
-{
-  typedef Gt_                                         Geom_traits;
-  typedef typename Geom_traits::FT                    FT;
-  typedef typename Geom_traits::Iso_cuboid_3          Iso_cuboid;
-
-  const FT px = p.x();
-  const FT py = p.y();
-  const FT pz = p.z();
-  FT sx = px, sy = py, sz = pz;
-
-  const Iso_cuboid& domain = gt.get_domain();
-  const FT dxm = domain.xmin();
-  const FT dym = domain.ymin();
-  const FT dzm = domain.zmin();
-  const FT dxM = domain.xmax();
-  const FT dyM = domain.ymax();
-  const FT dzM = domain.zmax();
-
-  // simply comparing to FT::epsilon() is probably not completely satisfactory
-  const FT eps = std::numeric_limits<FT>::epsilon();
-
-  if(CGAL::abs(px - dxm) < eps) sx = dxm;
-  if(CGAL::abs(px - dxM) < eps) sx = dxM;
-  if(CGAL::abs(py - dym) < eps) sy = dym;
-  if(CGAL::abs(py - dyM) < eps) sy = dyM;
-  if(CGAL::abs(pz - dzm) < eps) sz = dzm;
-  if(CGAL::abs(pz - dzM) < eps) sz = dzM;
-
-  return gt.construct_point_3_object()(sx, sy, sz);
-}
-
-template <typename Gt_>
-typename Gt_::Weighted_point_3
-snap_to_domain_border(const typename Gt_::Weighted_point_3& p,
-                      const Gt_& gt)
-{
-  typedef Gt_                                         Geom_traits;
-  typedef typename Geom_traits::Point_3               Bare_point;
-
-  typename Geom_traits::Compute_weight_3 cw = gt.compute_weight_3_object();
-
-  const Bare_point snapped_p = snap_to_domain_border(gt.construct_point_3_object()(p), gt);
-
-  return gt.construct_weighted_point_3_object()(snapped_p, cw(p));
+  return cp(x, y, z);
 }
 
 /// transform a bare point (living anywhere in space) into the canonical
 /// instance of the same bare point that lives inside the base domain
 template <typename Gt_>
 typename Gt_::Point_3
-robust_canonicalize_point(const typename Gt_::Point_3& p, const Gt_& gt)
+robust_canonicalize_point(const typename Gt_::Point_3& p,
+                          const Gt_& gt)
 {
   typedef Gt_                                         Geom_traits;
   typedef typename Geom_traits::Point_3               Bare_point;
   typedef typename Geom_traits::Periodic_3_offset_3   Offset;
   typedef typename Geom_traits::Iso_cuboid_3          Iso_cuboid;
 
+  typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
+
   const Iso_cuboid& domain = gt.get_domain();
   if(p.x() >= domain.xmin() && p.x() < domain.xmax() &&
      p.y() >= domain.ymin() && p.y() < domain.ymax() &&
      p.z() >= domain.zmin() && p.z() < domain.zmax())
     return p;
 
-  bool should_snap = false;
-  std::pair<Bare_point, Offset> pbp = construct_periodic_point(p, should_snap, gt);
-
-  if(!should_snap)
-  {
-    // Even if there is no issue while constructing the canonical point,
-    // snap the point if it's too close to a border of the domain
-    should_snap = is_point_too_close_to_border(pbp, gt);
-  }
-
-  if(should_snap)
-  {
-    Bare_point sp = snap_to_domain_border(p, gt);
-
-    // might have snapped to a 'max' of the domain, which is not in the domain
-    // note: we could snap to 'min' all the time in 'snap_to_domain_border'
-    // but this is clearer like that (and costs very little since we should
-    // not have to use exact computations too often)
-    return robust_canonicalize_point(sp, gt);
-  }
-
-  typename Geom_traits::Construct_point_3 cp = gt.construct_point_3_object();
-
+  bool encountered_issue = false;
+  std::pair<Bare_point, Offset> pbp = construct_periodic_point(p, encountered_issue, gt);
   Bare_point canonical_p = cp(pbp.first /*point*/, pbp.second /*offset*/);
-  CGAL_postcondition( !(canonical_p.x() < domain.xmin()) &&
-                       (canonical_p.x() < domain.xmax()) );
-  CGAL_postcondition( !(canonical_p.y() < domain.ymin()) &&
-                       (canonical_p.y() < domain.ymax()) );
-  CGAL_postcondition( !(canonical_p.z() < domain.zmin()) &&
-                       (canonical_p.z() < domain.zmax()) );
+
+  if(encountered_issue)
+  {
+    // If we encountered an issue, there's no guarantee that the double construction gives a point
+    // in the domain (even if we computed it exactly beforehand). So, forcefully put it into the domain.
+    canonical_p = constrain_to_canonical_domain(canonical_p, gt);
+  }
+
+  CGAL_postcondition( !(canonical_p.x() < domain.xmin()) && (canonical_p.x() < domain.xmax()));
+  CGAL_postcondition( !(canonical_p.y() < domain.ymin()) && (canonical_p.y() < domain.ymax()));
+  CGAL_postcondition( !(canonical_p.z() < domain.zmin()) && (canonical_p.z() < domain.zmax()));
 
   return canonical_p;
 }
@@ -352,7 +234,8 @@ robust_canonicalize_point(const typename Gt_::Point_3& p, const Gt_& gt)
 /// instance of the same weighted point that lives inside the base domain
 template <typename Gt_>
 typename Gt_::Weighted_point_3
-robust_canonicalize_point(const typename Gt_::Weighted_point_3& wp, const Gt_& gt)
+robust_canonicalize_point(const typename Gt_::Weighted_point_3& wp,
+                          const Gt_& gt)
 {
   typedef Gt_                                         Geom_traits;
   typedef typename Geom_traits::Point_3               Bare_point;

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_ds_cell_base_3.h

@@ -189,6 +189,11 @@ public:
 
   void set_offsets(int o0,int o1,int o2,int o3)
   {
+    CGAL_precondition(o0 >= 0 && o0 <= 7);
+    CGAL_precondition(o1 >= 0 && o1 <= 7);
+    CGAL_precondition(o2 >= 0 && o2 <= 7);
+    CGAL_precondition(o3 >= 0 && o3 <= 7);
+
     off = 0;
     // The following explicit cast are needed according to the Intel
     // Compiler version 12.

Periodic_3_triangulation_3/include/CGAL/Periodic_3_triangulation_traits_3.h

@@ -78,6 +78,12 @@ public:
       Construct_point_3;
 
   // Triangulation predicates
+  typedef Functor_with_offset_points_adaptor_3<Self, typename Kernel::Compare_x_3>
+      Compare_x_3;
+  typedef Functor_with_offset_points_adaptor_3<Self, typename Kernel::Compare_y_3>
+      Compare_y_3;
+  typedef Functor_with_offset_points_adaptor_3<Self, typename Kernel::Compare_z_3>
+      Compare_z_3;
   typedef Functor_with_offset_points_adaptor_3<Self, typename Kernel::Compare_xyz_3>
       Compare_xyz_3;
   typedef Functor_with_offset_points_adaptor_3<Self, typename Kernel::Orientation_3>
@@ -96,6 +102,15 @@ public:
     return Construct_point_3(&_domain, this->Kernel::construct_point_3_object());
   }
 
+  Compare_x_3 compare_x_3_object() const {
+    return Compare_x_3(this->Kernel::compare_x_3_object(), construct_point_3_object());
+  }
+  Compare_y_3 compare_y_3_object() const {
+    return Compare_y_3(this->Kernel::compare_y_3_object(), construct_point_3_object());
+  }
+  Compare_z_3 compare_z_3_object() const {
+    return Compare_z_3(this->Kernel::compare_z_3_object(), construct_point_3_object());
+  }
   Compare_xyz_3 compare_xyz_3_object() const {
     return Compare_xyz_3(this->Kernel::compare_xyz_3_object(), construct_point_3_object());
   }

Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/CMakeLists.txt

@@ -28,3 +28,4 @@ create_single_source_cgal_program(
   "test_periodic_3_triangulation_traits_SC_3.cpp")
 create_single_source_cgal_program(
   "test_periodic_3_triangulation_traits_SH_3.cpp")
+create_single_source_cgal_program("test_dummy_point_generation.cpp")

Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_dummy_point_generation.cpp

@@ -0,0 +1,148 @@
+#define CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY 1
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
+#include <CGAL/Periodic_3_regular_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_regular_triangulation_3.h>
+
+#include <CGAL/IO/OFF.h>
+#include <CGAL/Random.h>
+
+#include <array>
+#include <iostream>
+
+// - Use a lattice rather than a grid to have a Delaunay separation because constructions are inexact
+// - Use a reduced lattice as to use the 75 sufficient offsets giving the periodic star
+
+using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+using FT = K::FT;
+using Point_3 = K::Point_3;
+using Vector_3 = K::Vector_3;
+using Iso_cuboid = K::Iso_cuboid_3;
+
+using P3DTT3 = CGAL::Periodic_3_Delaunay_triangulation_traits_3<K>;
+using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3<P3DTT3>;
+
+using P3RTT3 = CGAL::Periodic_3_regular_triangulation_traits_3<K>;
+using P3RT3 = CGAL::Periodic_3_regular_triangulation_3<P3RTT3>;
+
+template <typename CellRangeIterator, typename Tr>
+void dump_cells(const std::string filename,
+                CellRangeIterator first,
+                CellRangeIterator beyond,
+                const Tr& tr,
+                const Vector_3& translation)
+{
+  using Cell_handle = typename Tr::Cell_handle;
+
+  auto cp = tr.geom_traits().construct_point_3_object();
+
+  std::vector<Point_3> points;
+  std::unordered_map<Point_3, std::size_t> points_to_ids;
+  std::vector<std::vector<std::size_t> > faces;
+
+  std::size_t pid = 0;
+  auto add_point = [&](const Point_3& p) -> std::size_t
+  {
+    auto insertion_result = points_to_ids.emplace(p, pid);
+    if(insertion_result.second)
+    {
+      points.push_back(p);
+      ++pid;
+    }
+
+    return insertion_result.first->second;
+  };
+
+  for(auto cell_it = first; cell_it != beyond; ++cell_it)
+  {
+    const Cell_handle ch = *cell_it;
+
+    const std::size_t v0_id = add_point(cp(tr.point(ch, 0)) + translation);
+    const std::size_t v1_id = add_point(cp(tr.point(ch, 1)) + translation);
+    const std::size_t v2_id = add_point(cp(tr.point(ch, 2)) + translation);
+    const std::size_t v3_id = add_point(cp(tr.point(ch, 3)) + translation);
+
+#if 1
+    faces.push_back(std::vector<std::size_t>{v0_id, v1_id, v2_id});
+    faces.push_back(std::vector<std::size_t>{v0_id, v2_id, v3_id});
+    faces.push_back(std::vector<std::size_t>{v0_id, v3_id, v1_id});
+    faces.push_back(std::vector<std::size_t>{v1_id, v2_id, v3_id});
+#else
+    faces.push_back(std::vector<std::size_t>{v0_id, v1_id, v2_id, v3_id});
+#endif
+  }
+
+  CGAL::IO::write_OFF(filename, points, faces, CGAL::parameters::stream_precision(17));
+}
+
+template <typename Tr>
+void dump(const std::string filename,
+          const Tr& tr)
+{
+  typedef typename Tr::Cell_handle Cell_handle;
+  std::vector<Cell_handle> chs;
+  for(auto cit = tr.finite_cells_begin(); cit != tr.finite_cells_end(); ++cit)
+    chs.push_back(cit);
+
+  dump_cells(filename, std::cbegin(chs), std::cend(chs), tr, CGAL::NULL_VECTOR);
+}
+
+template <typename Tr>
+void test(FT x_span, FT y_span, FT z_span)
+{
+  Tr p3t3;
+  p3t3.set_domain(Iso_cuboid(0,0,0, x_span,y_span,z_span));
+  p3t3.insert_generic_dummy_points();
+
+  // Abuse the multi cover function to check if cells have a small-enough orthoradius
+  p3t3.update_cover_data_after_converting_to_27_sheeted_covering();
+  if(!p3t3.can_be_converted_to_1_sheet())
+  {
+    std::cerr << "Error: dummy points do not create a 1-cover" << std::endl;
+    assert(false);
+  }
+
+//  dump("p3t3.off", p3t3);
+}
+
+int main (int argc, char** argv)
+{
+  CGAL::Random rnd = CGAL::get_default_random();
+  std::cout << "seed: " << rnd.get_seed() << std::endl;
+
+  if(argc != 1)
+  {
+    assert(argc == 4);
+    test<P3DT3>(std::stod(argv[1]), std::stod(argv[2]), std::stod(argv[3]));
+    std::cout << "OK" << std::endl;
+    return EXIT_SUCCESS;
+  }
+
+  for(int i=0; i<10; ++i)
+  {
+    test<P3DT3>(rnd.get_double(0.01, 1.),
+                rnd.get_double(0.01, 1.),
+                rnd.get_double(0.01, 1.));
+  }
+
+  for(int i=1; i<5; ++i)
+  {
+    for(int j=1; j<5; ++j)
+    {
+      for(int k=1; k<5; ++k)
+      {
+#if (CGAL_P3T3_DUMMY_GENERATION_DEBUG_VERBOSITY >= 1)
+        std::cout << "Test " << i << " " << j << " " << k << std::endl;
+#endif
+        test<P3DT3>(i,j,k);
+        test<P3RT3>(i,j,k);
+      }
+    }
+  }
+
+  std::cout << "OK" << std::endl;
+  return EXIT_SUCCESS;
+}

Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_periodic_3_delaunay_3.cpp

@@ -1,10 +1,9 @@
-#include <iostream>
-#include <fstream>
+#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
+#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
 
 #include <CGAL/Timer.h>
 
-#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
-#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
+#include <iostream>
 
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 typedef CGAL::Exact_predicates_inexact_constructions_kernel K1;
@@ -22,7 +21,8 @@ template class CGAL::Periodic_3_Delaunay_triangulation_3<PTT2>;
 #include <CGAL/Simple_homogeneous.h>
 typedef CGAL::Simple_homogeneous<CGAL::MP_Float> K3;
 typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K3> PTT3;
-// Explicit instantiation of the whole class :
+
+// Explicit instantiation of the whole class:
 template class CGAL::Periodic_3_Delaunay_triangulation_3<PTT3>;
 
 #include <CGAL/_test_cls_periodic_3_delaunay_3.h>

Periodic_3_triangulation_3/test/Periodic_3_triangulation_3/test_periodic_3_triangulation_3.cpp

@@ -1,13 +1,12 @@
-
-#include <iostream>
-
-#include <CGAL/Timer.h>
-
 #include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
 #include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
 #include <CGAL/Periodic_3_triangulation_3.h>
 #include <CGAL/Periodic_3_triangulation_traits_3.h>
 
+#include <CGAL/Timer.h>
+
+#include <iostream>
+
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 typedef CGAL::Exact_predicates_inexact_constructions_kernel     K1;
 typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K1>    PDTT1;
@@ -20,6 +19,7 @@ typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K2>    PDTT2;
 #include <CGAL/Simple_homogeneous.h>
 typedef CGAL::Simple_homogeneous<CGAL::MP_Float>                K3;
 typedef CGAL::Periodic_3_Delaunay_triangulation_traits_3<K3>    PDTT3;
+
 // Explicit instantiation of the whole class:
 template class CGAL::Periodic_3_triangulation_3<PDTT3>;
 

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Smoothing/mesh_smoothing_impl.h

@@ -468,7 +468,7 @@ public:
 //    std::cout << "y : " << initial_y << " -> " << y << "\n";
 //    std::cout << "z : " << initial_z << " -> " << z << "\n";
 
-    return Vector(FT(x - initial_x), FT(y - initial_y), FT(z - initial_z))l;
+    return Vector(FT(x - initial_x), FT(y - initial_y), FT(z - initial_z));
 #else
     CGAL_USE(v);
     return CGAL::NULL_VECTOR;

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap.h

@@ -152,11 +152,43 @@ void simplify_range(HalfedgeRange& halfedge_range,
             new_tolerance += CGAL::approximate_sqrt(CGAL::squared_distance(new_p, pt));
         }
 
-        if (!CGAL::Euler::does_satisfy_link_condition(edge(h, tm), tm))
+        // check that the collapse does not create a new degenerate face
+        bool do_collapse = true;
+        for(halfedge_descriptor he : halfedges_around_target(h, tm))
+        {
+          if(he != h &&
+             !is_border(he, tm) &&
+             collinear(get(vpm, source(he, tm)), new_p, get(vpm, target(next(he,tm),tm))))
+          {
+            do_collapse = false;
+            break;
+          }
+        }
+
+        if(!do_collapse)
           continue;
+
+        for(halfedge_descriptor he : halfedges_around_target(opposite(h,tm), tm))
+        {
+          if(he != opposite(h,tm) &&
+             !is_border(he, tm) &&
+             collinear(get(vpm, source(he, tm)), new_p, get(vpm, target(next(he,tm),tm))))
+          {
+            do_collapse = false;
+            break;
+          }
+        }
+
+        if(!do_collapse)
+          continue;
+
+        if(!CGAL::Euler::does_satisfy_link_condition(edge(h, tm), tm))
+          continue;
+
         const halfedge_descriptor opoh = opposite(prev(opposite(h, tm), tm), tm);
-        if (is_border(opoh, tm))
-          edges_to_test.erase( opoh );
+        if(is_border(opoh, tm))
+          edges_to_test.erase(opoh);
+
         vertex_descriptor v = Euler::collapse_edge(edge(h, tm), tm);
 
         put(vpm, v, new_p);
@@ -164,7 +196,7 @@ void simplify_range(HalfedgeRange& halfedge_range,
 
         if(get(range_halfedges, prev_h))
           edges_to_test.insert(prev_h);
-        if(next_h!=opoh && get(range_halfedges, next_h))
+        if(next_h != opoh && get(range_halfedges, next_h))
           edges_to_test.insert(next_h);
         ++collapsed_n;
       }
@@ -637,7 +669,6 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
 
   typedef typename boost::property_traits<VPMS>::value_type                       Point;
   typedef typename boost::property_traits<VPMT>::reference                        Point_ref;
-  typedef typename GeomTraits::Vector_3                                           Vector;
 
   typedef std::pair<halfedge_descriptor, Point>                                   Vertex_with_new_position;
   typedef std::vector<Vertex_with_new_position>                                   Vertices_with_new_position;
@@ -687,7 +718,12 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
 
       bool do_split = true;
 
-      // Some splits can create degenerate faces, avoid that
+      // In case of self-snapping, avoid degenerate caps
+      const bool is_same_mesh = (&tm_T == &tm_S);
+      if(is_same_mesh && target(next(opposite(h_to_split, tm_T), tm_T), tm_T) == splitter_v)
+        do_split = false;
+
+      // Do not split if it would create a degenerate needle
       if((new_position == get(vpm_T, target(h_to_split, tm_T))) ||
          (new_position == get(vpm_T, source(h_to_split, tm_T))))
         do_split = false;
@@ -695,11 +731,58 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
       if(!first_split && new_position == previous_split_position)
         do_split = false;
 
+      // check if the new faces after split will not be degenerate
+      const Point& p0 = new_position;
+      Point_ref p1 = get(vpm_T, source(h_to_split, tm_T));
+      Point_ref p2 = get(vpm_T, target(next(opposite(h_to_split, tm_T), tm_T), tm_T));
+      Point_ref p3 = get(vpm_T, target(h_to_split, tm_T));
+
+      /* Chooses the diagonal that will split the quad in two triangles that maximizes
+       * the scalar product of the un-normalized normals of the two triangles.
+       *
+       * The lengths of the un-normalized normals (computed using cross-products of two vectors)
+       * are proportional to the area of the triangles.
+       * Maximizing the scalar product of the two normals will avoid skinny triangles,
+       * and will also take into account the cosine of the angle between the two normals.
+       *
+       * In particular, if the two triangles are oriented in different directions,
+       * the scalar product will be negative.
+       */
+      auto p1p3 = CGAL::cross_product(p2-p1, p3-p2) * CGAL::cross_product(p0-p3, p1-p0);
+      auto p0p2 = CGAL::cross_product(p1-p0, p1-p2) * CGAL::cross_product(p3-p2, p3-p0);
+      bool first_split_face = (p0p2 > p1p3);
+
+      if(first_split_face)
+      {
+        if(p0p2 <= 0 || collinear(p0,p1,p2) || collinear(p0,p2,p3))
+          do_split = false;
+      }
+      else
+      {
+        if(p1p3 <= 0 || collinear(p0,p1,p3) || collinear(p1,p2,p3))
+          do_split = false;
+      }
+
+      if(do_split && !is_source_mesh_fixed)
+      {
+        for(halfedge_descriptor h : halfedges_around_target(splitter_v, tm_S))
+        {
+          if(!is_border(h,tm_S) && collinear(get(vpm_S, source(h,tm_S)), new_position, get(vpm_S, target(next(h,tm_S),tm_S))))
+          {
+            do_split = false;
+            break;
+          }
+        }
+
+        if(do_split)
+          put(vpm_S, splitter_v, new_position);
+      }
+
 #ifdef CGAL_PMP_SNAP_DEBUG_PP
       std::cout << " -.-.-. Splitting " << edge(h_to_split, tm_T) << " |||| "
                 << " Vs " << source(h_to_split, tm_T) << " (" << tm_T.point(source(h_to_split, tm_T)) << ")"
                 << " --- Vt " << target(h_to_split, tm_T) << " (" << tm_T.point(target(h_to_split, tm_T)) << ")" << std::endl;
-      std::cout << "With point: " << vnp.second << std::endl;
+      std::cout << "With point: " << new_position << " (init: " << vnp.second << ")" << std::endl;
       std::cout << "Actually split? " << do_split << std::endl;
 #endif
 
@@ -715,75 +798,27 @@ std::size_t split_edges(EdgesToSplitContainer& edges_to_split,
 
         visitor.after_vertex_edge_snap(new_v, tm_T);
       }
-
-      if(!is_source_mesh_fixed)
-        put(vpm_S, splitter_v, new_position);
+      else
+      {
+        continue;
+      }
 
       first_split = false;
       previous_split_position = new_position;
       ++snapped_n;
 
-      // Everything below is choosing the diagonal to triangulate the quad formed by the edge split
-      // So, it's only relevant if splitting has been performed
-      if(!do_split)
-        continue;
-
-      /*          new_p
-       *         /   \
-       *    res /     \ h_to_split
-       *       /       \
-       *      /         \
-       *    left       right
-       *     |         /
-       *     |        /
-       *     |       /
-       *     |      /
-       *     |     /
-       *     |    /
-       *      opp
-       */
-
-      const halfedge_descriptor res = prev(h_to_split, tm_T);
-      const Point_ref left_pt = get(vpm_T, source(res, tm_T));
-      const Point_ref right_pt = get(vpm_T, target(h_to_split, tm_T));
-      const Point_ref opp = get(vpm_T, target(next(opposite(res, tm_T), tm_T), tm_T));
-
-      // Check if 'p' is "visible" from 'opp' (i.e. its projection on the plane 'Pl(left, opp, right)'
-      // falls in the cone with apex 'opp' and sides given by 'left' and 'right')
-      const Vector n = gt.construct_orthogonal_vector_3_object()(right_pt, left_pt, opp);
-      const Point trans_left_pt = gt.construct_translated_point_3_object()(left_pt, n);
-      const Point trans_right_pt = gt.construct_translated_point_3_object()(right_pt, n);
-
-      const Point_ref new_p = get(vpm_T, new_v);
-      const bool left_of_left = (gt.orientation_3_object()(trans_left_pt, left_pt, opp, new_p) == CGAL::POSITIVE);
-      const bool right_of_right = (gt.orientation_3_object()(right_pt, trans_right_pt, opp, new_p) == CGAL::POSITIVE);
-
-      const bool is_visible = (!left_of_left && !right_of_right);
-
-#ifdef CGAL_PMP_SNAP_DEBUG_PP
-      std::cout << "Left/Right: " << left_of_left << " " << right_of_right << std::endl;
-      std::cout << "visible from " << opp << " ? " << is_visible << std::endl;
-#endif
-
-      // h_to_split is equal to 'next(res)' after splitting
-      const halfedge_descriptor h_to_split_opp = opposite(h_to_split, tm_T);
-
-      if(is_visible)
-      {
-        halfedge_descriptor h2 = prev(prev(h_to_split_opp, tm_T), tm_T);
-        halfedge_descriptor new_hd = CGAL::Euler::split_face(h_to_split_opp,
-                                                             h2, tm_T);
-        h_to_split = opposite(prev(new_hd, tm_T), tm_T);
-        visitor.after_split_face(h_to_split_opp, h2, tm_T);
-      }
+      halfedge_descriptor v0, v1, v2, v3;
+      v0 = opposite(h_to_split, tm_T);
+      v1 = next(v0, tm_T);
+      v2 = next(v1, tm_T);
+      v3 = next(v2, tm_T);
+      // halfedge_descriptor new_hd =
+      first_split_face ? CGAL::Euler::split_face(v0, v2, tm_T)
+                       : CGAL::Euler::split_face(v1, v3, tm_T);
+      if(first_split_face)
+        visitor.after_split_face(v0, v2, tm_T);
       else
-      {
-        halfedge_descriptor h2 = prev(h_to_split_opp, tm_T);
-        halfedge_descriptor new_hd = CGAL::Euler::split_face(opposite(res, tm_T),
-                                                             h2, tm_T);
-        h_to_split = opposite(next(new_hd, tm_T), tm_T);
-        visitor.after_split_face(opposite(res, tm_T), h2, tm_T);
-      }
+        visitor.after_split_face(v1, v3, tm_T);
     }
   }
 
@@ -1406,6 +1441,8 @@ std::size_t snap_borders(TriangleMesh& tm,
                                                       true /*self snapping*/, np, np);
 }
 
+//TODO:add an option to preserve orientation?
+
 } // end namespace experimental
 } // end namespace Polygon_mesh_processing
 } // end namespace CGAL

Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/internal/Snapping/snap_vertices.h

@@ -706,6 +706,7 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
   typedef typename GetGeomTraits<PolygonMesh, NamedParameters_A>::type         GT;
   typedef typename GT::FT                                                      FT;
   typedef typename boost::property_traits<VPM_B>::value_type                   Point;
+  typedef typename boost::property_traits<VPM_B>::reference                    Point_ref;
 
   typedef std::vector<halfedge_descriptor>                                     Vertex_container;
   typedef std::pair<Vertex_container, FT>                                      Unique_vertex;
@@ -955,8 +956,23 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
     {
       if(is_second_mesh_fixed)
       {
+        const Point_ref new_p = get(vpm_B, vb);
         for(const halfedge_descriptor ha : vs_a)
-          put(vpm_A, target(ha, tm_A), get(vpm_B, vb));
+        {
+          bool skip = false;
+          for(halfedge_descriptor haa : halfedges_around_target(ha, tm_A))
+          {
+            if(!is_border(haa,tm_A) &&
+               collinear(get(vpm_A, source(haa,tm_A)), new_p, get(vpm_A, target(next(haa,tm_A),tm_A))))
+            {
+              skip = true;
+              break;
+            }
+          }
+
+          if(!skip)
+            put(vpm_A, target(ha, tm_A), new_p);
+        }
       }
       else
       {
@@ -972,12 +988,40 @@ std::size_t snap_vertices_two_way(const HalfedgeRange_A& halfedge_range_A,
 #endif
 
         for(const halfedge_descriptor ha : vs_a)
-          put(vpm_A, target(ha, tm_A), new_p);
+        {
+          bool skip = false;
+          for(halfedge_descriptor haa : halfedges_around_target(ha, tm_A))
+          {
+            if(!is_border(haa,tm_A) &&
+               collinear(get(vpm_A, source(haa,tm_A)), new_p, get(vpm_A, target(next(haa,tm_A),tm_A))))
+            {
+              skip = true;
+              break;
+            }
+          }
+
+          if(!skip)
+            put(vpm_A, target(ha, tm_A), new_p);
+        }
 
         for(const halfedge_descriptor hb : vs_b)
-          put(vpm_B, target(hb, tm_B), new_p);
-      }
+        {
+          bool skip = false;
+          for(halfedge_descriptor hbb : halfedges_around_target(hb, tm_B))
+          {
+            if(!is_border(hbb,tm_B) &&
+               collinear(get(vpm_B, source(hbb,tm_B)), new_p, get(vpm_B, target(next(hbb,tm_B),tm_B))))
+            {
+              skip = true;
+              break;
+            }
+          }
 
+          if(!skip)
+            put(vpm_B, target(hb, tm_B), new_p);
+        }
+      }
+      //TODO: the counter shall depend on skip?
       ++counter;
     }
 

Polygonal_surface_reconstruction/examples/Polygonal_surface_reconstruction/polyfit_example_with_region_growing.cpp

@@ -3,7 +3,7 @@
 #include <CGAL/property_map.h>
 #include <CGAL/Surface_mesh.h>
 #include <CGAL/Shape_detection/Region_growing/Region_growing.h>
-#include <CGAL/Shape_detection/Region_growing/Region_growing_on_point_set.h>
+#include <CGAL/Shape_detection/Region_growing/Point_set.h>
 #include <CGAL/Polygonal_surface_reconstruction.h>
 
 #ifdef CGAL_USE_SCIP  // defined (or not) by CMake scripts, do not define by hand
@@ -22,6 +22,7 @@ typedef CGAL::GLPK_mixed_integer_program_traits<double>        MIP_Solver;
 
 #include <fstream>
 #include <CGAL/Timer.h>
+#include <boost/range/irange.hpp>
 
 typedef CGAL::Exact_predicates_inexact_constructions_kernel        Kernel;
 
@@ -37,45 +38,15 @@ typedef CGAL::Nth_of_tuple_property_map<0, PNI>        Point_map;
 typedef CGAL::Nth_of_tuple_property_map<1, PNI>        Normal_map;
 typedef CGAL::Nth_of_tuple_property_map<2, PNI>        Plane_index_map;
 
-typedef CGAL::Shape_detection::Point_set::
-Sphere_neighbor_query<Kernel, Point_vector, Point_map> Neighbor_query;
-typedef CGAL::Shape_detection::Point_set::
-Least_squares_plane_fit_region<Kernel, Point_vector, Point_map, Normal_map> Region_type;
-typedef CGAL::Shape_detection::
-Region_growing<Neighbor_query, Region_type> Region_growing;
+using Point_map_region_growing = CGAL::Compose_property_map<CGAL::Random_access_property_map<Point_vector>, Point_map >;
+using Normal_map_region_growing = CGAL::Compose_property_map<CGAL::Random_access_property_map<Point_vector>, Normal_map >;
+
+using Region_type = CGAL::Shape_detection::Point_set::Least_squares_plane_fit_region<Kernel, std::size_t, Point_map_region_growing, Normal_map_region_growing>;
+using Neighbor_query = CGAL::Shape_detection::Point_set::Sphere_neighbor_query<Kernel, std::size_t, Point_map_region_growing>;
+using Region_growing = CGAL::Shape_detection::Region_growing<Neighbor_query, Region_type>;
 
 typedef CGAL::Surface_mesh<Point>        Surface_mesh;
-typedef        CGAL::Polygonal_surface_reconstruction<Kernel> Polygonal_surface_reconstruction;
-
-class Index_map {
-
-public:
-  using key_type = std::size_t;
-  using value_type = int;
-  using reference = value_type;
-  using category = boost::readable_property_map_tag;
-
-  Index_map() { }
-  template<typename PointRange>
-  Index_map(const PointRange& points,
-            const std::vector< std::vector<std::size_t> >& regions)
-    : m_indices(new std::vector<int>(points.size(), -1))
-  {
-    for (std::size_t i = 0; i < regions.size(); ++i)
-      for (const std::size_t idx : regions[i])
-        (*m_indices)[idx] = static_cast<int>(i);
-  }
-
-  inline friend value_type get(const Index_map& index_map,
-                               const key_type key)
-  {
-    const auto& indices = *(index_map.m_indices);
-    return indices[key];
-  }
-
-private:
-  std::shared_ptr< std::vector<int> > m_indices;
-};
+typedef CGAL::Polygonal_surface_reconstruction<Kernel> Polygonal_surface_reconstruction;
 
 /*
 * This example first extracts planes from the input point cloud
@@ -100,7 +71,7 @@ int main(int argc, char* argv[])
   CGAL::Timer t;
   t.start();
   if (!CGAL::IO::read_points(input_file.c_str(), std::back_inserter(points),
-                             CGAL::parameters::point_map(Point_map()).normal_map(Normal_map()))) {
+    CGAL::parameters::point_map(Point_map()).normal_map(Normal_map()))) {
 
     std::cerr << "Error: cannot read file " << input_file << std::endl;
     return EXIT_FAILURE;
@@ -114,38 +85,41 @@ int main(int argc, char* argv[])
   // Shape detection.
 
   // Default parameter values for the data file cube.pwn.
-  const FT          search_sphere_radius  = FT(2) / FT(100);
+  const FT          search_sphere_radius = FT(2) / FT(100);
   const FT          max_distance_to_plane = FT(2) / FT(1000);
-  const FT          max_accepted_angle    = FT(25);
-  const std::size_t min_region_size       = 200;
+  const FT          max_accepted_angle = FT(25);
+  const std::size_t min_region_size = 200;
+
+  Point_map_region_growing point_map_rg(CGAL::make_random_access_property_map(points));
+  Normal_map_region_growing normal_map_rg(CGAL::make_random_access_property_map(points));
 
   // Create instances of the classes Neighbor_query and Region_type.
   Neighbor_query neighbor_query(
-    points,
-    search_sphere_radius);
+    boost::irange<std::size_t>(0, points.size()), CGAL::parameters::sphere_radius(search_sphere_radius).point_map(point_map_rg));
 
   Region_type region_type(
-    points,
-    max_distance_to_plane, max_accepted_angle, min_region_size);
+    CGAL::parameters::
+    maximum_distance(max_distance_to_plane).
+    maximum_angle(max_accepted_angle).
+    minimum_region_size(min_region_size).
+    point_map(point_map_rg).
+    normal_map(normal_map_rg));
 
   // Create an instance of the region growing class.
   Region_growing region_growing(
-    points, neighbor_query, region_type);
+    boost::irange<std::size_t>(0, points.size()), neighbor_query, region_type);
 
   std::cout << "Extracting planes...";
-  std::vector< std::vector<std::size_t> > regions;
+  std::vector<typename Region_growing::Primitive_and_region> regions;
   t.reset();
   region_growing.detect(std::back_inserter(regions));
   std::cout << " Done. " << regions.size() << " planes extracted. Time: "
-  << t.time() << " sec." << std::endl;
+    << t.time() << " sec." << std::endl;
 
   // Stores the plane index of each point as the third element of the tuple.
-  Index_map index_map(points, regions);
-  for (std::size_t i = 0; i < points.size(); ++i) {
+  for (std::size_t i = 0; i < points.size(); ++i)
     // Uses the get function from the property map that accesses the 3rd element of the tuple.
-    const int plane_index = get(index_map, i);
-    points[i].get<2>() = plane_index;
-  }
+    points[i].get<2>() = static_cast<int>(get(region_growing.region_map(), i));
 
   //////////////////////////////////////////////////////////////////////////
 
@@ -187,8 +161,8 @@ int main(int argc, char* argv[])
 
 int main(int, char**)
 {
-    std::cerr << "This test requires either GLPK or SCIP.\n";
-    return EXIT_SUCCESS;
+  std::cerr << "This test requires either GLPK or SCIP.\n";
+  return EXIT_SUCCESS;
 }
 
 #endif  // defined(CGAL_USE_GLPK) || defined(CGAL_USE_SCIP)

SMDS_3/examples/SMDS_3/c3t3_example.cpp

@@ -33,13 +33,20 @@ using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3<Triangulation>;
 
 int main(int argc, char* argv[])
 {
+  std::cout.precision(17);
+  std::cerr.precision(17);
+
   std::string filename = (argc > 1) ? std::string(argv[1])
                        : CGAL::data_file_path("meshes/elephant.mesh");
 
   Triangulation tr;
 
   std::ifstream is(filename, std::ios_base::in);
-  CGAL::IO::read_MEDIT(is, tr);
+  if(!CGAL::IO::read_MEDIT(is, tr))
+  {
+    std::cerr << "Failed to read" << std::endl;
+    return EXIT_FAILURE;
+  }
 
   // [call a remeshing algorithm]
 
@@ -49,8 +56,12 @@ int main(int argc, char* argv[])
 
   Triangulation tr2;
   std::ifstream is2("after_remeshing.mesh");
-  CGAL::IO::read_MEDIT(is2, tr2);
-  is2.close();
+  if(!CGAL::IO::read_MEDIT(is2, tr2))
+  {
+    std::cerr << "Failed to read (#2)" << std::endl;
+    return EXIT_FAILURE;
+  }
 
+  std::cout << "Done" << std::endl;
   return EXIT_SUCCESS;
 }

SMDS_3/include/CGAL/IO/File_medit.h

@@ -22,12 +22,11 @@
 
 #include <CGAL/utility.h>
 #include <CGAL/basic.h>
-
-#include <boost/unordered_map.hpp>
-
 #include <CGAL/Named_function_parameters.h>
 #include <CGAL/boost/graph/named_params_helper.h>
 
+#include <boost/unordered_map.hpp>
+
 #include <iomanip>
 #include <iostream>
 #include <map>
@@ -1090,9 +1089,26 @@ void write_MEDIT(std::ostream& os,
  * @tparam T3 can be instantiated with any 3D triangulation of \cgal provided that its
  *  vertex and cell base class are models of the concepts `MeshVertexBase_3` and `MeshCellBase_3`,
  *  respectively.
+ * @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters"
  *
  * @param in the input stream
  * @param t3 the triangulation
+ * @param np optional \ref bgl_namedparameters "Named Parameters" described below
+ *
+ * \cgalNamedParamsBegin
+ *   \cgalParamNBegin{verbose}
+ *     \cgalParamDescription{indicates whether output warnings and error messages should be printed or not.}
+ *     \cgalParamType{Boolean}
+ *     \cgalParamDefault{`false`}
+ *   \cgalParamNEnd
+ *   \cgalParamNBegin{allow_non_manifold}
+ *     \cgalParamDescription{allows the construction of a triangulation with non-manifold edges
+ *       and non manifold vertices. The triangulation is invalid if this situation is met,
+ *       so it should be used only in advanced cases, and the triangulation will be hardly usable.}
+ *     \cgalParamType{bool}
+ *     \cgalParamDefault{false}
+ *   \cgalParamNEnd
+ * \cgalNamedParamsEnd
  *
  * @returns `true` if the connectivity of the triangulation could be built consistently
  * from \p in,
@@ -1108,11 +1124,20 @@ void write_MEDIT(std::ostream& os,
  * positively oriented cells,
  * and cover the geometric convex hull of all points in `t3`.
  */
-template<typename T3>
-bool read_MEDIT(std::istream& in, T3& t3)
+template<typename T3, typename CGAL_NP_TEMPLATE_PARAMETERS>
+bool read_MEDIT(std::istream& in,
+                T3& t3,
+                const CGAL_NP_CLASS& np = parameters::default_values())
 {
-  CGAL_assertion(!(!in));
-  bool b = CGAL::SMDS_3::build_triangulation_from_file(in, t3);
+  using parameters::choose_parameter;
+  using parameters::get_parameter;
+
+  // Default non_manifold value is true if the triangulation periodic, false otherwise
+  const bool non_manifold = choose_parameter(get_parameter(np, internal_np::allow_non_manifold),
+                                             std::is_same<typename T3::Periodic_tag, Tag_true>::value);
+  const bool verbose = choose_parameter(get_parameter(np, internal_np::verbose), false);
+
+  bool b = CGAL::SMDS_3::build_triangulation_from_file(in, t3, verbose, false /*replace_domain_0*/, non_manifold);
   if(!b)
     t3.clear();
   return b;

SMDS_3/include/CGAL/SMDS_3/tet_soup_to_c3t3.h

@@ -19,21 +19,19 @@
 
 #include <CGAL/license/SMDS_3.h>
 
-#include <CGAL/disable_warnings.h>
 #include <CGAL/assertions.h>
 #include <CGAL/IO/File_medit.h>
 
-#include <array>
-#include <vector>
-#include <utility>
-#include <map>
 #include <boost/unordered_map.hpp>
 
+#include <array>
+#include <map>
+#include <utility>
+#include <vector>
+
+namespace CGAL {
+namespace SMDS_3 {
 
-namespace CGAL
-{
-namespace SMDS_3
-{
 template<typename Vh>
 std::array<Vh, 3> make_ordered_vertex_array(const Vh vh0, const Vh vh1, const Vh vh2)
 {
@@ -64,12 +62,12 @@ void build_vertices(Tr& tr,
   }
 }
 
-template<class Tr>
+template <class Tr>
 bool add_facet_to_incident_cells_map(const typename Tr::Cell_handle c, int i,
-    boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
-                         std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
-    const bool verbose,
-    const bool allow_non_manifold)
+                                     boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
+                                                          std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
+                                     const bool verbose,
+                                     const bool allow_non_manifold)
 {
   typedef typename Tr::Vertex_handle                                Vertex_handle;
   typedef typename Tr::Cell_handle                                  Cell_handle;
@@ -81,40 +79,49 @@ bool add_facet_to_incident_cells_map(const typename Tr::Cell_handle c, int i,
 
   // the opposite vertex of f in c is i
   Facet_vvv f = CGAL::SMDS_3::make_ordered_vertex_array(c->vertex((i + 1) % 4),
-                                          c->vertex((i + 2) % 4),
-                                          c->vertex((i + 3) % 4));
+                                                        c->vertex((i + 2) % 4),
+                                                        c->vertex((i + 3) % 4));
   CGAL_precondition(f[0] != f[1] && f[1] != f[2]);
 
   Incident_cell e = std::make_pair(c, i);
   std::vector<Incident_cell> vec;
   vec.push_back(e);
   std::pair<typename Incident_cells_map::iterator, bool> is_insert_successful =
-      incident_cells_map.insert(std::make_pair(f, vec));
+      incident_cells_map.emplace(f, vec);
   if(!is_insert_successful.second) // the entry already exists in the map
   {
-    // a facet must have exactly two incident cells
-    if (is_insert_successful.first->second.size() != 1)
+    // A finite facet must have exactly two incident cells
+    //
+    // If there is a non-manifold edge on the boundary, the infinite facet being
+    // that edge + the infinite vertex has (strictly) more than 2 incident cells
+    if(is_insert_successful.first->second.size() != 1)
     {
-      if(verbose)
-        std::cout << "Error in add_facet_to_incident_cells_map" << std::endl;
       if(!allow_non_manifold)
+      {
         success = false;
+        if(verbose)
+          std::cerr << "Error: " << is_insert_successful.first->second.size() << " previous incidences" << std::endl;
+      }
+      else if(verbose)
+      {
+        std::cerr << "Warning: " << is_insert_successful.first->second.size() << " previous incidences" << std::endl;
+      }
     }
     is_insert_successful.first->second.push_back(e);
   }
   return success;
 }
 
-template<class Tr, typename CellRange, typename SubdomainsRange, typename FacetPatchMap>
+template <class Tr, typename CellRange, typename SubdomainsRange, typename FacetPatchMap>
 bool build_finite_cells(Tr& tr,
-    const CellRange& finite_cells,
-    const SubdomainsRange& subdomains,
-    const std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
-    boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
-                        std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
-    const FacetPatchMap& border_facets,
-    const bool verbose,
-    const bool replace_domain_0)
+                        const CellRange& finite_cells,
+                        const SubdomainsRange& subdomains,
+                        const std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
+                        boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
+                                             std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
+                        const FacetPatchMap& border_facets,
+                        const bool verbose,
+                        const bool replace_domain_0)
 {
   typedef typename Tr::Vertex_handle                            Vertex_handle;
   typedef typename Tr::Cell_handle                              Cell_handle;
@@ -123,11 +130,10 @@ bool build_finite_cells(Tr& tr,
   bool success = true;
 
   CGAL_assertion_code(
-    typename Tr::Geom_traits::Construct_point_3 cp =
-      tr.geom_traits().construct_point_3_object();
-    typename Tr::Geom_traits::Orientation_3 orientation =
-      tr.geom_traits().orientation_3_object();
+    typename Tr::Geom_traits::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();
+    typename Tr::Geom_traits::Orientation_3 orientation = tr.geom_traits().orientation_3_object();
   )
+
   typename SubdomainsRange::value_type max_domain = 0;
   if(replace_domain_0)
   {
@@ -137,6 +143,7 @@ bool build_finite_cells(Tr& tr,
         max_domain = subdomains[i];
     }
   }
+
   // build the finite cells
   for(std::size_t i=0; i<finite_cells.size(); ++i)
   {
@@ -154,15 +161,12 @@ bool build_finite_cells(Tr& tr,
 
     // this assertion also tests for degeneracy
     CGAL_assertion(orientation(cp(tr.point(vs[0])), cp(tr.point(vs[1])),
-                               cp(tr.point(vs[2])), cp(tr.point(vs[3])))
-                     == POSITIVE);
+                               cp(tr.point(vs[2])), cp(tr.point(vs[3]))) == POSITIVE);
 
     Cell_handle c = tr.tds().create_cell(vs[0], vs[1], vs[2], vs[3]);
     c->set_subdomain_index(subdomains[i]); // the cell's info keeps the reference of the tetrahedron
     if(replace_domain_0 && subdomains[i] == 0)
-    {
       c->set_subdomain_index(max_domain+1); // the cell's info keeps the reference of the tetrahedron
-    }
 
     // assign cells to vertices
     for(int j=0; j<4; ++j)
@@ -174,39 +178,35 @@ bool build_finite_cells(Tr& tr,
     // build the map used for adjacency later
     for(int j=0; j<4; ++j)
     {
+      // do not allow non-manifoldness in the finite cells case
       if(!CGAL::SMDS_3::add_facet_to_incident_cells_map<Tr>(c, j, incident_cells_map, verbose, false))
-          //do not allow non-manifold in the finite cells case
         success = false;
+
       if(border_facets.size() != 0)
       {
         std::array<int,3> facet;
-        facet[0]=tet[(j+1) % 4];
-        facet[1]=tet[(j+2) % 4];
-        facet[2]=tet[(j+3) % 4];
-        //find the circular permutation that puts the smallest index in the first place.
-        int n0 = (std::min)((std::min)(facet[0], facet[1]), facet[2]);
-        int k=0;
-        std::array<int,3> f;
+        facet[0] = tet[(j+1) % 4];
+        facet[1] = tet[(j+2) % 4];
+        facet[2] = tet[(j+3) % 4];
+
+        // find the circular permutation that puts the smallest index in the first place.
+        int n0 = (std::min)({facet[0], facet[1], facet[2]});
         do
         {
-          f[0]=facet[(0+k)%3];
-          f[1]=facet[(1+k)%3];
-          f[2]=facet[(2+k)%3];
-          ++k;
-        } while(f[0] != n0);
+          std::rotate(std::begin(facet), std::next(std::begin(facet)), std::end(facet));
+        }
+        while(facet[0] != n0);
 
-        typename FacetPatchMap::const_iterator it = border_facets.find(f);
+        typename FacetPatchMap::const_iterator it = border_facets.find(facet);
         if(it != border_facets.end())
         {
           c->set_surface_patch_index(j, it->second);
         }
         else
         {
-          int temp = f[2];
-          f[2] = f[1];
-          f[1] = temp;
+          std::swap(facet[1], facet[2]); // facet[0] is still the smallest, no need to rotate again
 
-          it = border_facets.find(f);
+          it = border_facets.find(facet);
           if(it != border_facets.end())
             c->set_surface_patch_index(j, it->second);
           else
@@ -215,16 +215,17 @@ bool build_finite_cells(Tr& tr,
       }
     }
   }
+
   return success;
 }
 
 template<class Tr>
 bool add_infinite_facets_to_incident_cells_map(typename Tr::Cell_handle c,
-     int inf_vert_pos,
-     boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
-                          std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
-     const bool verbose,
-     const bool allow_non_manifold)
+                                               int inf_vert_pos,
+                                               boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
+                                                                    std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
+                                               const bool verbose,
+                                               const bool allow_non_manifold)
 {
   int l = (inf_vert_pos + 1) % 4;
   bool b1 = CGAL::SMDS_3::add_facet_to_incident_cells_map<Tr>(c, l, incident_cells_map, verbose, allow_non_manifold);
@@ -232,15 +233,16 @@ bool add_infinite_facets_to_incident_cells_map(typename Tr::Cell_handle c,
   bool b2 = CGAL::SMDS_3::add_facet_to_incident_cells_map<Tr>(c, l, incident_cells_map, verbose, allow_non_manifold);
   l = (inf_vert_pos + 3) % 4;
   bool b3 = CGAL::SMDS_3::add_facet_to_incident_cells_map<Tr>(c, l, incident_cells_map, verbose, allow_non_manifold);
+
   return b1 && b2 && b3;
 }
 
 template<class Tr>
 bool build_infinite_cells(Tr& tr,
-  boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
-                       std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
-  const bool verbose,
-  const bool allow_non_manifold)
+                          boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
+                                               std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
+                          const bool verbose,
+                          const bool allow_non_manifold)
 {
   typedef typename Tr::Vertex_handle                               Vertex_handle;
   typedef typename Tr::Cell_handle                                 Cell_handle;
@@ -260,13 +262,14 @@ bool build_infinite_cells(Tr& tr,
   {
     if(it->second.size() == 2) // facet already has both its incident cells
       continue;
+
     CGAL_assertion(it->second.size() == 1);
 
     Cell_handle c = it->second[0].first;
     int i = it->second[0].second;
 
-    Cell_handle opp_c;
     // the infinite cell that we are creating needs to be well oriented...
+    Cell_handle opp_c;
     if(i == 0 || i == 2)
       opp_c = tr.tds().create_cell(tr.infinite_vertex(),
                                    c->vertex((i + 2) % 4),
@@ -285,7 +288,7 @@ bool build_infinite_cells(Tr& tr,
       tr.infinite_vertex()->set_cell(opp_c);
 
     // the only finite facet
-    it->second.push_back(std::make_pair(opp_c, 0));
+    it->second.emplace_back(opp_c, 0);
     CGAL_assertion(it->second.size() == 2);
 
     opp_c->set_surface_patch_index(0, c->surface_patch_index(i));
@@ -301,23 +304,24 @@ bool build_infinite_cells(Tr& tr,
     for (int i = 1; i < 4; ++i)
     {
       std::array<Vertex_handle, 3> vs = CGAL::SMDS_3::make_ordered_vertex_array(c->vertex((i + 1) % 4),
-        c->vertex((i + 2) % 4),
-        c->vertex((i + 3) % 4));
+                                                                                c->vertex((i + 2) % 4),
+                                                                                c->vertex((i + 3) % 4));
       if (facets.find(vs) == facets.end())
-        facets.insert(std::make_pair(vs, 1));
+        facets.emplace(vs, 1);
       else
         facets[vs]++;
     }
   }
+
   for (auto fp : facets)
   {
     if (fp.second != 2)
     {
-      std::cout << "Warning : non manifold edge" << std::endl;
+      std::cout << "Warning: non manifold edge" << std::endl;
       std::cout << "fp.second = " << fp.second << std::endl;
-      std::cout << fp.first[0]->point() << " "
-        << fp.first[1]->point() << " "
-        << fp.first[2]->point() << std::endl;
+      std::cout << tr.point(fp.first[0]) << " "
+                << tr.point(fp.first[1]) << " "
+                << tr.point(fp.first[2]) << std::endl;
       success = false;
     }
 //    CGAL_assertion(fp.second == 2);
@@ -326,21 +330,23 @@ bool build_infinite_cells(Tr& tr,
 
   // add the facets to the incident cells map
   for (const Cell_handle& c : infinite_cells)
+  {
     if(!CGAL::SMDS_3::add_infinite_facets_to_incident_cells_map<Tr>(c,
-                                                                    0,
+                                                                    c->index(tr.infinite_vertex()),
                                                                     incident_cells_map,
                                                                     verbose,
                                                                     allow_non_manifold))
       success = false;
+  }
 
   return success;
 }
 
 template<typename Tr>
-bool has_infinite_vertex(const std::array<typename Tr::Vertex_handle, 3>& v,
-                         const Tr& tr)
+bool is_infinite(const std::array<typename Tr::Vertex_handle, 3>& f,
+                 const Tr& tr)
 {
-  for (auto vh : v)
+  for (auto vh : f)
   {
     if (tr.infinite_vertex() == vh)
       return true;
@@ -350,9 +356,9 @@ bool has_infinite_vertex(const std::array<typename Tr::Vertex_handle, 3>& v,
 
 template<class Tr>
 bool assign_neighbors(Tr& tr,
-  const boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
-                             std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
-  const bool allow_non_manifold)
+                      const boost::unordered_map<std::array<typename Tr::Vertex_handle, 3>,
+                                                 std::vector<std::pair<typename Tr::Cell_handle, int> > >& incident_cells_map,
+                      const bool allow_non_manifold)
 {
   typedef typename Tr::Cell_handle                                   Cell_handle;
   typedef std::pair<Cell_handle, int>                                Incident_cell;
@@ -374,10 +380,12 @@ bool assign_neighbors(Tr& tr,
 
       tr.tds().set_adjacency(c0, i0, c1, i1);
     }
-    else if(allow_non_manifold)// if (adjacent_cells.size() == 4)
+    else if(!allow_non_manifold)
     {
-      CGAL_assertion_code(const auto& v = icit->first);
-      CGAL_assertion(has_infinite_vertex(v, tr));
+      CGAL_assertion_code(const auto& f = icit->first);
+      CGAL_assertion(is_infinite(f, tr));
+
+      CGAL_assertion(adjacent_cells.size() % 2 == 0);
       success = false;
     }
   }
@@ -389,17 +397,18 @@ template<class Tr,
          typename CellRange,
          typename FacetPatchMap>
 bool build_triangulation_impl(Tr& tr,
-    const PointRange& points,
-    const CellRange& finite_cells,
-    const std::vector<typename Tr::Cell::Subdomain_index>& subdomains,
-    const FacetPatchMap& border_facets,
-    std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
-    const bool verbose,// = false,
-    const bool replace_domain_0,// = false,
-    const bool allow_non_manifold) // = false
+                              const PointRange& points,
+                              const CellRange& finite_cells,
+                              const std::vector<typename Tr::Cell::Subdomain_index>& subdomains,
+                              const FacetPatchMap& border_facets,
+                              std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
+                              const bool verbose,// = false,
+                              const bool replace_domain_0,// = false,
+                              const bool allow_non_manifold) // = false
 {
   if (verbose)
     std::cout << "build_triangulation_impl()..." << std::endl;
+
   typedef typename Tr::Vertex_handle            Vertex_handle;
   typedef typename Tr::Cell_handle              Cell_handle;
   typedef std::array<Vertex_handle, 3>          Facet_vvv;
@@ -408,64 +417,79 @@ bool build_triangulation_impl(Tr& tr,
   typedef std::pair<Cell_handle, int>                   Incident_cell;
   typedef boost::unordered_map<Facet_vvv, std::vector<Incident_cell> >  Incident_cells_map;
 
-  bool success = true;
-
-  Incident_cells_map incident_cells_map;
-  vertex_handle_vector.resize(points.size() + 1); // id to vertex_handle
-                                        //index 0 is for infinite vertex
-                                        // 1 to n for points in `points`
-
   CGAL_precondition(!points.empty());
 
+  bool success = true;
+  Incident_cells_map incident_cells_map;
+
+  // id to vertex_handle
+  // index 0 is for infinite vertex; 1 to n for points in `points`
+  vertex_handle_vector.resize(points.size() + 1);
+
   if(finite_cells.empty())
   {
     if (verbose)
-      std::cout << "WARNING: No finite cells were provided. Only the points will be loaded."<<std::endl;
+      std::cout << "WARNING: No finite cells were provided. Only the points will be loaded." << std::endl;
   }
 
-  tr.tds().clear(); // not tr.clear() since it calls tr.init() which we don't want
+  tr.tds().clear(); // not tr.clear() since it calls tr.init(), which we don't want
 
   build_vertices<Tr>(tr, points, vertex_handle_vector);
   for(Vertex_handle vh : vertex_handle_vector)
-  {
     vh->set_dimension(-1);
-  }
+
+  if(verbose)
+    std::cout << "build vertices done (" << tr.tds().number_of_vertices() << " vertices)" << std::endl;
+
   if (!finite_cells.empty())
   {
-    if (!CGAL::SMDS_3::build_finite_cells<Tr>(tr, finite_cells, subdomains, vertex_handle_vector, incident_cells_map,
-      border_facets, verbose, replace_domain_0))
+    if (!CGAL::SMDS_3::build_finite_cells<Tr>(tr, finite_cells, subdomains, vertex_handle_vector,
+                                              incident_cells_map, border_facets, verbose, replace_domain_0))
     {
-      if (verbose) std::cout << "build_finite_cells went wrong" << std::endl;
+      if (verbose)
+        std::cerr << "Error: build_finite_cells went wrong!" << std::endl;
       success = false;
     }
-    else
-      if (verbose) std::cout << "build finite cells done" << std::endl;
+    else if(verbose)
+    {
+      std::cout << "build finite cells done (" << tr.tds().cells().size() << " cells)" << std::endl;
+    }
+
     if (!CGAL::SMDS_3::build_infinite_cells<Tr>(tr, incident_cells_map, verbose, allow_non_manifold))
     {
-      if(verbose) std::cout << "build_infinite_cells went wrong" << std::endl;
+      if(verbose)
+        std::cerr << "Error: build_infinite_cells went wrong!" << std::endl;
       success = false;
     }
-    else
-      if (verbose) std::cout << "build infinite cells done" << std::endl;
+    else if(verbose)
+    {
+      std::cout << "build infinite cells done (" << tr.tds().cells().size() << " cells)" << std::endl;
+    }
+
     tr.tds().set_dimension(3);
+
     if (!CGAL::SMDS_3::assign_neighbors<Tr>(tr, incident_cells_map, allow_non_manifold))
     {
-      if(verbose) std::cout << "assign_neighbors went wrong" << std::endl;
+      if(verbose)
+        std::cerr << "Error: assign_neighbors went wrong!" << std::endl;
       success = false;
     }
-    else
-      if (verbose) std::cout << "assign neighbors done" << std::endl;
+    else if(verbose)
+    {
+      std::cout << "assign neighbors done" << std::endl;
+    }
+
     if (verbose)
     {
-      std::cout << "built triangulation : " << std::endl;
-      std::cout << tr.number_of_cells() << " cells" << std::endl;
+      std::cout << "built triangulation!" << std::endl;
     }
   }
-  if(verbose)
-    std::cout << tr.number_of_vertices() << " vertices" << std::endl;
 
-  return success;// tr.tds().is_valid();
-              //TDS not valid when cells do not cover the convex hull of vertices
+  // disabled because the TDS is not valid when cells do not cover the convex hull of vertices
+  // return tr.tds().is_valid();
+
+  return success;
+
 }
 
 template<class Tr,
@@ -473,14 +497,14 @@ template<class Tr,
          typename CellRange,
          typename FacetPatchMap>
 bool build_triangulation_one_subdomain(Tr& tr,
-    const PointRange& points,
-    const CellRange& finite_cells,
-    const typename Tr::Cell::Subdomain_index& subdomain,
-    const FacetPatchMap& border_facets,
-    std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
-    const bool verbose,// = false,
-    const bool replace_domain_0,// = false
-    const bool allow_non_manifold)// = false
+                                       const PointRange& points,
+                                       const CellRange& finite_cells,
+                                       const typename Tr::Cell::Subdomain_index& subdomain,
+                                       const FacetPatchMap& border_facets,
+                                       std::vector<typename Tr::Vertex_handle>& vertex_handle_vector,
+                                       const bool verbose,// = false,
+                                       const bool replace_domain_0,// = false
+                                       const bool allow_non_manifold)// = false
 {
   std::vector<typename Tr::Cell::Subdomain_index> subdomains(finite_cells.size(), subdomain);
   return build_triangulation_impl(tr, points, finite_cells, subdomains,
@@ -494,13 +518,13 @@ template<class Tr,
          typename CellRange,
          typename FacetPatchMap>
 bool build_triangulation_one_subdomain(Tr& tr,
-    const PointRange& points,
-    const CellRange& finite_cells,
-    const typename Tr::Cell::Subdomain_index& subdomain,
-    const FacetPatchMap& border_facets,
-    const bool verbose,// = false,
-    const bool replace_domain_0,// = false
-    const bool allow_non_manifold)//= false
+                                       const PointRange& points,
+                                       const CellRange& finite_cells,
+                                       const typename Tr::Cell::Subdomain_index& subdomain,
+                                       const FacetPatchMap& border_facets,
+                                       const bool verbose,// = false,
+                                       const bool replace_domain_0,// = false
+                                       const bool allow_non_manifold)//= false
 {
   std::vector<typename Tr::Cell::Subdomain_index> subdomains(finite_cells.size(), subdomain);
   std::vector<typename Tr::Vertex_handle> vertex_handle_vector;
@@ -516,13 +540,13 @@ template<class Tr,
          typename SubdomainsRange,
          typename FacetPatchMap>
 bool build_triangulation_with_subdomains_range(Tr& tr,
-    const PointRange& points,
-    const CellRange& finite_cells,
-    const SubdomainsRange& subdomains,
-    const FacetPatchMap& border_facets,
-    const bool verbose,// = false
-    const bool replace_domain_0,// = false,
-    const bool allow_non_manifold)
+                                               const PointRange& points,
+                                               const CellRange& finite_cells,
+                                               const SubdomainsRange& subdomains,
+                                               const FacetPatchMap& border_facets,
+                                               const bool verbose,// = false
+                                               const bool replace_domain_0,// = false,
+                                               const bool allow_non_manifold)
 {
   std::vector<typename Tr::Vertex_handle> vertex_handle_vector;
   std::vector<typename Tr::Cell::Subdomain_index> subdomains_vector(
@@ -546,12 +570,14 @@ bool build_triangulation_from_file(std::istream& is,
   using Facet        = std::array<int, 3>; // 3 = id
   using Tet_with_ref = std::array<int, 4>; // 4 = id
 
+  if(!is)
+    return false;
+
   std::vector<Tet_with_ref> finite_cells;
   std::vector<Subdomain_index> subdomains;
   std::vector<Point_3> points;
   boost::unordered_map<Facet, typename Tr::Cell::Surface_patch_index> border_facets;
 
-  // grab the vertices
   int dim;
   int nv, nf, ntet, ref;
   std::string word;
@@ -562,9 +588,13 @@ bool build_triangulation_from_file(std::istream& is,
   CGAL_assertion(dim == 3);
 
   if(verbose)
+  {
     std::cout << "Reading .mesh file..." << std::endl;
+    std::cout << "Replace domain #0 = " << replace_domain_0 << std::endl;
+    std::cout << "Allow non-manifoldness = " << allow_non_manifold << std::endl;
+  }
 
-  bool dont_replace_domain_0 = false;
+  bool is_CGAL_mesh = false;
 
   while(is >> word && word != "End")
   {
@@ -572,22 +602,30 @@ bool build_triangulation_from_file(std::istream& is,
     {
       is >> word;
       if (word == "End")
+      {
         break;
+      }
       else if (word == "CGAL::Mesh_complex_3_in_triangulation_3")
       {
-        dont_replace_domain_0 = true;//with CGAL meshes, domain 0 should be kept
+        is_CGAL_mesh = true; // with CGAL meshes, domain 0 should be kept
         continue;
       }
       //else skip other comments
     }
+
     if(word == "Vertices")
     {
       is >> nv;
       for(int i=0; i<nv; ++i)
       {
         typename Tr::Geom_traits::FT x,y,z;
-        is >> x >> y >> z >> ref;
-        points.push_back(Point_3(x,y,z));
+        if(!(is >> x >> y >> z >> ref))
+        {
+          if(verbose)
+            std::cerr << "Issue while reading vertices" << std::endl;
+          return false;
+        }
+        points.emplace_back(x,y,z);
       }
     }
 
@@ -596,25 +634,39 @@ bool build_triangulation_from_file(std::istream& is,
       is >> nf;
       for(int i=0; i<nf; ++i)
       {
-        int n1, n2, n3;
+        int n[3];
         typename Tr::Cell::Surface_patch_index surface_patch_id;
-        is >> n1 >> n2 >> n3 >> surface_patch_id;
+        if(!(is >> n[0] >> n[1] >> n[2] >> surface_patch_id))
+        {
+          if(verbose)
+            std::cerr << "Issue while reading triangles" << std::endl;
+          return false;
+        }
+
         Facet facet;
-        facet[0] = n1 - 1;
-        facet[1] = n2 - 1;
-        facet[2] = n3 - 1;
-        //find the circular permutation that puts the smallest index in the first place.
-        int n0 = (std::min)((std::min)(facet[0],facet[1]), facet[2]);
-        int k=0;
-        Facet f;
+        facet[0] = n[0] - 1;
+        facet[1] = n[1] - 1;
+        facet[2] = n[2] - 1;
+
+        if(verbose)
+          std::cout << "Looking at face #" << i << ": " << n[0] << " " << n[1] << " " << n[2] << std::endl;
+
+        CGAL_warning_code(
+        for(int j=0; j<3; ++j)
+          for(int k=0; k<3; ++k)
+            if(j != k)
+              CGAL_warning(n[j] != n[k]);
+        )
+
+        // find the circular permutation that puts the smallest index in the first place.
+        int n0 = (std::min)({facet[0],facet[1], facet[2]});
         do
         {
-          f[0] = facet[(0+k)%3];
-          f[1] = facet[(1+k)%3];
-          f[2] = facet[(2+k)%3];
-          ++k;
-        } while(f[0] != n0);
-        border_facets.insert(std::make_pair(f, surface_patch_id));
+          std::rotate(std::begin(facet), std::next(std::begin(facet)), std::end(facet));
+        }
+        while(facet[0] != n0);
+
+        border_facets.emplace(facet, surface_patch_id);
       }
     }
 
@@ -623,13 +675,32 @@ bool build_triangulation_from_file(std::istream& is,
       is >> ntet;
       for(int i=0; i<ntet; ++i)
       {
-        int n0, n1, n2, n3, reference;
-        is >> n0 >> n1 >> n2 >> n3 >> reference;
+        int n[4];
+        int reference;
+
+        if(!(is >> n[0] >> n[1] >> n[2] >> n[3] >> reference))
+        {
+          if(verbose)
+            std::cerr << "Issue while reading tetrahedra" << std::endl;
+          return false;
+        }
+
+        if(verbose)
+          std::cout << "Looking at tet #" << i << ": " << n[0] << " " << n[1] << " " << n[2] << " " << n[3] << std::endl;
+
+        CGAL_warning_code(
+        for(int j=0; j<4; ++j)
+          for(int k=0; k<4; ++k)
+            if(j != k)
+              CGAL_warning(n[j] != n[k]);
+        )
+
         Tet_with_ref t;
-        t[0] = n0 - 1;
-        t[1] = n1 - 1;
-        t[2] = n2 - 1;
-        t[3] = n3 - 1;
+        t[0] = n[0] - 1;
+        t[1] = n[1] - 1;
+        t[2] = n[2] - 1;
+        t[3] = n[3] - 1;
+
         finite_cells.push_back(t);
         subdomains.push_back(reference);
       }
@@ -645,18 +716,17 @@ bool build_triangulation_from_file(std::istream& is,
 
   if(finite_cells.empty())
     return false;
+
   CGAL_assertion(finite_cells.size() == subdomains.size());
 
   return build_triangulation_with_subdomains_range(tr,
-                      points, finite_cells, subdomains, border_facets,
-                      verbose,
-                      replace_domain_0 && !dont_replace_domain_0,
-                      allow_non_manifold);
+                                                   points, finite_cells, subdomains, border_facets,
+                                                   verbose,
+                                                   replace_domain_0 && !is_CGAL_mesh,
+                                                   allow_non_manifold);
 }
 
-}  // namespace SMDS_3
-}  // namespace CGAL
-
-#include <CGAL/enable_warnings.h>
+} // namespace SMDS_3
+} // namespace CGAL
 
 #endif // CGAL_SMDS_3_TET_SOUP_TO_C3T3_H

Scripts/developer_scripts/cgal_check_dependencies.sh

@@ -50,7 +50,7 @@ do
     fi
     PKG_DIFF=$(grep -Fxv -f "$pkg_path/package_info/$pkg/dependencies" "$pkg_path/package_info/$pkg/dependencies.old" || true)
     if [ -n "$PKG_DIFF" ]; then
-      TOTAL_RES="Differences in $pkg:\n$PKG_DIFF\mhave disappeared.\n$TOTAL_RES"
+      TOTAL_RES="Differences in $pkg:\n$PKG_DIFF\nhave disappeared.\n$TOTAL_RES"
     fi
     if [ -f $pkg_path/package_info/$pkg/dependencies.old ]; then
       rm $pkg_path/package_info/$pkg/dependencies.old

TDS_3/include/CGAL/Triangulation_data_structure_3.h

@@ -3935,21 +3935,21 @@ is_valid(Cell_handle c, bool verbose, int level) const
 
           int j1n=4,j2n=4,j3n=4;
           if ( ! n->has_vertex(c->vertex((i+1)&3),j1n) ) {
-            if (verbose) { std::cerr << "vertex " << ((i+1)&3)
+            if (verbose) { std::cerr << "vertex (+1) " << ((i+1)&3)
                                      << " not vertex of neighbor "
                                      << i << std::endl; }
             CGAL_assertion(false);
             return false;
           }
           if ( ! n->has_vertex(c->vertex((i+2)&3),j2n) ) {
-            if (verbose) { std::cerr << "vertex " << ((i+2)&3)
+            if (verbose) { std::cerr << "vertex (+2) " << ((i+2)&3)
                                      << " not vertex of neighbor "
                                      << i << std::endl; }
             CGAL_assertion(false);
             return false;
           }
           if ( ! n->has_vertex(c->vertex((i+3)&3),j3n) ) {
-            if (verbose) { std::cerr << "vertex " << ((i+3)&3)
+            if (verbose) { std::cerr << "vertex (+3) " << ((i+3)&3)
                                      << " not vertex of neighbor "
                                      << i << std::endl; }
             CGAL_assertion(false);

Triangulation_3/include/CGAL/Regular_triangulation_3.h

@@ -1003,9 +1003,10 @@ public:
 
   void dual_segment(Cell_handle c, int i, Bare_point& p, Bare_point&q) const;
   void dual_segment(const Facet& facet, Bare_point& p, Bare_point&q) const;
-  void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const;
   void dual_ray(Cell_handle c, int i, Ray& ray) const;
   void dual_ray(const Facet& facet, Ray& ray) const;
+  void dual_exact(const Facet& facet, const Weighted_point& p, Bare_point&q) const;
+  void dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const;
   void dual_ray_exact(const Facet& facet, Ray& ray) const;
 
   template < class Stream>
@@ -1825,14 +1826,42 @@ dual_ray(const Facet& facet, Ray& ray) const
   return dual_ray(facet.first, facet.second, ray);
 }
 
-// Exact versions of dual_segment() and dual_ray() for Mesh_3.
+// Exact versions of dual(), dual_segment(), and dual_ray() for Mesh_3.
 // These functions are really dirty: they assume that the point type is nice enough
 // such that EPECK can manipulate it (e.g. convert it to EPECK::Point_3) AND
 // that the result of these manipulations will make sense.
 template < class Gt, class Tds, class Lds >
 void
 Regular_triangulation_3<Gt,Tds,Lds>::
-dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point&q) const
+dual_exact(const Facet& f, const Weighted_point& s, Bare_point& cc) const
+{
+  typedef typename Kernel_traits<Bare_point>::Kernel           K;
+  typedef Exact_predicates_exact_constructions_kernel          EK;
+  typedef Cartesian_converter<K, EK>                           To_exact;
+  typedef Cartesian_converter<EK,K>                            Back_from_exact;
+
+  typedef EK                                                   Exact_Rt;
+
+  To_exact to_exact;
+  Back_from_exact back_from_exact;
+  Exact_Rt::Construct_weighted_circumcenter_3 exact_weighted_circumcenter =
+      Exact_Rt().construct_weighted_circumcenter_3_object();
+
+  const Cell_handle c = f.first;
+  const int i = f.second;
+
+  const typename Exact_Rt::Weighted_point_3& cp = to_exact(c->vertex((i+1)%4)->point());
+  const typename Exact_Rt::Weighted_point_3& cq = to_exact(c->vertex((i+2)%4)->point());
+  const typename Exact_Rt::Weighted_point_3& cr = to_exact(c->vertex((i+3)%4)->point());
+  const typename Exact_Rt::Weighted_point_3& cs = to_exact(s);
+
+  cc = back_from_exact(exact_weighted_circumcenter(cp, cq, cr, cs));
+}
+
+template < class Gt, class Tds, class Lds >
+void
+Regular_triangulation_3<Gt,Tds,Lds>::
+dual_segment_exact(const Facet& facet, Bare_point& p, Bare_point& q) const
 {
   typedef typename Kernel_traits<Bare_point>::Kernel           K;
   typedef Exact_predicates_exact_constructions_kernel          EK;