Tet remeshing for CCDT_3

Constrained_triangulation_3/include/CGAL/Conforming_Delaunay_triangulation_3.h

@@ -72,7 +72,10 @@ namespace CGAL {
 namespace CDT_3 {
 
 struct Debug_options {
-  enum class Flags {
+  enum class Flag;
+  bool get(Flag f) const { return flags[static_cast<int>(f)]; }
+  auto get(Flag f)  { return flags[static_cast<int>(f)]; }
+  enum class Flag {
     Steiner_points = 0,
     Steiner_points_construction,
     conforming,
@@ -91,71 +94,77 @@ struct Debug_options {
     debug_geometric_errors,
     debug_polygon_insertion,
     debug_restore_faces,
+    debug_move_Steiner_vertices,
     display_statistics,
     use_epeck_for_normals,
     use_epeck_for_Steiner_points,
     nb_of_flags
   };
-  bool Steiner_points_construction() const { return flags[static_cast<int>(Flags::Steiner_points_construction)]; }
-  void Steiner_points_construction(bool b) { flags.set(static_cast<int>(Flags::Steiner_points_construction), b); }
 
-  bool Steiner_points() const { return flags[static_cast<int>(Flags::Steiner_points)]; }
-  void Steiner_points(bool b) { flags.set(static_cast<int>(Flags::Steiner_points), b); }
+  using Flags = std::bitset<static_cast<int>(Flag::nb_of_flags)>;
 
-  bool input_faces() const { return flags[static_cast<int>(Flags::input_faces)]; }
-  void input_faces(bool b) { flags.set(static_cast<int>(Flags::input_faces), b); }
+  bool Steiner_points_construction() const { return get(Flag::Steiner_points_construction); }
+  void Steiner_points_construction(bool b) { get(Flag::Steiner_points_construction) = b; }
 
-  bool missing_region() const { return flags[static_cast<int>(Flags::missing_region)]; }
-  void missing_region(bool b) { flags.set(static_cast<int>(Flags::missing_region), b); }
+  bool Steiner_points() const { return get(Flag::Steiner_points); }
+  void Steiner_points(bool b) { get(Flag::Steiner_points) = b; }
 
-  bool regions() const { return flags[static_cast<int>(Flags::regions)]; }
-  void regions(bool b) { flags.set(static_cast<int>(Flags::regions), b); }
+  bool input_faces() const { return get(Flag::input_faces); }
+  void input_faces(bool b) { get(Flag::input_faces) = b; }
 
-  bool copy_triangulation_into_hole() const { return flags[static_cast<int>(Flags::copy_triangulation_into_hole)]; }
-  void copy_triangulation_into_hole(bool b) { flags.set(static_cast<int>(Flags::copy_triangulation_into_hole), b); }
+  bool missing_region() const { return get(Flag::missing_region); }
+  void missing_region(bool b) { get(Flag::missing_region) = b; }
 
-  bool validity() const { return flags[static_cast<int>(Flags::validity)]; }
-  void validity(bool b) { flags.set(static_cast<int>(Flags::validity), b); }
+  bool regions() const { return get(Flag::regions); }
+  void regions(bool b) { get(Flag::regions) = b; }
 
+  bool copy_triangulation_into_hole() const { return get(Flag::copy_triangulation_into_hole); }
+  void copy_triangulation_into_hole(bool b) { get(Flag::copy_triangulation_into_hole) = b; }
 
-  bool finite_edges_map() const { return flags[static_cast<int>(Flags::debug_finite_edges_map)]; }
-  void finite_edges_map(bool b) { flags.set(static_cast<int>(Flags::debug_finite_edges_map), b); }
+  bool validity() const { return get(Flag::validity); }
+  void validity(bool b) { get(Flag::validity) = b; }
 
-  bool subconstraints_to_conform() const { return flags[static_cast<int>(Flags::debug_subconstraints_to_conform)]; }
-  void subconstraints_to_conform(bool b) { flags.set(static_cast<int>(Flags::debug_subconstraints_to_conform), b); }
+  bool finite_edges_map() const { return get(Flag::debug_finite_edges_map); }
+  void finite_edges_map(bool b) { get(Flag::debug_finite_edges_map) = b; }
 
-  bool use_finite_edges_map_flag() const { return flags[static_cast<int>(Flags::use_finite_edges_map)]; }
-  void use_finite_edges_map(bool b) { flags.set(static_cast<int>(Flags::use_finite_edges_map), b); }
+  bool subconstraints_to_conform() const { return get(Flag::debug_subconstraints_to_conform); }
+  void subconstraints_to_conform(bool b) { get(Flag::debug_subconstraints_to_conform) = b; }
 
-  bool verbose_special_cases() const { return flags[static_cast<int>(Flags::verbose_special_cases)]; }
-  void verbose_special_cases(bool b) { flags.set(static_cast<int>(Flags::verbose_special_cases), b); }
+  bool use_finite_edges_map_flag() const { return get(Flag::use_finite_edges_map); }
+  void use_finite_edges_map(bool b) { get(Flag::use_finite_edges_map) = b; }
 
-  bool encroaching_vertices() const { return flags[static_cast<int>(Flags::debug_encroaching_vertices)]; }
-  void encroaching_vertices(bool b) { flags.set(static_cast<int>(Flags::debug_encroaching_vertices), b); }
+  bool verbose_special_cases() const { return get(Flag::verbose_special_cases); }
+  void verbose_special_cases(bool b) { get(Flag::verbose_special_cases) = b; }
 
-  bool conforming_validation() const { return flags[static_cast<int>(Flags::debug_conforming_validation)]; }
-  void conforming_validation(bool b) { flags.set(static_cast<int>(Flags::debug_conforming_validation), b); }
+  bool encroaching_vertices() const { return get(Flag::debug_encroaching_vertices); }
+  void encroaching_vertices(bool b) { get(Flag::debug_encroaching_vertices) = b; }
 
-  bool constraint_hierarchy() const { return flags[static_cast<int>(Flags::debug_constraint_hierarchy)]; }
-  void constraint_hierarchy(bool b) { flags.set(static_cast<int>(Flags::debug_constraint_hierarchy), b); }
+  bool conforming_validation() const { return get(Flag::debug_conforming_validation); }
+  void conforming_validation(bool b) { get(Flag::debug_conforming_validation) = b; }
 
-  bool geometric_errors() const { return flags[static_cast<int>(Flags::debug_geometric_errors)]; }
-  void geometric_errors(bool b) { flags.set(static_cast<int>(Flags::debug_geometric_errors), b); }
+  bool constraint_hierarchy() const { return get(Flag::debug_constraint_hierarchy); }
+  void constraint_hierarchy(bool b) { get(Flag::debug_constraint_hierarchy) = b; }
 
-  bool polygon_insertion() const { return flags[static_cast<int>(Flags::debug_polygon_insertion)]; }
-  void polygon_insertion(bool b) { flags.set(static_cast<int>(Flags::debug_polygon_insertion), b); }
+  bool geometric_errors() const { return get(Flag::debug_geometric_errors); }
+  void geometric_errors(bool b) { get(Flag::debug_geometric_errors) = b; }
 
-  bool restore_faces() const { return flags[static_cast<int>(Flags::debug_restore_faces)]; }
-  void restore_faces(bool b) { flags.set(static_cast<int>(Flags::debug_restore_faces), b); }
+  bool polygon_insertion() const { return get(Flag::debug_polygon_insertion); }
+  void polygon_insertion(bool b) { get(Flag::debug_polygon_insertion) = b; }
 
-  bool display_statistics() const { return flags[static_cast<int>(Flags::display_statistics)]; }
-  void display_statistics(bool b) { flags.set(static_cast<int>(Flags::display_statistics), b); }
+  bool restore_faces() const { return get(Flag::debug_restore_faces); }
+  void restore_faces(bool b) { get(Flag::debug_restore_faces) = b; }
 
-  bool use_epeck_for_normals() const { return flags[static_cast<int>(Flags::use_epeck_for_normals)]; }
-  void use_epeck_for_normals(bool b) { flags.set(static_cast<int>(Flags::use_epeck_for_normals), b); }
+  bool move_Steiner_vertices() const { return get(Flag::debug_move_Steiner_vertices); }
+  void move_Steiner_vertices(bool b) { get(Flag::debug_move_Steiner_vertices) = b; }
 
-  bool use_epeck_for_Steiner_points() const { return flags[static_cast<int>(Flags::use_epeck_for_Steiner_points)]; }
-  void use_epeck_for_Steiner_points(bool b) { flags.set(static_cast<int>(Flags::use_epeck_for_Steiner_points), b); }
+  bool display_statistics() const { return get(Flag::display_statistics); }
+  void display_statistics(bool b) { get(Flag::display_statistics) = b; }
+
+  bool use_epeck_for_normals() const { return get(Flag::use_epeck_for_normals); }
+  void use_epeck_for_normals(bool b) { get(Flag::use_epeck_for_normals) = b; }
+
+  bool use_epeck_for_Steiner_points() const { return get(Flag::use_epeck_for_Steiner_points); }
+  void use_epeck_for_Steiner_points(bool b) { get(Flag::use_epeck_for_Steiner_points) = b; }
 
   double segment_vertex_epsilon() const { return segment_vertex_epsilon_; }
   void set_segment_vertex_epsilon(double eps) { segment_vertex_epsilon_ = eps; }
@@ -164,7 +173,7 @@ struct Debug_options {
   void set_vertex_vertex_epsilon(double eps) { vertex_vertex_epsilon_ = eps; }
 
 private:
-  std::bitset<static_cast<int>(Flags::nb_of_flags)> flags{};
+  Flags flags{};
   double segment_vertex_epsilon_ = 0.0;
   double vertex_vertex_epsilon_  = 0.0;
 }; // end struct Debug_options
@@ -173,32 +182,34 @@ namespace internal {
 
 inline auto& tasks_manager() {
   struct Tasks_manager {
-    enum {
+    enum class Tm {
       READ_INPUT = 0,
       MERGE_FACETS,
       INSERT_VERTICES,
       COMPUTE_DISTANCES,
       CONFORMING,
       CDT,
+      MOVE_STEINER_VERTICES,
       OUTPUT,
       VALIDATION,
       NB_TASKS
     };
 
   #if CGAL_USE_ITT
     __itt_domain* cdt_3_domain = __itt_domain_create("org.cgal.CDT_3");
-    const std::array<__itt_string_handle*, NB_TASKS> task_handles = {
+    const std::array<__itt_string_handle*, static_cast<std::size_t>(Tm::NB_TASKS)> task_handles = {
       __itt_string_handle_create("CDT_3: read input file"),
       __itt_string_handle_create("CDT_3: merge facets"),
       __itt_string_handle_create("CDT_3: insert vertices"),
       __itt_string_handle_create("CDT_3: compute distances"),
       __itt_string_handle_create("CDT_3: conforming"),
       __itt_string_handle_create("CDT_3: cdt"),
+      __itt_string_handle_create("CDT_3: move Steiner vertices"),
       __itt_string_handle_create("CDT_3: outputs"),
       __itt_string_handle_create("CDT_3: validation")
     };
   #endif
-    std::array<CGAL::Real_timer, NB_TASKS> timers{};
+    std::array<CGAL::Real_timer, static_cast<std::size_t>(Tm::NB_TASKS)> timers{};
     struct Scope_guard {
       Tasks_manager *instance = nullptr;
       int task_id;
@@ -222,18 +233,19 @@ inline auto& tasks_manager() {
       }
     };
 
-    Scope_guard make_task_scope_guard(int task_id) {
-      return Scope_guard(this, task_id);
+    Scope_guard make_task_scope_guard(Tm task_id) {
+      return Scope_guard(this, static_cast<int>(task_id));
     }
 
-    Scope_guard READ_INPUT_TASK_guard() { return make_task_scope_guard(READ_INPUT); }
-    Scope_guard MERGE_FACETS_TASK_guard() { return make_task_scope_guard(MERGE_FACETS); }
-    Scope_guard INSERT_VERTICES_TASK_guard() { return make_task_scope_guard(INSERT_VERTICES); }
-    Scope_guard COMPUTE_DISTANCES_TASK_guard() { return make_task_scope_guard(COMPUTE_DISTANCES); }
-    Scope_guard CONFORMING_TASK_guard() { return make_task_scope_guard(CONFORMING); }
-    Scope_guard CDT_TASK_guard() { return make_task_scope_guard(CDT); }
-    Scope_guard OUTPUT_TASK_guard() { return make_task_scope_guard(OUTPUT); }
-    Scope_guard VALIDATION_TASK_guard() { return make_task_scope_guard(VALIDATION); }
+    Scope_guard READ_INPUT_TASK_guard() { return make_task_scope_guard(Tm::READ_INPUT); }
+    Scope_guard MERGE_FACETS_TASK_guard() { return make_task_scope_guard(Tm::MERGE_FACETS); }
+    Scope_guard INSERT_VERTICES_TASK_guard() { return make_task_scope_guard(Tm::INSERT_VERTICES); }
+    Scope_guard COMPUTE_DISTANCES_TASK_guard() { return make_task_scope_guard(Tm::COMPUTE_DISTANCES); }
+    Scope_guard CONFORMING_TASK_guard() { return make_task_scope_guard(Tm::CONFORMING); }
+    Scope_guard CDT_TASK_guard() { return make_task_scope_guard(Tm::CDT); }
+    Scope_guard MOVE_STEINER_VERTICES_TASK_guard() { return make_task_scope_guard(Tm::MOVE_STEINER_VERTICES); }
+    Scope_guard OUTPUT_TASK_guard() { return make_task_scope_guard(Tm::OUTPUT); }
+    Scope_guard VALIDATION_TASK_guard() { return make_task_scope_guard(Tm::VALIDATION); }
 
   }; // end struct Intel_OneAPI_ITT_API
 
@@ -273,6 +285,10 @@ inline auto CDT_3_OUTPUT_TASK_guard() {
   return CDT_3::internal::tasks_manager().OUTPUT_TASK_guard();
 }
 
+inline auto CDT_3_MOVE_STEINER_VERTICES_TASK_guard() {
+  return CDT_3::internal::tasks_manager().MOVE_STEINER_VERTICES_TASK_guard();
+}
+
 inline auto CDT_3_VALIDATION_TASK_guard() {
   return CDT_3::internal::tasks_manager().VALIDATION_TASK_guard();
 }
@@ -318,10 +334,9 @@ class Conforming_Delaunay_triangulation_3 : public T_3 {
   using Subconstraint = typename Constraint_hierarchy::Subconstraint;
 
   auto display_vert(Vertex_handle v) const{
-    std::stringstream os;
-    os.precision(17);
-    os << IO::oformat(v, with_point);
-    return os.str();
+    return IO::oformat([&](std::ostream& os) -> auto& {
+      return os << IO::oformat(v, with_point);
+    }, IO_manip_tag{});
   }
 
   auto display_subcstr(Subconstraint subconstraint) const {
@@ -863,6 +878,7 @@ class Conforming_Delaunay_triangulation_3 : public T_3 {
     return false;
   }
 
+public:
   Constrained_polyline_id constraint_from_extremities(Vertex_handle va, Vertex_handle vb) const {
     if(va->ccdt_3_data().number_of_incident_constraints() == 0 || vb->ccdt_3_data().number_of_incident_constraints() == 0)
     {
@@ -882,15 +898,16 @@ class Conforming_Delaunay_triangulation_3 : public T_3 {
     return c_id;
   }
 
+protected:
   auto constraint_extremities(Constrained_polyline_id c_id) const {
-      CGAL_assertion(this->constraint_hierarchy.vertices_in_constraint_begin(c_id) !=
-                     this->constraint_hierarchy.vertices_in_constraint_end(c_id));
+      CGAL_assertion(this->constraint_hierarchy.is_valid_constraint(c_id));
       if(debug().constraint_hierarchy()) {
         std::cerr << "constraint " << static_cast<void*>(c_id.vl_ptr()) << " has "
                   << c_id.vl_ptr()->skip_size() << " vertices\n";
       }
       const auto begin = this->constraint_hierarchy.vertices_in_constraint_begin(c_id);
       const auto end = this->constraint_hierarchy.vertices_in_constraint_end(c_id);
+      CGAL_assertion(begin != end);
       const auto c_va = *begin;
       const auto c_vb = *std::prev(end);
     return std::make_pair(c_va, c_vb);