diff --git a/.gitignore b/.gitignore index 97a3235e0..de097e94b 100644 --- a/.gitignore +++ b/.gitignore @@ -11,3 +11,4 @@ LocalPreferences.toml docs/src/applications/cell_simulation.mp4 varying_weight_power.mp4 varying_weight.mp4 +*settings.local.json diff --git a/NEWS.md b/NEWS.md index 87a33c704..82c8730b7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # Changelog +## 1.6.6 + +- The `Triangulation` struct now stores a `has_ghosts::Bool` field and a `boundary_vertex_to_ghost` map for efficient lookup. See [#240](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/240). + ## 1.6.5 - Clarified the counter-clockwise requirement for the `clip_polygon` argument in `voronoi`. See [#230](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/230). diff --git a/Project.toml b/Project.toml index 504126d89..7dd676f0b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DelaunayTriangulation" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" authors = ["Daniel VandenHeuvel "] -version = "1.6.6" +version = "1.6.7" [deps] AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7" diff --git a/src/algorithms/triangulation/basic_operations/add_boundary_information.jl b/src/algorithms/triangulation/basic_operations/add_boundary_information.jl index 07a7c83da..5b8a03433 100644 --- a/src/algorithms/triangulation/basic_operations/add_boundary_information.jl +++ b/src/algorithms/triangulation/basic_operations/add_boundary_information.jl @@ -7,6 +7,8 @@ function add_boundary_information!(tri::Triangulation) I = integer_type(tri) ghost_vertex = I(𝒢) bn = get_boundary_nodes(tri) + # Clear any existing boundary_vertex_to_ghost entries before repopulating + empty!(get_boundary_vertex_to_ghost(tri)) if has_multiple_curves(tri) add_boundary_curve_information!(tri, bn, ghost_vertex) elseif has_multiple_sections(tri) @@ -18,14 +20,18 @@ function add_boundary_information!(tri::Triangulation) end function add_boundary_node_information!(tri::Triangulation, bn, ghost_vertex) n_edge = num_boundary_edges(bn) + n_edge == 0 && return ghost_vertex u = get_boundary_nodes(bn, n_edge + 1) for j in n_edge:-1:1 v = get_boundary_nodes(bn, j) add_adjacent!(tri, u, v, ghost_vertex) add_adjacent2vertex!(tri, ghost_vertex, u, v) add_neighbour!(tri, ghost_vertex, u, v) + add_boundary_vertex_to_ghost!(tri, u, ghost_vertex) u = v end + # Add the last vertex + add_boundary_vertex_to_ghost!(tri, u, ghost_vertex) return ghost_vertex end function add_boundary_segment_information!(tri::Triangulation, bn, ghost_vertex) diff --git a/src/algorithms/triangulation/basic_operations/add_ghost_triangles.jl b/src/algorithms/triangulation/basic_operations/add_ghost_triangles.jl index ff1cbff47..7607c8ab8 100644 --- a/src/algorithms/triangulation/basic_operations/add_ghost_triangles.jl +++ b/src/algorithms/triangulation/basic_operations/add_ghost_triangles.jl @@ -5,15 +5,31 @@ Adds all the ghost triangles to `tri`. """ function add_ghost_triangles!(tri::Triangulation) T = get_triangles(tri) - for g in each_ghost_vertex(tri) + # If we have boundary nodes but the boundary_vertex_to_ghost map is empty, + # we need to populate it. This can happen when ghost triangles were deleted + # and are now being re-added. + if has_boundary_nodes(tri) && isempty(get_boundary_vertex_to_ghost(tri)) + add_boundary_information!(tri) + end + # Use all_ghost_vertices instead of each_ghost_vertex because each_ghost_vertex + # depends on has_ghost_vertices(graph) which may be false if ghost vertices + # were removed from the graph by clear_empty_vertices!. + adj2v_dict = get_adjacent2vertex(get_adjacent2vertex(tri)) + for g in all_ghost_vertices(tri) + # Skip ghost vertices that don't have any edges in adjacent2vertex. + # This can happen when ghost_vertex_ranges has a default entry but no + # actual boundary data was set up. + haskey(adj2v_dict, g) || continue for e in each_edge(get_adjacent2vertex(tri, g)) u, v = edge_vertices(e) add_adjacent!(tri, v, g, u) add_adjacent!(tri, g, u, v) add_adjacent2vertex!(tri, u, v, g) add_adjacent2vertex!(tri, v, g, u) + add_neighbour!(tri, g, u, v) add_triangle!(T, u, v, g) end end + set_has_ghosts!(tri, true) return tri end diff --git a/src/algorithms/triangulation/basic_operations/delete_ghost_triangles.jl b/src/algorithms/triangulation/basic_operations/delete_ghost_triangles.jl index cb32b441b..415719b55 100644 --- a/src/algorithms/triangulation/basic_operations/delete_ghost_triangles.jl +++ b/src/algorithms/triangulation/basic_operations/delete_ghost_triangles.jl @@ -5,13 +5,13 @@ Deletes all the ghost triangles from `tri`. !!! warning "Ghost vertices" - Ghost vertices are still used in the `keys` of the [`Adjacent2Vertex`](@ref) - of `tri`, and are still present in the [`Graph`](@ref). If you want to delete the - ghost vertex `keys` from the [`Adjacent2Vertex`](@ref), you need to use - [`delete_adjacent2vertex!`](@ref). For deleting the ghost vertices from the - [`Graph`](@ref), you need [`delete_ghost_vertices_from_graph!`](@ref). Additionally, - edges in [`Adjacent`](@ref) can still map to ghost vertices. If you also want to delete - those, you need to filter through the `values` of the [`Adjacent`](@ref) map + Ghost vertices are still used in the `keys` of the [`Adjacent2Vertex`](@ref) + of `tri`, and are still present in the [`Graph`](@ref). If you want to delete the + ghost vertex `keys` from the [`Adjacent2Vertex`](@ref), you need to use + [`delete_adjacent2vertex!`](@ref). For deleting the ghost vertices from the + [`Graph`](@ref), you need [`delete_ghost_vertices_from_graph!`](@ref). Additionally, + edges in [`Adjacent`](@ref) can still map to ghost vertices. If you also want to delete + those, you need to filter through the `values` of the [`Adjacent`](@ref) map that are ghost vertices, and use [`delete_adjacent!`](@ref). """ function delete_ghost_triangles!(tri::Triangulation) @@ -26,5 +26,7 @@ function delete_ghost_triangles!(tri::Triangulation) delete_triangle!(T, u, v, g) end end + empty!(get_boundary_vertex_to_ghost(tri)) + set_has_ghosts!(tri, false) return tri end diff --git a/src/algorithms/triangulation/basic_operations/lock_convex_hull.jl b/src/algorithms/triangulation/basic_operations/lock_convex_hull.jl index f33b8948e..244408769 100644 --- a/src/algorithms/triangulation/basic_operations/lock_convex_hull.jl +++ b/src/algorithms/triangulation/basic_operations/lock_convex_hull.jl @@ -30,16 +30,24 @@ function lock_convex_hull!(tri::Triangulation; rng::Random.AbstractRNG = Random. interior_segments = get_interior_segments(tri) interior_segments_on_hull = get_interior_segments_on_hull(get_cache(tri)) empty!(interior_segments_on_hull) + ghost_vertex = I(𝒢) for i in 1:ne u = get_boundary_nodes(bn, i) v = get_boundary_nodes(bn, i + 1) e = construct_edge(E, u, v) bnn_map[e] = (bn, i) + # Populate the boundary_vertex_to_ghost map to maintain consistency + # with the boundary_nodes. This prevents add_ghost_triangles! from + # needing to repopulate it later. + add_boundary_vertex_to_ghost!(tri, u, ghost_vertex) if contains_unoriented_edge(e, interior_segments) delete_unoriented_edge!(interior_segments, e) add_edge!(interior_segments_on_hull, e) end end + # Add the last vertex (closed boundary) + u = get_boundary_nodes(bn, ne + 1) + add_boundary_vertex_to_ghost!(tri, u, ghost_vertex) for e in keys(bnn_map) add_segment!(tri, e; rng, predicates) end diff --git a/src/algorithms/triangulation/basic_operations/unlock_convex_hull.jl b/src/algorithms/triangulation/basic_operations/unlock_convex_hull.jl index a8ec74930..85f8259be 100644 --- a/src/algorithms/triangulation/basic_operations/unlock_convex_hull.jl +++ b/src/algorithms/triangulation/basic_operations/unlock_convex_hull.jl @@ -25,6 +25,9 @@ function unlock_convex_hull!(tri::Triangulation; reconstruct = false) bnn_map = get_boundary_edge_map(tri) empty!(bn) bn_map[I(𝒢)] = bn + # Clear the boundary_vertex_to_ghost map since we're unlocking the convex hull + # and these are no longer constrained boundary nodes + empty!(get_boundary_vertex_to_ghost(tri)) segments = get_interior_segments(tri) all_segments = get_all_segments(tri) for e in keys(bnn_map) diff --git a/src/algorithms/triangulation/constrained_triangulation.jl b/src/algorithms/triangulation/constrained_triangulation.jl index 441b17b7c..a721c52ca 100644 --- a/src/algorithms/triangulation/constrained_triangulation.jl +++ b/src/algorithms/triangulation/constrained_triangulation.jl @@ -332,6 +332,8 @@ function remake_triangulation_with_constraints(tri::Triangulation, segments, bou polygon_hierarchy = get_polygon_hierarchy(tri) boundary_enricher = get_boundary_enricher(tri) cache = get_cache(tri) + has_ghosts_flag = has_ghosts(tri) + boundary_vertex_to_ghost_map = get_boundary_vertex_to_ghost(tri) return new_ghost_vertex_map, new_ghost_vertex_ranges, Triangulation( points, triangles, @@ -351,6 +353,8 @@ function remake_triangulation_with_constraints(tri::Triangulation, segments, bou polygon_hierarchy, boundary_enricher, cache, + has_ghosts_flag, + boundary_vertex_to_ghost_map, ) end @@ -385,6 +389,8 @@ function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map, representative_point_list = get_representative_point_list(tri) boundary_enricher = get_boundary_enricher(tri) cache = get_cache(tri) + has_ghosts_flag = has_ghosts(tri) + boundary_vertex_to_ghost_map = get_boundary_vertex_to_ghost(tri) return Triangulation( points, triangles, @@ -404,6 +410,8 @@ function replace_ghost_vertex_information(tri::Triangulation, ghost_vertex_map, polygon_hierarchy, boundary_enricher, cache, + has_ghosts_flag, + boundary_vertex_to_ghost_map, ) end diff --git a/src/algorithms/triangulation/triangulate_convex.jl b/src/algorithms/triangulation/triangulate_convex.jl index 97b9e48ea..ffb37c0d4 100644 --- a/src/algorithms/triangulation/triangulate_convex.jl +++ b/src/algorithms/triangulation/triangulate_convex.jl @@ -159,27 +159,36 @@ function postprocess_triangulate_convex!(tri::Triangulation, S; delete_ghosts, d append!(hull, S) push!(hull, S[begin]) I = integer_type(tri) + ghost_vertex = I(𝒢) + # Only populate boundary_vertex_to_ghost map for triangulations with actual constrained boundary nodes. + # Convex hull vertices (from unconstrained triangulations) should not be in this map. + # This matches the behavior of add_ghost_triangles! + populate_map = has_boundary_nodes(tri) && !delete_ghosts for i in firstindex(S):(lastindex(S) - 1) u = S[i] v = S[i + 1] if !delete_ghosts - add_triangle!(tri, v, u, I(𝒢); protect_boundary = true, update_ghost_edges = false) + add_triangle!(tri, v, u, ghost_vertex; protect_boundary = true, update_ghost_edges = false) else # Still want the ghost vertex in the graph, adjacent2vertex map, and the adjacent map. - add_neighbour!(tri, I(𝒢), u) - add_adjacent2vertex!(tri, I(𝒢), v, u) - add_adjacent!(tri, v, u, I(𝒢)) + add_neighbour!(tri, ghost_vertex, u) + add_adjacent2vertex!(tri, ghost_vertex, v, u) + add_adjacent!(tri, v, u, ghost_vertex) end + # Add boundary vertex to ghost mapping only for constrained boundaries + populate_map && add_boundary_vertex_to_ghost!(tri, u, ghost_vertex) end u = S[end] v = S[begin] if !delete_ghosts - add_triangle!(tri, v, u, I(𝒢); protect_boundary = true, update_ghost_edges = false) + add_triangle!(tri, v, u, ghost_vertex; protect_boundary = true, update_ghost_edges = false) else - add_neighbour!(tri, I(𝒢), u) - add_adjacent2vertex!(tri, I(𝒢), v, u) - add_adjacent!(tri, v, u, I(𝒢)) + add_neighbour!(tri, ghost_vertex, u) + add_adjacent2vertex!(tri, ghost_vertex, v, u) + add_adjacent!(tri, v, u, ghost_vertex) end + populate_map && add_boundary_vertex_to_ghost!(tri, u, ghost_vertex) + set_has_ghosts!(tri, !delete_ghosts) delete_empty_features && clear_empty_features!(tri) empty_representative_points!(tri) cx, cy = mean_points(get_points(tri), S) diff --git a/src/algorithms/triangulation/unconstrained_triangulation.jl b/src/algorithms/triangulation/unconstrained_triangulation.jl index d9382621c..8f8002947 100644 --- a/src/algorithms/triangulation/unconstrained_triangulation.jl +++ b/src/algorithms/triangulation/unconstrained_triangulation.jl @@ -83,6 +83,7 @@ function initialise_bowyer_watson!(tri::Triangulation, insertion_order, predicat add_triangle!(tri, v, u, g; protect_boundary = true, update_ghost_edges = false) add_triangle!(tri, w, v, g; protect_boundary = true, update_ghost_edges = false) add_triangle!(tri, u, w, g; protect_boundary = true, update_ghost_edges = false) + set_has_ghosts!(tri, true) new_representative_point!(tri, I(1)) for i in triangle_vertices(initial_triangle) p = get_point(tri, i) diff --git a/src/data_structures/triangulation/methods/boundary_nodes.jl b/src/data_structures/triangulation/methods/boundary_nodes.jl index 42981d25e..be55eb871 100644 --- a/src/data_structures/triangulation/methods/boundary_nodes.jl +++ b/src/data_structures/triangulation/methods/boundary_nodes.jl @@ -126,11 +126,22 @@ function split_boundary_edge!(tri::Triangulation, i, j, node) delete_unoriented_edge!(segments, construct_edge(E, i, j)) !contains_edge(construct_edge(E, node, i), segments) && add_edge!(segments, construct_edge(E, i, node)) !contains_edge(construct_edge(E, j, node), segments) && add_edge!(segments, construct_edge(E, node, j)) - # Sometimes, the user might accidentally also put an interior segment in that forms part of the boundary. Let's remove it. + # Sometimes, the user might accidentally also put an interior segment in that forms part of the boundary. Let's remove it. interior_segments = get_interior_segments(tri) if contains_unoriented_edge(construct_edge(E, i, j), interior_segments) delete_unoriented_edge!(interior_segments, construct_edge(E, i, j)) end + # Update boundary_vertex_to_ghost map for the new boundary node + # We need to find the ghost vertex for the section being split, not just use + # vertex i's ghost (which may be from a different section in multiply-connected domains) + section_path = pos[1] + ghost_vertex_map = get_ghost_vertex_map(tri) + for (g, path) in ghost_vertex_map + if path == section_path + add_boundary_vertex_to_ghost!(tri, node, g) + break + end + end return tri end @@ -149,6 +160,7 @@ function merge_boundary_edge!(tri::Triangulation, i, j, node) node_pos = (pos[1], pos[2] + 1) bnn = get_boundary_edge_map(tri) delete_boundary_node!(tri, node_pos) + delete_boundary_vertex_from_ghost_map!(tri, node) E = edge_type(tri) delete!(bnn, construct_edge(E, i, node)) delete!(bnn, construct_edge(E, node, j)) diff --git a/src/data_structures/triangulation/methods/checks.jl b/src/data_structures/triangulation/methods/checks.jl index f36116455..1fba91147 100644 --- a/src/data_structures/triangulation/methods/checks.jl +++ b/src/data_structures/triangulation/methods/checks.jl @@ -23,15 +23,10 @@ unoriented_edge_exists(tri::Triangulation, i, j) = unoriented_edge_exists(tri, c has_ghost_triangles(tri::Triangulation) -> Bool Returns `true` if `tri` has ghost triangles, and `false` otherwise. + +This is an alias for [`has_ghosts`](@ref). """ -function has_ghost_triangles(tri::Triangulation) - num_ghost_vertices(tri) == 0 && return false - I = integer_type(tri) - some_outer_boundary_edges = get_adjacent2vertex(tri, I(𝒢)) - isempty(some_outer_boundary_edges) && return false - e = (first ∘ each_edge)(some_outer_boundary_edges) - return edge_exists(tri, terminal(e), I(𝒢)) -end +has_ghost_triangles(tri::Triangulation) = has_ghosts(tri) """ is_constrained(tri::Triangulation) -> Bool diff --git a/src/data_structures/triangulation/triangulation.jl b/src/data_structures/triangulation/triangulation.jl index 06e0d9916..415bfa1ce 100644 --- a/src/data_structures/triangulation/triangulation.jl +++ b/src/data_structures/triangulation/triangulation.jl @@ -83,30 +83,40 @@ The [`BoundaryEnricher`](@ref) used for triangulating a curve-bounded domain. If A [`TriangulationCache`](@ref) used as a cache for [`add_segment!`](@ref) which requires a separate `Triangulation` structure for use. This will not contain any segments or boundary nodes. Also stores segments useful for [`lock_convex_hull!`](@ref) and [`unlock_convex_hull!`](@ref). +- `has_ghosts::Bool` + +A flag indicating whether the triangulation currently has ghost triangles. This is toggled by [`add_ghost_triangles!`](@ref) and [`delete_ghost_triangles!`](@ref). +- `boundary_vertex_to_ghost::Dict{I, I}` + +A `Dict` that maps each boundary vertex to its corresponding ghost vertex. This provides O(1) lookup for [`is_boundary_node`](@ref) +instead of iterating through all ghost vertices. This map is maintained automatically when boundary nodes are added or removed from the triangulation. +For unconstrained triangulations with no boundary nodes, this will be an empty `Dict`. """ -struct Triangulation{P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,BE} +mutable struct Triangulation{P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,BE} # Geometry - points::P - triangles::T - boundary_nodes::BN - interior_segments::Es - all_segments::Es - weights::W + const points::P + const triangles::T + const boundary_nodes::BN + const interior_segments::Es + const all_segments::Es + const weights::W # Topology - adjacent::Adjacent{I,E} - adjacent2vertex::Adjacent2Vertex{I,Es} - graph::Graph{I} + const adjacent::Adjacent{I,E} + const adjacent2vertex::Adjacent2Vertex{I,Es} + const graph::Graph{I} # Boundary handling - boundary_curves::BC - boundary_edge_map::BEM - ghost_vertex_map::GVM - ghost_vertex_ranges::GVR - # Other - convex_hull::ConvexHull{P,I} - representative_point_list::BPL - polygon_hierarchy::PolygonHierarchy{I} - boundary_enricher::BE - cache::C + const boundary_curves::BC + const boundary_edge_map::BEM + const ghost_vertex_map::GVM + const ghost_vertex_ranges::GVR + # Other + const convex_hull::ConvexHull{P,I} + const representative_point_list::BPL + const polygon_hierarchy::PolygonHierarchy{I} + const boundary_enricher::BE + const cache::C + has_ghosts::Bool + const boundary_vertex_to_ghost::Dict{I,I} end function Base.copy(tri::Triangulation) @@ -126,15 +136,17 @@ function Base.copy(tri::Triangulation) convex_hull = copy(get_convex_hull(tri)) representative_point_list = copy(get_representative_point_list(tri)) polygon_hierarchy = copy(get_polygon_hierarchy(tri)) - boundary_enricher = enrcopy(get_boundary_enricher(tri); points, + boundary_enricher = enrcopy(get_boundary_enricher(tri); points, boundary_nodes, segments=interior_segments,boundary_curves, - polygon_hierarchy,boundary_edge_map) - cache = _copy_cache(get_cache(tri); weights) + polygon_hierarchy,boundary_edge_map) + cache = _copy_cache(get_cache(tri); weights) + has_ghosts_flag = has_ghosts(tri) + boundary_vertex_to_ghost_map = copy(get_boundary_vertex_to_ghost(tri)) return Triangulation( points, triangles, boundary_nodes, interior_segments, all_segments, weights, adjacent, adjacent2vertex, graph, boundary_curves, boundary_edge_map, ghost_vertex_map, ghost_vertex_ranges, convex_hull, representative_point_list, - polygon_hierarchy, boundary_enricher, cache + polygon_hierarchy, boundary_enricher, cache, has_ghosts_flag, boundary_vertex_to_ghost_map ) end @@ -151,6 +163,8 @@ function Base.show(io::IO, ::MIME"text/plain", tri::Triangulation) end function Base.:(==)(tri1::Triangulation, tri2::Triangulation) + has_ghosts(tri1) ≠ has_ghosts(tri2) && return false + get_boundary_vertex_to_ghost(tri1) ≠ get_boundary_vertex_to_ghost(tri2) && return false !has_ghost_triangles(tri1) && has_ghost_triangles(tri2) && return false !has_ghost_triangles(tri2) && has_ghost_triangles(tri1) && return false get_points(tri1) ≠ get_points(tri2) && return false @@ -358,6 +372,53 @@ Returns the cache of `tri`. This is a [`TriangulationCache`](@ref) used as a cac """ get_cache(tri::Triangulation) = tri.cache +""" + has_ghosts(tri::Triangulation) -> Bool + +Returns `true` if the triangulation currently has ghost triangles, and `false` otherwise. +This flag is toggled by [`add_ghost_triangles!`](@ref) and [`delete_ghost_triangles!`](@ref). +""" +has_ghosts(tri::Triangulation) = tri.has_ghosts + +""" + set_has_ghosts!(tri::Triangulation, flag::Bool) + +Sets the `has_ghosts` flag of the triangulation to `flag`. +""" +function set_has_ghosts!(tri::Triangulation, flag::Bool) + tri.has_ghosts = flag + return tri +end + +""" + get_boundary_vertex_to_ghost(tri::Triangulation) -> Dict{I, I} + +Returns the `Dict` that maps each boundary vertex to its corresponding ghost vertex. +This provides O(1) lookup for [`is_boundary_node`](@ref) instead of iterating through +all ghost vertices. +""" +get_boundary_vertex_to_ghost(tri::Triangulation) = tri.boundary_vertex_to_ghost + +""" + add_boundary_vertex_to_ghost!(tri::Triangulation, vertex, ghost_vertex) + +Adds a mapping from `vertex` to `ghost_vertex` in the `boundary_vertex_to_ghost` map. +""" +function add_boundary_vertex_to_ghost!(tri::Triangulation, vertex, ghost_vertex) + get_boundary_vertex_to_ghost(tri)[vertex] = ghost_vertex + return tri +end + +""" + delete_boundary_vertex_from_ghost_map!(tri::Triangulation, vertex) + +Deletes the mapping for `vertex` from the `boundary_vertex_to_ghost` map. +""" +function delete_boundary_vertex_from_ghost_map!(tri::Triangulation, vertex) + delete!(get_boundary_vertex_to_ghost(tri), vertex) + return tri +end + """ get_incircle_cache(tri::Triangulation) -> Tuple @@ -511,10 +572,13 @@ end ) where {I,E,V,Es,Ts,B} T, adj, adj2v, graph, all_segments, boundary_edge_map, ghost_vertex_map, ghost_vertex_ranges, ch, polygon_hierarchy = __Triangulation(points, boundary_nodes, I, E, V, Es, Ts) cache = _build_cache(points, I, E, V, Es, Ts, weights, build_cache) + has_ghosts_flag = false + boundary_vertex_to_ghost_map = Dict{I,I}() tri = Triangulation( points, T, boundary_nodes, segments, all_segments, weights, adj, adj2v, graph, boundary_curves, boundary_edge_map, ghost_vertex_map, ghost_vertex_ranges, ch, representative_point_list, polygon_hierarchy, boundary_enricher, cache, + has_ghosts_flag, boundary_vertex_to_ghost_map, ) return tri end diff --git a/src/predicates/boundaries_and_ghosts.jl b/src/predicates/boundaries_and_ghosts.jl index 182736533..1d4fc267a 100644 --- a/src/predicates/boundaries_and_ghosts.jl +++ b/src/predicates/boundaries_and_ghosts.jl @@ -107,18 +107,35 @@ end Tests if the vertex `i` is a boundary node of `tri`. -# Arguments +This function uses an O(1) lookup into the `boundary_vertex_to_ghost` map when available, +providing significant performance improvement for triangulations with many boundary sections. +If the map is not yet populated (e.g., during triangulation), it falls back to O(k) iteration +over ghost vertices. + +# Arguments - `tri::Triangulation`: The [`Triangulation`](@ref). -- `i`: The vertex to test. +- `i`: The vertex to test. -# Outputs +# Outputs - `flag`: `true` if `i` is a boundary node, and `false` otherwise. - `g`: Either the ghost vertex corresponding with the section that `i` lives on if `flag` is true, or $∅ otherwise. """ function is_boundary_node(tri::Triangulation, i) - for g in each_ghost_vertex(tri) - i ∈ get_neighbours(tri, g) && return (true, g) + bv_map = get_boundary_vertex_to_ghost(tri) + # Use the map if it's populated, otherwise fall back to iteration + if !isempty(bv_map) + if haskey(bv_map, i) + return (true, bv_map[i]) + else + I = integer_type(tri) + return (false, I(∅)) + end + else + # Fallback to O(k) iteration when map is not populated + for g in each_ghost_vertex(tri) + i ∈ get_neighbours(tri, g) && return (true, g) + end + I = integer_type(tri) + return (false, I(∅)) end - I = integer_type(tri) - return (false, I(∅)) end diff --git a/test/data_structures/triangulation.jl b/test/data_structures/triangulation.jl index 4e99d23f3..d2edc459c 100644 --- a/test/data_structures/triangulation.jl +++ b/test/data_structures/triangulation.jl @@ -53,6 +53,8 @@ global tri = Triangulation(pts; IntegerType=Int32) DT.construct_polygon_hierarchy(pts; IntegerType=Int32), nothing, # boundary_enricher DT.TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), + false, # has_ghosts + Dict{Int32,Int32}(), # boundary_vertex_to_ghost ), DT.Triangulation( pts, @@ -73,6 +75,8 @@ global tri = Triangulation(pts; IntegerType=Int32) DT.construct_polygon_hierarchy(pts; IntegerType=Int32), nothing, # boundary_enricher DT.TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), + false, # has_ghosts + Dict{Int32,Int32}(), # boundary_vertex_to_ghost ), Int32[], Set{NTuple{2,Int32}}(), @@ -82,6 +86,8 @@ global tri = Triangulation(pts; IntegerType=Int32) AdaptivePredicates.orient3adapt_cache(Float64), AdaptivePredicates.insphereexact_cache(Float64), ), + false, # has_ghosts + Dict{Int32,Int32}(), # boundary_vertex_to_ghost ) end @@ -1498,7 +1504,7 @@ end for f in fieldnames(typeof(tri)) f in (:weights, :boundary_enricher) && continue @test getfield(tri, f) == getfield(tri2, f) - f == :boundary_curves && continue + f in (:boundary_curves, :has_ghosts) && continue @test !(getfield(tri, f) === getfield(tri2, f)) end end @@ -1515,7 +1521,7 @@ end for f in fieldnames(typeof(tri)) f in (:weights, :boundary_enricher, :cache) && continue @test getfield(tri, f) == getfield(tri2, f) - f == :boundary_curves && continue + f in (:boundary_curves, :has_ghosts) && continue @test !(getfield(tri, f) === getfield(tri2, f)) end bem = DT.get_boundary_edge_map(tri) @@ -1563,6 +1569,7 @@ end for f in fieldnames(typeof(tri)) f in (:weights, :cache) && continue @test getfield(tri, f) == getfield(tri2, f) + f in (:has_ghosts,) && continue @test !(getfield(tri, f) === getfield(tri2, f)) end enricher = DT.get_boundary_enricher(tri) @@ -1594,7 +1601,7 @@ end for f in fieldnames(typeof(tri)) f in (:boundary_enricher, :cache) && continue @test getfield(tri, f) == getfield(tri2, f) - f in (:boundary_curves,) && continue + f in (:boundary_curves, :has_ghosts) && continue @test !(getfield(tri, f) === getfield(tri2, f)) end @test get_weights(tri) == get_weights(tri2) && !(get_weights(tri) === get_weights(tri2)) @@ -1603,4 +1610,119 @@ end @test DT.get_weights(DT.get_triangulation(cache)) === DT.get_weights(DT.get_triangulation_2(cache)) === DT.get_weights(tri) @test DT.get_weights(DT.get_triangulation(_cache)) === DT.get_weights(DT.get_triangulation_2(_cache)) === DT.get_weights(tri2) end +end + +@testset "has_ghosts getter and setter" begin + # Test basic getter functionality + points = [(0.0, 0.0), (1.0, 0.0), (0.5, 0.5)] + tri = DT.Triangulation(points) + + # Initial state should be false + @test !DT.has_ghosts(tri) + @test DT.has_ghosts(tri) == false + + # Test setter + DT.set_has_ghosts!(tri, true) + @test DT.has_ghosts(tri) + @test DT.has_ghosts(tri) == true + + DT.set_has_ghosts!(tri, false) + @test !DT.has_ghosts(tri) + @test DT.has_ghosts(tri) == false + + # Test that setter returns the triangulation (for chaining) + result = DT.set_has_ghosts!(tri, true) + @test result === tri + @test DT.has_ghosts(tri) + + # Test has_ghost_triangles is an alias for has_ghosts + @test DT.has_ghost_triangles(tri) == DT.has_ghosts(tri) + DT.set_has_ghosts!(tri, false) + @test DT.has_ghost_triangles(tri) == DT.has_ghosts(tri) + DT.set_has_ghosts!(tri, true) + @test DT.has_ghost_triangles(tri) == DT.has_ghosts(tri) +end + +@testset "has_ghosts integration with add/delete ghost triangles" begin + points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)] + tri = DT.Triangulation(points) + + # Initial state + @test !DT.has_ghosts(tri) + + # Add boundary information first + DT.add_boundary_information!(tri) + @test !DT.has_ghosts(tri) # Still false before adding ghost triangles + + # Add ghost triangles + DT.add_ghost_triangles!(tri) + @test DT.has_ghosts(tri) + @test DT.has_ghost_triangles(tri) + + # Delete ghost triangles + DT.delete_ghost_triangles!(tri) + @test !DT.has_ghosts(tri) + @test !DT.has_ghost_triangles(tri) + + # Add them again to test toggle + DT.add_ghost_triangles!(tri) + @test DT.has_ghosts(tri) +end + +@testset "Mutating has_ghosts" begin + points = [(0.0, 0.0), (1.0, 0.0), (0.5, 0.5)] + tri = DT.Triangulation(points) + @test !DT.has_ghosts(tri) + DT.set_has_ghosts!(tri, true) + @test DT.has_ghosts(tri) + @test ismutable(tri) + @test tri.has_ghosts == DT.has_ghosts(tri) +end + +@testset "has_ghosts with copy" begin + points = [(0.0, 0.0), (1.0, 0.0), (0.5, 0.5)] + tri = DT.Triangulation(points) + + DT.set_has_ghosts!(tri, true) + @test DT.has_ghosts(tri) + + tri_copy = Base.copy(tri) + + @test DT.has_ghosts(tri_copy) == DT.has_ghosts(tri) + @test DT.has_ghosts(tri_copy) == true + + DT.set_has_ghosts!(tri, false) + @test !DT.has_ghosts(tri) + @test DT.has_ghosts(tri_copy) + + tri2 = DT.Triangulation(points) + @test !DT.has_ghosts(tri2) + tri2_copy = Base.copy(tri2) + @test !DT.has_ghosts(tri2_copy) +end + +@testset "has_ghosts with Base.==" begin + points = [(0.0, 0.0), (1.0, 0.0), (0.5, 0.5)] + tri1 = DT.Triangulation(points) + tri2 = DT.Triangulation(points) + + @test !DT.has_ghosts(tri1) + @test !DT.has_ghosts(tri2) + @test tri1 == tri2 + + DT.set_has_ghosts!(tri1, true) + @test DT.has_ghosts(tri1) + @test !DT.has_ghosts(tri2) + @test tri1 != tri2 + + DT.set_has_ghosts!(tri2, true) + @test DT.has_ghosts(tri1) + @test DT.has_ghosts(tri2) + @test tri1 == tri2 + + DT.set_has_ghosts!(tri1, false) + DT.set_has_ghosts!(tri2, false) + @test !DT.has_ghosts(tri1) + @test !DT.has_ghosts(tri2) + @test tri1 == tri2 end \ No newline at end of file diff --git a/test/helper_functions.jl b/test/helper_functions.jl index 20f10e889..aa9f9e7c9 100644 --- a/test/helper_functions.jl +++ b/test/helper_functions.jl @@ -239,6 +239,8 @@ function example_triangulation() pts, Int, NTuple{2,Int}, NTuple{3,Int}, Set{NTuple{2,Int}}, Set{NTuple{3,Int}}, DT.ZeroWeight(), Val(true), ), + false, # has_ghosts + Dict{Int,Int}(), # boundary_vertex_to_ghost ) DT.compute_representative_points!(tri) return tri @@ -264,6 +266,8 @@ function example_empty_triangulation() pts, Int, NTuple{2,Int}, NTuple{3,Int}, Set{NTuple{2,Int}}, Set{NTuple{3,Int}}, DT.ZeroWeight(), Val(true), ), + false, # has_ghosts + Dict{Int,Int}(), # boundary_vertex_to_ghost ) DT.compute_representative_points!(tri) return tri @@ -1568,6 +1572,8 @@ function why_not_equal(tri1, tri2) DT.get_polygon_hierarchy(tri1) ≠ DT.get_polygon_hierarchy(tri2) && println("get_polygon_hierarchy(tri1) ≠ get_polygon_hierarchy(tri2)") DT.get_boundary_nodes(tri1) ≠ DT.get_boundary_nodes(tri2) && println("get_boundary_nodes(tri1) ≠ get_boundary_nodes(tri2)") DT.get_weights(tri1) ≠ DT.get_weights(tri2) && println("get_weights(tri1) ≠ get_weights(tri2)") + DT.has_ghosts(tri1) ≠ DT.has_ghosts(tri2) && println("has_ghosts(tri1) ≠ has_ghosts(tri2)") + DT.get_boundary_vertex_to_ghost(tri1) ≠ DT.get_boundary_vertex_to_ghost(tri2) && println("get_boundary_vertex_to_ghost(tri1) ≠ get_boundary_vertex_to_ghost(tri2)") end #= diff --git a/test/operations/add_ghost_triangles.jl b/test/operations/add_ghost_triangles.jl index 4d57282e3..e108344be 100644 --- a/test/operations/add_ghost_triangles.jl +++ b/test/operations/add_ghost_triangles.jl @@ -6,7 +6,13 @@ using StaticArrays @testset "Adding ghost triangles" begin tri, label_map, index_map = simple_geometry() + @test !DT.has_ghosts(tri) + @test !DT.has_ghost_triangles(tri) + DT.add_ghost_triangles!(tri) + + @test DT.has_ghosts(tri) + @test DT.has_ghost_triangles(tri) outer_edges = [ ("a", "b") => DT.𝒢, ("b", "c") => DT.𝒢, diff --git a/test/operations/delete_ghost_triangles.jl b/test/operations/delete_ghost_triangles.jl index ac8890b3e..55651f71d 100644 --- a/test/operations/delete_ghost_triangles.jl +++ b/test/operations/delete_ghost_triangles.jl @@ -8,9 +8,20 @@ using StaticArrays @testset "Deleting ghost triangles" begin tri, label_map, index_map = simple_geometry() _tri = deepcopy(tri) + + @test !DT.has_ghosts(tri) + @test !DT.has_ghost_triangles(tri) + DT.add_ghost_triangles!(tri) + + @test DT.has_ghosts(tri) + @test DT.has_ghost_triangles(tri) + DT.delete_ghost_triangles!(tri) + @test !DT.has_ghosts(tri) + @test !DT.has_ghost_triangles(tri) + @test tri == _tri end diff --git a/test/operations/lock_convex_hull.jl b/test/operations/lock_convex_hull.jl index fb178d61c..aabe21aef 100644 --- a/test/operations/lock_convex_hull.jl +++ b/test/operations/lock_convex_hull.jl @@ -67,6 +67,9 @@ end @test tri.ghost_vertex_map == tri2.ghost_vertex_map @test tri.boundary_edge_map == tri2.boundary_edge_map @test tri.convex_hull == tri2.convex_hull + # Test that unlock_convex_hull! clears the boundary_vertex_to_ghost map + @test isempty(DT.get_boundary_vertex_to_ghost(tri)) + @test isempty(DT.get_boundary_vertex_to_ghost(tri2)) @test_throws ArgumentError("Cannot unlock the convex hull of a triangulation without boundary nodes.") unlock_convex_hull!(tri2) lock_convex_hull!(tri) push!(tri.convex_hull.vertices, 17) diff --git a/test/predicates/boundaries_and_ghosts.jl b/test/predicates/boundaries_and_ghosts.jl index b406e8708..81c5ece77 100644 --- a/test/predicates/boundaries_and_ghosts.jl +++ b/test/predicates/boundaries_and_ghosts.jl @@ -15,6 +15,20 @@ A = get_area(tri) refine!(tri; max_area = 1.0e-2A, rng, use_circumcenter = true) tri2, label_map, index_map = simple_geometry() +let bn = get_boundary_nodes(tri2) + bv_map = DT.get_boundary_vertex_to_ghost(tri2) + ghost = DT.𝒢 + for curve_idx in 1:DT.num_curves(bn) + curve = DT.get_boundary_nodes(bn, curve_idx) + for section_idx in 1:DT.num_sections(curve) + section = DT.get_boundary_nodes(curve, section_idx) + for node in section + bv_map[node] = ghost + end + ghost -= 1 + end + end +end add_ghost_triangles!(tri2) DT.compute_representative_points!(tri2) pts = get_points(tri2) @@ -188,10 +202,16 @@ end @testset "has_ghost_triangles" begin @test DT.has_ghost_triangles(tri) + @test DT.has_ghosts(tri) DT.delete_ghost_triangles!(tri) @test !DT.has_ghost_triangles(tri) + @test !DT.has_ghosts(tri) DT.add_ghost_triangles!(tri) @test DT.has_ghost_triangles(tri) + @test DT.has_ghosts(tri) + + # Test that has_ghost_triangles is an alias for has_ghosts + @test DT.has_ghost_triangles(tri) == DT.has_ghosts(tri) end @testset "has_boundary_nodes and is_constrained" begin @@ -240,6 +260,21 @@ end end end tri2, label_map, index_map = simple_geometry() + # Manually populate boundary_vertex_to_ghost for tri2 since simple_geometry uses delete_ghosts=true + let bn = get_boundary_nodes(tri2) + bv_map = DT.get_boundary_vertex_to_ghost(tri2) + ghost = DT.𝒢 + for curve_idx in 1:DT.num_curves(bn) + curve = DT.get_boundary_nodes(bn, curve_idx) + for section_idx in 1:DT.num_sections(curve) + section = DT.get_boundary_nodes(curve, section_idx) + for node in section + bv_map[node] = ghost + end + ghost -= 1 + end + end + end for (ghost_vertex, segment_index) in get_ghost_vertex_map(tri2) nodes = get_boundary_nodes(tri2, segment_index) for node in nodes @@ -249,9 +284,9 @@ end @test res1 ∈ DT.get_ghost_vertex_range(tri2, ghost_vertex) end end - reduced_bn = reduce(vcat, reduce(vcat, get_boundary_nodes(tri2))) - for node in each_vertex(tri) - if node ∉ reduced_bn + reduced_bn2 = reduce(vcat, reduce(vcat, get_boundary_nodes(tri2))) + for node in each_vertex(tri2) + if node ∉ reduced_bn2 flag, res = DT.is_boundary_node(tri2, node) @test !flag && res == DT.∅ end @@ -270,4 +305,100 @@ end @test !flag2 && res2 == DT.∅ end end + + @testset "boundary_vertex_to_ghost consistency" begin + # Test with complicated geometry (multiple curves and sections) + # Note: For junction nodes, the map will only + # contain one ghost vertex, whichever was set last during processing. + # So we check that each boundary node is in the map and maps to a valid + # ghost vertex (one that is in the ghost vertex range for that curve). + bv_map = DT.get_boundary_vertex_to_ghost(tri) + for (ghost_vertex, segment_index) in get_ghost_vertex_map(tri) + nodes = get_boundary_nodes(tri, segment_index) + for node in nodes + @test haskey(bv_map, node) + # The mapped ghost vertex should be in the valid range for this curve + mapped_ghost = bv_map[node] + @test mapped_ghost ∈ DT.get_ghost_vertex_range(tri, ghost_vertex) + end + end + + # Test with simple geometry + bv_map2 = DT.get_boundary_vertex_to_ghost(tri2) + for (ghost_vertex, segment_index) in get_ghost_vertex_map(tri2) + nodes = get_boundary_nodes(tri2, segment_index) + for node in nodes + @test haskey(bv_map2, node) + mapped_ghost = bv_map2[node] + @test mapped_ghost ∈ DT.get_ghost_vertex_range(tri2, ghost_vertex) + end + end + + # Test that non-boundary nodes are not in the map + for node in each_vertex(tri) + if node ∉ reduced_bn + @test !haskey(bv_map, node) + end + end + end +end + +@testset "boundary_vertex_to_ghost getter and setters" begin + tri, label_map, index_map = simple_geometry() + + bv_map = DT.get_boundary_vertex_to_ghost(tri) + @test bv_map isa Dict + @test isempty(bv_map) + + I = DT.integer_type(tri) + test_vertex = I(5) + test_ghost = I(-2) + + DT.add_boundary_vertex_to_ghost!(tri, test_vertex, test_ghost) + @test haskey(DT.get_boundary_vertex_to_ghost(tri), test_vertex) + @test DT.get_boundary_vertex_to_ghost(tri)[test_vertex] == test_ghost + + DT.delete_boundary_vertex_from_ghost_map!(tri, test_vertex) + @test !haskey(DT.get_boundary_vertex_to_ghost(tri), test_vertex) + + x, y = complicated_geometry() + rng = StableRNG(99988) + boundary_nodes, points = convert_boundary_points_to_indices(x, y) + tri_with_boundary = triangulate(points; rng, boundary_nodes, delete_ghosts = false) + + bv_map = DT.get_boundary_vertex_to_ghost(tri_with_boundary) + @test !isempty(bv_map) + + all_boundary_nodes = reduce(vcat, reduce(vcat, get_boundary_nodes(tri_with_boundary))) + for bn in all_boundary_nodes + @test haskey(bv_map, bn) + end + + for v in each_vertex(tri_with_boundary) + if v ∉ all_boundary_nodes + @test !haskey(bv_map, v) + end + end +end + +@testset "boundary_vertex_to_ghost with copy and ==" begin + x, y = complicated_geometry() + rng = StableRNG(99988) + boundary_nodes, points = convert_boundary_points_to_indices(x, y) + tri = triangulate(points; rng, boundary_nodes, delete_ghosts = false) + + tri_copy = copy(tri) + @test DT.get_boundary_vertex_to_ghost(tri) == DT.get_boundary_vertex_to_ghost(tri_copy) + @test DT.get_boundary_vertex_to_ghost(tri) !== DT.get_boundary_vertex_to_ghost(tri_copy) # Different objects + + @test tri == tri_copy + + I = DT.integer_type(tri_copy) + fake_vertex = I(999999) + fake_ghost = I(-999) + DT.add_boundary_vertex_to_ghost!(tri_copy, fake_vertex, fake_ghost) + @test tri != tri_copy + + DT.delete_boundary_vertex_from_ghost_map!(tri_copy, fake_vertex) + @test tri == tri_copy end diff --git a/test/predicates/index_and_ghost_handling.jl b/test/predicates/index_and_ghost_handling.jl index f8466919b..5abba49c9 100644 --- a/test/predicates/index_and_ghost_handling.jl +++ b/test/predicates/index_and_ghost_handling.jl @@ -19,6 +19,7 @@ for PT in subtypes(DT.AbstractPredicateKernel) _tri.all_segments, _tri.weights, _tri.adjacent, _tri.adjacent2vertex, _tri.graph, _tri.boundary_curves, _tri.boundary_edge_map, _tri.ghost_vertex_map, _tri.ghost_vertex_ranges, DT.ConvexHull(_pts, _tri.convex_hull.vertices), _tri.representative_point_list, _tri.polygon_hierarchy, _tri.boundary_enricher, _tri.cache, + _tri.has_ghosts, copy(_tri.boundary_vertex_to_ghost), ) DT.compute_representative_points!(tri) global rep = DT.get_representative_point_list(tri) diff --git a/test/runtests.jl b/test/runtests.jl index 96fcc941c..e9a2c44b9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -163,4 +163,4 @@ end end @test isempty(missing_set) end -end +end \ No newline at end of file diff --git a/test/triangulation/triangulate.jl b/test/triangulation/triangulate.jl index 426cad229..d507ce3d3 100644 --- a/test/triangulation/triangulate.jl +++ b/test/triangulation/triangulate.jl @@ -9,12 +9,12 @@ using CairoMakie for PT in subtypes(DT.AbstractPredicateKernel) for _ in 1:100 pts = rand(2, 38) - tri = triangulate(pts; predicates = PT()) - @test validate_triangulation(tri; predicates = PT()) - _tri = DT.triangulate(pts; predicates = PT()) + tri = triangulate(pts; predicates=PT()) + @test validate_triangulation(tri; predicates=PT()) + _tri = DT.triangulate(pts; predicates=PT()) @test tri == _tri - __tri = retriangulate(_tri; predicates = PT()) - @inferred retriangulate(_tri; predicates = PT()) + __tri = retriangulate(_tri; predicates=PT()) + @inferred retriangulate(_tri; predicates=PT()) @test __tri == _tri end end @@ -22,18 +22,18 @@ end @testset "Retriangulate should ignore deleted points" begin points = [(0.0, 0.0), (0.87, 0.0), (1.0006, 0.7766), (0.0, 1.0), (0.5, 0.5)] - tri = triangulate(points; skip_points = 5) + tri = triangulate(points; skip_points=5) _tri = retriangulate(tri) @test tri == _tri && !DelaunayTriangulation.has_vertex(_tri, 5) && validate_triangulation(_tri) end @testset "Lots of collinearity" begin for PT in (DT.AdaptiveKernel, DT.ExactKernel) - _tri = triangulate_rectangle(-3.0, 2.0, 5.0, 17.3, 23, 57; single_boundary = true, predicates = PT()) - @test validate_triangulation(_tri; predicates = PT()) + _tri = triangulate_rectangle(-3.0, 2.0, 5.0, 17.3, 23, 57; single_boundary=true, predicates=PT()) + @test validate_triangulation(_tri; predicates=PT()) for _ in 1:10 - tri = triangulate(_tri.points; predicates = PT()) - @test validate_triangulation(tri; predicates = PT()) + tri = triangulate(_tri.points; predicates=PT()) + @test validate_triangulation(tri; predicates=PT()) end end end \ No newline at end of file