Merge pull request #8072 from janetournois/Mesh_3-use_one_step_loop-j…

…tournois

Mesh_3 - fix mesher using the `one_step()` loop

Data/data/meshes/halfcube.off

@@ -0,0 +1,13 @@
+OFF
+5 4 0
+
+-1 -1 1
+-1 1 1
+1 -1 -1
+-1 1 -1
+-1 -1 -1
+3  0 4 2
+3  2 4 3
+3  4 0 1
+3  1 3 4
+

Mesh_3/include/CGAL/Mesh_3/Mesher_3.h

@@ -404,6 +404,13 @@ Mesher_3<C3T3,MC,MD>::Mesher_3(C3T3& c3t3,
   cells_mesher_.set_stop_pointer(stop_ptr);
   facets_mesher_.set_stop_pointer(stop_ptr);
 #endif
+
+  // First surface mesh could modify c3t3 without notifying cells_mesher
+  // So we have to ensure that no old cell will be left in c3t3
+  // Second, the c3t3 object could have been corrupted since the last call
+  // to `refine_mesh`, for example by inserting new vertices in the
+  // triangulation.
+  r_c3t3_.clear_cells_and_facets_from_c3t3();
 }
 
 
@@ -424,13 +431,6 @@ refine_mesh(std::string dump_after_refine_surface_prefix)
   auto refine_volume_mesh_task_handle = __itt_string_handle_create("Mesher_3 refine volume mesh");
 #endif // CGAL_MESH_3_USE_INTEL_ITT
 
-  // First surface mesh could modify c3t3 without notifying cells_mesher
-  // So we have to ensure that no old cell will be left in c3t3
-  // Second, the c3t3 object could have been corrupted since the last call
-  // to `refine_mesh`, for example by inserting new vertices in the
-  // triangulation.
-  r_c3t3_.clear_cells_and_facets_from_c3t3();
-
   const Triangulation& r_tr = r_c3t3_.triangulation();
   CGAL_USE(r_tr);
 

Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h

@@ -289,7 +289,6 @@ class Refine_facets_3_base
     return Container_::no_longer_element_to_refine_impl();
   }
 
-  // Gets the point to insert from the element to refine
   /// Gets the point to insert from the element to refine
   Bare_point refinement_point_impl(const Facet& facet) const
   {

Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h

@@ -382,11 +382,16 @@ inside_protecting_balls(const Tr& tr,
                         const Vertex_handle v,
                         const Bare_point& p) const
 {
+  if(tr.number_of_vertices() == 0)
+    return false;
+
   typename Gt::Compare_weighted_squared_radius_3 cwsr =
     tr.geom_traits().compare_weighted_squared_radius_3_object();
 
-  Vertex_handle nv = tr.nearest_power_vertex(p, v->cell());
+  Cell_handle hint = (v == Vertex_handle()) ? Cell_handle() : v->cell();
+  Vertex_handle nv = tr.nearest_power_vertex(p, hint);
   const Point& nvwp = tr.point(nv);
+
   if(cwsr(nvwp, FT(0)) == CGAL::SMALLER)
   {
     typename Tr::Geom_traits::Construct_point_3 cp = tr.geom_traits().construct_point_3_object();

Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h

@@ -296,7 +296,7 @@ class Polyhedral_mesh_domain_3
       }
       tree_.build();
     }
-    bounding_tree_ = 0;
+    set_surface_only();
   }
 
   /// Destructor

Mesh_3/include/CGAL/make_mesh_3.h

