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 a64a44fa331e..0a75dbae932f 100644 --- a/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp +++ b/Polyhedron/demo/Polyhedron/Plugins/Tetrahedral_remeshing/Tetrahedral_remeshing_plugin.cpp @@ -3,6 +3,7 @@ //#define CGAL_TETRAHEDRAL_REMESHING_DEBUG //#define CGAL_TETRAHEDRAL_REMESHING_VERBOSE_PROGRESS //#define CGAL_TETRAHEDRAL_REMESHING_PROFILE +//#define CGAL_FLIP_ON_SURFACE_DISABLE_44_FLIP #include @@ -120,7 +121,7 @@ public Q_SLOTS: 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); + Vertex_pair evv = CGAL::Tetrahedral_remeshing::make_vertex_pair(e); constraints.insert(evv); } } diff --git a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp index b5a10784e929..6416f6fd3341 100644 --- a/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp +++ b/Tetrahedral_remeshing/examples/Tetrahedral_remeshing/mesh_and_remesh_c3t3.cpp @@ -12,27 +12,27 @@ #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; +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Polyhedron = CGAL::Mesh_polyhedron_3::type; +using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3; #ifdef CGAL_CONCURRENT_MESH_3 -typedef CGAL::Parallel_tag Concurrency_tag; +using Concurrency_tag = CGAL::Parallel_tag; #else -typedef CGAL::Sequential_tag Concurrency_tag; +using Concurrency_tag = CGAL::Sequential_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; +using Tr = CGAL::Mesh_triangulation_3::type; +using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3< + Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>; // Criteria -typedef CGAL::Mesh_criteria_3 Mesh_criteria; +using Mesh_criteria = CGAL::Mesh_criteria_3; // Triangulation for Remeshing -typedef CGAL::Triangulation_3 Triangulation_3; +using Triangulation_3 = CGAL::Triangulation_3; using Vertex_handle = Triangulation_3::Vertex_handle; using Vertex_pair = std::pair; @@ -68,8 +68,7 @@ int main(int argc, char* argv[]) // Mesh generation C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); - CGAL::dump_c3t3(c3t3, "out"); - + // Property map of constraints Constraints_set constraints; Constraints_pmap constraints_pmap(constraints); 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 5880052b43c7..2014827dd51b 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 @@ -4,42 +4,52 @@ #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; +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using Polyhedron = CGAL::Surface_mesh; +using Mesh_domain = CGAL::Polyhedral_mesh_domain_with_features_3; #ifdef CGAL_CONCURRENT_MESH_3 -typedef CGAL::Parallel_tag Concurrency_tag; +using Concurrency_tag = CGAL::Parallel_tag; #else -typedef CGAL::Sequential_tag Concurrency_tag; +using Concurrency_tag = CGAL::Sequential_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; +// Triangulation +using Tr = CGAL::Mesh_triangulation_3::type; +using C3t3 = CGAL::Mesh_complex_3_in_triangulation_3< + Tr, Mesh_domain::Corner_index, Mesh_domain::Curve_index>; // Criteria -typedef CGAL::Mesh_criteria_3 Mesh_criteria; +using Mesh_criteria = CGAL::Mesh_criteria_3; // Triangulation for Remeshing -typedef CGAL::Triangulation_3 Triangulation_3; +using Triangulation_3 = CGAL::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"); + const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/anchor.off"); + const int nb_iter = (argc > 2) ? atoi(argv[2]) : 5; + std::ifstream input(fname); Polyhedron polyhedron; + + std::string filename(fname); input >> polyhedron; if (input.fail()) { std::cerr << "Error: Cannot read file " << fname << std::endl; @@ -58,24 +68,29 @@ int main(int argc, char* argv[]) 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); + const double size = 0.072; + Mesh_criteria criteria(edge_size = size, + facet_angle = 25, + facet_size = size, + cell_radius_edge_ratio = 2, + cell_size = size); // Mesh generation - C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria); - - Triangulation_3 tr = CGAL::convert_to_triangulation_3(std::move(c3t3)); - //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)); + C3t3 c3t3 = CGAL::make_mesh_3(domain, criteria, no_perturb(), no_exude()); + + 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)); + + // Remeshing + CGAL::tetrahedral_isotropic_remeshing(tr, size, + CGAL::parameters::number_of_iterations(nb_iter)); + + std::ofstream out("out_remeshed.mesh"); + CGAL::IO::write_MEDIT(out, tr); + out.close(); return EXIT_SUCCESS; } diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h index d2ababd2779e..c36603287bce 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/collapse_short_edges.h @@ -44,6 +44,7 @@ enum Edge_type { FEATURE, BOUNDARY, INSIDE, MIXTE, enum Collapse_type { TO_MIDPOINT, TO_V0, TO_V1, IMPOSSIBLE }; enum Result_type { VALID, V_PROBLEM, C_PROBLEM, E_PROBLEM, + ANGLE_PROBLEM, TOPOLOGICAL_PROBLEM, ORIENTATION_PROBLEM, SHARED_NEIGHBOR_PROBLEM }; template @@ -214,6 +215,15 @@ class CollapseTriangulation } while (++circ != done); +#ifdef PROTECT_ANGLES_FROM_COLLAPSE + Dihedral_angle_cosine curr_max_cos + = max_cos_dihedral_angle(triangulation, cells_to_remove[0]); + for (std::size_t i = 1; i < cells_to_remove.size(); ++i) + { + curr_max_cos = (std::max)(curr_max_cos, + max_cos_dihedral_angle(triangulation, cells_to_remove[i])); + } +#endif vh0->set_point(Point_3(v0_new_pos.x(), v0_new_pos.y(), v0_new_pos.z())); vh1->set_point(Point_3(v0_new_pos.x(), v0_new_pos.y(), v0_new_pos.z())); @@ -269,6 +279,10 @@ class CollapseTriangulation return ORIENTATION_PROBLEM; if (!triangulation.tds().is_valid(cit, true)) return C_PROBLEM; +#ifdef PROTECT_ANGLES_FROM_COLLAPSE + if (curr_max_cos < max_cos_dihedral_angle(triangulation, cit)) + return ANGLE_PROBLEM; +#endif } for (Vertex_handle vit : triangulation.finite_vertex_handles()) @@ -791,6 +805,9 @@ collapse(const typename C3t3::Cell_handle ch, } while (++circ != done); + if(c3t3.is_in_complex(ch->vertex(from), ch->vertex(to))) + c3t3.remove_from_complex(ch->vertex(from), ch->vertex(to)); + bool valid = true; std::vector cells_to_remove; std::unordered_set invalid_cells; @@ -1109,11 +1126,9 @@ typename C3t3::Vertex_handle collapse_edge(typename C3t3::Edge& edge, if (are_edge_lengths_valid(edge, c3t3, new_pos, sqhigh, cell_selector/*, adaptive = false*/) && collapse_preserves_surface_star(edge, c3t3, new_pos, cell_selector)) { - CGAL_assertion_code(typename Tr::Cell_handle dc); - CGAL_assertion_code(int di); - CGAL_assertion_code(int dj); - CGAL_assertion(c3t3.triangulation().is_edge(edge.first->vertex(edge.second), - edge.first->vertex(edge.third), dc, di, dj)); + CGAL_assertion(c3t3.triangulation().tds().is_edge( + edge.first->vertex(edge.second), + edge.first->vertex(edge.third))); Vertex_handle v0_init = edge.first->vertex(edge.second); Vertex_handle v1_init = edge.first->vertex(edge.third); @@ -1223,7 +1238,7 @@ void collapse_short_edges(C3T3& c3t3, FT sqlen = sql(tr.segment(e)); if (sqlen < sq_low) - short_edges.insert(short_edge(make_vertex_pair(e), sqlen)); + short_edges.insert(short_edge(make_vertex_pair(e), sqlen)); } #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG @@ -1285,7 +1300,7 @@ void collapse_short_edges(C3T3& c3t3, const FT sqlen = sql(tr.segment(eshort)); if (sqlen < sq_low) - short_edges.insert(short_edge(make_vertex_pair(eshort), sqlen)); + short_edges.insert(short_edge(make_vertex_pair(eshort), sqlen)); } //debug::dump_c3t3(c3t3, "dump_after_collapse"); diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h index 798284fa04ea..a80a883fa17d 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/flip_edges.h @@ -21,7 +21,6 @@ #include #include -#include #include #include @@ -38,12 +37,6 @@ namespace internal enum Flip_Criterion{ MIN_ANGLE_BASED, AVERAGE_ANGLE_BASED, VALENCE_BASED, VALENCE_MIN_DH_BASED }; -//template -//void flip_inside_edges(std::vector&) -//{ -// //TODO -//} - //outer_mirror_facets contains the set of facets of the outer hull //of the set of cells modified by the flip operation, //"seen from" outside @@ -88,13 +81,13 @@ void update_c3t3_facets(C3t3& c3t3, } } -template +template Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge, C3t3& c3t3, const std::vector& vertices_around_edge, const Flip_Criterion& criterion, IncCellsVectorMap& inc_cells, - Cell_selector& cell_selector) + CellSelector& cell_selector) { typedef typename C3t3::Triangulation Tr; typedef typename C3t3::Facet Facet; @@ -317,7 +310,7 @@ Sliver_removal_result flip_3_to_2(typename C3t3::Edge& edge, } ch->vertex(v)->set_cell(ch); - inc_cells[ch->vertex(v)] = std::nullopt; + inc_cells[ch->vertex(v)].clear(); ch->reset_cache_validity(); } } @@ -385,7 +378,6 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, { typedef typename C3t3::Triangulation Tr; typedef typename C3t3::Vertex_handle Vertex_handle; - typedef typename C3t3::Cell_handle Cell_handle; typedef typename C3t3::Facet Facet; typedef typename Tr::Facet_circulator Facet_circulator; typedef typename Tr::Cell_circulator Cell_circulator; @@ -431,9 +423,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, indices(facet_circulator->second, j)); if (curr != vh0 && curr != vh1) { - Cell_handle ch; - int i0, i1; - if (tr.is_edge(curr, vh, ch, i0, i1)) + if (tr.tds().is_edge(curr, vh)) is_edge = true; } } @@ -456,12 +446,12 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, Cell_circulator cell_circulator = tr.incident_cells(edge); Cell_circulator done = cell_circulator; - boost::container::small_vector facets; for (std::size_t i = 0; i < opposite_vertices.size(); ++i) { Vertex_handle vh = opposite_vertices[i]; bool keep = true; + boost::container::small_vector facets; do { //Store it if it do not have vh @@ -482,7 +472,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, Dihedral_angle_cosine max_flip_cos_dh(CGAL::NEGATIVE, 1., 1.); for (const Facet& fi : facets) { - if (!tr.is_infinite(fi.first)) + if (!tr.is_infinite(fi.first) && c3t3.is_in_complex(fi.first)) { if (is_well_oriented(tr, vh, fi.first->vertex(indices(fi.second, 0)), fi.first->vertex(indices(fi.second, 1)), @@ -596,13 +586,9 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, if(tr.is_infinite(vh)) continue; - std::optional>& o_inc_vh = inc_cells[vh]; - if (o_inc_vh == std::nullopt) - { - boost::container::small_vector inc_vec; - tr.incident_cells(vh, std::back_inserter(inc_vec)); - o_inc_vh = inc_vec; - } + boost::container::small_vector& o_inc_vh = inc_cells[vh]; + if (o_inc_vh.empty()) + tr.incident_cells(vh, std::back_inserter(o_inc_vh)); Facet_circulator facet_circulator = curr_fcirc; Facet_circulator facet_done = curr_fcirc; @@ -620,7 +606,7 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, indices(facet_circulator->second, i)); if (curr_vertex != vh0 && curr_vertex != vh1) { - if (is_edge_uv(vh, curr_vertex, o_inc_vh.value())) + if (is_edge_uv(vh, curr_vertex, o_inc_vh)) { is_edge = true; break; @@ -707,13 +693,13 @@ void find_best_flip_to_improve_dh(C3t3& c3t3, template Sliver_removal_result flip_n_to_m(C3t3& c3t3, typename C3t3::Edge& edge, typename C3t3::Vertex_handle vh, IncCellsVectorMap& inc_cells, - Cell_selector& cell_selector, + CellSelector& cell_selector, Visitor& visitor, bool check_validity = false) { @@ -728,19 +714,26 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, Tr& tr = c3t3.triangulation(); - Vertex_handle vh0 = edge.first->vertex(edge.second); - Vertex_handle vh1 = edge.first->vertex(edge.third); + const Vertex_handle vh0 = edge.first->vertex(edge.second); + const Vertex_handle vh1 = edge.first->vertex(edge.third); //This vertex will have its valence augmenting a lot, //TODO take the best one //TODO!!!! Check that the created edges do not exist!!! + boost::container::small_vector facets_in_complex; + Facet_circulator facet_circulator = tr.incident_facets(edge); Facet_circulator done_facet_circulator = facet_circulator; bool look_for_vh_iterator = true; do { + if (c3t3.is_in_complex(*facet_circulator)) + { + facets_in_complex.push_back(*facet_circulator); + } + facet_circulator++; //Get the ids of the opposite vertices @@ -749,6 +742,7 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, if (facet_circulator->first->vertex(indices(facet_circulator->second, i)) == vh) look_for_vh_iterator = false; } + } while (facet_circulator != done_facet_circulator && look_for_vh_iterator); if (look_for_vh_iterator) { @@ -761,13 +755,9 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, facet_circulator++; facet_circulator++; - std::optional>& o_inc_vh = inc_cells[vh]; - if (o_inc_vh == std::nullopt) - { - boost::container::small_vector inc_vec; - tr.incident_cells(vh, std::back_inserter(inc_vec)); - o_inc_vh = inc_vec; - } + boost::container::small_vector& o_inc_vh = inc_cells[vh]; + if (o_inc_vh.empty()) + tr.incident_cells(vh, std::back_inserter(o_inc_vh)); do { @@ -778,7 +768,7 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, indices(facet_circulator->second, i)); if (curr_vertex != vh0 && curr_vertex != vh1) { - if (is_edge_uv(vh, curr_vertex, o_inc_vh.value())) + if (is_edge_uv(vh, curr_vertex, o_inc_vh)) return NOT_FLIPPABLE; } } @@ -866,16 +856,16 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, //} ///*************************************************************/ - //Subdomain index? - typename C3t3::Subdomain_index subdomain = to_remove[0]->subdomain_index(); + //Surface + for (const Facet& f : facets_in_complex) + c3t3.remove_from_complex(f); + + //Subdomain index + typedef typename C3t3::Subdomain_index Subdomain_index; + const Subdomain_index subdomain = to_remove[0]->subdomain_index(); bool selected = get(cell_selector, to_remove[0]); visitor.before_flip(to_remove[0]); -#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG - for (std::size_t i = 1; i < to_remove.size(); ++i) - CGAL_assertion(subdomain == to_remove[i]->subdomain_index()); -#endif - std::vector cells_to_update; //Create new cells @@ -948,9 +938,9 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, } ch->vertex(v)->set_cell(ch); - inc_cells[ch->vertex(v)] = std::nullopt; - ch->reset_cache_validity(); + inc_cells[ch->vertex(v)].clear(); } + ch->reset_cache_validity(); } // Update c3t3 @@ -960,6 +950,7 @@ Sliver_removal_result flip_n_to_m(C3t3& c3t3, for (Cell_handle ch : to_remove) { treat_before_delete(ch, cell_selector, c3t3); + ch->reset_cache_validity(); tr.tds().delete_cell(ch); } @@ -1045,16 +1036,16 @@ Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge, Cell_circulator done = circ; Dihedral_angle_cosine curr_max_cosdh = max_cos_dihedral_angle(tr, circ++); - while (circ != done) + do { - curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, circ++)); - } + curr_max_cosdh = (std::max)(curr_max_cosdh, max_cos_dihedral_angle(tr, circ)); + } while (++circ != done); + if (boundary_vertices.size() == 2) find_best_flip_to_improve_dh(c3t3, edge, boundary_vertices[0], boundary_vertices[1], candidates, curr_max_cosdh); else - find_best_flip_to_improve_dh(c3t3, edge, candidates, curr_max_cosdh, - inc_cells); + find_best_flip_to_improve_dh(c3t3, edge, candidates, curr_max_cosdh, inc_cells); bool flip_performed = false; while (!candidates.empty() && !flip_performed) @@ -1062,10 +1053,10 @@ Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge, CosAngle_and_vertex curr_cost_vpair = candidates.top(); candidates.pop(); -// std::cout << "\tcurrent cos = " << curr_max_cosdh -// << "\t angle = " << std::acos(curr_max_cosdh) * 180./CGAL_PI << std::endl; -// std::cout << "\tcandidate cos = " << curr_cost_vpair.first -// << "\t angle = " << std::acos(curr_cost_vpair.first) * 180./CGAL_PI << std::endl; +// std::cout << "\tcurrent cos = " << curr_max_cosdh.value() +// << "\t angle = " << std::acos(curr_max_cosdh.value()) * 180./CGAL_PI << std::endl; +// std::cout << "\tcandidate cos = " << curr_cost_vpair.first.value() +// << "\t angle = " << std::acos(curr_cost_vpair.first.value()) * 180./CGAL_PI << std::endl; // std::cout << std::endl; if (curr_max_cosdh <= curr_cost_vpair.first) @@ -1082,12 +1073,12 @@ Sliver_removal_result flip_n_to_m(typename C3t3::Edge& edge, return result; } -template +template Sliver_removal_result find_best_flip(typename C3t3::Edge& edge, C3t3& c3t3, const Flip_Criterion& criterion, IncCellsVectorMap& inc_cells, - Cell_selector& cell_selector, + CellSelector& cell_selector, Visitor& visitor) { typedef typename C3t3::Triangulation Tr; @@ -1171,38 +1162,31 @@ Sliver_removal_result find_best_flip(typename C3t3::Edge& edge, } -template +template std::size_t flip_all_edges(const std::vector& edges, C3t3& c3t3, + IncidentCellsVectorMap& inc_cells, const Flip_Criterion& criterion, - Cell_selector& cell_selector, + CellSelector& cell_selector, Visitor& visitor) { typedef typename C3t3::Triangulation Tr; typedef typename Tr::Cell_handle Cell_handle; - typedef typename Tr::Vertex_handle Vertex_handle; typedef typename Tr::Edge Edge; Tr& tr = c3t3.triangulation(); - std::unordered_map > > inc_cells; - std::size_t count = 0; for (const VertexPair& vp : edges) { - std::optional>& - o_inc_vh = inc_cells[vp.first]; - if (o_inc_vh == std::nullopt) - { - boost::container::small_vector inc_vec; - tr.incident_cells(vp.first, std::back_inserter(inc_vec)); - o_inc_vh = inc_vec; - } + boost::container::small_vector& o_inc_vh = inc_cells[vp.first]; + if (o_inc_vh.empty()) + tr.incident_cells(vp.first, std::back_inserter(o_inc_vh)); Cell_handle ch; int i0, i1; - if (is_edge_uv(vp.first, vp.second, o_inc_vh.value(), ch, i0, i1)) + if (is_edge_uv(vp.first, vp.second, o_inc_vh, ch, i0, i1)) { Edge edge(ch, i0, i1); @@ -1225,10 +1209,684 @@ std::size_t flip_all_edges(const std::vector& edges, } } - for(Cell_handle c : c3t3.triangulation().finite_cell_handles()) + return count; +} + +template +void collect_subdomains_on_boundary(const C3t3& c3t3, + boost::unordered_map >& vertices_subdomain_indices) +{ + for (auto c : c3t3.triangulation().all_cell_handles()) + { + for (auto v : c3t3.triangulation().vertices(c)) + { + const int dim = v->in_dimension(); + if(dim >= 0 && dim < 3) + vertices_subdomain_indices[v].insert(c->subdomain_index()); + } + } +} + +template +void collectBoundaryEdgesAndComputeVerticesValences( + const C3T3& c3t3, + const CellSelector& cell_selector, + std::vector& boundary_edges, + boost::unordered_map >& + boundary_vertices_valences, + boost::unordered_map >& + vertices_subdomain_indices) +{ + typedef typename C3T3::Surface_patch_index Surface_patch_index; + typedef typename C3T3::Vertex_handle Vertex_handle; + typedef typename C3T3::Edge Edge; + typedef typename C3T3::Triangulation::Facet_circulator Facet_circulator; + + const typename C3T3::Triangulation& tr = c3t3.triangulation(); + + boundary_edges.clear(); + boundary_vertices_valences.clear(); + + for (const Edge& e : tr.finite_edges()) + { + if (is_boundary(c3t3, e, cell_selector)) + boundary_edges.push_back(e); + } + +#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG + CGAL::Tetrahedral_remeshing::debug::dump_edges(boundary_edges, + "boundary_edges.polylines.txt"); +#endif + + // collect incident subdomain indices at vertices + collect_subdomains_on_boundary(c3t3, vertices_subdomain_indices); + + for (const Edge& e : boundary_edges) + { + const Vertex_handle v0 = e.first->vertex(e.second); + const Vertex_handle v1 = e.first->vertex(e.third); + + //In case of feature edge + if (vertices_subdomain_indices[v0].size() > 2 + && vertices_subdomain_indices[v1].size() > 2) + { + Facet_circulator facet_circulator = tr.incident_facets(e); + Facet_circulator done(facet_circulator); + do + { + if (c3t3.is_in_complex(*facet_circulator)) + { + Surface_patch_index surfi = c3t3.surface_patch_index(*facet_circulator); + boundary_vertices_valences[v0][surfi]++; + boundary_vertices_valences[v1][surfi]++; + } + } while (++facet_circulator != done); + } + else if (vertices_subdomain_indices[v0].size() == 2 + || vertices_subdomain_indices[v1].size() == 2) + { + Facet_circulator facet_circulator = tr.incident_facets(e); + Facet_circulator done(facet_circulator); + Surface_patch_index first_patch = Surface_patch_index(); + do + { + if (c3t3.is_in_complex(*facet_circulator)) + { + Surface_patch_index surfi = c3t3.surface_patch_index(*facet_circulator); + if (first_patch == Surface_patch_index()) + first_patch = surfi; + else if (first_patch == surfi) + continue; + + boundary_vertices_valences[v0][surfi]++; + boundary_vertices_valences[v1][surfi]++; + } + } while (++facet_circulator != done); + } + } +} + +template +Sliver_removal_result flip_n_to_m_on_surface(typename C3T3::Edge& edge, + C3T3& c3t3, + typename C3T3::Vertex_handle v0i,//v0 of new edge that will replace edge + typename C3T3::Vertex_handle v1i,//v1 of new edge that will replace edge + const IncCellsVector& cells_around_edge, + Flip_Criterion /*flip_criterion*/, + Visitor& /*visitor*/) +{ + typedef typename C3T3::Vertex_handle Vertex_handle; + typedef typename C3T3::Cell_handle Cell_handle; + + typename C3T3::Triangulation& tr = c3t3.triangulation(); + + const Vertex_handle u = edge.first->vertex(edge.second); + const Vertex_handle v = edge.first->vertex(edge.third); + + typedef std::pair IndInCell; + std::map indices; + for (Cell_handle c : cells_around_edge) + { + indices[c] = std::make_pair(c->index(u), c->index(v)); + } + + for (Cell_handle c : cells_around_edge) + { + int i = indices[c].first; + int j = indices[c].second; + c->set_vertex(i, v0i); + c->set_vertex(j, v1i); + + if (!is_well_oriented(tr, c)) + { + c->set_vertex(j, v0i); + c->set_vertex(i, v1i); + if (!is_well_oriented(tr, c)) + { + //rollback all changes + for (Cell_handle cc : cells_around_edge) + { + const int ii = indices[cc].first; + if (cc->vertex(ii) != u) + { + cc->set_vertex(ii, u); + cc->set_vertex(indices[cc].second, v); + } + } + return NOT_FLIPPABLE; + } + } + } + + for (Cell_handle c : cells_around_edge) c->reset_cache_validity(); - return count; + return VALID_FLIP; +} + +//v0i and v1i are the vertices opposite to `edge` +//on facets of the surface +template +Sliver_removal_result flip_on_surface(C3T3& c3t3, + typename C3T3::Edge& edge, + typename C3T3::Vertex_handle v0i,//v0 of new edge that will replace edge + typename C3T3::Vertex_handle v1i,//v1 of new edge that will replace edge + IncCellsVectorMap& inc_cells, + Flip_Criterion flip_criterion, + Visitor& visitor) +{ + typedef typename C3T3::Triangulation Tr; + typedef typename Tr::Cell_handle Cell_handle; + typedef typename Tr::Vertex_handle Vertex_handle; + typedef typename Tr::Cell_circulator Cell_circulator; + + Tr& tr = c3t3.triangulation(); + Cell_circulator circ = tr.incident_cells(edge); + Cell_circulator done(circ); + + std::vector cells_around_edge; + do + { + cells_around_edge.push_back(circ); + } while (++circ != done); + + if (cells_around_edge.size() != 4) + { + if (cells_around_edge.size() > 4){ +////// if (flip_criterion == VALENCE_BASED){ +////// return find_best_n_m_flip(edge, vh0_index, vh1_index); +////// } +////// else { +// std::vector boundary_vertices; +// boundary_vertices.push_back(v0i); +// boundary_vertices.push_back(v1i); +#ifdef CGAL_FLIP_ON_SURFACE_DISABLE_NM_FLIP + return NOT_FLIPPABLE; +#else + return flip_n_to_m_on_surface(edge, c3t3, v0i, v1i, + cells_around_edge, flip_criterion, + visitor); +#endif +////// } + } + else + return NOT_FLIPPABLE; + } + +#ifdef CGAL_FLIP_ON_SURFACE_DISABLE_44_FLIP + return NOT_FLIPPABLE; +#endif + + inc_cells[edge.first->vertex(edge.second)].clear(); + inc_cells[edge.first->vertex(edge.third)].clear(); + + Cell_handle ch0, ch1, ch2, ch3; + ch0 = cells_around_edge[0]; + ch1 = cells_around_edge[1]; + ch2 = cells_around_edge[2]; + ch3 = cells_around_edge[3]; + + Dihedral_angle_cosine curr_max_cosdh = max_cos_dihedral_angle(tr, ch0); + for (int i = 1; i < 4; ++i) + curr_max_cosdh = (std::max)(curr_max_cosdh, + max_cos_dihedral_angle(tr, cells_around_edge[i])); + + Vertex_handle vh0, vh1, vh2, vh3, vh4, vh5; + + int ivh4_in_ch0 = ch0->index(ch1); + vh4 = ch0->vertex(ivh4_in_ch0); + + int ivh2_in_ch0 = ch0->index(ch3); + vh2 = ch0->vertex(ivh2_in_ch0); + + vh5 = ch1->vertex(ch1->index(ch0)); + vh0 = ch2->vertex(ch2->index(ch1)); + + for (int j = 0; j < 3; j++){ + if (indices(ivh4_in_ch0, j) == ivh2_in_ch0){ + int j1 = (j + 1) % 3; + int j2 = (j + 2) % 3; + vh1 = ch0->vertex(indices(ivh4_in_ch0, j1)); + vh3 = ch0->vertex(indices(ivh4_in_ch0, j2)); + break; + } + } + + bool planar_flip; + if ((vh0 == v0i && vh2 == v1i) || (vh2 == v0i && vh0 == v1i)) + planar_flip = true; + else if ((vh4 == v0i && vh5 == v1i) || (vh5 == v0i && vh4 == v1i)) + planar_flip = false; + else + return NOT_FLIPPABLE; + + typedef typename C3T3::Facet Facet; + typedef typename C3T3::Surface_patch_index Surface_patch_index; + + if (planar_flip) + { +#ifdef CGAL_FLIP_ON_SURFACE_DISABLE_PLANAR_44_FLIP + return NOT_FLIPPABLE; +#endif + Surface_patch_index patch = c3t3.surface_patch_index(ch0, ch0->index(vh4)); + CGAL_assertion(patch != Surface_patch_index()); + CGAL_assertion(c3t3.is_in_complex(ch0, ch0->index(vh4))); + c3t3.remove_from_complex(ch0, ch0->index(vh4)); + CGAL_assertion(c3t3.is_in_complex(ch3, ch3->index(vh4))); + c3t3.remove_from_complex(ch3, ch3->index(vh4)); + + boost::unordered_map opposite_facet_in_complex; + for (Cell_handle chi : cells_around_edge) + { + Facet f1(chi, chi->index(vh1)); + Facet f2(chi, chi->index(vh3)); + + if (c3t3.is_in_complex(f1)) + { + Surface_patch_index spi = c3t3.surface_patch_index(f1); + opposite_facet_in_complex[c3t3.triangulation().mirror_facet(f1)] = spi; + c3t3.remove_from_complex(f1); + } + if (c3t3.is_in_complex(f2)) + { + Surface_patch_index spi = c3t3.surface_patch_index(f2); + opposite_facet_in_complex[c3t3.triangulation().mirror_facet(f2)] = spi; + c3t3.remove_from_complex(f2); + } + } + + Cell_handle n_ch3_vh1 = ch3->neighbor(ch3->index(vh1)); + Cell_handle n_ch0_vh3 = ch0->neighbor(ch0->index(vh3)); + + Cell_handle n_ch2_vh1 = ch2->neighbor(ch2->index(vh1)); + Cell_handle n_ch1_vh3 = ch1->neighbor(ch1->index(vh3)); + + ch3->set_vertex(ch3->index(vh3), vh2); + ch0->set_vertex(ch0->index(vh1), vh0); + ch2->set_vertex(ch2->index(vh3), vh2); + ch1->set_vertex(ch1->index(vh1), vh0); + + Sliver_removal_result db = VALID_FLIP; + if (!is_well_oriented(tr, ch0) + || !is_well_oriented(tr, ch1) + || !is_well_oriented(tr, ch2) + || !is_well_oriented(tr, ch3)) + db = NOT_FLIPPABLE; + else if (curr_max_cosdh < max_cos_dihedral_angle(tr, ch0, false) + || curr_max_cosdh < max_cos_dihedral_angle(tr, ch1, false) + || curr_max_cosdh < max_cos_dihedral_angle(tr, ch2, false) + || curr_max_cosdh < max_cos_dihedral_angle(tr, ch3, false)) + db = NO_BEST_CONFIGURATION; + + if(db != VALID_FLIP) + { + ch3->set_vertex(ch3->index(vh2), vh3); + ch0->set_vertex(ch0->index(vh0), vh1); + ch2->set_vertex(ch2->index(vh2), vh3); + ch1->set_vertex(ch1->index(vh0), vh1); + + c3t3.add_to_complex(ch0, ch0->index(vh4), patch); + c3t3.add_to_complex(ch3, ch3->index(vh4), patch); + + for (Cell_handle chi : cells_around_edge) + { + Facet f1(chi, chi->index(vh1)); + Facet f2(chi, chi->index(vh3)); + + auto it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f1)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f1, it->second); + + it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f2)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f2, it->second); + } + + return db; + } + + //Top cells 2-2 flip + ch3->set_neighbor(ch3->index(vh1), ch0); + ch3->set_neighbor(ch3->index(vh0), n_ch0_vh3); + n_ch0_vh3->set_neighbor(n_ch0_vh3->index(ch0), ch3); + + ch0->set_neighbor(ch0->index(vh3), ch3); + ch0->set_neighbor(ch0->index(vh2), n_ch3_vh1); + n_ch3_vh1->set_neighbor(n_ch3_vh1->index(ch3), ch0); + + //Bottom cells 2-2 flip + ch2->set_neighbor(ch2->index(vh1), ch1); + ch2->set_neighbor(ch2->index(vh0), n_ch1_vh3); + n_ch1_vh3->set_neighbor(n_ch1_vh3->index(ch1), ch2); + + ch1->set_neighbor(ch1->index(vh3), ch2); + ch1->set_neighbor(ch1->index(vh2), n_ch2_vh1); + n_ch2_vh1->set_neighbor(n_ch2_vh1->index(ch2), ch1); + + for (Cell_handle ci : cells_around_edge) + { + for (int j = 0; j < 4; j++) + ci->vertex(j)->set_cell(ci); + } + + c3t3.add_to_complex(ch0, ch0->index(vh4), patch); + c3t3.add_to_complex(ch3, ch3->index(vh4), patch); + + for (Cell_handle chi : cells_around_edge) + { + Facet f1(chi, chi->index(vh0)); + Facet f2(chi, chi->index(vh2)); + + auto it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f1)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f1, it->second); + + it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f2)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f2, it->second); + } + + for(Cell_handle c : cells_around_edge) + c->reset_cache_validity(); + + return db; + } + else //Non planar flip + { +#ifdef CGAL_FLIP_ON_SURFACE_DISABLE_NON_PLANAR_44_FLIP + return NOT_FLIPPABLE; +#endif + typename C3T3::Surface_patch_index patch = c3t3.surface_patch_index(ch0, ch0->index(vh2)); + CGAL_assertion(patch != typename C3T3::Surface_patch_index()); + + CGAL_assertion(c3t3.is_in_complex(ch0, ch0->index(vh2))); + c3t3.remove_from_complex(ch0, ch0->index(vh2)); + CGAL_assertion(c3t3.is_in_complex(ch1, ch1->index(vh2))); + c3t3.remove_from_complex(ch1, ch1->index(vh2)); + + boost::unordered_map opposite_facet_in_complex; + for (Cell_handle chi : cells_around_edge) + { + Facet f1(chi, chi->index(vh1)); + Facet f2(chi, chi->index(vh3)); + + if (c3t3.is_in_complex(f1)) + { + Surface_patch_index spi = c3t3.surface_patch_index(f1); + opposite_facet_in_complex[c3t3.triangulation().mirror_facet(f1)] = spi; + c3t3.remove_from_complex(f1); + } + if (c3t3.is_in_complex(f2)) + { + Surface_patch_index spi = c3t3.surface_patch_index(f2); + opposite_facet_in_complex[c3t3.triangulation().mirror_facet(f2)] = spi; + c3t3.remove_from_complex(f2); + } + } + + // Top Flip + ch3->set_vertex(ch3->index(vh1), vh5); + ch2->set_vertex(ch2->index(vh3), vh4); + ch0->set_vertex(ch0->index(vh1), vh5); + ch1->set_vertex(ch1->index(vh3), vh4); + + Sliver_removal_result db = VALID_FLIP; + if (!is_well_oriented(tr, ch0) + || !is_well_oriented(tr, ch1) + || !is_well_oriented(tr, ch2) + || !is_well_oriented(tr, ch3)) + db = NOT_FLIPPABLE; + else if (curr_max_cosdh < max_cos_dihedral_angle(tr, ch0, false) + || curr_max_cosdh < max_cos_dihedral_angle(tr, ch1, false) + || curr_max_cosdh < max_cos_dihedral_angle(tr, ch2, false) + || curr_max_cosdh < max_cos_dihedral_angle(tr, ch3, false)) + db = NO_BEST_CONFIGURATION; + + if (db == NOT_FLIPPABLE || db == NO_BEST_CONFIGURATION) + { + ch3->set_vertex(ch3->index(vh5), vh1); + ch2->set_vertex(ch2->index(vh4), vh3); + ch0->set_vertex(ch0->index(vh5), vh1); + ch1->set_vertex(ch1->index(vh4), vh3); + + c3t3.add_to_complex(ch0, ch0->index(vh2), patch); + c3t3.add_to_complex(ch1, ch1->index(vh2), patch); + + for (Cell_handle chi : cells_around_edge) + { + Facet f1(chi, chi->index(vh1)); + Facet f2(chi, chi->index(vh3)); + + auto it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f1)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f1, it->second); + + it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f2)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f2, it->second); + } + + return db; + } + + //Left cells 2-2 flip + Cell_handle n_ch3_vh3 = ch3->neighbor(ch3->index(vh3)); + Cell_handle n_ch2_vh1 = ch2->neighbor(ch2->index(vh1)); + + ch3->set_neighbor(ch3->index(vh3), ch2); + ch3->set_neighbor(ch3->index(vh4), n_ch2_vh1); + n_ch2_vh1->set_neighbor(n_ch2_vh1->index(ch2), ch3); + + ch2->set_neighbor(ch2->index(vh1), ch3); + ch2->set_neighbor(ch2->index(vh5), n_ch3_vh3); + n_ch3_vh3->set_neighbor(n_ch3_vh3->index(ch3), ch2); + + //Right cells 2-2 flip + Cell_handle n_ch0_vh3 = ch0->neighbor(ch0->index(vh3)); + Cell_handle n_ch1_vh1 = ch1->neighbor(ch1->index(vh1)); + + ch0->set_neighbor(ch0->index(vh3), ch1); + ch0->set_neighbor(ch0->index(vh4), n_ch1_vh1); + n_ch1_vh1->set_neighbor(n_ch1_vh1->index(ch1), ch0); + + ch1->set_neighbor(ch1->index(vh1), ch0); + ch1->set_neighbor(ch1->index(vh5), n_ch0_vh3); + n_ch0_vh3->set_neighbor(n_ch0_vh3->index(ch0), ch1); + + for (const Cell_handle ce : cells_around_edge) + { + for (int j = 0; j < 4; j++) + ce->vertex(j)->set_cell(ce); + } + + c3t3.add_to_complex(ch0, ch0->index(vh2), patch); + c3t3.add_to_complex(ch1, ch1->index(vh2), patch); + + for (Cell_handle chi : cells_around_edge) + { + Facet f1(chi, chi->index(vh4)); + Facet f2(chi, chi->index(vh5)); + + auto it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f1)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f1, it->second); + + it = opposite_facet_in_complex.find(c3t3.triangulation().mirror_facet(f2)); + if (it != opposite_facet_in_complex.end()) + c3t3.add_to_complex(f2, it->second); + } + + for (Cell_handle c : cells_around_edge) + c->reset_cache_validity(); + + return VALID_FLIP; + } + + return NOT_FLIPPABLE; +} + +template +std::size_t flipBoundaryEdges( + C3T3& c3t3, + const std::vector& boundary_edges, + SurfaceIndexMapMap& boundary_vertices_valences, + IncidentCellsVectorMap& inc_cells, + const Flip_Criterion& flip_criterion, + CellSelector& cell_selector, + Visitor& visitor) +{ + typedef typename C3T3::Vertex_handle Vertex_handle; + typedef typename C3T3::Cell_handle Cell_handle; + typedef typename C3T3::Facet Facet; + typedef typename C3T3::Edge Edge; + typedef typename C3T3::Surface_patch_index Surface_patch_index; + typedef typename C3T3::Triangulation Tr; + typedef std::pair Edge_vv; + + std::size_t nb_success = 0; + + Tr& tr = c3t3.triangulation(); + + std::vector candidate_edges_for_flip; + for (const Edge& e : boundary_edges) + { + if (!c3t3.is_in_complex(e)) + candidate_edges_for_flip.push_back(make_vertex_pair(e)); + } + + for (const auto& [vh0, vh1] : candidate_edges_for_flip) + { + boost::container::small_vector& inc_vh0 = inc_cells[vh0]; + if (inc_vh0.empty()) + tr.incident_cells(vh0, std::back_inserter(inc_vh0)); + + Cell_handle c; + int i, j; + if (!is_edge_uv(vh0, vh1, inc_vh0, c, i, j)) + continue; + + Edge edge(c, i, j); + std::vector boundary_facets; + const bool on_boundary = is_boundary_edge(edge, c3t3, cell_selector, boundary_facets); + +// if (on_boundary && boundary_facets.empty()) +// { +// std::cerr << vh0->point().point() << "\t " << vh1->point().point() << std::endl; +// bool b = is_boundary_edge(vh0, vh1, c3t3, cell_selector); +// CGAL::Tetrahedral_remeshing::debug::dump_c3t3(c3t3, "dump_c3t3_about_boundary_"); +// CGAL::Tetrahedral_remeshing::debug::dump_facets_in_complex(c3t3, "dump_facets_about_boundary_.off"); +// CGAL::Tetrahedral_remeshing::debug::dump_facets_from_selection( +// c3t3, cell_selector, "dump_facets_from_selection_.off"); +// std::cerr << "valid = " << tr.tds().is_valid(true) << std::endl; +// std::cerr << "boundary = " << b << std::endl; +// CGAL_assertion(on_boundary); +// } +// else if (on_boundary && boundary_facets.size() != 2) +// { +// std::cerr << vh0->point().point() << "\t " << vh1->point().point() << std::endl; +// CGAL::Tetrahedral_remeshing::debug::dump_c3t3(c3t3, "dump_c3t3_about_boundary_"); +// CGAL::Tetrahedral_remeshing::debug::dump_facets(boundary_facets, "dump_boundary_facets.polylines.txt"); +// std::vector dummy_facets; +// bool b = is_boundary_edge(edge, c3t3, cell_selector, dummy_facets, true/**/); +// std::cerr << "boundary = " << b << std::endl; +// } + + if (!on_boundary) + continue; + CGAL_assertion(boundary_facets.size() == 2); + + const Facet& f0 = boundary_facets[0]; + const Facet& f1 = boundary_facets[1]; + + // find 3rd and 4th vertices to flip on surface + const Vertex_handle vh2 = third_vertex(f0, vh0, vh1, tr); + const Vertex_handle vh3 = third_vertex(f1, vh0, vh1, tr); + + CGAL_assertion_code(debug::check_facets(vh0, vh1, vh2, vh3, c3t3)); + + if (!tr.tds().is_edge(vh2, vh3)) // most-likely to happen early exit + { + const Surface_patch_index surfi = c3t3.surface_patch_index(boundary_facets[0]); + + int v0 = boundary_vertices_valences.at(vh0)[surfi]; + int v1 = boundary_vertices_valences.at(vh1)[surfi]; + int v2 = boundary_vertices_valences.at(vh2)[surfi]; + int v3 = boundary_vertices_valences.at(vh3)[surfi]; + int m0 = (boundary_vertices_valences.at(vh0).size() > 1 ? 4 : 6); + int m1 = (boundary_vertices_valences.at(vh1).size() > 1 ? 4 : 6); + int m2 = (boundary_vertices_valences.at(vh2).size() > 1 ? 4 : 6); + int m3 = (boundary_vertices_valences.at(vh3).size() > 1 ? 4 : 6); + + int initial_cost = (v0 - m0)*(v0 - m0) + + (v1 - m1)*(v1 - m1) + + (v2 - m2)*(v2 - m2) + + (v3 - m3)*(v3 - m3); + v0--; + v1--; + v2++; + v3++; + + int final_cost = (v0 - m0)*(v0 - m0) + + (v1 - m1)*(v1 - m1) + + (v2 - m2)*(v2 - m2) + + (v3 - m3)*(v3 - m3); + if (initial_cost > final_cost) + { + CGAL_assertion_code(std::size_t nbf = + std::distance(c3t3.facets_in_complex_begin(), + c3t3.facets_in_complex_end())); + CGAL_assertion_code(std::size_t nbe = + std::distance(c3t3.edges_in_complex_begin(), + c3t3.edges_in_complex_end())); + + Sliver_removal_result db = flip_on_surface(c3t3, edge, vh2, vh3, + inc_cells, + flip_criterion, + visitor); + if (db == VALID_FLIP) + { + CGAL_assertion(tr.tds().is_edge(vh2, vh3)); + Cell_handle c; + int li, lj, lk; + CGAL_assertion_code(bool b =) + tr.tds().is_facet(vh2, vh3, vh0, c, li, lj, lk); + CGAL_assertion(b); + c3t3.add_to_complex(c, (6 - li - lj - lk), surfi); + + CGAL_assertion_code(b = ) + tr.tds().is_facet(vh2, vh3, vh1, c, li, lj, lk); + CGAL_assertion(b); + c3t3.add_to_complex(c, (6 - li - lj - lk), surfi); + + CGAL_assertion_code(std::size_t nbf_post = + std::distance(c3t3.facets_in_complex_begin(), + c3t3.facets_in_complex_end())); + CGAL_assertion(nbf == nbf_post); + CGAL_assertion_code(std::size_t nbe_post = + std::distance(c3t3.edges_in_complex_begin(), + c3t3.edges_in_complex_end())); + CGAL_assertion(nbe == nbe_post); + + boundary_vertices_valences[vh0][surfi]--; + boundary_vertices_valences[vh1][surfi]--; + boundary_vertices_valences[vh2][surfi]++; + boundary_vertices_valences[vh3][surfi]++; + + nb_success++; + } + else + continue; + } + } + } + CGAL_assertion(tr.tds().is_valid()); + + return nb_success; } template @@ -1240,54 +1898,82 @@ void flip_edges(C3T3& c3t3, CGAL_USE(protect_boundaries); typedef typename C3T3::Triangulation T3; typedef typename T3::Vertex_handle Vertex_handle; + typedef typename T3::Cell_handle Cell_handle; + typedef typename T3::Edge Edge; typedef typename std::pair Edge_vv; + typedef typename C3T3::Subdomain_index Subdomain_index; #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE std::cout << "Flip edges..."; std::cout.flush(); - std::size_t nb_flips = 0; + std::size_t nb_flips_in_volume = 0; + std::size_t nb_flips_on_surface = 0; #endif + for (auto c : c3t3.cells_in_complex()) + c->reset_cache_validity();//we will use sliver_value + //to store the cos_dihedral_angle + //const Flip_Criterion criterion = VALENCE_MIN_DH_BASED; - //compute vertices normals map? + std::vector inside_edges; + get_internal_edges(c3t3, + cell_selector, + std::back_inserter(inside_edges)); - // typedef typename C3T3::Surface_patch_index Surface_patch_index; - // typedef std::unordered_map Spi_map; - //if (!protect_boundaries) - //{ - // std::cout << "\tBoundary flips" << std::endl; - // //Boundary flip - // std::vector boundary_vertices_valences; - // std::vector boundary_edges; +// //if (criterion == VALENCE_BASED) +// // flip_inside_edges(inside_edges); +// //else +// //{ +//#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE +// nb_flips += +//#endif +// flip_all_edges(inside_edges, c3t3, MIN_ANGLE_BASED, visitor); +// //} - // collectBoundaryEdges(boundary_edges); + std::unordered_map > inc_cells; - // computeVerticesValences(boundary_vertices_valences); +#ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE + nb_flips_in_volume += +#endif + flip_all_edges(inside_edges, c3t3, inc_cells, MIN_ANGLE_BASED, cell_selector, visitor); + if (!protect_boundaries) + { + typedef typename C3T3::Surface_patch_index Surface_patch_index; + typedef boost::unordered_map Spi_map; + + //Boundary flip + std::vector boundary_edges; + boost::unordered_map boundary_vertices_valences; + boost::unordered_map > vertices_subdomain_indices; + collectBoundaryEdgesAndComputeVerticesValences(c3t3, + cell_selector, + boundary_edges, + boundary_vertices_valences, + vertices_subdomain_indices); + +#ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG + if(!debug::are_cell_orientations_valid(c3t3.triangulation())) + std::cerr << "ERROR in ORIENTATION" << std::endl; +#endif // if (criterion == VALENCE_BASED) // flipBoundaryEdges(boundary_edges, boundary_vertices_valences, VALENCE_BASED); // else - // flipBoundaryEdges(boundary_edges, boundary_vertices_valences, MIN_ANGLE_BASED); - //} - - std::vector inside_edges; - get_internal_edges(c3t3, - cell_selector, - std::back_inserter(inside_edges)); - - //if (criterion == VALENCE_BASED) - // flip_inside_edges(inside_edges); - //else - //{ #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - nb_flips = + nb_flips_on_surface += #endif - flip_all_edges(inside_edges, c3t3, MIN_ANGLE_BASED, cell_selector, visitor); - //} + flipBoundaryEdges(c3t3, boundary_edges, boundary_vertices_valences, + inc_cells, + MIN_ANGLE_BASED, + cell_selector, visitor); + } #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE - std::cout << " done (" << nb_flips << " flips)." << std::endl; + std::cout << "\rFlip edges... done (" + << nb_flips_on_surface << "/" + << nb_flips_in_volume << " surface/volume flips done)." << std::endl; #endif } diff --git a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h index 5106fb1835bf..197753c9784b 100644 --- a/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h +++ b/Tetrahedral_remeshing/include/CGAL/Tetrahedral_remeshing/internal/split_long_edges.h @@ -41,6 +41,7 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e, typedef typename C3t3::Triangulation Tr; typedef typename C3t3::Subdomain_index Subdomain_index; typedef typename C3t3::Surface_patch_index Surface_patch_index; + typedef typename C3t3::Curve_index Curve_index; typedef typename Tr::Geom_traits::Point_3 Point; typedef typename Tr::Facet Facet; typedef typename Tr::Vertex_handle Vertex_handle; @@ -70,6 +71,11 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e, } CGAL_assertion(dimension > 0); + // remove complex edge before splitting + const Curve_index curve_index = (dimension == 1) ? c3t3.curve_index(e) : Curve_index(); + if (dimension == 1) + c3t3.remove_from_complex(e); + struct Cell_info { Subdomain_index subdomain_index_; bool selected_; @@ -194,6 +200,13 @@ typename C3t3::Vertex_handle split_edge(const typename C3t3::Edge& e, // will have its patch tagged from the other side, if needed } + // re-insert complex sub-edges + if (dimension == 1) + { + c3t3.add_to_complex(new_v, v1, curve_index); + c3t3.add_to_complex(new_v, v2, curve_index); + } + set_index(new_v, c3t3); return new_v; @@ -273,7 +286,7 @@ void split_long_edges(C3T3& c3t3, = tr.geom_traits().compute_squared_length_3_object(); FT sqlen = sql(tr.segment(e)); if (sqlen > sq_high) - long_edges.insert(long_edge(make_vertex_pair(e), sqlen)); + long_edges.insert(long_edge(make_vertex_pair(e), sqlen)); } #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG 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 7a9050928706..8118105559a9 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 @@ -32,6 +32,7 @@ #include #include +#include namespace CGAL { @@ -192,6 +193,8 @@ class Adaptive_remesher CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( m_c3t3.triangulation(), "1-c3t3_vertices_after_split"); CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3); + const double mdh = CGAL::Tetrahedral_remeshing::min_dihedral_angle(m_c3t3); + std::cout << "Min dihedral angle = " << mdh << std::endl; #endif #ifdef CGAL_DUMP_REMESHING_STEPS CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "1-split"); @@ -214,6 +217,8 @@ class Adaptive_remesher CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( m_c3t3.triangulation(), "2-c3t3_vertices_after_collapse"); CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3); + const double mdh = CGAL::Tetrahedral_remeshing::min_dihedral_angle(m_c3t3); + std::cout << "Min dihedral angle = " << mdh << std::endl; #endif #ifdef CGAL_DUMP_REMESHING_STEPS CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "2-collapse"); @@ -231,6 +236,8 @@ class Adaptive_remesher CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( m_c3t3.triangulation(), "3-c3t3_vertices_after_flip"); CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3); + const double mdh = CGAL::Tetrahedral_remeshing::min_dihedral_angle(m_c3t3); + std::cout << "Min dihedral angle = " << mdh << std::endl; #endif #ifdef CGAL_DUMP_REMESHING_STEPS CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "3-flip"); @@ -247,6 +254,8 @@ class Adaptive_remesher CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( m_c3t3.triangulation(), "4-c3t3_vertices_after_smooth"); CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3); + const double mdh = CGAL::Tetrahedral_remeshing::min_dihedral_angle(m_c3t3); + std::cout << "Min dihedral angle = " << mdh << std::endl; #endif #ifdef CGAL_DUMP_REMESHING_STEPS CGAL::Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, "4-smooth"); @@ -400,10 +409,11 @@ class Adaptive_remesher const Curve_index default_curve_id = default_curve_index(); for (const Edge& e : tr().finite_edges()) { - if (get(ecmap, CGAL::Tetrahedral_remeshing::make_vertex_pair(e)) + 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) + || nb_incident_surface_patches(e, m_c3t3) > 1 + || nb_incident_complex_facets(e, m_c3t3) > 2)//non-manifold edges { const bool in_complex = m_c3t3.is_in_complex(e); typename C3t3::Curve_index curve_id = in_complex @@ -430,8 +440,12 @@ class Adaptive_remesher Corner_index corner_id = 0; for (Vertex_handle vit : tr().finite_vertex_handles()) { - if ( vit->in_dimension() == 0 - || nb_incident_complex_edges(vit, m_c3t3) > 2) + boost::container::small_vector incident_edges; + incident_complex_edges(vit, m_c3t3, std::back_inserter(incident_edges)); + if ( incident_edges.size() == 1 //tip or endpoint + || incident_edges.size() > 2 //corner + || (incident_edges.size() == 2 + && edges_form_a_sharp_angle(incident_edges, 60, m_c3t3))) { if (!m_c3t3.is_in_complex(vit)) m_c3t3.add_to_complex(vit, ++corner_id); @@ -455,6 +469,8 @@ class Adaptive_remesher std::cout << "\t facets = " << nbf << std::endl; std::cout << "\t edges = " << nbe << std::endl; std::cout << "\t vertices = " << nbv << std::endl; + const double mdh = CGAL::Tetrahedral_remeshing::min_dihedral_angle(m_c3t3); + std::cout << "\t Min dihedral angle = " << mdh << std::endl; CGAL::Tetrahedral_remeshing::debug::dump_vertices_by_dimension( m_c3t3.triangulation(), "0-c3t3_vertices_after_init_"); @@ -600,6 +616,9 @@ class Adaptive_remesher ossi << "statistics_" << it_nb << ".txt"; Tetrahedral_remeshing::internal::compute_statistics( tr(), m_cell_selector, ossi.str().c_str()); + std::ostringstream oss_it; + oss_it << "iteration_" << it_nb; + Tetrahedral_remeshing::debug::dump_c3t3(m_c3t3, oss_it.str().c_str()); #endif #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG CGAL::Tetrahedral_remeshing::debug::check_surface_patch_indices(m_c3t3); 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 9ea06ce7d4ee..cef3b46ed4a1 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 @@ -37,7 +37,7 @@ namespace Tetrahedral_remeshing { enum Subdomain_relation { EQUAL, DIFFERENT, INCLUDED, INCLUDES }; -enum Sliver_removal_result { INVALID_ORIENTATION, INVALID_CELL, INVALID_VERTEX, +enum Sliver_removal_result { INVALID_ORIENTATION = 1, INVALID_CELL, INVALID_VERTEX, NOT_FLIPPABLE, EDGE_PROBLEM, VALID_FLIP, NO_BEST_CONFIGURATION, EXISTING_EDGE }; template @@ -79,6 +79,23 @@ inline int indices(const int& i, const int& j) return 0; } +template +typename Tr::Vertex_handle +third_vertex(const typename Tr::Facet& f, + const typename Tr::Vertex_handle v0, + const typename Tr::Vertex_handle v1, + const Tr& tr) +{ + for(auto v : tr.vertices(f)) + if(v != v0 && v != v1) + return v; + + CGAL_assertion(false); + return + typename Tr::Vertex_handle(); +} + +// returns angle in degrees template std::array vertices(const typename Tr::Edge& e , const Tr&) @@ -164,6 +181,8 @@ template typename Tr::Geom_traits::FT min_dihedral_angle(const Tr& tr, const typename Tr::Cell_handle c) { + if (tr.is_infinite(c)) + return 180.; return min_dihedral_angle(tr, c->vertex(0), c->vertex(1), @@ -171,6 +190,18 @@ typename Tr::Geom_traits::FT min_dihedral_angle(const Tr& tr, c->vertex(3)); } +template +typename C3t3::Triangulation::Geom_traits::FT min_dihedral_angle(const C3t3& c3t3) +{ + typename C3t3::Triangulation::Geom_traits::FT min_dh = 180.; + + for (const auto c : c3t3.cells_in_complex()) + min_dh = (std::min)(min_dh, min_dihedral_angle(c3t3.triangulation(), c)); + + return min_dh; +} + + struct Dihedral_angle_cosine { CGAL::Sign m_sgn; @@ -201,6 +232,20 @@ struct Dihedral_angle_cosine }; } + double value() const + { + switch (m_sgn) + { + case CGAL::POSITIVE: + return CGAL::approximate_sqrt(m_sq_num / m_sq_den); + case CGAL::ZERO: + return 0.; + default: + CGAL_assertion(m_sgn == CGAL::NEGATIVE); + return -1. * CGAL::approximate_sqrt(m_sq_num / m_sq_den); + }; + } + friend bool operator<(const Dihedral_angle_cosine& l, const Dihedral_angle_cosine& r) { @@ -276,8 +321,12 @@ Dihedral_angle_cosine cos_dihedral_angle(const typename Gt::Point_3& i, const double sqden = CGAL::to_double( scalar_product(ijik, ijik) * scalar_product(ilij, ilij)); + const double sqnum = CGAL::square(CGAL::to_double(num)); - return Dihedral_angle_cosine(CGAL::sign(num), CGAL::square(num), sqden); + if(sqnum > sqden) + return Dihedral_angle_cosine(CGAL::sign(num), 1., 1.);//snap to 1 or -1 + else + return Dihedral_angle_cosine(CGAL::sign(num), CGAL::square(num), sqden); } template @@ -310,6 +359,9 @@ Dihedral_angle_cosine max_cos_dihedral_angle(const Point& p, a = cos_dihedral_angle(r, s, p, q, gt); if (max_cos_dh < a) max_cos_dh = a; + CGAL_assertion(max_cos_dh.value() <= 1.); + CGAL_assertion(max_cos_dh.value() >= -1.); + return max_cos_dh; } @@ -329,21 +381,84 @@ Dihedral_angle_cosine max_cos_dihedral_angle(const Tr& tr, template Dihedral_angle_cosine max_cos_dihedral_angle(const Tr& tr, - const typename Tr::Cell_handle c) + const typename Tr::Cell_handle c, + const bool use_cache = true) { - if (c->is_cache_valid()) + if (use_cache && c->is_cache_valid()) return Dihedral_angle_cosine(CGAL::sign(c->sliver_value()), CGAL::abs(c->sliver_value()), 1.); + else if(tr.is_infinite(c)) + return Dihedral_angle_cosine(CGAL::ZERO, 0., 1.); + + typedef typename Tr::Triangulation_data_structure::Cell::Subdomain_index Subdomain_index; + if(c->subdomain_index() == Subdomain_index()) + return Dihedral_angle_cosine(CGAL::ZERO, 0., 1.); Dihedral_angle_cosine cos_dh = max_cos_dihedral_angle(tr, c->vertex(0), c->vertex(1), c->vertex(2), c->vertex(3)); - c->set_sliver_value(cos_dh.signed_square_value()); + if(use_cache) + c->set_sliver_value(cos_dh.signed_square_value()); return cos_dh; } +// incident_edges must share a vertex +template +bool edges_form_a_sharp_angle(const EdgesVector& incident_edges, + const double angle_bound, + const C3t3& c3t3) +{ + CGAL_assertion(incident_edges.size() == 2); + + typedef typename C3t3::Vertex_handle Vertex_handle; + + const std::array ev0 + = c3t3.triangulation().vertices(incident_edges[0]); + const std::array ev1 + = c3t3.triangulation().vertices(incident_edges[1]); + + Vertex_handle shared_vertex, v1, v2; + if (ev0[0] == ev1[0]) + { + shared_vertex = ev0[0]; + v1 = ev0[1]; + v2 = ev1[1]; + } + else if(ev0[0] == ev1[1]) + { + shared_vertex = ev0[0]; + v1 = ev0[1]; + v2 = ev1[0]; + } + else if(ev0[1] == ev1[0]) + { + shared_vertex = ev0[1]; + v1 = ev0[0]; + v2 = ev1[1]; + } + else if(ev0[1] == ev1[1]) + { + shared_vertex = ev0[1]; + v1 = ev0[0]; + v2 = ev1[0]; + } + else + { + CGAL_assertion(false); + return false; + } + CGAL_assertion(shared_vertex != v1); + CGAL_assertion(shared_vertex != v2); + CGAL_assertion(v1 != v2); + + const double angle = CGAL::approximate_angle(point(v1->point()), + point(shared_vertex->point()), + point(v2->point())); + return angle < angle_bound; +} + template bool is_peelable(const C3t3& c3t3, const typename C3t3::Cell_handle ch, @@ -398,14 +513,10 @@ std::pair make_vertex_pair(const Vh v1, const Vh v2) else return std::make_pair(v1, v2); } -template -std::pair -make_vertex_pair(const typename Tr::Edge& e) +template +auto make_vertex_pair(const Edge& e) { - typedef typename Tr::Vertex_handle Vertex_handle; - Vertex_handle v1 = e.first->vertex(e.second); - Vertex_handle v2 = e.first->vertex(e.third); - return make_vertex_pair(v1, v2); + return make_vertex_pair(e.first->vertex(e.second), e.first->vertex(e.third)); } template @@ -460,23 +571,45 @@ bool is_well_oriented(const Tr& tr, template bool is_well_oriented(const Tr& tr, const typename Tr::Cell_handle ch) { - return is_well_oriented(tr, ch->vertex(0), ch->vertex(1), - ch->vertex(2), ch->vertex(3)); + return tr.is_infinite(ch) + || is_well_oriented(tr, + ch->vertex(0), + ch->vertex(1), + ch->vertex(2), + ch->vertex(3)); } template bool is_boundary(const C3T3& c3t3, const typename C3T3::Facet& f, - const CellSelector& cell_selector) + const CellSelector& cell_selector, + const bool verbose = false) { - return c3t3.is_in_complex(f) - || get(cell_selector, f.first) != get(cell_selector, f.first->neighbor(f.second)); + if (c3t3.triangulation().is_infinite(f)) + return false; + + const auto& mf = c3t3.triangulation().mirror_facet(f); + const bool res = c3t3.is_in_complex(f) + || get(cell_selector, f.first) != get(cell_selector, mf.first); + + if (verbose && res) + { + std::cout << "is_boundary(f) :" + << "\n\t in_complex = " << c3t3.is_in_complex(f) + << "\n\t selector(f.first) = " << get(cell_selector, f.first) + << "\n\t selector(mirror ) = " << get(cell_selector, mf.first) + << "\n\t subdomain(f.first)= " << f.first->subdomain_index() + << "\n\t subdomain(mirror) = " << mf.first->subdomain_index() + << std::endl; + } + + return res; } template bool is_boundary(const C3T3& c3t3, const typename C3T3::Triangulation::Edge& e, - CellSelector cell_selector) + const CellSelector& cell_selector) { typedef typename C3T3::Triangulation Tr; typedef typename Tr::Facet_circulator Facet_circulator; @@ -496,6 +629,23 @@ bool is_boundary(const C3T3& c3t3, return false; } +template +typename C3T3::Edge get_edge(const typename C3T3::Vertex_handle v0, + const typename C3T3::Vertex_handle v1, + const C3T3& c3t3) +{ + typedef typename C3T3::Edge Edge; + typedef typename C3T3::Cell_handle Cell_handle; + + Cell_handle cell; + int i0, i1; + if (c3t3.triangulation().tds().is_edge(v0, v1, cell, i0, i1)) + return Edge(cell, i0, i1); + else + CGAL_assertion(false); + return Edge(); +} + template bool is_boundary_edge(const typename C3t3::Vertex_handle& v0, const typename C3t3::Vertex_handle& v1, @@ -513,6 +663,33 @@ bool is_boundary_edge(const typename C3t3::Vertex_handle& v0, return false; } +template +bool is_boundary_edge(const typename C3t3::Edge& e, + const C3t3& c3t3, + const CellSelector& cell_selector, + std::vector& boundary_facets, + const bool verbose = false) +{ + using Facet = typename C3t3::Facet; + + auto fcirc = c3t3.triangulation().incident_facets(e); + auto fend = fcirc; + bool result = false; + + do + { + const Facet& f = *fcirc; + if (is_boundary(c3t3, f, cell_selector, verbose)) + { + boundary_facets.push_back(f); + result = true; + } + } + while (++fcirc != fend); + + return result; +} + template bool is_boundary_vertex(const typename C3t3::Vertex_handle& v, const C3t3& c3t3, @@ -665,6 +842,27 @@ OutputIterator incident_surface_patches(const typename C3t3::Edge& e, return oit; } +template +OutputIterator incident_surface_patches(const typename C3t3::Vertex_handle& v, + const C3t3& c3t3, + OutputIterator oit) +{ + typedef typename C3t3::Triangulation::Facet Facet; + boost::unordered_set facets; + c3t3.triangulation().incident_facets(v, std::inserter(facets, facets.begin())); + + for (typename boost::unordered_set::iterator fit = facets.begin(); + fit != facets.end(); + ++fit) + { + const Facet& f = *fit; + if (c3t3.is_in_complex(f)) + *oit++ = c3t3.surface_patch_index(f); + } + + return oit; +} + template std::size_t nb_incident_subdomains(const typename C3t3::Vertex_handle v, const C3t3& c3t3) @@ -701,6 +899,23 @@ std::size_t nb_incident_surface_patches(const typename C3t3::Edge& e, return indices.size(); } +template +OutputIterator incident_complex_edges(const typename C3t3::Vertex_handle v, + const C3t3& c3t3, + OutputIterator oit) +{ + typedef typename C3t3::Edge Edge; + std::unordered_set edges; + c3t3.triangulation().finite_incident_edges(v, std::inserter(edges, edges.begin())); + + for (const Edge& e : edges) + { + if (c3t3.is_in_complex(e)) + *oit++ = e; + } + return oit; +} + template std::size_t nb_incident_complex_edges(const typename C3t3::Vertex_handle v, const C3t3& c3t3) @@ -718,6 +933,23 @@ std::size_t nb_incident_complex_edges(const typename C3t3::Vertex_handle v, return count; } +template +std::size_t nb_incident_complex_facets(const typename C3t3::Edge& e, + const C3t3& c3t3) +{ + typedef typename C3t3::Triangulation::Facet_circulator Facet_circulator; + Facet_circulator circ = c3t3.triangulation().incident_facets(e); + Facet_circulator end = circ; + std::size_t count = 0; + do + { + const typename C3t3::Facet& f = *circ; + if (c3t3.is_in_complex(f)) + ++count; + } + while (++circ != end); + return count; +} template bool is_feature(const typename C3t3::Vertex_handle v, @@ -944,7 +1176,7 @@ OutputIterator get_internal_edges(const C3t3& c3t3, const typename C3t3::Edge& e = *eit; if (is_internal(e, c3t3, cell_selector)) { - *oit++ = make_vertex_pair(e); + *oit++ = make_vertex_pair(e); } } return oit; @@ -1189,6 +1421,20 @@ void get_edge_info(const typename C3t3::Edge& edge, } } +template +std::array +cell_edges(const typename Tr::Cell_handle c, const Tr&) +{ + using Edge = typename Tr::Edge; + std::array edges_array = { { Edge(c, 0, 1), + Edge(c, 0, 2), + Edge(c, 0, 3), + Edge(c, 1, 2), + Edge(c, 1, 3), + Edge(c, 2, 3) } }; + return edges_array; +} + namespace internal { template @@ -1211,6 +1457,9 @@ namespace internal { //update C3t3 using Subdomain_index = typename C3t3::Subdomain_index; + if(c3t3.is_in_complex(c)) + c3t3.remove_from_complex(c); + if (Subdomain_index() != subdomain) c3t3.add_to_complex(c, subdomain); else @@ -1224,6 +1473,24 @@ namespace internal namespace debug { +template +bool check_facets(const typename C3t3::Vertex_handle vh0, + const typename C3t3::Vertex_handle vh1, + const typename C3t3::Vertex_handle vh2, + const typename C3t3::Vertex_handle vh3, + const C3t3& c3t3) +{ + int li, lj, lk; + typename C3t3::Cell_handle c; + + bool b1 = c3t3.triangulation().is_facet(vh0, vh1, vh2, c, li, lj, lk); + bool b2 = c3t3.is_in_complex(c, (6 - li - lj - lk)); + bool b3 = c3t3.triangulation().is_facet(vh0, vh1, vh3, c, li, lj, lk); + bool b4 = c3t3.is_in_complex(c, (6 - li - lj - lk)); + + return b1 && b2 && b3 && b4; +} + // forward-declaration template void dump_cells(const CellRange& cells, const char* filename); @@ -1242,6 +1509,20 @@ void dump_edges(const Bimap& edges, const char* filename) ofs.close(); } +template +void dump_edges(const std::vector& edges, const char* filename) +{ + std::ofstream ofs(filename); + ofs.precision(17); + + for (const Edge& e : edges) + { + ofs << "2 " << point(e.first->vertex(e.second)->point()) + << " " << point(e.first->vertex(e.third)->point()) << std::endl; + } + ofs.close(); +} + template void dump_facet(const Facet& f, OutputStream& os) { @@ -1264,6 +1545,20 @@ void dump_facets(const FacetRange& facets, const char* filename) os.close(); } +template +void dump_facets_from_selection(const C3t3& c3t3, + const CellSelector& selector, + const char* filename) +{ + std::vector facets; + for (const auto& f : c3t3.triangulation().finite_facets()) + { + if (get(selector, f.first) != get(selector, f.first->neighbor(f.second))) + facets.push_back(f); + } + dump_facets(facets, filename); +} + template void dump_polylines(const CellRange& cells, const char* filename) { diff --git a/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h b/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h index d11d72973b6a..3eac9a222bdd 100644 --- a/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h +++ b/Tetrahedral_remeshing/include/CGAL/tetrahedral_remeshing.h @@ -280,6 +280,10 @@ void tetrahedral_isotropic_remeshing( // perform remeshing std::size_t nb_extra_iterations = 3; +#ifdef CGAL_TETRAHEDRAL_REMESHING_NO_EXTRA_ITERATIONS + nb_extra_iterations = 0; +#endif + remesher.remesh(max_it, nb_extra_iterations); #ifdef CGAL_TETRAHEDRAL_REMESHING_DEBUG @@ -363,7 +367,7 @@ convert_to_triangulation_3( for (auto e : c3t3.edges_in_complex()) { const Edge_vv evv - = CGAL::Tetrahedral_remeshing::make_vertex_pair(e);//ordered pair + = CGAL::Tetrahedral_remeshing::make_vertex_pair(e);//ordered pair put(ecmap, evv, true); } } @@ -501,6 +505,10 @@ void tetrahedral_isotropic_remeshing( // perform remeshing std::size_t nb_extra_iterations = 3; +#ifdef CGAL_TETRAHEDRAL_REMESHING_NO_EXTRA_ITERATIONS + nb_extra_iterations = 0; +#endif + remesher.remesh(max_it, nb_extra_iterations); #ifdef CGAL_TETRAHEDRAL_REMESHING_VERBOSE