diff --git a/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h b/Mesh_3/include/CGAL/Mesh_3/C3T3_helpers.h index 49217739942..ae33317599b 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,45 +3441,47 @@ 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; - typename Facet_vector::iterator fit = facets.begin(), fend = facets.end(); - for ( ; fit != fend; ++fit ) + std::vector triangles; + typename C3T3::Facet ref_facet; + + 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; // 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); + + 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(), + + CGAL::linear_least_squares_fitting_3(triangles.begin(), + triangles.end(), plane, point, - Dimension_tag<0>(), + Dimension_tag<2>(), 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_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); 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_; } diff --git a/Mesh_3/include/CGAL/Mesh_triangulation_3.h b/Mesh_3/include/CGAL/Mesh_triangulation_3.h index 7bc643f4444..8d97e0696fe 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; @@ -77,11 +78,16 @@ 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; } + const 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/Mesh_3/test/Mesh_3/CMakeLists.txt b/Mesh_3/test/Mesh_3/CMakeLists.txt index 90275373611..6dd24aad08b 100644 --- a/Mesh_3/test/Mesh_3/CMakeLists.txt +++ b/Mesh_3/test/Mesh_3/CMakeLists.txt @@ -62,6 +62,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 @@ -93,6 +94,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) @@ -113,6 +115,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; +} 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..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) @@ -471,6 +471,48 @@ public: return cwp(get_closest_point(cp(wp), cp(wq)), cw(wq)); } + Triangle get_closest_triangle(const Bare_point& p, const Triangle& t) const + { + 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(); + + // 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 = -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, 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); + + if(sq_dist == FT(0)) + return rt; + + if(sq_dist < min_sq_dist) { + rt = tt; + min_sq_dist = sq_dist; + } + } + } + } + + 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. diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt index cc0ba60e6e2..f797bcd0657 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/CMakeLists.txt @@ -19,6 +19,16 @@ 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) + +# 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 # ########################################################## @@ -26,7 +36,9 @@ 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 PRIVATE 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;