@@ -45,9 +45,10 @@ init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&,
 {
   typedef typename MeshDomain::Point_3 Point_3;
   typedef typename MeshDomain::Index Index;
-  typedef std::vector<std::pair<Point_3, Index> > Initial_points_vector;
-  typedef typename Initial_points_vector::iterator Ipv_iterator;
+  typedef typename std::pair<Point_3, Index> PI;
+  typedef std::vector<PI> Initial_points_vector;
   typedef typename C3T3::Vertex_handle Vertex_handle;
+  typedef CGAL::Mesh_3::Triangulation_helpers<typename C3T3::Triangulation> Th;
 
   // Mesh initialization : get some points and add them to the mesh
   Initial_points_vector initial_points;
@@ -61,17 +62,18 @@ init_c3t3(C3T3& c3t3, const MeshDomain& domain, const MeshCriteria&,
       c3t3.triangulation().geom_traits().construct_weighted_point_3_object();
 
   // Insert points and set their index and dimension
-  for ( Ipv_iterator it = initial_points.begin() ;
-       it != initial_points.end() ;
-       ++it )
+  for (const PI& pi : initial_points)
   {
-    Vertex_handle v = c3t3.triangulation().insert(cwp(it->first));
+    if(Th().inside_protecting_balls(c3t3.triangulation(), Vertex_handle(), pi.first))
+      continue;
+
+    Vertex_handle v = c3t3.triangulation().insert(cwp(pi.first));
 
     // v could be null if point is hidden
     if ( v != Vertex_handle() )
     {
       c3t3.set_dimension(v,2); // by construction, points are on surface
-      c3t3.set_index(v,it->second);
+      c3t3.set_index(v, pi.second);
     }
   }
 }

Mesh_3/test/Mesh_3/CMakeLists.txt

@@ -53,6 +53,7 @@ create_single_source_cgal_program( "test_meshing_without_features_determinism.cp
 create_single_source_cgal_program( "test_mesh_3_issue_1554.cpp" )
 create_single_source_cgal_program( "test_mesh_polyhedral_domain_with_features_deprecated.cpp" )
 create_single_source_cgal_program( "test_meshing_with_one_step.cpp" )
+create_single_source_cgal_program( "test_meshing_with_one_step_with_features.cpp" )
 create_single_source_cgal_program( "test_mesh_cell_base_3.cpp")
 create_single_source_cgal_program( "test_min_edge_length.cpp")
 
@@ -85,6 +86,7 @@ foreach(target
     test_mesh_polyhedral_domain_with_features_deprecated
     test_mesh_cell_base_3
     test_meshing_with_one_step
+    test_meshing_with_one_step_with_features
     test_min_edge_length
     test_min_size_criteria
     )

Mesh_3/test/Mesh_3/test_meshing_with_one_step_with_features.cpp

@@ -0,0 +1,86 @@
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+#include <CGAL/Mesh_triangulation_3.h>
+#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
+#include <CGAL/Mesh_criteria_3.h>
+#include <CGAL/Surface_mesh.h>
+#include <CGAL/boost/graph/helpers.h>
+#include <CGAL/Polyhedral_mesh_domain_with_features_3.h>
+#include <CGAL/make_mesh_3.h>
+#include <CGAL/refine_mesh_3.h>
+
+// Domain
+typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
+typedef CGAL::Surface_mesh<K::Point_3> Polyhedron;
+typedef CGAL::Polyhedral_mesh_domain_with_features_3<K, Polyhedron> Mesh_domain;
+
+typedef CGAL::Sequential_tag Concurrency_tag;
+
+// Triangulation
+typedef CGAL::Mesh_triangulation_3<Mesh_domain,
+                                   CGAL::Default,
+                                   Concurrency_tag>::type Tr;
+
+typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
+
+// Criteria
+typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
+
+typedef CGAL::Mesh_3::Mesher_3<C3t3,
+                               Mesh_criteria,
+                               Mesh_domain>   Mesher;
+
+// To avoid verbose function and named parameters call
+using namespace CGAL::parameters;
+
+int main()
+{
+  const std::string fname = CGAL::data_file_path("meshes/halfcube.off");
+  // Create input polyhedron
+  Polyhedron polyhedron;
+  std::ifstream input(fname);
+  input >> polyhedron;
+  if(input.fail()){
+    std::cerr << "Error: Cannot read file " <<  fname << std::endl;
+    return EXIT_FAILURE;
+  }
+  input.close();
+
+  if (!CGAL::is_triangle_mesh(polyhedron)){
+    std::cerr << "Input geometry is not triangulated." << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  // Create domain
+
+  std::vector<Polyhedron*> polyhedra;
+  polyhedra.push_back(&polyhedron);
+
+  Mesh_domain domain(polyhedra.begin(), polyhedra.end());
+  domain.detect_features();
+
+  // Mesh criteria (no cell_size set)
+  Mesh_criteria criteria(facet_angle = 25,
+                         facet_size = 0.17,
+                         facet_distance = 0.017,
+                         edge_size = 0.17);
+
+  // Mesh generation
+  namespace p = CGAL::parameters;
+  C3t3 c3t3;
+  CGAL::Mesh_3::internal::C3t3_initializer<C3t3,
+    Mesh_domain,
+    Mesh_criteria,
+    true>()(c3t3, domain, criteria, true);
+
+  Mesher mesher(c3t3, domain, criteria, CGAL::FACET_VERTICES_ON_SURFACE);
+  mesher.initialize();
+  mesher.display_number_of_bad_elements();
+  while ( ! mesher.is_algorithm_done() ) mesher.one_step();
+  assert(c3t3.triangulation().number_of_vertices() > 200);
+
+  // Output
+  mesher.display_number_of_bad_elements();
+
+  return EXIT_SUCCESS;
+}

SMDS_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h

@@ -1423,28 +1423,32 @@ class Mesh_complex_3_in_triangulation_3
     }
   }
 
-  void clear_cells_and_facets_from_c3t3() {
-    for (typename Tr::Finite_cells_iterator
-      cit = this->triangulation().finite_cells_begin(),
-      end = this->triangulation().finite_cells_end();
-      cit != end; ++cit)
+  void clear_cells_and_facets_from_c3t3()
+  {
+    //clear cells
+    for (typename Tr::All_cells_iterator cit = this->triangulation().all_cells_begin();
+         cit != this->triangulation().all_cells_end();
+         ++cit)
     {
       set_subdomain_index(cit, Subdomain_index());
     }
     this->number_of_cells_ = 0;
-    for (typename Tr::Finite_facets_iterator
-      fit = this->triangulation().finite_facets_begin(),
-      end = this->triangulation().finite_facets_end();
-      fit != end; ++fit)
+
+    //clear facets
+    for (typename Tr::All_facets_iterator fit = this->triangulation().all_facets_begin();
+         fit != this->triangulation().all_facets_end();
+         ++fit)
     {
-      Facet facet = *fit;
+      const auto& facet = *fit;
       set_surface_patch_index(facet.first, facet.second, Surface_patch_index());
       if (this->triangulation().dimension() > 2) {
-        Facet mirror = tr_.mirror_facet(facet);
+        const Facet& mirror = tr_.mirror_facet(facet);
         set_surface_patch_index(mirror.first, mirror.second, Surface_patch_index());
       }
     }
     this->number_of_facets_ = 0;
+
+    //clear manifold info
     clear_manifold_info();
   }