diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h index fd9981c1426..cf73ab88891 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/smooth_vertices.h @@ -91,6 +91,13 @@ private: std::vector m_free_vertices{}; bool m_flip_smooth_steps{false}; + struct Move + { + Vector_3 move; + int neighbors; + FT mass; + }; + public: Tetrahedral_remeshing_smoother(const SizingFunction& sizing, const CellSelector& cell_selector, @@ -627,8 +634,8 @@ private: const C3t3& c3t3, const bool boundary_edge = false) const { - const auto mwi = midpoint_with_info(e, boundary_edge, c3t3); - const FT s = sizing_at_midpoint(e, mwi.dim, mwi.index, m_sizing, c3t3, m_cell_selector); + const auto [pt, dim, index] = midpoint_with_info(e, boundary_edge, c3t3); + const FT s = sizing_at_midpoint(e, pt, dim, index, m_sizing, c3t3, m_cell_selector); const FT density = 1. / s; //density = 1 / size^(dimension) //edge dimension is 1, so density = 1 / size //to have mass = length * density with no dimension @@ -655,9 +662,8 @@ private: auto& tr = c3t3.triangulation(); const std::size_t nbv = tr.number_of_vertices(); - std::vector moves(nbv, CGAL::NULL_VECTOR); - std::vector neighbors(nbv, 0); - std::vector masses(nbv, 0.); + const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; + std::vector moves(nbv, default_move); //collect neighbors for (const Edge& e : c3t3.edges_in_complex()) @@ -683,33 +689,35 @@ private: if (vh0_moving) { - moves[i0] += density * Vector_3(p0, p1); - neighbors[i0]++; - masses[i0] += density; + moves[i0].move += density * Vector_3(p0, p1); + moves[i0].mass += density; + ++moves[i0].neighbors; } if (vh1_moving) { - moves[i1] += density * Vector_3(p1, p0); - neighbors[i1]++; - masses[i1] += density; + moves[i1].move += density * Vector_3(p1, p0); + moves[i1].mass += density; + ++moves[i1].neighbors; } } - // iterate over map of - for (auto [v, vid] : m_vertex_id) + // iterate over vertices and move + for(Vertex_handle v : tr.finite_vertex_handles()) { + const std::size_t vid = vertex_id(v); + if (!is_free(vid) || !is_on_feature(v)) continue; const Point_3 current_pos = point(v->point()); - const std::size_t nb_neighbors = neighbors[vid]; + const std::size_t nb_neighbors = moves[vid].neighbors; if(nb_neighbors == 0) continue; - CGAL_assertion(masses[vid] > 0); + CGAL_assertion(moves[vid].mass > 0); const Vector_3 move = (nb_neighbors > 0) - ? moves[vid] / masses[vid] + ? moves[vid].move / moves[vid].mass : CGAL::NULL_VECTOR; const Point_3 smoothed_position = current_pos + move; @@ -771,9 +779,8 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3, auto& tr = c3t3.triangulation(); const std::size_t nbv = tr.number_of_vertices(); - std::vector moves(nbv, CGAL::NULL_VECTOR); - std::vector neighbors(nbv, 0); - std::vector masses(nbv, 0.); + const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; + std::vector moves(nbv, default_move); for (const Edge& e : tr.finite_edges()) { @@ -797,26 +804,28 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3, if (vh0_moving) { - moves[i0] += density * Vector_3(p0, p1); - neighbors[i0]++; - masses[i0] += density; + moves[i0].move += density * Vector_3(p0, p1); + moves[i0].mass += density; + ++moves[i0].neighbors; } if (vh1_moving) { - moves[i1] += density * Vector_3(p1, p0); - neighbors[i1]++; - masses[i1] += density; + moves[i1].move += density * Vector_3(p1, p0); + moves[i1].mass += density; + ++moves[i1].neighbors; } } } - // iterate over map of - for (auto [v, vid] : m_vertex_id) + // iterate over vertices and move + for(Vertex_handle v : tr.finite_vertex_handles()) { + const std::size_t vid = vertex_id(v); + if (!is_free(vid) || v->in_dimension() != 2) continue; - const std::size_t nb_neighbors = neighbors[vid]; + const std::size_t nb_neighbors = moves[vid].neighbors; const Point_3 current_pos = point(v->point()); const auto& incident_surface_patches = vertices_surface_indices.at(v); @@ -830,7 +839,7 @@ std::size_t smooth_vertices_on_surfaces(C3t3& c3t3, if (nb_neighbors > 1) { - const Vector_3 move = moves[vid] / masses[vid]; + const Vector_3 move = moves[vid].move / moves[vid].mass; const Point_3 smoothed_position = point(v->point()) + move; #ifdef CGAL_TET_REMESHING_SMOOTHING_WITH_MLS @@ -969,9 +978,9 @@ std::size_t smooth_internal_vertices(C3t3& c3t3, auto& tr = c3t3.triangulation(); const std::size_t nbv = tr.number_of_vertices(); - std::vector moves(nbv, CGAL::NULL_VECTOR); - std::vector neighbors(nbv, 0);/*for dim 3 vertices, start counting directly from 0*/ - std::vector masses(nbv, 0.); + const Move default_move{CGAL::NULL_VECTOR, 0 /*neighbors*/, 0. /*mass*/}; + std::vector moves(nbv, default_move); + /*for dim 3 vertices, start counting neighbors directly from 0*/ for (const Edge& e : tr.finite_edges()) { @@ -979,8 +988,7 @@ std::size_t smooth_internal_vertices(C3t3& c3t3, continue; else { - const Vertex_handle vh0 = e.first->vertex(e.second); - const Vertex_handle vh1 = e.first->vertex(e.third); + const auto [vh0, vh1] = make_vertex_pair(e); const std::size_t& i0 = vertex_id(vh0); const std::size_t& i1 = vertex_id(vh1); @@ -997,31 +1005,32 @@ std::size_t smooth_internal_vertices(C3t3& c3t3, if (vh0_moving) { - moves[i0] += density * Vector_3(p0, p1); - neighbors[i0]++; - masses[i0] += density; + moves[i0].move += density * Vector_3(p0, p1); + moves[i0].mass += density; + ++moves[i0].neighbors; } if (vh1_moving) { - moves[i1] += density * Vector_3(p1, p0); - neighbors[i1]++; - masses[i1] += density; + moves[i1].move += density * Vector_3(p1, p0); + moves[i1].mass += density; + ++moves[i1].neighbors; } } } - // iterate over map of - for (auto [v, vid] : m_vertex_id) + // iterate over vertices and move + for(Vertex_handle v : tr.finite_vertex_handles()) { + const std::size_t vid = vertex_id(v); if (!is_free(vid)) continue; - if (c3t3.in_dimension(v) == 3 && neighbors[vid] > 1) + if (c3t3.in_dimension(v) == 3 && moves[vid].neighbors > 1) { #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG os_vol << "2 " << point(v->point()); #endif - const Vector_3 move = moves[vid] / masses[vid];// static_cast(neighbors[vid]); + const Vector_3 move = moves[vid].move / moves[vid].mass;// static_cast(neighbors[vid]); Point_3 new_pos = point(v->point()) + move; if (check_inversion_and_move(v, new_pos, inc_cells[vid], tr, total_move)){ nb_done_3d++; 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 c667cfe9d12..2c3cfe46cd4 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 @@ -1439,6 +1439,7 @@ auto sizing_at_vertex(const Vertex_handle v, template auto sizing_at_midpoint(const typename C3t3::Edge& e, + const typename C3t3::Triangulation::Geom_traits::Point_3& m, // edge midpoint const int dim, const typename C3t3::Index& index, const Sizing& sizing, @@ -1446,17 +1447,13 @@ auto sizing_at_midpoint(const typename C3t3::Edge& e, const Cell_selector& cell_selector) { using FT = typename C3t3::Triangulation::Geom_traits::FT; - using Point_3 = typename C3t3::Triangulation::Geom_traits::Point_3; - auto cp = c3t3.triangulation().geom_traits().construct_point_3_object(); - const Point_3 m = CGAL::midpoint(cp(e.first->vertex(e.second)->point()), - cp(e.first->vertex(e.third)->point())); + const FT size = sizing(m, dim, index); if (dim < 3 && size == 0) { - const auto u = e.first->vertex(e.second); - const auto v = e.first->vertex(e.third); + const auto [u, v] = make_vertex_pair(e); const FT size_at_u = sizing(cp(u->point()), u->in_dimension(), u->index()); const FT size_at_v = sizing(cp(v->point()), v->in_dimension(), v->index()); @@ -1572,7 +1569,6 @@ auto midpoint_with_info(const typename C3t3::Edge& e, const C3t3& c3t3) { using Tr = typename C3t3::Triangulation; - using Vertex_handle = typename Tr::Vertex_handle; using Gt = typename Tr::Geom_traits; using Point_3 = typename Gt::Point_3; using Index = typename C3t3::Index; @@ -1584,9 +1580,7 @@ auto midpoint_with_info(const typename C3t3::Edge& e, Index index; }; - const auto vs = c3t3.triangulation().vertices(e); - const Vertex_handle u = vs[0]; - const Vertex_handle v = vs[1]; + const auto [u, v] = make_vertex_pair(e); const auto& gt = c3t3.triangulation().geom_traits(); auto cp = gt.construct_point_3_object(); diff --git a/Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt b/Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt index 63bc0c06023..3747ae19125 100644 --- a/Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt +++ b/Tetrahedral_remeshing/test/Tetrahedral_remeshing/CMakeLists.txt @@ -16,6 +16,7 @@ create_single_source_cgal_program("test_tetrahedral_remeshing_with_features.cpp" create_single_source_cgal_program("test_tetrahedral_remeshing_of_one_subdomain.cpp") create_single_source_cgal_program("test_tetrahedral_remeshing_io.cpp") create_single_source_cgal_program("test_tetrahedral_remeshing_from_mesh_file.cpp") +create_single_source_cgal_program("test_tetrahedral_remeshing_determinism.cpp") # Test MLS projection add_executable(test_tetrahedral_remeshing_mls diff --git a/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_determinism.cpp b/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_determinism.cpp new file mode 100644 index 00000000000..7bf3bc26a7f --- /dev/null +++ b/Tetrahedral_remeshing/test/Tetrahedral_remeshing/test_tetrahedral_remeshing_determinism.cpp @@ -0,0 +1,104 @@ +//#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE +//#define CGAL_TETRAHEDRAL_REMESHING_DEBUG + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation_3 Remeshing_triangulation; + +template +void generate_input_one_subdomain(const std::size_t nbv, T3& tr) +{ + CGAL::Random rng; + + typedef typename T3::Point Point; + std::vector pts; + while (pts.size() < nbv) + { + const float x = rng.uniform_real(-10.f, 10.f); + const float y = rng.uniform_real(-10.f, 10.f); + const float z = rng.uniform_real(-10.f, 10.f); + + pts.push_back(Point(x, y, z)); + } + tr.insert(pts.begin(), pts.end()); + + for (typename T3::Cell_handle c : tr.finite_cell_handles()) + c->set_subdomain_index(1); + + assert(tr.is_valid(true)); + +#ifdef CGAL_TETRAHEDRAL_REMESHING_GENERATE_INPUT_FILES + std::ofstream out("data/triangulation_one_subdomain.binary.cgal", + std::ios_base::out | std::ios_base::binary); + CGAL::save_binary_triangulation(out, tr); + out.close(); +#endif +} + +int main(int argc, char* argv[]) +{ + std::cout << "CGAL Random seed = " << CGAL::get_default_random().get_seed() << std::endl; + + // Collect options + std::size_t nb_runs = 5; + + Remeshing_triangulation tr; + generate_input_one_subdomain(1000, tr); + + const int target_edge_length = (argc > 1) ? atoi(argv[1]) : 1; + + std::ofstream ofs0("in.mesh"); + CGAL::IO::write_MEDIT(ofs0, tr); + ofs0.close(); + + const Remeshing_triangulation tr_ref = tr; // Keep a reference for comparison + + std::vector output_tr; + output_tr.reserve(nb_runs); + + for(std::size_t i = 0; i < nb_runs; ++i) + { + std::cout << "Run " << i << " of " << nb_runs << std::endl; + + tr = tr_ref; // Reset triangulation to reference + CGAL::tetrahedral_isotropic_remeshing(tr, target_edge_length); + + std::ostringstream oss; + CGAL::IO::write_MEDIT(oss, tr); + output_tr.push_back(oss.str()); + oss.clear(); + + if(i == 0) + continue; // skip first run, it is the reference + else + { + if(0 != output_tr[i-1].compare(output_tr[i])) + { + std::cerr << "******************************************" << std::endl; + std::cerr << "*** Run " << i << " differs from run " << i-1 << std::endl; + assert(false); // This should not happen + } + else + { + std::cout << "******************************************" << std::endl; + std::cout << "*** Run " << i << " is identical to run " << i - 1 << std::endl; + std::cout << "******************************************" << std::endl; + } + } + } + + return EXIT_SUCCESS; +}