From 9666b3cb7316517be7c4cea09448e17a8d835718 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Thu, 26 Nov 2020 16:25:55 +0100 Subject: [PATCH 01/13] replace PCA of points by PCA of triangles to make projection more precise and avoid moving a point inside a protecting ball --- Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 2c162e7cc61..6aa627628a2 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -3441,6 +3441,8 @@ get_least_square_surface_plane(const Vertex_handle& v, // Get adjacent surface points std::vector surface_point_vector; + std::vector triangles; + typename Facet_vector::iterator fit = facets.begin(), fend = facets.end(); for ( ; fit != fend; ++fit ) { @@ -3460,6 +3462,8 @@ get_least_square_surface_plane(const Vertex_handle& v, const Bare_point& bp = tr_.get_closest_point(cp(position), cell->get_facet_surface_center(i)); surface_point_vector.push_back(bp); + + triangles.push_back(c3t3_.triangulation().triangle(*fit)); } } @@ -3470,11 +3474,19 @@ 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(surface_point_vector.begin(), - surface_point_vector.end(), +// CGAL::linear_least_squares_fitting_3(surface_point_vector.begin(), +// surface_point_vector.end(), +// plane, +// point, +// Dimension_tag<0>(), +// tr_.geom_traits(), +// Default_diagonalize_traits()); + + CGAL::linear_least_squares_fitting_3(triangles.begin(), + triangles.end(), plane, point, - Dimension_tag<0>(), + Dimension_tag<2>(), tr_.geom_traits(), Default_diagonalize_traits()); From 0ed6aca5a43bcb953dea35db8d62ca65f1e09c27 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 27 Nov 2020 13:04:06 +0100 Subject: [PATCH 02/13] replace get_closest_point() by get_closest_triangle() to use this function with Periodic_3_mesh_3 --- Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h | 38 ++++++++----------- Mesh_3/include/CGAL/Mesh_triangulation_3.h | 6 +++ .../CGAL/Periodic_3_mesh_triangulation_3.h | 36 ++++++++++++++++++ 3 files changed, 58 insertions(+), 22 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 6aa627628a2..3191a31f85b 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -3420,6 +3420,7 @@ get_least_square_surface_plane(const Vertex_handle& v, Bare_point& reference_point, 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 @@ -3440,47 +3441,40 @@ get_least_square_surface_plane(const Vertex_handle& v, const Weighted_point& position = tr_.point(v); // Get adjacent surface points - std::vector surface_point_vector; - std::vector triangles; + std::vector triangles; + typename C3T3::Facet ref_facet; - typename Facet_vector::iterator fit = facets.begin(), fend = facets.end(); - for ( ; fit != fend; ++fit ) + for (typename C3T3::Facet f : facets) { - if ( c3t3_.is_in_complex(*fit) && + if ( c3t3_.is_in_complex(f) && (patch_index == Surface_patch_index() || - c3t3_.surface_patch_index(*fit) == patch_index) ) + c3t3_.surface_patch_index(f) == patch_index) ) { - const Cell_handle cell = fit->first; - const int i = fit->second; + ref_facet = f; + + const Cell_handle cell = f.first; + const int i = f.second; // 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_point(p, q)' that simply returns q + // 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 - const Bare_point& bp = tr_.get_closest_point(cp(position), - cell->get_facet_surface_center(i)); - surface_point_vector.push_back(bp); - triangles.push_back(c3t3_.triangulation().triangle(*fit)); + Triangle t = c3t3_.triangulation().triangle(f); + Triangle ct = tr_.get_closest_triangle(cp(position), t); + triangles.push_back(ct); } } // In some cases point is not a real surface point - if ( surface_point_vector.empty() ) + if ( triangles.empty() ) return boost::none; // Compute least square fitting plane Plane_3 plane; Bare_point point; -// CGAL::linear_least_squares_fitting_3(surface_point_vector.begin(), -// surface_point_vector.end(), -// plane, -// point, -// Dimension_tag<0>(), -// tr_.geom_traits(), -// Default_diagonalize_traits()); CGAL::linear_least_squares_fitting_3(triangles.begin(), triangles.end(), @@ -3490,7 +3484,7 @@ get_least_square_surface_plane(const Vertex_handle& v, tr_.geom_traits(), Default_diagonalize_traits()); - reference_point = surface_point_vector.front(); + reference_point = ref_facet.first->get_facet_surface_center(ref_facet.second); return plane; } diff --git a/Mesh_3/include/CGAL/Mesh_triangulation_3.h b/Mesh_3/include/CGAL/Mesh_triangulation_3.h index 8d5373b1bd7..67c8dad87e0 100644 --- a/Mesh_3/include/CGAL/Mesh_triangulation_3.h +++ b/Mesh_3/include/CGAL/Mesh_triangulation_3.h @@ -59,6 +59,7 @@ public: typedef typename Geom_traits::FT FT; typedef typename Base::Bare_point Bare_point; typedef typename Base::Weighted_point Weighted_point; + typedef typename Base::Triangle Triangle; typedef typename Base::Vertex_handle Vertex_handle; typedef typename Base::Cell_handle Cell_handle; @@ -82,6 +83,11 @@ public: return q; } + Triangle get_closest_triangle(const Bare_point& /*p*/, const Triangle& t) const + { + return t; + } + void set_point(const Vertex_handle v, const Vector& /*move*/, const Weighted_point& new_position) diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 4bef87373f5..68b8b84059c 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -471,6 +471,42 @@ public: return cwp(get_closest_point(cp(wp), cp(wq)), cw(wq)); } + Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const + { + Triangle rt; + const std::array ct = { canonicalize_point(t[0]), + canonicalize_point(t[1]), + canonicalize_point(t[2]) }; + FT min_sq_dist = std::numeric_limits::infinity(); + + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + for (int k = 0; k < 3; ++k) { + + const Triangle tt( + construct_point(std::make_pair(ct[0], Offset(i - 1, j - 1, k - 1))), + construct_point(std::make_pair(ct[1], Offset(i - 1, j - 1, k - 1))), + construct_point(std::make_pair(ct[2], Offset(i - 1, j - 1, k - 1)))); + + for (int v = 0; v < 3; ++v) {//vertices + const Bare_point ttv = tt[v]; + FT sq_dist = geom_traits().compute_squared_distance_3_object()(p, ttv); + + if (sq_dist < min_sq_dist) + { + rt = tt; + min_sq_dist = sq_dist; + } + if(min_sq_dist == 0.) + return rt; + } + } + } + } + + return rt; + } + // Warning: This is a periodic version that computes the smallest possible // distances between p and q, and between p and r FOR ALL POSSIBLE OFFSETS // before comparing these distances. From edafc8c8d57d707cf53c260a08b63d9fd05cfe57 Mon Sep 17 00:00:00 2001 From: Maxime Gimeno Date: Mon, 30 Nov 2020 09:11:50 +0100 Subject: [PATCH 03/13] Add TBB to Tet_remesh example --- .../examples/Tetrahedral_remeshing/CMakeLists.txt | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt index e14159f08e0..71408a3cecb 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt @@ -22,6 +22,8 @@ endif() # Use Eigen for Mesh_3 find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater) include(CGAL_Eigen_support) +find_package(TBB QUIET) +include(CGAL_TBB_support) # Creating entries for all C++ files with "main" routine # ########################################################## @@ -32,3 +34,6 @@ create_single_source_cgal_program( "tetrahedral_remeshing_from_mesh.cpp") 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::Eigen_support) +if(TARGET CGAL::TBB_support) +target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features PRIVATE CGAL::TBB_support) +endif() From 39b9a0fd262cf3c907c34dca8c6586fd3d3fb424 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Mon, 30 Nov 2020 09:19:48 +0100 Subject: [PATCH 04/13] add possibility to mesh with concurrent mesh_3 before remeshing --- .../Tetrahedral_remeshing/CMakeLists.txt | 18 ++++++++++++++++-- ..._remesh_polyhedral_domain_with_features.cpp | 8 +++++++- 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt index e14159f08e0..7c96ee4aeb0 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt @@ -23,12 +23,26 @@ endif() find_package(Eigen3 3.1.0 REQUIRED) #(3.1.0 or greater) include(CGAL_Eigen_support) +# Concurrent Mesh_3 +option(CGAL_ACTIVATE_CONCURRENT_MESH_3 "Activate parallelism in Mesh_3" OFF) +if(CGAL_ACTIVATE_CONCURRENT_MESH_3 OR ENV{CGAL_ACTIVATE_CONCURRENT_MESH_3}) + add_definitions(-DCGAL_CONCURRENT_MESH_3) + find_package(TBB REQUIRED) + include(CGAL_TBB_support) +endif() + # Creating entries for all C++ files with "main" routine # ########################################################## create_single_source_cgal_program( "tetrahedral_remeshing_example.cpp" ) create_single_source_cgal_program( "tetrahedral_remeshing_with_features.cpp") create_single_source_cgal_program( "tetrahedral_remeshing_of_one_subdomain.cpp") create_single_source_cgal_program( "tetrahedral_remeshing_from_mesh.cpp") +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::Eigen_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::Eigen_support) +if(CGAL_ACTIVATE_CONCURRENT_MESH_3 AND TARGET CGAL::TBB_support) + target_link_libraries(mesh_and_remesh_polyhedral_domain_with_features + PUBLIC CGAL::TBB_support) +endif() diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp index d8ae3c79401..bcfe750e891 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_polyhedral_domain_with_features.cpp @@ -14,8 +14,14 @@ 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_triangulation_3::type Tr; typedef CGAL::Mesh_complex_3_in_triangulation_3< Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index> C3t3; From 5ee03985bf9f596c45f26067e1c6e81286990efd Mon Sep 17 00:00:00 2001 From: Maxime Gimeno Date: Mon, 30 Nov 2020 10:09:48 +0100 Subject: [PATCH 05/13] Fix atomic wrong function --- Mesh_3/include/CGAL/Mesh_3/Mesh_surface_cell_base_3.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/Mesh_surface_cell_base_3.h b/Mesh_3/include/CGAL/Mesh_3/Mesh_surface_cell_base_3.h index f4b279a0fdc..90a0e1ec909 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Mesh_surface_cell_base_3.h +++ b/Mesh_3/include/CGAL/Mesh_3/Mesh_surface_cell_base_3.h @@ -97,7 +97,7 @@ public: { CGAL_precondition(facet>=0 && facet<4); char current_bits = bits_; - while (bits_.compare_and_swap(current_bits | (1 << facet), current_bits) != current_bits) + while (bits_.compare_exchange_weak(current_bits, current_bits | (1 << facet)) ) { current_bits = bits_; } @@ -108,7 +108,7 @@ public: { CGAL_precondition(facet>=0 && facet<4); char current_bits = bits_; - while (bits_.compare_and_swap(current_bits & (15 & ~(1 << facet)), current_bits) != current_bits) + while (bits_.compare_exchange_weak(current_bits, current_bits & (15 & ~(1 << facet)))) { current_bits = bits_; } From 0ed99d0182915b49321e064e9a84dd9bbfb9488e Mon Sep 17 00:00:00 2001 From: Maxime Gimeno Date: Mon, 30 Nov 2020 12:54:59 +0100 Subject: [PATCH 06/13] Add a test for old Mesh_cell_base_3 --- Mesh_3/test/Mesh_3/CMakeLists.txt | 3 + Mesh_3/test/Mesh_3/test_mesh_cell_base_3.cpp | 100 +++++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 Mesh_3/test/Mesh_3/test_mesh_cell_base_3.cpp diff --git a/Mesh_3/test/Mesh_3/CMakeLists.txt b/Mesh_3/test/Mesh_3/CMakeLists.txt index d731e12d4e4..1c92b091c2e 100644 --- a/Mesh_3/test/Mesh_3/CMakeLists.txt +++ b/Mesh_3/test/Mesh_3/CMakeLists.txt @@ -60,6 +60,7 @@ if ( CGAL_FOUND ) 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_mesh_cell_base_3.cpp") foreach(target test_boost_has_xxx @@ -91,6 +92,7 @@ if ( CGAL_FOUND ) test_c3t3_extract_subdomains_boundaries test_mesh_3_issue_1554 test_mesh_polyhedral_domain_with_features_deprecated + test_mesh_cell_base_3 test_meshing_with_one_step.cpp) if(TARGET ${target}) target_link_libraries(${target} PUBLIC CGAL::Eigen_support) @@ -111,6 +113,7 @@ if ( CGAL_FOUND ) test_mesh_capsule_var_distance_bound test_mesh_3_issue_1554 test_mesh_polyhedral_domain_with_features_deprecated + test_mesh_cell_base_3 ) if(TARGET ${target}) target_link_libraries(${target} PUBLIC CGAL::TBB_support) diff --git a/Mesh_3/test/Mesh_3/test_mesh_cell_base_3.cpp b/Mesh_3/test/Mesh_3/test_mesh_cell_base_3.cpp new file mode 100644 index 00000000000..0e433bb0922 --- /dev/null +++ b/Mesh_3/test/Mesh_3/test_mesh_cell_base_3.cpp @@ -0,0 +1,100 @@ +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +int main (int argc, char** argv){ + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + typedef CGAL::Polyhedral_mesh_domain_with_features_3 Polyhedral_mesh_domain; + + // Traits classes + typedef CGAL::Kernel_traits::Kernel Robust_intersections_traits; + typedef CGAL::details::Mesh_geom_traits_generator< + Robust_intersections_traits>::type Robust_K; +#ifdef CGAL_LINKED_WITH_TBB + typedef CGAL::Parallel_tag Concurrency_tag; +#else + typedef CGAL::Sequential_tag Concurrency_tag; +#endif + + // Triangulation + typedef CGAL::Mesh_cell_base_3 Cell_base; + typedef CGAL::Triangulation_cell_base_with_info_3 Cell_base_with_info; + typedef CGAL::Mesh_triangulation_3::type Tr; + + typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + + + // Open file + std::ifstream in (argc > 1 ? argv[1] : "data/elephant.mesh", + std::ios_base::in); + if(!in) { + std::cerr << "Error! Cannot open file " << argv[1] << std::endl; + return 1; + } + C3t3 c3t3; + if(CGAL::build_triangulation_from_file(in, c3t3.triangulation())) + { + for( C3t3::Triangulation::Finite_cells_iterator + cit = c3t3.triangulation().finite_cells_begin(); + cit != c3t3.triangulation().finite_cells_end(); + ++cit) + { + CGAL_assertion(cit->subdomain_index() >= 0); + c3t3.add_to_complex(cit, cit->subdomain_index()); + for(int i=0; i < 4; ++i) + { + if(cit->surface_patch_index(i)>0) + { + c3t3.add_to_complex(cit, i, cit->surface_patch_index(i)); + } + } + } + + //if there is no facet in the complex, we add the border facets. + if(c3t3.number_of_facets_in_complex() == 0) + { + for( C3t3::Triangulation::All_cells_iterator + cit = c3t3.triangulation().all_cells_begin(); + cit != c3t3.triangulation().all_cells_end(); + ++cit) + { + if(c3t3.triangulation().is_infinite(cit)) + { + for(int i=0; i<4; ++i) + { + if(!c3t3.triangulation().is_infinite(cit, i)) + c3t3.add_to_complex(cit, i, 1); + } + } + } + } + } + + CGAL::Polyhedron_3 poly; + CGAL::facets_in_complex_3_to_triangle_mesh(c3t3, poly); + + std::cout << "Graph has " << num_faces(poly) << " faces" << std::endl; + std::cout << "Graph has " << num_vertices(poly) << " vertices" << std::endl; + + std::ofstream out("graph.off"); + out << poly; + + CGAL_assertion(is_valid(poly)); + return EXIT_SUCCESS; +} From cea9cd509b584125936f8e6133859a7512dd7eb4 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Tue, 1 Dec 2020 15:10:07 +0100 Subject: [PATCH 07/13] fix get_closest_triangle with canonicalize_point() used for the 3 pts, the resulting triangle may not be a facet but only a triplet of points --- .../CGAL/Periodic_3_mesh_triangulation_3.h | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 68b8b84059c..c733a277a76 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -473,12 +473,22 @@ public: Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const { - Triangle rt; - const std::array ct = { canonicalize_point(t[0]), - canonicalize_point(t[1]), - canonicalize_point(t[2]) }; - FT min_sq_dist = std::numeric_limits::infinity(); + //canonicalize t + Bare_point min_p = t[0]; + if (t[1] < min_p) { + min_p = t[1]; + } + if (t[2] < min_p) { + min_p = t[2]; + } + Vector_3 move_to_canonical(min_p, canonicalize_point(min_p)); + const std::array ct = { t[0] + move_to_canonical, + t[1] + move_to_canonical, + t[2] + move_to_canonical }; + + FT min_sq_dist = std::numeric_limits::infinity(); + Triangle rt; for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) { for (int k = 0; k < 3; ++k) { From 35fa50213af87f7e51393e7b2c6bbea36e2597f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 3 Dec 2020 12:25:37 +0100 Subject: [PATCH 08/13] Avoid copies in Mesh_triangulation_3's trivial functions These exist because of P3M3 --- Mesh_3/include/CGAL/Mesh_triangulation_3.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_triangulation_3.h b/Mesh_3/include/CGAL/Mesh_triangulation_3.h index 67c8dad87e0..34ffb55fbc2 100644 --- a/Mesh_3/include/CGAL/Mesh_triangulation_3.h +++ b/Mesh_3/include/CGAL/Mesh_triangulation_3.h @@ -78,12 +78,12 @@ public: // possibilities). To allow Periodic_3_mesh_3 to use Mesh_3's files, // each mesh triangulation implements its own version. - Bare_point get_closest_point(const Bare_point& /*p*/, const Bare_point& q) const + const Bare_point& get_closest_point(const Bare_point& /*p*/, const Bare_point& q) const { return q; } - Triangle get_closest_triangle(const Bare_point& /*p*/, const Triangle& t) const + const Triangle& get_closest_triangle(const Bare_point& /*p*/, const Triangle& t) const { return t; } From 8bb4b5884875b0e1136837385a6c334abe1bbc74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 3 Dec 2020 12:26:26 +0100 Subject: [PATCH 09/13] Misc minor improvements for speed (avoid some copies and constructions) --- .../CGAL/Periodic_3_mesh_triangulation_3.h | 48 ++++++++++--------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index c733a277a76..6f1f40ca40a 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -473,42 +473,46 @@ public: Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const { - //canonicalize t - Bare_point min_p = t[0]; - if (t[1] < min_p) { - min_p = t[1]; - } - if (t[2] < min_p) { - min_p = t[2]; - } + typename Geom_traits::Less_xyz_3 less = geom_traits().less_xyz_3_object(); + typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_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(); - Vector_3 move_to_canonical(min_p, canonicalize_point(min_p)); - const std::array ct = { t[0] + move_to_canonical, - t[1] + move_to_canonical, - t[2] + move_to_canonical }; + //canonicalize t + std::size_t min_p_id = 0; + if(less(t[1], t[min_p_id])) + min_p_id = 1; + if(less(t[2], t[min_p_id])) + min_p_id = 2; + + Bare_point canon_min_p = canonicalize_point(t[min_p_id]); + Vector_3 move_to_canonical = cv(t[min_p_id], canon_min_p); + const std::array ct = { (min_p_id == 0) ? canon_min_p : tr(t[0], move_to_canonical), + (min_p_id == 1) ? canon_min_p : tr(t[1], move_to_canonical), + (min_p_id == 2) ? canon_min_p : tr(t[2], move_to_canonical) }; FT min_sq_dist = std::numeric_limits::infinity(); Triangle rt; - for (int i = 0; i < 3; ++i) { - for (int j = 0; j < 3; ++j) { - for (int k = 0; k < 3; ++k) { + for(int i = 0; i < 3; ++i) { + for(int j = 0; j < 3; ++j) { + for(int k = 0; k < 3; ++k) { const Triangle tt( construct_point(std::make_pair(ct[0], Offset(i - 1, j - 1, k - 1))), construct_point(std::make_pair(ct[1], Offset(i - 1, j - 1, k - 1))), construct_point(std::make_pair(ct[2], Offset(i - 1, j - 1, k - 1)))); - for (int v = 0; v < 3; ++v) {//vertices - const Bare_point ttv = tt[v]; - FT sq_dist = geom_traits().compute_squared_distance_3_object()(p, ttv); + for(int v = 0; v < 3; ++v) {//vertices + const Bare_point& ttv = tt[v]; + const FT sq_dist = csd(p, ttv); - if (sq_dist < min_sq_dist) - { + if(sq_dist == FT(0)) + return rt; + + if(sq_dist < min_sq_dist) { rt = tt; min_sq_dist = sq_dist; } - if(min_sq_dist == 0.) - return rt; } } } From a8f368d43756712e77e0eb3897164fd86cca76fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Thu, 3 Dec 2020 12:33:37 +0100 Subject: [PATCH 10/13] Use compute_squared_distance_3(Point, Triangle) --- .../CGAL/Periodic_3_mesh_triangulation_3.h | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index 6f1f40ca40a..f894433fc3a 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -502,17 +502,14 @@ public: construct_point(std::make_pair(ct[1], Offset(i - 1, j - 1, k - 1))), construct_point(std::make_pair(ct[2], Offset(i - 1, j - 1, k - 1)))); - for(int v = 0; v < 3; ++v) {//vertices - const Bare_point& ttv = tt[v]; - const FT sq_dist = csd(p, ttv); + const FT sq_dist = csd(p, tt); - if(sq_dist == FT(0)) - return rt; + if(sq_dist == FT(0)) + return rt; - if(sq_dist < min_sq_dist) { - rt = tt; - min_sq_dist = sq_dist; - } + if(sq_dist < min_sq_dist) { + rt = tt; + min_sq_dist = sq_dist; } } } From 232cf10af8fdb78a2e183ef1446261f434654773 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 4 Dec 2020 10:40:42 +0100 Subject: [PATCH 11/13] remove unused code --- Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 3191a31f85b..0d22f06d9b7 100644 --- a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h +++ b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h @@ -3452,9 +3452,6 @@ get_least_square_surface_plane(const Vertex_handle& v, { ref_facet = f; - const Cell_handle cell = f.first; - const int i = f.second; - // 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 From 2d73ef361eefc0744361e5ef32cf6baad54abc73 Mon Sep 17 00:00:00 2001 From: Jane Tournois Date: Fri, 4 Dec 2020 17:18:04 +0100 Subject: [PATCH 12/13] remove wrong assertions when v has dimension 3, its incident cells can be either all inside or all outside the complex --- Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h | 3 --- 1 file changed, 3 deletions(-) diff --git a/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h b/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h index 387f1e9c1f1..d4a638c78c6 100644 --- a/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h +++ b/Mesh_3/include/CGAL/Mesh_3/Lloyd_move.h @@ -551,7 +551,6 @@ private: Cell_circulator current_cell = tr.incident_cells(edge); Cell_circulator done = current_cell; - CGAL_assertion(c3t3.is_in_complex(current_cell)); // a & b are fixed points const Weighted_point& wa = tr.point(v); @@ -561,7 +560,6 @@ private: const Weighted_point& a_b = tr.point(current_cell, current_cell->index(v)); Vector_3 ba = Vector_3(cp(a_b), b); ++current_cell; - CGAL_assertion(c3t3.is_in_complex(current_cell)); CGAL_assertion(current_cell != done); // c & d are moving points @@ -573,7 +571,6 @@ private: while ( current_cell != done ) { - CGAL_assertion(c3t3.is_in_complex(current_cell)); Bare_point d = tr.dual(current_cell); const Weighted_point& a_d = tr.point(current_cell, current_cell->index(v)); Vector_3 da = Vector_3(cp(a_d), d); From 80806019704e93ac6dda12ba4b03672bc3a35a98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mael=20Rouxel-Labb=C3=A9?= Date: Sat, 5 Dec 2020 17:58:52 +0100 Subject: [PATCH 13/13] Improve P3M3's get_closest_triangle(): any point can be used to canonicalize --- .../CGAL/Periodic_3_mesh_triangulation_3.h | 41 ++++++++----------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h index f894433fc3a..72f731760ec 100644 --- a/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h +++ b/Periodic_3_mesh_3/include/CGAL/Periodic_3_mesh_triangulation_3.h @@ -444,10 +444,10 @@ public: const Bare_point cq = canonicalize_point(q); FT min_sq_dist = std::numeric_limits::infinity(); - for(int i = 0; i < 3; ++i) { - for(int j = 0; j < 3; ++j) { - for(int k = 0; k < 3; ++k) { - const Bare_point tcq = construct_point(std::make_pair(cq, Offset(i-1, j-1, k-1))); + 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) @@ -473,34 +473,29 @@ public: Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const { - typename Geom_traits::Less_xyz_3 less = geom_traits().less_xyz_3_object(); typename Geom_traits::Construct_vector_3 cv = geom_traits().construct_vector_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(); - //canonicalize t - std::size_t min_p_id = 0; - if(less(t[1], t[min_p_id])) - min_p_id = 1; - if(less(t[2], t[min_p_id])) - min_p_id = 2; - - Bare_point canon_min_p = canonicalize_point(t[min_p_id]); - Vector_3 move_to_canonical = cv(t[min_p_id], canon_min_p); - const std::array ct = { (min_p_id == 0) ? canon_min_p : tr(t[0], move_to_canonical), - (min_p_id == 1) ? canon_min_p : tr(t[1], move_to_canonical), - (min_p_id == 2) ? canon_min_p : tr(t[2], move_to_canonical) }; + // 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 ct = { canon_p0, + tr(t[1], move_to_canonical), + tr(t[2], move_to_canonical) }; FT min_sq_dist = std::numeric_limits::infinity(); Triangle rt; - for(int i = 0; i < 3; ++i) { - for(int j = 0; j < 3; ++j) { - for(int k = 0; k < 3; ++k) { + for(int i = -1; i < 2; ++i) { + for(int j = -1; j < 2; ++j) { + for(int k = -1; k < 2; ++k) { const Triangle tt( - construct_point(std::make_pair(ct[0], Offset(i - 1, j - 1, k - 1))), - construct_point(std::make_pair(ct[1], Offset(i - 1, j - 1, k - 1))), - construct_point(std::make_pair(ct[2], Offset(i - 1, j - 1, k - 1)))); + 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);