diff --git a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h index f2bfe4dcbcd..50ff93d7356 100644 --- a/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h +++ b/Constrained_triangulation_3/include/CGAL/Conforming_constrained_Delaunay_triangulation_3.h @@ -1914,8 +1914,8 @@ private: }; template - int does_edge_intersect_region(Cell_handle cell, int index_vc, int index_vd, - const CDT_2& cdt_2, const Fh_region& fh_region) + int does_edge_interior_intersect_region(Cell_handle cell, int index_vc, int index_vd, + const CDT_2& cdt_2, const Fh_region& fh_region) { const auto vc = cell->vertex(index_vd); const auto vd = cell->vertex(index_vc); @@ -1938,7 +1938,7 @@ private: if(opc == CGAL::COPLANAR || opd == CGAL::COPLANAR || opc == opd) { continue; } else { - // otherwise the segment intersects the plane of the triangle + // otherwise the segment interior intersects the plane of the triangle if(CGAL::orientation(pc, pd, t0, t1) != opc && CGAL::orientation(pc, pd, t1, t2) != opc && CGAL::orientation(pc, pd, t2, t0) != opc) @@ -1987,7 +1987,7 @@ private: //write_segment(dump_edges_around, cell_circ->vertex(index_vc), cell_circ->vertex(index_vd)); if(is_marked(cell_circ->vertex(index_vc), Vertex_marker::REGION_BORDER)) continue; if(is_marked(cell_circ->vertex(index_vd), Vertex_marker::REGION_BORDER)) continue; - int cd_intersects_region = does_edge_intersect_region(cell_circ, index_vc, index_vd, cdt_2, fh_region); + int cd_intersects_region = does_edge_interior_intersect_region(cell_circ, index_vc, index_vd, cdt_2, fh_region); if(cd_intersects_region == 1) { return Search_first_intersection_result_type{ Edge{cell_circ, index_vc, index_vd}, border_edge }; } @@ -2044,7 +2044,7 @@ private: std::set> non_intersecting_edges_set; if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - debug_dump_tetrahedra_intersect_region(face_index, region_index, cdt_2, fh_region); + expensive_debug_dump_tetrahedra_intersect_region(face_index, region_index, cdt_2, fh_region); } detect_edges_and_cells_intersecting_region(face_index, region_index, cdt_2, fh_region, region_border_vertices, @@ -2085,7 +2085,9 @@ private: } } - // use edges of the border facets to unify sets + if(this->debug_regions()) { + std::cerr << "...use edges of the border facets to unify sets\n"; + } for(auto facet: facets_of_border) { auto vertices = tr().vertices(facet); for(int i = 0; i < 3; ++i) { @@ -2098,8 +2100,10 @@ private: } } - // use non-intersecting edges to unify sets, until we have at most 2 sets if(vertices_of_cavity_union_find.number_of_sets() > 2) { + if(this->debug_regions()) { + std::cerr << "...use non-intersecting edges to unify sets, until we have at most 2 sets\n"; + } for(auto c : cr_intersecting_cells) { for(int i = 0; i < 4; ++i) { @@ -2147,10 +2151,15 @@ private: return std::make_pair(Vertex_handle{}, Facet{}); }); CGAL_assume(vertex_above != Vertex_handle{}); - + if(this->debug_regions()) { + std::cerr << "The vertex above the region is " << IO::oformat(vertex_above, with_point_and_info) << "\n"; + } // if there are still more than 2 sets, we need to propagate the information // using a DFS on the border facets if(vertices_of_cavity_union_find.number_of_sets() > 2) { + if(this->debug_regions()) { + std::cerr << "...propagate the information using a DFS on the border facets\n"; + } const auto border_edges_set = std::invoke([&] { using Ordered_pair = std::pair; using Hash = boost::hash; @@ -2527,8 +2536,12 @@ private: // std::ofstream dump(std::string("dump_no_segment_found_") + std::to_string(face_index) + "_" + // std::to_string(region_index) + ".binary.cgal"); // CGAL::IO::save_binary_file(dump, *this); - dump_region(face_index, region_index, cdt_2); // } + if(this->debug_regions()) { + std::cerr << "ERROR: No first segment found intersecting region " << region_index + << " of face #" << face_index << "\n"; + dump_region(face_index, region_index, cdt_2); + } throw Next_region{"No segment found", fh_region[0]}; } CGAL_assertion(found_edge_opt != std::nullopt); @@ -2571,7 +2584,7 @@ private: }; if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - std::cerr << cdt_3_format("Cavity has {} cells and {} edges, " + std::cerr << cdt_3_format("Cavity has {} piercing cells and {} piercing edges, " "{} vertices in upper cavity and {} in lower, " "{} facets in upper cavity and {} in lower\n", original_intersecting_cells.size(), @@ -2580,20 +2593,6 @@ private: original_vertices_of_lower_cavity.size(), original_facets_of_upper_cavity.size(), original_facets_of_lower_cavity.size()); - if(original_intersecting_cells.size() > 3 || intersecting_edges.size() > 1) { - std::cerr << "!! Interesting case !!\n"; - // dump_region(face_index, region_index, cdt_2); - // { - // std::ofstream out(std::string("dump_intersecting_edges_") + std::to_string(face_index) + "_" + - // std::to_string(region_index) + ".polylines.txt"); - // out.precision(17); - // for(auto edge: intersecting_edges) { - // write_segment(out, edge); - // } - // } - // dump_facets_of_cavity(face_index, region_index, "lower", original_facets_of_lower_cavity); - // dump_facets_of_cavity(face_index, region_index, "upper", original_facets_of_upper_cavity); - } } auto register_internal_constrained_facet = [this](Facet f) { this->register_facet_to_be_constrained(f); }; @@ -3234,9 +3233,10 @@ private: bool restore_face(CDT_3_signed_index face_index) { CDT_2& non_const_cdt_2 = face_cdt_2[face_index]; const CDT_2& cdt_2 = non_const_cdt_2; - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_copy_triangulation_into_hole()) { - std::cerr << cdt_3_format("restore_face({}): CDT_2 has {} vertices\n", face_index, cdt_2.number_of_vertices()); - } + if constexpr (cdt_3_can_use_cxx20_format()) + if(this->debug_copy_triangulation_into_hole() || this->debug_verbose_special_cases()) { + std::cerr << cdt_3_format("restore_face({}): CDT_2 has {} vertices\n", face_index, cdt_2.number_of_vertices()); + } for(const auto& edge : cdt_2.finite_edges()) { const auto fh = edge.first; const auto i = edge.second; @@ -3280,7 +3280,7 @@ private: processed_faces.insert(fh_region.begin(), fh_region.end()); auto handle_error_with_region = [&](const char* what, CDT_2_face_handle fh_2d) { - if(this->debug_regions()) { + if(this->debug_regions() || this->debug_verbose_special_cases()) { std::cerr << "NOTE: " << what << " in sub-region " << (region_index - 1) << " of face #F" << face_index << '\n'; } @@ -3490,7 +3490,9 @@ public: i = face_constraint_misses_subfaces_find_next(i); if(i == npos) { std::cerr << "ERROR: No more missing face to restore after a PLC error\n"; - dump_region(e.face_index, e.region_index); + if(this->debug_regions() || this->debug_verbose_special_cases()) { + dump_region(e.face_index, e.region_index); + } throw; } std::cerr << "Next face is face #F " << i << '\n'; @@ -3625,9 +3627,6 @@ public: return visited_edges.emplace(CGAL::make_sorted_pair(v0, v1), does_intersect); }; - // Alias for C++20 extension warning workaround - auto& intersecting_edges_ = intersecting_edges; - intersecting_edges.push_back(first_intersecting_edge); const auto [v0, v1] = tr().vertices(first_intersecting_edge); (void)new_edge(v0, v1, true); @@ -3640,24 +3639,25 @@ public: } auto test_edge = [&](Cell_handle cell, Vertex_handle v0, int index_v0, Vertex_handle v1, int index_v1, - int expected) { - auto value_returned = [this](bool b) { + int expected) + { + auto value_returned = [this, v0, v1](bool b, bool not_visited) { CGAL_USE(this); if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - std::cerr << cdt_3_format(" return {}\n", b); + std::cerr << cdt_3_format(" test_edge {} {} return {} {}\n", + IO::oformat(v0, with_point_and_info), + IO::oformat(v1, with_point_and_info), + b, + not_visited ? "(new)" : "(cached)"); } return b; }; - if constexpr (cdt_3_can_use_cxx20_format()) if(this->debug_regions()) { - std::cerr << cdt_3_format(" test_edge {} {} ", IO::oformat(v0, with_point_and_info), - IO::oformat(v1, with_point_and_info)); - } auto [cached_value_it, not_visited] = new_edge(v0, v1, false); - if(!not_visited) return value_returned(cached_value_it->second); - int v0v1_intersects_region = (is_marked(v0, Vertex_marker::REGION_INSIDE) || - is_marked(v1, Vertex_marker::REGION_INSIDE)) - ? expected - : does_edge_intersect_region(cell, index_v0, index_v1, cdt_2, fh_region); + if(!not_visited) return value_returned(cached_value_it->second, not_visited); + int v0v1_intersects_region = + (is_marked(v0, Vertex_marker::REGION_INSIDE) || is_marked(v1, Vertex_marker::REGION_INSIDE)) + ? expected + : does_edge_interior_intersect_region(cell, index_v0, index_v1, cdt_2, fh_region); if(v0v1_intersects_region != 0) { if(this->use_older_cavity_algorithm()) { if(v0v1_intersects_region != expected) { @@ -3669,13 +3669,13 @@ public: if(v0v1_intersects_region < 0) { std::swap(index_v0, index_v1); } - intersecting_edges_.emplace_back(cell, index_v0, index_v1); + intersecting_edges.emplace_back(cell, index_v0, index_v1); cached_value_it->second = true; - return value_returned(true); + return value_returned(true, not_visited); } else { non_intersecting_edges_set.insert(make_sorted_pair(v0, v1)); cached_value_it->second = false; - return value_returned(false); + return value_returned(false, not_visited); } }; @@ -3741,10 +3741,10 @@ public: { intersecting_cells.insert(n_ch); if(this->debug_regions()) { - std::cerr << "tetrahedron #" << n_ch->time_stamp() << " intersects the region\n"; + std::cerr << "new tetrahedron #" << n_ch->time_stamp() << " intersects the region\n"; } } else if(this->debug_regions()) { - std::cerr << "NO, tetrahedron #" << n_ch->time_stamp() << " does not intersect the region\n"; + std::cerr << "NO, new tetrahedron #" << n_ch->time_stamp() << " does not intersect the region\n"; } for(int i = 0; i < 4; ++i) { for(int j = i + 1; j < 4; ++j) { @@ -3820,12 +3820,13 @@ public: { using EK = CGAL::Exact_predicates_exact_constructions_kernel; const auto to_exact = CGAL::Cartesian_converter(); + const auto& cdt_2 = this->face_cdt_2[face_index]; std::cerr << cdt_3_format("restore_subface_region face index: {}, region #{}, intersecting edge #{}: ({} {})\n", face_index, region_index, edge_index, IO::oformat(v_above, with_point_and_info), IO::oformat(v_below, with_point_and_info)); - dump_region(face_index, region_index, this->face_cdt_2[face_index]); + dump_region(face_index, region_index, cdt_2); const auto p_above = this->point(v_above); const auto p_below = this->point(v_below); @@ -3851,30 +3852,17 @@ public: intersect_out.close(); auto cells_around_intersecting_edge = Container_from_circulator{this->incident_cells(intersecting_edge)}; - for(const auto& cell: cells_around_intersecting_edge) { - CGAL_assertion(!cell.has_vertex(tr().infinite_vertex())); - auto tetrahedron = - typename Geom_traits::Tetrahedron_3{tr().point(cell.vertex(0)), tr().point(cell.vertex(1)), - tr().point(cell.vertex(2)), tr().point(cell.vertex(3))}; - for(auto fh: fh_region) { - auto v0 = fh->vertex(0)->info().vertex_handle_3d; - auto v1 = fh->vertex(1)->info().vertex_handle_3d; - auto v2 = fh->vertex(2)->info().vertex_handle_3d; - auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; - auto exact_triangle = to_exact(triangle); - } - + for(const auto& ch: make_prevent_deref_range(cells_around_intersecting_edge)) { + CGAL_assertion(!ch->has_vertex(tr().infinite_vertex())); + auto tetrahedron = tr().tetrahedron(ch.current_circulator()); std::cerr << cdt_3_format("Test tetrahedron (#{}):\n {}\n {}\n {}\n {}\n", - cell.time_stamp(), - IO::oformat(cell.vertex(0), with_point_and_info), - IO::oformat(cell.vertex(1), with_point_and_info), - IO::oformat(cell.vertex(2), with_point_and_info), - IO::oformat(cell.vertex(3), with_point_and_info)); + ch->time_stamp(), + IO::oformat(ch->vertex(0), with_point_and_info), + IO::oformat(ch->vertex(1), with_point_and_info), + IO::oformat(ch->vertex(2), with_point_and_info), + IO::oformat(ch->vertex(3), with_point_and_info)); if(!std::any_of(fh_region.begin(), fh_region.end(), [&](const auto fh) { - auto v0 = fh->vertex(0)->info().vertex_handle_3d; - auto v1 = fh->vertex(1)->info().vertex_handle_3d; - auto v2 = fh->vertex(2)->info().vertex_handle_3d; - auto triangle = typename Geom_traits::Triangle_3{tr().point(v0), tr().point(v1), tr().point(v2)}; + auto triangle = cdt_2.triangle(fh); bool b = does_tetrahedron_intersect_triangle_interior(tetrahedron, triangle, tr().geom_traits()); if(b) { std::cerr << " intersects the region\n"; @@ -3884,18 +3872,18 @@ public: { std::cerr << cdt_3_format( "ERROR: The following tetrahedron (#{}) does not intersect the region:\n {}\n {}\n {}\n {}\n", - cell.time_stamp(), - IO::oformat(cell.vertex(0), with_point_and_info), IO::oformat(cell.vertex(1), with_point_and_info), - IO::oformat(cell.vertex(2), with_point_and_info), IO::oformat(cell.vertex(3), with_point_and_info)); + ch->time_stamp(), + IO::oformat(ch->vertex(0), with_point_and_info), IO::oformat(ch->vertex(1), with_point_and_info), + IO::oformat(ch->vertex(2), with_point_and_info), IO::oformat(ch->vertex(3), with_point_and_info)); } } } template - void debug_dump_tetrahedra_intersect_region(CDT_3_signed_index face_index, - int region_index, - const CDT_2& cdt_2, - const Fh_region& fh_region) + void expensive_debug_dump_tetrahedra_intersect_region(CDT_3_signed_index face_index, + int region_index, + const CDT_2& cdt_2, + const Fh_region& fh_region) { using Mesh = Surface_mesh; using Face_index = typename Mesh::Face_index; @@ -3922,7 +3910,7 @@ public: bool intersects = false; for(int i = 0; i < 4; ++i) { for(int j = i + 1; j < 4; ++j) { - int intersects_region = does_edge_intersect_region(ch, i, j, cdt_2, fh_region); + int intersects_region = does_edge_interior_intersect_region(ch, i, j, cdt_2, fh_region); if(intersects_region != 0) { intersects = true; }