diff --git a/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp b/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp index 05875d4e6ca..a64a44fa331 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp @@ -12,6 +12,7 @@ #include "C3t3_type.h" #include +#include #include #include @@ -80,6 +81,7 @@ public Q_SLOTS: Scene_c3t3_item* c3t3_item = qobject_cast(scene->item(index)); + const auto& c3t3 = c3t3_item->c3t3(); if (c3t3_item) { @@ -100,18 +102,37 @@ public Q_SLOTS: bool protect = ui.protect_checkbox->isChecked(); bool smooth_edges = ui.smoothEdges_checkBox->isChecked(); + // collect constraints + using Vertex_handle = Tr::Vertex_handle; + using Vertex_pair = std::pair; + using Constraints_set = std::unordered_set>; + using Constraints_pmap = CGAL::Boolean_property_map; + // wait cursor QApplication::setOverrideCursor(Qt::WaitCursor); QElapsedTimer time; time.start(); + Constraints_set constraints; + for (const auto e : c3t3_item->c3t3().triangulation().finite_edges()) + { + if (c3t3_item->c3t3().is_in_complex(e) + || CGAL::Tetrahedral_remeshing::protecting_balls_intersect(e, c3t3)) + { + Vertex_pair evv = CGAL::Tetrahedral_remeshing::make_vertex_pair(e); + constraints.insert(evv); + } + } + CGAL::tetrahedral_isotropic_remeshing( - c3t3_item->c3t3(), - target_length, - CGAL::parameters::remesh_boundaries(!protect) - .number_of_iterations(nb_iter) - .smooth_constrained_edges(smooth_edges)); + c3t3_item->c3t3(), + target_length, + CGAL::parameters::remesh_boundaries(!protect) + .number_of_iterations(nb_iter) + .smooth_constrained_edges(smooth_edges) + .edge_is_constrained_map(Constraints_pmap(constraints)) + ); std::cout << "Remeshing done (" << time.elapsed() << " ms)" << std::endl; diff --git a/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt b/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt index 0500556f074..a2e1c0fb136 100644 --- a/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt +++ b/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/Tetrahedral_remeshing.txt @@ -103,6 +103,12 @@ setting the named parameter `remesh_boundaries` to `false`. \cgalExample{Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp } +It is also possible to define the polyline features as the ones +stored as complex edges in a `Mesh_complex_3_in_triangulation_3` +(e.g. generated by the \ref PkgMesh3 package mesh generation algorithms). + +\cgalExample{Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp } + \subsection ssecEx4 Tetrahedral Remeshing After Mesh Generation diff --git a/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt b/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt index 018dde26767..17d90d1f20d 100644 --- a/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt +++ b/Tetrahedral_remeshing/doc/Tetrahedral_remeshing/examples.txt @@ -4,5 +4,6 @@ \example Tetrahedral_remeshing/tetrahedral_remeshing_of_one_subdomain.cpp \example Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp \example Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp +\example Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp \example Tetrahedral_remeshing/tetrahedral_remeshing_from_mesh.cpp */ diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt index b5c4c1da2ab..94981f0550b 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt @@ -31,9 +31,13 @@ if(TARGET CGAL::Eigen3_support) create_single_source_cgal_program( "mesh_and_remesh_polyhedral_domain_with_features.cpp" ) target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PUBLIC CGAL::Eigen3_support) + create_single_source_cgal_program("mesh_and_remesh_c3t3.cpp") + target_link_libraries(mesh_and_remesh_c3t3 PUBLIC CGAL::Eigen3_support) + if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support) message(STATUS "Found TBB") target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PRIVATE CGAL::TBB_support) + target_link_libraries(mesh_and_remesh_c3t3 PRIVATE CGAL::TBB_support) endif() else() message(STATUS "NOTICE: Some examples require Eigen 3.1 (or greater), and will not be compiled.") diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp new file mode 100644 index 00000000000..b5a10784e92 --- /dev/null +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp @@ -0,0 +1,96 @@ +#include + +#include +#include +#include + +#include +#include + +#include + +#include + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Mesh_polyhedron_3::type Polyhedron; +typedef CGAL::Polyhedral_mesh_domain_with_features_3 Mesh_domain; + +#ifdef CGAL_CONCURRENT_MESH_3 +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Triangulation for Meshing +typedef CGAL::Mesh_triangulation_3::type Tr; +typedef CGAL::Mesh_complex_3_in_triangulation_3< + Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// Triangulation for Remeshing +typedef CGAL::Triangulation_3 Triangulation_3; +using Vertex_handle = Triangulation_3::Vertex_handle; + +using Vertex_pair = std::pair; +using Constraints_set = std::unordered_set>; +using Constraints_pmap = CGAL::Boolean_property_map; + + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +int main(int argc, char* argv[]) +{ + const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/fandisk.off"); + std::ifstream input(fname); + Polyhedron polyhedron; + input >> polyhedron; + if (input.fail() || !CGAL::is_triangle_mesh(polyhedron)) { + std::cerr << "Error: Input invalid " << fname << std::endl; + return EXIT_FAILURE; + } + + // Create domain + Mesh_domain domain(polyhedron); + + // Get sharp features + domain.detect_features(); + + // Mesh criteria + Mesh_criteria criteria(edge_size = 0.025, + facet_angle = 25, facet_size = 0.05, facet_distance = 0.005, + cell_radius_edge_ratio = 3, cell_size = 0.05); + + // Mesh generation + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); + + CGAL::dump_c3t3(c3t3, "out"); + + Constraints_set constraints; + Constraints_pmap constraints_pmap(constraints); + + Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3), + CGAL::parameters::edge_is_constrained_map(constraints_pmap)); + + //note we use the move semantic, with std::move(c3t3), + // to avoid a copy of the triangulation by the function + // `CGAL::convert_to_triangulation_3()` + // After the call to this function, c3t3 is an empty and valid C3t3. + //It is possible to use : CGAL::convert_to_triangulation_3(c3t3), + // Then the triangulation is copied and duplicated, and c3t3 remains as is. + + const double target_edge_length = 0.1;//coarsen the mesh + CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length, + CGAL::parameters::number_of_iterations(3) + .edge_is_constrained_map(constraints_pmap)); + + std::ofstream out("out_remeshed.mesh"); + CGAL::IO::write_MEDIT(out, tr); + out.close(); + + return EXIT_SUCCESS; +} diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp index 032cca7e1d3..4cb187d96f2 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/tetrahedral_remeshing_with_features.cpp @@ -12,59 +12,16 @@ #include "tetrahedral_remeshing_generate_input.h" -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +using K = CGAL::Exact_predicates_inexact_constructions_kernel; -typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3 Remeshing_triangulation; +using Remeshing_triangulation = CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3; +using Point = Remeshing_triangulation::Point; +using Vertex_handle = Remeshing_triangulation::Vertex_handle; -typedef Remeshing_triangulation::Point Point; -typedef Remeshing_triangulation::Vertex_handle Vertex_handle; -typedef Remeshing_triangulation::Cell_handle Cell_handle; -typedef Remeshing_triangulation::Edge Edge; +using Vertex_pair = std::pair; +using Constraints_set = std::unordered_set>; +using Constraints_pmap = CGAL::Boolean_property_map; -class Constrained_edges_property_map -{ -public: - typedef bool value_type; - typedef bool reference; - typedef std::pair key_type; - typedef boost::read_write_property_map_tag category; - -private: - std::unordered_set>* m_set_ptr; - -public: - Constrained_edges_property_map() - : m_set_ptr(nullptr) - {} - Constrained_edges_property_map(std::unordered_set>* set_) - : m_set_ptr(set_) - {} - -public: - friend void put(Constrained_edges_property_map& map, - const key_type& k, - const bool b) - { - assert(map.m_set_ptr != nullptr); - assert(k.first < k.second); - if (b) map.m_set_ptr->insert(k); - else map.m_set_ptr->erase(k); - } - - friend value_type get(const Constrained_edges_property_map& map, - const key_type& k) - { - assert(map.m_set_ptr != nullptr); - assert(k.first < k.second); - return (map.m_set_ptr->count(k) > 0); - } -}; - -void set_subdomain(Remeshing_triangulation& tr, const int index) -{ - for (auto cit : tr.finite_cell_handles()) - cit->set_subdomain_index(index); -} int main(int argc, char* argv[]) { @@ -73,16 +30,14 @@ int main(int argc, char* argv[]) const int nbv = (argc > 3) ? atoi(argv[3]) : 500; Remeshing_triangulation t3; - typedef std::pair Vertex_pair; - std::unordered_set> constraints; + Constraints_set constraints; CGAL::Tetrahedral_remeshing::generate_input_cube(nbv, t3, constraints); make_constraints_from_cube_edges(t3, constraints); CGAL::tetrahedral_isotropic_remeshing(t3, target_edge_length, - CGAL::parameters::edge_is_constrained_map( - Constrained_edges_property_map(&constraints)) - .number_of_iterations(nb_iter)); + CGAL::parameters::edge_is_constrained_map(Constraints_pmap(constraints)) + .number_of_iterations(nb_iter)); return EXIT_SUCCESS; } diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h index be3dfe8f276..f07a0a2d793 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_adaptive_remeshing_impl.h @@ -471,17 +471,8 @@ private: const Curve_index default_curve_id = default_curve_index(); for (const Edge& e : tr().finite_edges()) { - if (m_c3t3.is_in_complex(e)) - { - CGAL_assertion(m_c3t3.in_dimension(e.first->vertex(e.second)) <= 1); - CGAL_assertion(m_c3t3.in_dimension(e.first->vertex(e.third)) <= 1); -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - ++nbe; -#endif - continue; - } - if (get(ecmap, CGAL::Tetrahedral_remeshing::make_vertex_pair(e)) + || m_c3t3.is_in_complex(e) || nb_incident_subdomains(e, m_c3t3) > 2 || nb_incident_surface_patches(e, m_c3t3) > 1) { @@ -539,6 +530,7 @@ private: CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( m_c3t3.triangulation(), "0-c3t3_vertices_after_init_"); CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3); + CGAL::Tetrahedral_remeshing::debug::count_far_points(m_c3t3); #endif } diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h index f566ab1890f..9ea06ce7d4e 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/tetrahedral_remeshing_helpers.h @@ -568,6 +568,8 @@ void set_index(typename C3t3::Vertex_handle v, const C3t3& c3t3) case 0: v->set_index(Mesh_3::internal::get_index(v->index())); break; + case -1://far points from concurrent Mesh_3 + break; default: CGAL_assertion(false); } @@ -589,6 +591,27 @@ bool is_edge_in_complex(const typename C3t3::Vertex_handle& v0, return false; } +template +bool protecting_balls_intersect(const typename C3t3::Edge& e, + const C3t3& c3t3) +{ + const auto vv = c3t3.triangulation().vertices(e); + if( c3t3.in_dimension(vv[0]) > 1 + || c3t3.in_dimension(vv[1]) > 1) + return false; + + const auto& p0 = vv[0]->point(); + const auto& p1 = vv[1]->point(); + if(p0.weight() == 0 || p1.weight() == 0) + return false; + + const auto r0 = CGAL::approximate_sqrt(p0.weight()); + const auto r1 = CGAL::approximate_sqrt(p1.weight()); + const auto d = CGAL::approximate_sqrt(CGAL::squared_distance(p0, p1)); + + return d < r0 + r1; +} + template OutputIterator incident_subdomains(const typename C3t3::Vertex_handle v, const C3t3& c3t3, @@ -1268,6 +1291,18 @@ void check_surface_patch_indices(const C3t3& c3t3) } } +template +void count_far_points(const C3t3& c3t3) +{ + std::size_t count = 0; + for (auto v : c3t3.triangulation().finite_vertex_handles()) + { + if(c3t3.in_dimension(v) == -1) + ++count; + } + std::cout << "Nb far points : " << count << std::endl; +} + template bool are_cell_orientations_valid(const Tr& tr) { @@ -1678,13 +1713,17 @@ void dump_vertices_by_dimension(const Tr& tr, const char* prefix) typedef typename Tr::Vertex_handle Vertex_handle; std::vector< std::vector > vertices_per_dimension(4); + std::size_t nb_far_points = 0; for (typename Tr::Finite_vertices_iterator vit = tr.finite_vertices_begin(); vit != tr.finite_vertices_end(); ++vit) { if (vit->in_dimension() == -1) + { + ++nb_far_points; continue;//far point + } CGAL_assertion(vit->in_dimension() >= 0 && vit->in_dimension() < 4); vertices_per_dimension[vit->in_dimension()].push_back(vit); @@ -1712,6 +1751,7 @@ void dump_vertices_by_dimension(const Tr& tr, const char* prefix) ofs.close(); } + std::cout << "Nb far points : " << nb_far_points << std::endl; } template diff --git a/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h b/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h index 8a8ebe925b3..d11d72973b6 100644 --- a/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h +++ b/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h @@ -310,21 +310,64 @@ void tetrahedral_isotropic_remeshing( * @tparam CurveIndex is the type of the indices for feature curves. * If `c3t3` has been generated using `CGAL::make_mesh_3()`, it must match * `MeshDomainWithFeatures_3::Curve_index`. +* @tparam NamedParameters a sequence of \ref bgl_namedparameters "Named Parameters" * * @param c3t3 the complex containing the triangulation to be remeshed. +* @param np optional sequence of \ref bgl_namedparameters "Named Parameters" +* among the ones listed below +* +* \cgalNamedParamsBegin +* \cgalParamNBegin{edge_is_constrained_map} +* \cgalParamDescription{a property map containing the constrained-or-not status of each edge of +* `c3t3.triangulation()`. +* For each edge `e` for which `c3t3.is_in_complex(e)` returns `true`, +* the constrained status of `e` is set to `true`.} +* \cgalParamType{a class model of `ReadWritePropertyMap` +* with `std::pair` +* as key type and `bool` as value type. It must be default constructible.} +* \cgalParamDefault{a default property map where no edge is constrained} +* \cgalParamNEnd +* +* \cgalNamedParamsEnd */ template + typename CurveIndex, + typename NamedParameters = parameters::Default_named_parameters> CGAL::Triangulation_3 convert_to_triangulation_3( - CGAL::Mesh_complex_3_in_triangulation_3 c3t3) + CGAL::Mesh_complex_3_in_triangulation_3 c3t3, + const NamedParameters& np = parameters::default_values()) { + using parameters::get_parameter; + using parameters::choose_parameter; + using GT = typename Tr::Geom_traits; using TDS = typename Tr::Triangulation_data_structure; + using Vertex_handle = typename Tr::Vertex_handle; + using Edge_vv = std::pair; + using Default_pmap = Constant_property_map; + using ECMap = typename internal_np::Lookup_named_param_def < + internal_np::edge_is_constrained_t, + NamedParameters, + Default_pmap + >::type; + ECMap ecmap = choose_parameter(get_parameter(np, internal_np::edge_is_constrained), + Default_pmap(false)); + + if (!std::is_same_v) + { + for (auto e : c3t3.edges_in_complex()) + { + const Edge_vv evv + = CGAL::Tetrahedral_remeshing::make_vertex_pair(e);//ordered pair + put(ecmap, evv, true); + } + } + CGAL::Triangulation_3 tr; tr.swap(c3t3.triangulation()); return tr;