diff --git a/Data/data/meshes/halfcube.off b/Data/data/meshes/halfcube.off new file mode 100644 index 00000000000..4c9a166bc42 --- /dev/null +++ b/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 + diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h index 185efe6229c..a4374383b5b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesher_3.h @@ -404,6 +404,13 @@ Mesher_3::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); diff --git a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h index 8374a106b62..b7d489a096b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Refine_facets_3.h @@ -289,7 +289,6 @@ public: 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 { diff --git a/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h b/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h index 6459aaeb047..c3b0b3ac034 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Triangulation_helpers.h +++ b/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(); diff --git a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h index 8c7e57e1dd9..efb5c975470 100644 --- a/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h +++ b/Mesh_3/include/CGAL/Polyhedral_mesh_domain_3.h @@ -296,7 +296,7 @@ public: } tree_.build(); } - bounding_tree_ = 0; + set_surface_only(); } /// Destructor diff --git a/Mesh_3/include/CGAL/make_mesh_3.h b/Mesh_3/include/CGAL/make_mesh_3.h index 79277f1643f..68d6198cbbb 100644 --- a/Mesh_3/include/CGAL/make_mesh_3.h +++ b/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 > Initial_points_vector; - typedef typename Initial_points_vector::iterator Ipv_iterator; + typedef typename std::pair PI; + typedef std::vector Initial_points_vector; typedef typename C3T3::Vertex_handle Vertex_handle; + typedef CGAL::Mesh_3::Triangulation_helpers 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); } } } diff --git a/Mesh_3/test/Mesh_3/CMakeLists.txt b/Mesh_3/test/Mesh_3/CMakeLists.txt index 2a8e19e1c38..37e41ec079a 100644 --- a/Mesh_3/test/Mesh_3/CMakeLists.txt +++ b/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 ) diff --git a/Mesh_3/test/Mesh_3/test_meshing_with_one_step_with_features.cpp b/Mesh_3/test/Mesh_3/test_meshing_with_one_step_with_features.cpp new file mode 100644 index 00000000000..c3b01f084e7 --- /dev/null +++ b/Mesh_3/test/Mesh_3/test_meshing_with_one_step_with_features.cpp @@ -0,0 +1,86 @@ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Surface_mesh Polyhedron; +typedef CGAL::Polyhedral_mesh_domain_with_features_3 Mesh_domain; + +typedef CGAL::Sequential_tag Concurrency_tag; + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; + +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +typedef CGAL::Mesh_3::Mesher_3 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 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, 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; +} diff --git a/SMDS_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h b/SMDS_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h index 643de3ff7da..34d8b73f7c2 100644 --- a/SMDS_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h +++ b/SMDS_3/include/CGAL/Mesh_complex_3_in_triangulation_3.h @@ -1423,28 +1423,32 @@ public: } } - 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(); }