From 55a9848852bf66d12b81010802c44bf3461e6f4f Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 15 Sep 2024 16:16:38 +0100 Subject: [PATCH 1/3] Upgrade --- NEWS.md | 4 + Project.toml | 4 +- src/DelaunayTriangulation.jl | 1 + src/algorithms/point_location/brute_force.jl | 12 +- .../basic_operations/add_point.jl | 2 +- .../basic_operations/legalise_edge.jl | 2 +- .../constrained_triangulation.jl | 2 +- .../triangulation/triangulate_convex.jl | 2 +- .../unconstrained_triangulation.jl | 6 +- .../triangulation/methods/weights.jl | 36 ++- .../triangulation/triangulation.jl | 142 +++++---- .../triangulation/triangulation_cache.jl | 52 +++- src/predicates/predicate_definitions.jl | 34 ++- src/predicates/predicates.jl | 108 +++++-- src/public.jl | 4 +- test/data_structures/triangulation.jl | 278 +++++++++--------- test/data_structures/triangulation_cache.jl | 18 +- test/helper_functions.jl | 16 + test/predicates/general.jl | 2 +- test/triangulation/weighted.jl | 54 ++-- test/utils.jl | 35 +++ 21 files changed, 531 insertions(+), 283 deletions(-) diff --git a/NEWS.md b/NEWS.md index 868363c6a..8babe42c9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # Changelog +## 1.4.0 + +- Updated to AdaptivePredicates.jl v1.2, now allowing caches to be passed to the predicates involving `incircle` and `orient3`. These are only useful when using the `AdaptiveKernel()` kernel. Outside of triangulating, these caches are not passed by default, but can be provided. The functions `get_incircle_cache` and `get_orient3_cache` can be used for this purpose on a triangulation (without a triangulation, refer to AdaptivePredicate.jl's `incircleadapt_cache` and `orient3adapt_cache`). + ## 1.3.1 - Fix an issue with a weighted triangulation where the lifted points' convex hull was entirely coplanar. See [#184](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/184) diff --git a/Project.toml b/Project.toml index b01604a23..c9b9d7807 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DelaunayTriangulation" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" authors = ["Daniel VandenHeuvel "] -version = "1.3.1" +version = "1.4.0" [deps] AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7" @@ -10,7 +10,7 @@ ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -AdaptivePredicates = "1" +AdaptivePredicates = "1.2" EnumX = "1.0" ExactPredicates = "2.2" Random = "1" diff --git a/src/DelaunayTriangulation.jl b/src/DelaunayTriangulation.jl index ab6c7c0ab..89f8c0750 100644 --- a/src/DelaunayTriangulation.jl +++ b/src/DelaunayTriangulation.jl @@ -25,6 +25,7 @@ import EnumX import Random abstract type AbstractPredicateKernel end # needs to be defined early for use in data_structures.jl +const PredicateCacheType = Union{Nothing, <:Tuple} include("data_structures.jl") include("predicates.jl") diff --git a/src/algorithms/point_location/brute_force.jl b/src/algorithms/point_location/brute_force.jl index dc0a309ab..e70d2aa33 100644 --- a/src/algorithms/point_location/brute_force.jl +++ b/src/algorithms/point_location/brute_force.jl @@ -38,16 +38,22 @@ function brute_force_search(tri::Triangulation, q; itr = each_triangle(tri), pre end """ - brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel=AdaptiveKernel()) -> Triangle + brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel=AdaptiveKernel(); cache = nothing) -> Triangle Searches for a triangle in `tri` containing the vertex `i` in its circumcircle using brute force. If `tri` is a weighted Delaunay triangulation, the triangle returned instead has the lifted vertex `i` below its witness plane. If no such triangle exists, `($∅, $∅, $∅)` is returned. You can control the method used for computing predicates via the `predicates` argument. + +The `cache` argument is passed to [`point_position_relative_to_circumcircle`] and should be one of +- `nothing`: No cache is used. +- `get_incircle_cache(tri)`: The cache stored inside `tri`. +- `AdaptivePredicates.incircleadapt_cache(number_type(tri))`: Compute a new cache. +The cache is only needed if an `AdaptiveKernel()` is used. """ -function brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel = AdaptiveKernel()) +function brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel = AdaptiveKernel(); cache::PredicateCacheType = nothing) for V in each_triangle(tri) - cert = point_position_relative_to_circumcircle(predicates, tri, V, i) + cert = point_position_relative_to_circumcircle(predicates, tri, V, i; cache) !is_outside(cert) && return V end tri_type = triangle_type(tri) diff --git a/src/algorithms/triangulation/basic_operations/add_point.jl b/src/algorithms/triangulation/basic_operations/add_point.jl index 425ba4ca1..b12cf9327 100644 --- a/src/algorithms/triangulation/basic_operations/add_point.jl +++ b/src/algorithms/triangulation/basic_operations/add_point.jl @@ -67,7 +67,7 @@ function add_point!( end q = get_point(tri, new_point) if is_weighted(tri) - cert = point_position_relative_to_circumcircle(predicates, tri, V, new_point) # redirects to point_position_relative_to_witness_plane + cert = point_position_relative_to_circumcircle(predicates, tri, V, new_point; cache = get_orient3_cache(tri)) # redirects to point_position_relative_to_witness_plane is_outside(cert) && return V # If the point is submerged, then we don't add it end flag = point_position_relative_to_triangle(predicates, tri, V, q) diff --git a/src/algorithms/triangulation/basic_operations/legalise_edge.jl b/src/algorithms/triangulation/basic_operations/legalise_edge.jl index 8a07f4f0c..4e90c3b5a 100644 --- a/src/algorithms/triangulation/basic_operations/legalise_edge.jl +++ b/src/algorithms/triangulation/basic_operations/legalise_edge.jl @@ -27,7 +27,7 @@ There is no output, as `tri` is updated in-place. so that at a triangle that might have appeared in both will only appear in one. """ function legalise_edge!(tri::Triangulation, i, j, r, store_event_history = Val(false), event_history = nothing; predicates::AbstractPredicateKernel = AdaptiveKernel()) - cert = is_legal(predicates, tri, i, j) + cert = is_legal(predicates, tri, i, j; cache = get_incircle_cache(tri)) if is_illegal(cert) e = get_adjacent(tri, j, i) flip_edge!(tri, i, j, e, r, store_event_history, event_history) diff --git a/src/algorithms/triangulation/constrained_triangulation.jl b/src/algorithms/triangulation/constrained_triangulation.jl index 7effaae74..4cc653888 100644 --- a/src/algorithms/triangulation/constrained_triangulation.jl +++ b/src/algorithms/triangulation/constrained_triangulation.jl @@ -652,7 +652,7 @@ function add_point_cavity_cdt!(tri::Triangulation, u, v, w, marked_vertices, pre insert_flag = true else p, q, r, s = get_point(tri, w, v, x, u) # Don't want to deal with boundary handling here - incircle_test = point_position_relative_to_circle(predicates, p, q, r, s) + incircle_test = point_position_relative_to_circle(predicates, p, q, r, s; cache = get_incircle_cache(tri)) orient_test = triangle_orientation(predicates, tri, u, v, w) insert_flag = !is_inside(incircle_test) && is_positively_oriented(orient_test) end diff --git a/src/algorithms/triangulation/triangulate_convex.jl b/src/algorithms/triangulation/triangulate_convex.jl index 3fa2d8fc3..97b9e48ea 100644 --- a/src/algorithms/triangulation/triangulate_convex.jl +++ b/src/algorithms/triangulation/triangulate_convex.jl @@ -121,7 +121,7 @@ to make. """ function add_point_convex_triangulation!(tri::Triangulation, u, v, w, S, predicates::AbstractPredicateKernel = AdaptiveKernel()) x = get_adjacent(tri, w, v) - if edge_exists(x) && is_inside(point_position_relative_to_circumcircle(predicates, tri, u, v, w, x)) + if edge_exists(x) && is_inside(point_position_relative_to_circumcircle(predicates, tri, u, v, w, x; cache = get_incircle_cache(tri))) # uvw and wvx are not Delaunay delete_triangle!(tri, w, v, x; protect_boundary = true, update_ghost_edges = false) add_point_convex_triangulation!(tri, u, v, x, S, predicates) diff --git a/src/algorithms/triangulation/unconstrained_triangulation.jl b/src/algorithms/triangulation/unconstrained_triangulation.jl index 7b08eb782..d9382621c 100644 --- a/src/algorithms/triangulation/unconstrained_triangulation.jl +++ b/src/algorithms/triangulation/unconstrained_triangulation.jl @@ -153,7 +153,7 @@ function add_point_bowyer_watson!(tri::Triangulation, new_point, initial_search_ q = get_point(tri, _new_point) V = find_triangle(tri, q; predicates, m = nothing, point_indices = nothing, try_points = nothing, k = initial_search_point, rng) if is_weighted(tri) - cert = point_position_relative_to_circumcircle(predicates, tri, V, _new_point) # redirects to point_position_relative_to_witness_plane + cert = point_position_relative_to_circumcircle(predicates, tri, V, _new_point; cache = get_orient3_cache(tri)) # redirects to point_position_relative_to_witness_plane is_outside(cert) && return V # If the point is submerged, then we don't add it end flag = point_position_relative_to_triangle(predicates, tri, V, q) @@ -298,9 +298,9 @@ Determines whether to enter the cavity in `tri` through the edge `(i, j)` when i function enter_cavity(tri::Triangulation, r, i, j, ℓ, predicates::AbstractPredicateKernel = AdaptiveKernel()) contains_segment(tri, i, j) && return false if is_ghost_vertex(ℓ) - cert = point_position_relative_to_circumcircle(tri, j, i, ℓ, r) + cert = point_position_relative_to_circumcircle(predicates, tri, j, i, ℓ, r; cache = get_incircle_cache(tri)) else - cert = point_position_relative_to_circumcircle(tri, r, i, j, ℓ) + cert = point_position_relative_to_circumcircle(predicates, tri, r, i, j, ℓ; cache = get_incircle_cache(tri)) end return is_weighted(tri) ? !is_outside(cert) : is_inside(cert) end diff --git a/src/data_structures/triangulation/methods/weights.jl b/src/data_structures/triangulation/methods/weights.jl index cb418fef3..38909e14e 100644 --- a/src/data_structures/triangulation/methods/weights.jl +++ b/src/data_structures/triangulation/methods/weights.jl @@ -98,7 +98,7 @@ function get_power_distance(tri::Triangulation, i, j) end """" - get_distance_to_witness_plane(tri::Triangulation, i, V) -> Number + get_distance_to_witness_plane([kernel::AbstractPredicateKernel = AdaptiveKernel(), ] tri::Triangulation, i, V; cache = nothing) -> Number Computes the distance between the lifted companion of the vertex `i` and the witness plane to the triangle `V`. If `V` is a ghost triangle and `i` is not on its solid edge, then the distance is `-Inf` if it is below the ghost triangle's witness plane and `Inf` if it is above. If `V` is a ghost triangle and `i` @@ -106,10 +106,16 @@ is on its solid edge, then the distance returned is the distance associated with In general, the distance is positive if the lifted vertex is above the witness plane, negative if it is below, and zero if it is on the plane. + +The `kernel` argument determines how this result is computed, and should be +one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). +See the documentation for more information about these choices. + +The `cache` keyword argument is passed to [`point_position_relative_to_circumcircle`]. Please see the documentation for that function for more information. See also [`point_position_relative_to_witness_plane`](@ref) and [`get_distance_to_plane`](@ref). """ -function get_distance_to_witness_plane(tri::Triangulation, i, V) +function get_distance_to_witness_plane(kernel::AbstractPredicateKernel, tri::Triangulation, i, V; cache::PredicateCacheType = nothing) if !is_ghost_triangle(V) u, v, w = triangle_vertices(V) p⁺ = get_lifted_point(tri, u) @@ -118,17 +124,19 @@ function get_distance_to_witness_plane(tri::Triangulation, i, V) s⁺ = get_lifted_point(tri, i) return get_distance_to_plane(p⁺, q⁺, r⁺, s⁺) else - cert = point_position_relative_to_circumcircle(tri, V, i) + cert = point_position_relative_to_circumcircle(kernel, tri, V, i; cache) if is_inside(cert) return -Inf elseif is_outside(cert) return Inf else # is_on(cert) V′ = replace_ghost_triangle_with_boundary_triangle(tri, V) - return get_distance_to_witness_plane(tri, i, V′) + return get_distance_to_witness_plane(kernel, tri, i, V′; cache) end end end +get_distance_to_witness_plane(tri::Triangulation, i, V; cache::PredicateCacheType = nothing) = get_distance_to_witness_plane(AdaptiveKernel(), tri, i, V; cache) + """ get_weighted_nearest_neighbour(tri::Triangulation, i, j = rand(each_solid_vertex(tri))) -> Vertex @@ -164,24 +172,32 @@ function _get_weighted_nearest_neighbour(tri, i, j) end @doc """ - is_submerged(tri::Triangulation, i) -> Bool - is_submerged(tri::Triangulation, i, V) -> Bool + is_submerged([kernel::AbstractPredicateKernel = AdaptiveKernel(), ] tri::Triangulation, i; cache = nothing) -> Bool + is_submerged([kernel::AbstractPredicateKernel = AdaptiveKernel(), ] tri::Triangulation, i, V; cache = nothing) -> Bool Returns `true` if the vertex `i` is submerged in `tri` and `false` otherwise. In the second method, `V` is a triangle containing `tri`. + +The `kernel` argument determines how this result is computed, and should be +one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). +See the documentation for more information about these choices. + +The `cache` keyword argument is passed to [`point_position_relative_to_circumcircle`]. Please see the documentation for that function for more information. """ is_submerged -function is_submerged(tri::Triangulation, i) +function is_submerged(kernel::AbstractPredicateKernel, tri::Triangulation, i; cache::PredicateCacheType = nothing) # A source that mentions that testing if `i` is submerged only needs to consider the triangle that contains it # is given in https://otik.uk.zcu.cz/bitstream/11025/21574/1/Zemek.pdf on p.17. # (If the link dies, it is the PhD thesis of Michal Zemek, "Regular Triangulation in 3D and Its Applications".) is_ghost_vertex(i) && return false q = get_point(tri, i) V = find_triangle(tri, q) - return is_submerged(tri, i, V) + return is_submerged(kernel, tri, i, V; cache) end -function is_submerged(tri::Triangulation, i, V) +function is_submerged(kernel::AbstractPredicateKernel, tri::Triangulation, i, V; cache::PredicateCacheType = nothing) is_ghost_vertex(i) && return false - cert = point_position_relative_to_circumcircle(tri, V, i) + cert = point_position_relative_to_circumcircle(kernel, tri, V, i; cache) return is_outside(cert) end +is_submerged(tri::Triangulation, i; cache::PredicateCacheType = nothing) = is_submerged(AdaptiveKernel(), tri, i; cache) +is_submerged(tri::Triangulation, i, V; cache::PredicateCacheType = nothing) = is_submerged(AdaptiveKernel(), tri, i, V; cache) \ No newline at end of file diff --git a/src/data_structures/triangulation/triangulation.jl b/src/data_structures/triangulation/triangulation.jl index 35a48aa70..1bb1d21a2 100644 --- a/src/data_structures/triangulation/triangulation.jl +++ b/src/data_structures/triangulation/triangulation.jl @@ -84,7 +84,7 @@ 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). """ -struct Triangulation{P, T, BN, W, I, E, Es, BC, BEM, GVM, GVR, BPL, C, BE} +struct Triangulation{P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,BE} # Geometry points::P triangles::T @@ -93,8 +93,8 @@ struct Triangulation{P, T, BN, W, I, E, Es, BC, BEM, GVM, GVR, BPL, C, BE} all_segments::Es weights::W # Topology - adjacent::Adjacent{I, E} - adjacent2vertex::Adjacent2Vertex{I, Es} + adjacent::Adjacent{I,E} + adjacent2vertex::Adjacent2Vertex{I,Es} graph::Graph{I} # Boundary handling boundary_curves::BC @@ -102,7 +102,7 @@ struct Triangulation{P, T, BN, W, I, E, Es, BC, BEM, GVM, GVR, BPL, C, BE} ghost_vertex_map::GVM ghost_vertex_ranges::GVR # Other - convex_hull::ConvexHull{P, I} + convex_hull::ConvexHull{P,I} representative_point_list::BPL polygon_hierarchy::PolygonHierarchy{I} boundary_enricher::BE @@ -320,6 +320,27 @@ Returns the cache of `tri`. This is a [`TriangulationCache`](@ref) used as a cac """ get_cache(tri::Triangulation) = tri.cache +""" + get_incircle_cache(tri::Triangulation) -> Tuple + +Returns the incircle cache stored in `tri`. +""" +get_incircle_cache(tri::Triangulation) = get_incircle_cache(get_cache(tri)) + +""" + get_orient3_cache(tri::Triangulation) -> Tuple + +Returns the orient3 cache stored in `tri`. +""" +get_orient3_cache(tri::Triangulation) = get_orient3_cache(get_cache(tri)) + +""" + get_insphere_cache(tri::Triangulation) -> Tuple + +Returns the insphere cache stored in `tri`. +""" +get_insphere_cache(tri::Triangulation) = get_insphere_cache(get_cache(tri)) + """ number_type(tri::Triangulation) -> DataType @@ -332,28 +353,28 @@ number_type(::Triangulation{P}) where {P} = number_type(P) Returns the type used for representing vertices in `tri`. """ -integer_type(::Triangulation{P, T, BN, W, I}) where {P, T, BN, W, I} = I +integer_type(::Triangulation{P,T,BN,W,I}) where {P,T,BN,W,I} = I """ edges_type(tri::Triangulation) -> DataType Returns the type used for representing collections of edges in `tri`. """ -edges_type(::Triangulation{P, T, BN, W, I, E, Es}) where {P, T, BN, W, I, E, Es} = Es +edges_type(::Triangulation{P,T,BN,W,I,E,Es}) where {P,T,BN,W,I,E,Es} = Es """ edge_type(tri::Triangulation) -> DataType Returns the type used for representing individual edges in `tri`. """ -edge_type(::Triangulation{P, T, BN, W, I, E}) where {P, T, BN, W, I, E} = E +edge_type(::Triangulation{P,T,BN,W,I,E}) where {P,T,BN,W,I,E} = E """ triangle_type(tri::Triangulation) -> DataType Returns the type used for representing collections of triangles in `tri`. """ -triangles_type(::Triangulation{P, T}) where {P, T} = T +triangles_type(::Triangulation{P,T}) where {P,T} = T """ triangle_type(tri::Triangulation) -> DataType @@ -363,12 +384,12 @@ Returns the type used for representing individual triangles in `tri`. triangle_type(tri::Triangulation) = triangle_type(triangles_type(tri)) @inline function __Triangulation( - points::P, boundary_nodes, IntegerType::Type{I}, EdgeType::Type{E}, - TriangleType::Type{V}, EdgesType::Type{Es}, TrianglesType::Type{Ts}, - ) where {P, I, E, V, Es, Ts} + points::P, boundary_nodes, IntegerType::Type{I}, EdgeType::Type{E}, + TriangleType::Type{V}, EdgesType::Type{Es}, TrianglesType::Type{Ts}, +) where {P,I,E,V,Es,Ts} T = TrianglesType() - adj = Adjacent{IntegerType, EdgeType}() - adj2v = Adjacent2Vertex{IntegerType, EdgesType}() + adj = Adjacent{IntegerType,EdgeType}() + adj2v = Adjacent2Vertex{IntegerType,EdgesType}() graph = Graph{IntegerType}() all_segments = EdgesType() boundary_edge_map = construct_boundary_edge_map(boundary_nodes, IntegerType, EdgeType) @@ -386,26 +407,33 @@ triangle_type(tri::Triangulation) = triangle_type(triangles_type(tri)) end @inline function _build_cache( - points::P, IntegerType::Type{I}, EdgeType::Type{E}, - TriangleType::Type{V}, EdgesType::Type{Es}, TrianglesType::Type{T}, weights::W, build_cache::Val{B}, - ) where {P, I, E, V, Es, T, W, B} + points::P, IntegerType::Type{I}, EdgeType::Type{E}, + TriangleType::Type{V}, EdgesType::Type{Es}, TrianglesType::Type{T}, weights::W, build_cache::Val{B}, +) where {P,I,E,V,Es,T,W,B} if B + NT = number_type(points) BN = Vector{IntegerType} BC = Tuple{} - BEM = Dict{EdgeType, Tuple{Vector{IntegerType}, IntegerType}} - GVM = Dict{IntegerType, Vector{IntegerType}} - GVR = Dict{IntegerType, UnitRange{IntegerType}} - BPL = Dict{IntegerType, RepresentativeCoordinates{IntegerType, number_type(points)}} - C = TriangulationCache{Nothing, Nothing, Nothing, Nothing, Nothing} + BEM = Dict{EdgeType,Tuple{Vector{IntegerType},IntegerType}} + GVM = Dict{IntegerType,Vector{IntegerType}} + GVR = Dict{IntegerType,UnitRange{IntegerType}} + BPL = Dict{IntegerType,RepresentativeCoordinates{IntegerType,NT}} + incircle_cache = AP.incircleadapt_cache(NT) + orient3_cache = AP.orient3adapt_cache(NT) + insphere_cache = AP.insphereexact_cache(NT) + IC = typeof(incircle_cache) + OC = typeof(orient3_cache) + IS = typeof(insphere_cache) + C = EmptyTriangulationCache cache = TriangulationCache( - Triangulation(points; weights, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType, build_cache = Val(false)), - Triangulation(points; weights, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType, build_cache = Val(false)), - I[], Es(), I[], T(), + Triangulation(points; weights, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType, build_cache=Val(false)), + Triangulation(points; weights, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType, build_cache=Val(false)), + I[], Es(), I[], T(), incircle_cache, orient3_cache, insphere_cache )::TriangulationCache{ - Triangulation{P, T, BN, W, I, E, Es, BC, BEM, GVM, GVR, BPL, C, Nothing}, Vector{I}, Es, Vector{I}, T, + Triangulation{P,T,BN,W,I,E,Es,BC,BEM,GVM,GVR,BPL,C,Nothing},Vector{I},Es,Vector{I},T,IC,OC,IS } else - cache = TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing)::TriangulationCache{Nothing, Nothing, Nothing, Nothing, Nothing} + cache = TriangulationCache()::EmptyTriangulationCache end return cache end @@ -417,20 +445,20 @@ Initialises an empty `Triangulation` for triangulating `points`. The keyword arg `kwargs...` match those of [`triangulate`](@ref). """ @inline function Triangulation( - points::P; - IntegerType::Type{I} = Int, - EdgeType::Type{E} = NTuple{2, IntegerType}, - TriangleType::Type{T} = NTuple{3, IntegerType}, - EdgesType::Type{Es} = Set{EdgeType}, - TrianglesType::Type{Vs} = Set{TriangleType}, - boundary_nodes = IntegerType[], - segments = EdgesType(), - weights = ZeroWeight(), - representative_point_list = Dict{IntegerType, RepresentativeCoordinates{IntegerType, number_type(points)}}(), - boundary_curves = (), - build_cache = true, - boundary_enricher = nothing, - ) where {P, I, E, T, Es, Vs} + points::P; + IntegerType::Type{I}=Int, + EdgeType::Type{E}=NTuple{2,IntegerType}, + TriangleType::Type{T}=NTuple{3,IntegerType}, + EdgesType::Type{Es}=Set{EdgeType}, + TrianglesType::Type{Vs}=Set{TriangleType}, + boundary_nodes=IntegerType[], + segments=EdgesType(), + weights=ZeroWeight(), + representative_point_list=Dict{IntegerType,RepresentativeCoordinates{IntegerType,number_type(points)}}(), + boundary_curves=(), + build_cache=true, + boundary_enricher=nothing, +) where {P,I,E,T,Es,Vs} return _Triangulation( points, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType, @@ -439,10 +467,10 @@ Initialises an empty `Triangulation` for triangulating `points`. The keyword arg ) end @inline function _Triangulation( - points, ::Type{I}, ::Type{E}, ::Type{V}, ::Type{Es}, ::Type{Ts}, - boundary_nodes, segments, weights, representative_point_list, - boundary_curves, build_cache::Val{B}, boundary_enricher, - ) where {I, E, V, Es, Ts, B} + points, ::Type{I}, ::Type{E}, ::Type{V}, ::Type{Es}, ::Type{Ts}, + boundary_nodes, segments, weights, representative_point_list, + boundary_curves, build_cache::Val{B}, boundary_enricher, +) 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) tri = Triangulation( @@ -477,18 +505,18 @@ Returns the `Triangulation` corresponding to the triangulation of `points` with - `tri`: The [`Triangulation`](@ref). """ @inline function Triangulation( - points::P, triangles::T, boundary_nodes::BN; - IntegerType::Type{I} = Int, - EdgeType::Type{E} = NTuple{2, IntegerType}, - TriangleType::Type{V} = NTuple{3, IntegerType}, - EdgesType::Type{Es} = Set{EdgeType}, - TrianglesType::Type{Ts} = Set{TriangleType}, - weights = ZeroWeight(), - delete_ghosts = false, - predicates::AbstractPredicateKernel = AdaptiveKernel(), - ) where {P, T, BN, I, E, V, Es, Ts} + points::P, triangles::T, boundary_nodes::BN; + IntegerType::Type{I}=Int, + EdgeType::Type{E}=NTuple{2,IntegerType}, + TriangleType::Type{V}=NTuple{3,IntegerType}, + EdgesType::Type{Es}=Set{EdgeType}, + TrianglesType::Type{Ts}=Set{TriangleType}, + weights=ZeroWeight(), + delete_ghosts=false, + predicates::AbstractPredicateKernel=AdaptiveKernel(), +) where {P,T,BN,I,E,V,Es,Ts} _bn = copy(boundary_nodes) - tri = Triangulation(points; boundary_nodes = _bn, weights, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType) + tri = Triangulation(points; boundary_nodes=_bn, weights, IntegerType, EdgeType, TriangleType, EdgesType, TrianglesType) return build_triangulation_from_data!(tri, triangles, _bn, delete_ghosts, predicates) end @@ -502,7 +530,7 @@ The `kernel` argument determines how predicates are computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. """ -@inline function build_triangulation_from_data!(tri::Triangulation, triangles, boundary_nodes, delete_ghosts, predicates::AbstractPredicateKernel = AdaptiveKernel()) +@inline function build_triangulation_from_data!(tri::Triangulation, triangles, boundary_nodes, delete_ghosts, predicates::AbstractPredicateKernel=AdaptiveKernel()) Es = edges_type(tri) points = get_points(tri) polygon_hierarchy = get_polygon_hierarchy(tri) @@ -519,7 +547,7 @@ See the documentation for more information about these choices. end add_boundary_information!(tri) add_ghost_triangles!(tri) - convex_hull!(tri; predicates, reconstruct = true) + convex_hull!(tri; predicates, reconstruct=true) segments = get_all_segments(tri) ghost_vertex_map = get_ghost_vertex_map(tri) all_segments = merge_segments(ghost_vertex_map, boundary_nodes, Es()) diff --git a/src/data_structures/triangulation/triangulation_cache.jl b/src/data_structures/triangulation/triangulation_cache.jl index 3917d2a07..bcb9f2737 100644 --- a/src/data_structures/triangulation/triangulation_cache.jl +++ b/src/data_structures/triangulation/triangulation_cache.jl @@ -1,5 +1,5 @@ """ - TriangulationCache{T,M,I,S} + TriangulationCache{T,M,I,S,IC,OC,IS} A cache to be used as a field in [`Triangulation`](@ref). @@ -14,6 +14,9 @@ A cache to be used as a field in [`Triangulation`](@ref). interior segments. - `surrounding_polygon::S`: The polygon surrounding the triangulation. This is needed for [`delete_point!`](@ref). - `fan_triangles::F`: Triangles in a fan. This is needed for sorting fans for constrained triangulations. +- `incircle_cache::IC`: Cache for incircle tests. +- `orient3_cache::OC`: Cache for orient3 tests. +- `insphere_cache::IS`: Cache for insphere tests. !!! note "Caches of caches" @@ -24,13 +27,16 @@ A cache to be used as a field in [`Triangulation`](@ref). The `points` of the cache's `triangulation` will be aliased to the `points` of the parent triangulation. """ -struct TriangulationCache{T, M, I, S, F} +struct TriangulationCache{T, M, I, S, F, IC, OC, IS} triangulation::T triangulation_2::T marked_vertices::M interior_segments_on_hull::I surrounding_polygon::S - fan_triangles::F + fan_triangles::F + incircle_cache::IC + orient3_cache::OC + insphere_cache::IS end function Base.show(io::IO, ::MIME"text/plain", cache::TriangulationCache) if isnothing(get_triangulation(cache)) # all fields are nothing, or none are @@ -40,6 +46,12 @@ function Base.show(io::IO, ::MIME"text/plain", cache::TriangulationCache) end end +const EmptyTriangulationCache = TriangulationCache{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing} +TriangulationCache{Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}() = TriangulationCache() +function TriangulationCache() + return TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing) +end + """ get_triangulation(cache::TriangulationCache) -> Triangulation @@ -82,6 +94,27 @@ Returns the triangles in a fan stored in `cache`. """ get_fan_triangles(cache::TriangulationCache) = cache.fan_triangles +""" + get_incircle_cache(cache::TriangulationCache) -> Tuple + +Returns the incircle cache stored in `cache`. +""" +get_incircle_cache(cache::TriangulationCache) = cache.incircle_cache + +""" + get_orient3_cache(cache::TriangulationCache) -> Tuple + +Returns the orient3 cache stored in `cache`. +""" +get_orient3_cache(cache::TriangulationCache) = cache.orient3_cache + +""" + get_insphere_cache(cache::TriangulationCache) -> Tuple + +Returns the insphere cache stored in `cache`. +""" +get_insphere_cache(cache::TriangulationCache) = cache.insphere_cache + """ empty!(cache::TriangulationCache) @@ -100,7 +133,18 @@ function Base.empty!(cache::TriangulationCache) empty!(surrounding_polygon) fan_triangles = get_fan_triangles(cache) empty!(fan_triangles) - return cache + #= + incircle_cache = get_incircle_cache(cache) + incircle_parent = parent(incircle_cache[1]) + fill!(zero(eltype(incircle_parent)), incircle_parent) + orient3_cache = get_orient3_cache(cache) + orient3_parent = parent(orient3_cache[1]) + fill!(zero(eltype(orient3_parent)), orient3_parent) + insphere_cache = get_insphere_cache(cache) + insphere_parent = parent(insphere_cache[1]) + fill!(zero(eltype(insphere_parent)), insphere_parent) + =# # No point clearing these caches. + return cache end """ diff --git a/src/predicates/predicate_definitions.jl b/src/predicates/predicate_definitions.jl index 0e97343ec..0f1c35d15 100644 --- a/src/predicates/predicate_definitions.jl +++ b/src/predicates/predicate_definitions.jl @@ -96,7 +96,7 @@ end @inline orient(::AdaptiveKernel, p, q, r) = AP.orient2p(p, q, r) """ - orient_predicate([kernel::AbstractPredicateKernel,] p, q, r) -> Integer + orient_predicate([kernel::AbstractPredicateKernel,] p, q, r; cache = nothing) -> Integer Returns @@ -126,18 +126,21 @@ Here, a positively oriented tetrahedron `(p, q, r, s)` takes the form The `kernel` argument determines how this result is computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. + +The optional `cache` keyword argument can be used for preallocating memory for intermediate results, passing the argument from `AdaptivePredicates.orient3adapt_cache(T)`, +where `T` is the number type of the input points. If `nothing` is passed, no cache is used. This is only needed if an `AdaptiveKernel()` is used. """ -@inline function orient_predicate(kernel::AbstractPredicateKernel, p, q, r, s) - return orient(kernel, getxyz(p), getxyz(q), getxyz(r), getxyz(s)) +@inline function orient_predicate(kernel::AbstractPredicateKernel, p, q, r, s; cache::PredicateCacheType = nothing) + return orient(kernel, getxyz(p), getxyz(q), getxyz(r), getxyz(s), cache) end -@inline orient_predicate(p, q, r, s) = orient_predicate(AdaptiveKernel(), p, q, r, s) +@inline orient_predicate(p, q, r, s; cache::PredicateCacheType = nothing) = orient_predicate(AdaptiveKernel(), p, q, r, s; cache) -@inline orient(::FastKernel, p, q, r, s) = sgn(AP.orient3fast(p, q, r, s)) -@inline orient(::ExactKernel, p, q, r, s) = EP.orient(_getxyz(p), _getxyz(q), _getxyz(r), _getxyz(s)) -@inline orient(::AdaptiveKernel, p, q, r, s) = AP.orient3p(p, q, r, s) +@inline orient(::FastKernel, p, q, r, s, _::PredicateCacheType = nothing) = sgn(AP.orient3fast(p, q, r, s)) +@inline orient(::ExactKernel, p, q, r, s, _::PredicateCacheType = nothing) = EP.orient(_getxyz(p), _getxyz(q), _getxyz(r), _getxyz(s)) +@inline orient(::AdaptiveKernel, p, q, r, s, cache::PredicateCacheType = nothing) = AP.orient3p(p, q, r, s, cache) """ - incircle_predicate([kernel::AbstractPredicateKernel,] a, b, c, p) -> Integer + incircle_predicate([kernel::AbstractPredicateKernel,] a, b, c, p; cache = nothing) -> Integer Assuming that `(a, b, c)` is a positively oriented triangle, returns @@ -148,15 +151,18 @@ Assuming that `(a, b, c)` is a positively oriented triangle, returns The `kernel` argument determines how this result is computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. + +The optional `cache` keyword argument can be used for preallocating memory for intermediate results, passing the argument from `AdaptivePredicates.incircleadapt_cache(T)`, +where `T` is the number type of the input points. If `nothing` is passed, no cache is used. This is only needed if an `AdaptiveKernel()` is used. """ -@inline function incircle_predicate(kernel::AbstractPredicateKernel, a, b, c, p) - return incircle(kernel, getxy(a), getxy(b), getxy(c), getxy(p)) +@inline function incircle_predicate(kernel::AbstractPredicateKernel, a, b, c, p; cache::PredicateCacheType = nothing) + return incircle(kernel, getxy(a), getxy(b), getxy(c), getxy(p), cache) end -@inline incircle_predicate(a, b, c, p) = incircle_predicate(AdaptiveKernel(), a, b, c, p) +@inline incircle_predicate(a, b, c, p; cache::PredicateCacheType = nothing) = incircle_predicate(AdaptiveKernel(), a, b, c, p; cache) -@inline incircle(::FastKernel, a, b, c, p) = sgn(AP.incirclefast(a, b, c, p)) -@inline incircle(::ExactKernel, a, b, c, p) = EP.incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p)) -@inline incircle(::AdaptiveKernel, a, b, c, p) = AP.incirclep(a, b, c, p) +@inline incircle(::FastKernel, a, b, c, p, _::PredicateCacheType = nothing) = sgn(AP.incirclefast(a, b, c, p)) +@inline incircle(::ExactKernel, a, b, c, p, _::PredicateCacheType = nothing) = EP.incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p)) +@inline incircle(::AdaptiveKernel, a, b, c, p, cache::PredicateCacheType = nothing) = AP.incirclep(a, b, c, p, cache) """ parallelorder_predicate([kernel::AbstractPredicateKernel,] a, b, p, q) -> Integer diff --git a/src/predicates/predicates.jl b/src/predicates/predicates.jl index 1bc4f8c8b..7c7b497cf 100644 --- a/src/predicates/predicates.jl +++ b/src/predicates/predicates.jl @@ -34,7 +34,7 @@ triangle_orientation(kernel::AbstractPredicateKernel, tri::Triangulation, T) = t triangle_orientation(tri::Triangulation, T) = triangle_orientation(AdaptiveKernel(), tri, T) """ - point_position_relative_to_circle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] a, b, c, p) -> Certificate + point_position_relative_to_circle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] a, b, c, p; cache = nothing) -> Certificate Given a circle through the coordinates `(a, b, c)`, assumed to be positively oriented, computes the position of `p` relative to the circle. Returns a [`Certificate`](@ref), which is one of: @@ -46,12 +46,17 @@ computes the position of `p` relative to the circle. Returns a [`Certificate`](@ The `kernel` argument determines how this result is computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. + +The `cache` keyword argument is passed to [`incircle_predicate`] and should be one of +- `nothing`: No cache is used. +- `AdaptivePredicates.incircleadapt_cache(T)`, where `T` is the number type used (`Float64` or `Float32`). +The cache is only needed if an `AdaptiveKernel()` is used. """ -function point_position_relative_to_circle(kernel::AbstractPredicateKernel, a, b, c, p) - cert = incircle_predicate(kernel, a, b, c, p) +function point_position_relative_to_circle(kernel::AbstractPredicateKernel, a, b, c, p; cache::PredicateCacheType = nothing) + cert = incircle_predicate(kernel, a, b, c, p; cache) return convert_certificate(cert, Cert.Outside, Cert.On, Cert.Inside) end -point_position_relative_to_circle(a, b, c, p) = point_position_relative_to_circle(AdaptiveKernel(), a, b, c, p) +point_position_relative_to_circle(a, b, c, p; cache::PredicateCacheType = nothing) = point_position_relative_to_circle(AdaptiveKernel(), a, b, c, p; cache) @doc """ point_position_relative_to_line([kernel::AbstractPredicateKernel=AdaptiveKernel(),] a, b, p) -> Certificate @@ -307,8 +312,8 @@ end point_position_relative_to_oriented_outer_halfplane(a, b, p) = point_position_relative_to_oriented_outer_halfplane(AdaptiveKernel(), a, b, p) @doc """ - is_legal([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, i, j) -> Certificate - is_legal([kernel::AbstractPredicateKernel=AdaptiveKernel(),] p, q, r, s) -> Certificate + is_legal([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, i, j; cache = nothing) -> Certificate + is_legal([kernel::AbstractPredicateKernel=AdaptiveKernel(),] p, q, r, s[, cache = nothing]) -> Certificate Tests if the edge `(p, q)` (or the edge `(i, j)` of `tri`) is legal, where the edge `(p, q)` is incident to two triangles `(p, q, r)` and `(q, p, s)`. In partiuclar, tests that `s` is not inside @@ -322,19 +327,24 @@ If the edge `(i, j)` is a segment of `tri` or is a ghost edge, then the edge is The `kernel` argument determines how this result is computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. + +The `cache` keyword argument is passed to [`point_position_relative_to_circumcircle`] and should be one of +- `nothing`: No cache is used. +- `AdaptivePredicates.incircleadapt_cache(T)`, where `T` is the number type used (`Float64` or `Float32`). +The cache is only needed if an `AdaptiveKernel()` is used. """ is_legal -function is_legal(kernel::AbstractPredicateKernel, p, q, r, s) - incirc = point_position_relative_to_circle(kernel, p, q, r, s) +function is_legal(kernel::AbstractPredicateKernel, p, q, r, s; cache::PredicateCacheType = nothing) + incirc = point_position_relative_to_circle(kernel, p, q, r, s; cache) if is_inside(incirc) return Cert.Illegal else return Cert.Legal end end -is_legal(p, q, r, s) = is_legal(AdaptiveKernel(), p, q, r, s) +is_legal(p, q, r, s; cache::PredicateCacheType = nothing) = is_legal(AdaptiveKernel(), p, q, r, s; cache) -function is_legal(kernel::AbstractPredicateKernel, tri::Triangulation, i, j) +function is_legal(kernel::AbstractPredicateKernel, tri::Triangulation, i, j; cache::PredicateCacheType = nothing) if contains_segment(tri, i, j) || is_boundary_edge(tri, j, i) || is_boundary_edge(tri, i, j) || !edge_exists(tri, i, j) || !edge_exists(tri, j, i) || @@ -344,11 +354,11 @@ function is_legal(kernel::AbstractPredicateKernel, tri::Triangulation, i, j) k = get_adjacent(tri, i, j) ℓ = get_adjacent(tri, j, i) p, q, r, s = get_point(tri, i, j, k, ℓ) - cert = is_legal(kernel, p, q, r, s) + cert = is_legal(kernel, p, q, r, s; cache) return cert end end -is_legal(tri::Triangulation, i, j) = is_legal(AdaptiveKernel(), tri, i, j) +is_legal(tri::Triangulation, i, j; cache::PredicateCacheType = nothing) = is_legal(AdaptiveKernel(), tri, i, j; cache) @doc """ triangle_line_segment_intersection([kernel::AbstractPredicateKernel=AdaptiveKernel(),] p, q, r, a, b) -> Certificate @@ -533,8 +543,8 @@ end find_edge(tri::Triangulation, T, ℓ) = find_edge(AdaptiveKernel(), tri, T, ℓ) @doc """ - point_position_relative_to_circumcircle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, i, j, k, ℓ) -> Certificate - point_position_relative_to_circumcircle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, T, ℓ) -> Certificate + point_position_relative_to_circumcircle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, i, j, k, ℓ; cache = nothing) -> Certificate + point_position_relative_to_circumcircle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, T, ℓ; cache = nothing) -> Certificate Tests the position of the vertex `ℓ` of `tri` relative to the circumcircle of the triangle `T = (i, j, k)`. The returned value is a [`Certificate`](@ref), which is one of: @@ -546,6 +556,12 @@ The `kernel` argument determines how this result is computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. +The `cache` keyword argument is useful when an `AdaptiveKernel()` is used. Depending on whether the triangulation is weighted or not, the cache should be of the following: +- `nothing`: No cache is used. +- `AdaptivePredicates.incircleadapt_cache(T)`, where `T` is the number type used (`Float64` or `Float32`). This is used for unweighted triangulations. +- `AdaptivePredicates.orient3adapt_cache(T)`, where `T` is the number type used (`Float64` or `Float32`). This is used for weighted triangulations. +In case the wrong cache type is used, it is replaced with `nothing`. + !!! note "Ghost triangles" The circumcircle of a ghost triangle is defined as the oriented outer halfplane of the solid edge of the triangle. See [`point_position_relative_to_oriented_outer_halfplane`](@ref). @@ -557,7 +573,7 @@ See the documentation for more information about these choices. in addition to checking [`point_position_relative_to_oriented_outer_halfplane`](@ref), we must check if the new vertex is not submerged by the adjoining solid triangle. """ point_position_relative_to_circumcircle -function point_position_relative_to_circumcircle(kernel::AbstractPredicateKernel, tri::Triangulation, i, j, k, ℓ) +function point_position_relative_to_circumcircle(kernel::AbstractPredicateKernel, tri::Triangulation, i, j, k, ℓ; cache::PredicateCacheType = nothing) u, v, w = sort_triangle(i, j, k) a, b, c, p = get_point(tri, u, v, w, ℓ) if is_ghost_vertex(w) @@ -565,27 +581,65 @@ function point_position_relative_to_circumcircle(kernel::AbstractPredicateKernel if is_on(cert) && is_weighted(tri) u′, v′, w′ = replace_ghost_triangle_with_boundary_triangle(tri, (u, v, w)) !edge_exists(w′) && return Cert.Inside # needed for the case tri = triangulate(get_points(triangulate_rectangle(0, 10, 0, 10, 3, 3)), weights = zeros(9), insertion_order = [9, 7, 6, 8, 4, 5, 3, 1, 2]) - sub_cert = point_position_relative_to_witness_plane(kernel, tri, u′, v′, w′, ℓ) + sub_cert = point_position_relative_to_witness_plane(kernel, tri, u′, v′, w′, ℓ; cache = fix_orient3_cache(tri, cache)) is_above(sub_cert) && return Cert.Outside return cert else return cert end elseif !is_weighted(tri) - return point_position_relative_to_circle(kernel, a, b, c, p) + return point_position_relative_to_circle(kernel, a, b, c, p; cache = fix_incircle_cache(tri, cache)) else - cert = point_position_relative_to_witness_plane(kernel, tri, i, j, k, ℓ) + cert = point_position_relative_to_witness_plane(kernel, tri, i, j, k, ℓ; cache = fix_orient3_cache(tri, cache)) return is_above(cert) ? Cert.Outside : is_below(cert) ? Cert.Inside : Cert.On end end -point_position_relative_to_circumcircle(tri::Triangulation, i, j, k, ℓ) = point_position_relative_to_circumcircle(AdaptiveKernel(), tri, i, j, k, ℓ) -point_position_relative_to_circumcircle(kernel::AbstractPredicateKernel, tri::Triangulation, T, ℓ) = point_position_relative_to_circumcircle(kernel, tri, geti(T), getj(T), getk(T), ℓ) -point_position_relative_to_circumcircle(tri::Triangulation, T, ℓ) = point_position_relative_to_circumcircle(AdaptiveKernel(), tri, geti(T), getj(T), getk(T), ℓ) +point_position_relative_to_circumcircle(tri::Triangulation, i, j, k, ℓ; cache::PredicateCacheType = nothing) = point_position_relative_to_circumcircle(AdaptiveKernel(), tri, i, j, k, ℓ; cache) +point_position_relative_to_circumcircle(kernel::AbstractPredicateKernel, tri::Triangulation, T, ℓ; cache::PredicateCacheType = nothing) = point_position_relative_to_circumcircle(kernel, tri, geti(T), getj(T), getk(T), ℓ; cache) +point_position_relative_to_circumcircle(tri::Triangulation, T, ℓ; cache::PredicateCacheType = nothing) = point_position_relative_to_circumcircle(AdaptiveKernel(), tri, geti(T), getj(T), getk(T), ℓ; cache) + +""" + validate_incircle_cache(tri::Triangulation, cache) -> Bool + +Checks if the cache `cache` is of the correct type for computing incircle predicates for the triangulation `tri`. If `isnothing(cache)` or +the cache is of the correct type, then `true` is returned. Otherwise, `false` is returned. +""" +@inline function validate_incircle_cache(tri::Triangulation, cache) + isnothing(cache) && return true + incircle_cache = get_incircle_cache(tri) + return typeof(incircle_cache) === typeof(cache) +end + +""" + validate_orient3_cache(tri::Triangulation, cache) -> Bool + +Checks if the cache `cache` is of the correct type for cmoputing orient3 predicates for the triangulation `tri`. If `isnothing(cache)` or +the cache is of the correct type, then `true` is returned. Otherwise, `false` is returned. +""" +@inline function validate_orient3_cache(tri::Triangulation, cache) + isnothing(cache) && return true + orient3_cache = get_orient3_cache(tri) + return typeof(orient3_cache) === typeof(cache) +end + +""" + fix_incircle_cache(tri::Triangulation, cache) + +Returns `cache` if `validate_incircle_cache(tri, cache)` is `true`, otherwise returns `nothing`. +""" +@inline fix_incircle_cache(tri::Triangulation, cache) = validate_incircle_cache(tri, cache) ? cache : nothing """ - point_position_relative_to_witness_plane([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, i, j, k, ℓ) -> Certificate + fix_orient3_cache(tri::Triangulation, cache) + +Returns `cache` if `validate_orient3_cache(tri, cache)` is `true`, otherwise returns `nothing`. +""" +@inline fix_orient3_cache(tri::Triangulation, cache) = validate_orient3_cache(tri, cache) ? cache : nothing + +""" + point_position_relative_to_witness_plane([kernel::AbstractPredicateKernel=AdaptiveKernel(),] tri::Triangulation, i, j, k, ℓ; cache = nothing) -> Certificate Given a positively oriented triangle `T = (i, j, k)` of `tri` and a vertex `ℓ` of `tri`, returns the position of `ℓ` relative to the witness plane of `T`. The returned value is a [`Certificate`](@ref), which is one of: @@ -597,19 +651,23 @@ The `kernel` argument determines how this result is computed, and should be one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default). See the documentation for more information about these choices. +The `cache` keyword argument is useful when an `AdaptiveKernel()` is used, and should be one of: +- `nothing`: No cache is used. +- `AdaptivePredicates.orient3adapt_cache(T)`, where `T` is the number type used (`Float64` or `Float32`). + # Extended help The witness plane of a triangle is defined as the plane through the triangle `(i⁺, j⁺, k⁺)` where, for example, `pᵢ⁺ = (x, y, x^2 + y^2 - wᵢ)` is the lifted companion of `i` and `(x, y)` are the coordinates of the `i`th vertex. Moreover, by the orientation of `ℓ` relative to this witness plane we are referring to `ℓ⁺`'s position, not the plane point `ℓ`. """ -function point_position_relative_to_witness_plane(kernel::AbstractPredicateKernel, tri::Triangulation, i, j, k, ℓ) +function point_position_relative_to_witness_plane(kernel::AbstractPredicateKernel, tri::Triangulation, i, j, k, ℓ; cache::PredicateCacheType = nothing) p⁺ = get_lifted_point(tri, i) q⁺ = get_lifted_point(tri, j) r⁺ = get_lifted_point(tri, k) a⁺ = get_lifted_point(tri, ℓ) - cert = orient_predicate(kernel, p⁺, q⁺, r⁺, a⁺) + cert = orient_predicate(kernel, p⁺, q⁺, r⁺, a⁺; cache) return convert_certificate(cert, Cert.Above, Cert.On, Cert.Below) end -point_position_relative_to_witness_plane(tri::Triangulation, i, j, k, ℓ) = point_position_relative_to_witness_plane(AdaptiveKernel(), tri, i, j, k, ℓ) +point_position_relative_to_witness_plane(tri::Triangulation, i, j, k, ℓ; cache::PredicateCacheType = nothing) = point_position_relative_to_witness_plane(AdaptiveKernel(), tri, i, j, k, ℓ; cache) """ opposite_angle([kernel::AbstractPredicateKernel=AdaptiveKernel(),] p, q, r) -> Certificate diff --git a/src/public.jl b/src/public.jl index a6d764fa0..82637b7b1 100644 --- a/src/public.jl +++ b/src/public.jl @@ -185,7 +185,9 @@ meet_predicate, sort_triangle, triangle_orthocenter, - triangle_orthoradius_squared + triangle_orthoradius_squared, + get_incircle_cache, + get_orient3_cache """, ), ) diff --git a/test/data_structures/triangulation.jl b/test/data_structures/triangulation.jl index edf2224a0..73c3a6194 100644 --- a/test/data_structures/triangulation.jl +++ b/test/data_structures/triangulation.jl @@ -6,80 +6,81 @@ using Setfield using Random using StatsBase using StableRNGs +using AdaptivePredicates using ..DelaunayTriangulation: add_weight!, get_weight, get_weights -@struct_equal DT.TriangulationCache - - using ..DelaunayTriangulation: Triangulation global pts = rand(2, 500) -global tri = Triangulation(pts; IntegerType = Int32) +global tri = Triangulation(pts; IntegerType=Int32) @testset "Initialising a triangulation" begin @test tri ⊢ Triangulation( pts, # points - Set{NTuple{3, Int32}}(), # triangles + Set{NTuple{3,Int32}}(), # triangles Int32[], # boundary_nodes - Set{NTuple{2, Int32}}(), # interior_segments - Set{NTuple{2, Int32}}(), # all_segments + Set{NTuple{2,Int32}}(), # interior_segments + Set{NTuple{2,Int32}}(), # all_segments DT.ZeroWeight(), # weights - DT.Adjacent{Int32, NTuple{2, Int32}}(), # adjacent - DT.Adjacent2Vertex{Int32, Set{NTuple{2, Int32}}}(), # adjacent2vertex + DT.Adjacent{Int32,NTuple{2,Int32}}(), # adjacent + DT.Adjacent2Vertex{Int32,Set{NTuple{2,Int32}}}(), # adjacent2vertex DT.Graph{Int32}(), # graph (), # boundary_curves - Dict{NTuple{2, Int32}, Tuple{Vector{Int32}, Int32}}(), # boundary_edge_map - Dict{Int32, Vector{Int32}}(DT.𝒢 => Int32[]), # ghost_vertex_map - Dict{Int32, UnitRange{Int32}}(-1 => -1:-1), # ghost_vertex_ranges + Dict{NTuple{2,Int32},Tuple{Vector{Int32},Int32}}(), # boundary_edge_map + Dict{Int32,Vector{Int32}}(DT.𝒢 => Int32[]), # ghost_vertex_map + Dict{Int32,UnitRange{Int32}}(-1 => -1:-1), # ghost_vertex_ranges DT.ConvexHull(pts, Int32[]), # convex_hull - Dict{Int32, DT.RepresentativeCoordinates{Int32, Float64}}(), # representative_point_coordinate - DT.construct_polygon_hierarchy(pts; IntegerType = Int32), + Dict{Int32,DT.RepresentativeCoordinates{Int32,Float64}}(), # representative_point_coordinate + DT.construct_polygon_hierarchy(pts; IntegerType=Int32), nothing, DT.TriangulationCache( DT.Triangulation( pts, - Set{NTuple{3, Int32}}(), # triangles + Set{NTuple{3,Int32}}(), # triangles Int32[], # boundary_nodes - Set{NTuple{2, Int32}}(), # interior_segments - Set{NTuple{2, Int32}}(), # all_segments + Set{NTuple{2,Int32}}(), # interior_segments + Set{NTuple{2,Int32}}(), # all_segments DT.ZeroWeight(), # weights - DT.Adjacent{Int32, NTuple{2, Int32}}(), # adjacent - DT.Adjacent2Vertex{Int32, Set{NTuple{2, Int32}}}(), # adjacent2vertex + DT.Adjacent{Int32,NTuple{2,Int32}}(), # adjacent + DT.Adjacent2Vertex{Int32,Set{NTuple{2,Int32}}}(), # adjacent2vertex DT.Graph{Int32}(), # graph (), # boundary_curves - Dict{NTuple{2, Int32}, Tuple{Vector{Int32}, Int32}}(), # boundary_edge_map - Dict{Int32, Vector{Int32}}(DT.𝒢 => Int32[]), # ghost_vertex_map - Dict{Int32, UnitRange{Int32}}(-1 => -1:-1), # ghost_vertex_ranges + Dict{NTuple{2,Int32},Tuple{Vector{Int32},Int32}}(), # boundary_edge_map + Dict{Int32,Vector{Int32}}(DT.𝒢 => Int32[]), # ghost_vertex_map + Dict{Int32,UnitRange{Int32}}(-1 => -1:-1), # ghost_vertex_ranges DT.ConvexHull(pts, Int32[]), # convex_hull - Dict{Int32, DT.RepresentativeCoordinates{Int32, Float64}}(), # representative_point_coordinate - DT.construct_polygon_hierarchy(pts; IntegerType = Int32), + Dict{Int32,DT.RepresentativeCoordinates{Int32,Float64}}(), # representative_point_coordinate + DT.construct_polygon_hierarchy(pts; IntegerType=Int32), nothing, # boundary_enricher - DT.TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing), + DT.TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), ), DT.Triangulation( pts, - Set{NTuple{3, Int32}}(), # triangles + Set{NTuple{3,Int32}}(), # triangles Int32[], # boundary_nodes - Set{NTuple{2, Int32}}(), # interior_segments - Set{NTuple{2, Int32}}(), # all_segments + Set{NTuple{2,Int32}}(), # interior_segments + Set{NTuple{2,Int32}}(), # all_segments DT.ZeroWeight(), # weights - DT.Adjacent{Int32, NTuple{2, Int32}}(), # adjacent - DT.Adjacent2Vertex{Int32, Set{NTuple{2, Int32}}}(), # adjacent2vertex + DT.Adjacent{Int32,NTuple{2,Int32}}(), # adjacent + DT.Adjacent2Vertex{Int32,Set{NTuple{2,Int32}}}(), # adjacent2vertex DT.Graph{Int32}(), # graph (), # boundary_curves - Dict{NTuple{2, Int32}, Tuple{Vector{Int32}, Int32}}(), # boundary_edge_map - Dict{Int32, Vector{Int32}}(DT.𝒢 => Int32[]), # ghost_vertex_map - Dict{Int32, UnitRange{Int32}}(-1 => -1:-1), # ghost_vertex_ranges + Dict{NTuple{2,Int32},Tuple{Vector{Int32},Int32}}(), # boundary_edge_map + Dict{Int32,Vector{Int32}}(DT.𝒢 => Int32[]), # ghost_vertex_map + Dict{Int32,UnitRange{Int32}}(-1 => -1:-1), # ghost_vertex_ranges DT.ConvexHull(pts, Int32[]), # convex_hull - Dict{Int32, DT.RepresentativeCoordinates{Int32, Float64}}(), # representative_point_coordinate - DT.construct_polygon_hierarchy(pts; IntegerType = Int32), + Dict{Int32,DT.RepresentativeCoordinates{Int32,Float64}}(), # representative_point_coordinate + DT.construct_polygon_hierarchy(pts; IntegerType=Int32), nothing, # boundary_enricher - DT.TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing), + DT.TriangulationCache(nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), ), Int32[], - Set{NTuple{2, Int32}}(), + Set{NTuple{2,Int32}}(), Vector{Int32}(), - Set{NTuple{3, Int32}}(), + Set{NTuple{3,Int32}}(), + AdaptivePredicates.incircleadapt_cache(Float64), + AdaptivePredicates.orient3adapt_cache(Float64), + AdaptivePredicates.insphereexact_cache(Float64), ), ) end @@ -88,21 +89,21 @@ _x, _y = complicated_geometry() global x = _x global y = _y boundary_nodes, points = convert_boundary_points_to_indices(x, y) -global tri = triangulate(points; boundary_nodes, delete_ghosts = false) +global tri = triangulate(points; boundary_nodes, delete_ghosts=false) A = get_area(tri) -refine!(tri; max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) +refine!(tri; max_area=1.0e-1A, use_circumcenter=true, use_lens=false) boundary_nodes, points = convert_boundary_points_to_indices(x[1], y[1]) -global tri_2 = triangulate(points; boundary_nodes, delete_ghosts = false) +global tri_2 = triangulate(points; boundary_nodes, delete_ghosts=false) A = get_area(tri_2) -refine!(tri_2; max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) +refine!(tri_2; max_area=1.0e-1A, use_circumcenter=true, use_lens=false) boundary_nodes, points = convert_boundary_points_to_indices([0.0, 2.0, 2.0, 0.0, 0.0], [0.0, 0.0, 2.0, 2.0, 0.0]) -global tri_3 = triangulate(points; boundary_nodes, delete_ghosts = false) +global tri_3 = triangulate(points; boundary_nodes, delete_ghosts=false) A = get_area(tri_3) -refine!(tri_3; max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) +refine!(tri_3; max_area=1.0e-1A, use_circumcenter=true, use_lens=false) boundary_nodes, points = convert_boundary_points_to_indices(reverse(reverse.(x[2])), reverse(reverse.(y[2]))) -global tri_4 = triangulate(points; boundary_nodes, delete_ghosts = false) +global tri_4 = triangulate(points; boundary_nodes, delete_ghosts=false) A = get_area(tri_4) -refine!(tri_4; max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) +refine!(tri_4; max_area=1.0e-1A, use_circumcenter=true, use_lens=false) @testset "Triangulation getters" begin @test DT.get_points(tri) == tri.points @@ -129,6 +130,9 @@ refine!(tri_4; max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) @test DT.get_ghost_vertex_ranges(tri) == tri.ghost_vertex_ranges == DT.construct_ghost_vertex_ranges(get_boundary_nodes(tri)) @test DT.get_boundary_edge_map(tri) == tri.boundary_edge_map == DT.construct_boundary_edge_map(get_boundary_nodes(tri)) @test DT.get_cache(tri) == tri.cache + @test DT.get_incircle_cache(tri) == tri.cache.incircle_cache + @test DT.get_orient3_cache(tri) == tri.cache.orient3_cache + @test DT.get_insphere_cache(tri) == tri.cache.insphere_cache @inferred DT.get_points(tri) @inferred DT.get_triangles(tri) @inferred DT.get_boundary_nodes(tri) @@ -145,6 +149,12 @@ refine!(tri_4; max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) @inferred DT.get_convex_hull(tri) @inferred DT.get_representative_point_list(tri) @inferred DT.get_exterior_curve_indices(tri) + @inferred DT.get_boundary_enricher(tri) + @inferred DT.get_polygon_hierarchy(tri) + @inferred DT.get_cache(tri) + @inferred DT.get_incircle_cache(tri) + @inferred DT.get_orient3_cache(tri) + @inferred DT.get_insphere_cache(tri) end @testset "Forwarded methods" begin @@ -210,14 +220,14 @@ end @testset "ConvexHull" begin @test DT.get_convex_hull(tri) == tri.convex_hull @test DT.get_convex_hull_vertices(tri) == tri.convex_hull.vertices - _tri = triangulate_rectangle(0.0, 2.0, 5.0, 7.3, 5, 15; single_boundary = true) + _tri = triangulate_rectangle(0.0, 2.0, 5.0, 7.3, 5, 15; single_boundary=true) points = get_points(_tri) ch = convex_hull(points) @test ch.vertices == _tri.boundary_nodes _indices = deepcopy(ch.vertices) __indices = DT.get_convex_hull_vertices(_tri) empty!(__indices) - DT.convex_hull!(_tri; reconstruct = false) + DT.convex_hull!(_tri; reconstruct=false) @test length(__indices) == length(_indices) unique!(__indices) unique!(ch.vertices) @@ -231,7 +241,7 @@ end ch.vertices .= circshift(ch.vertices, 1 - shift) @test __indices == ch.vertices delete_ghost_triangles!(_tri) - DT.convex_hull!(_tri; reconstruct = false) + DT.convex_hull!(_tri; reconstruct=false) @test length(__indices) == length(_indices) unique!(__indices) unique!(ch.vertices) @@ -261,7 +271,7 @@ end c, d = 3.0, 7.0 nx = 12 ny = 15 - @test DT.num_curves(triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 10, 10; delete_ghosts = false, single_boundary = false)) == 1 + @test DT.num_curves(triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 10, 10; delete_ghosts=false, single_boundary=false)) == 1 @test DT.num_sections(tri_2) == 4 @inferred DT.num_sections(tri_2) @test DT.get_boundary_nodes(tri, 1) == tri.boundary_nodes[1] @@ -284,10 +294,10 @@ end @testset "Triangles" begin rng = StableRNG(9882881) boundary_nodes, points = convert_boundary_points_to_indices(x, y) - tri = triangulate(points; rng, boundary_nodes, delete_ghosts = false) + tri = triangulate(points; rng, boundary_nodes, delete_ghosts=false) A = get_area(tri) - refine!(tri; max_area = 1.0e-1A, rng, use_circumcenter = true) - @test DT.triangle_type(tri) == NTuple{3, Int} + refine!(tri; max_area=1.0e-1A, rng, use_circumcenter=true) + @test DT.triangle_type(tri) == NTuple{3,Int} @inferred DT.triangle_type(tri) @test DT.num_triangles(tri) == length(tri.triangles) @test DT.each_triangle(tri) == tri.triangles @@ -303,7 +313,7 @@ end @test DelaunayTriangulation.each_triangle(_solid_itr) == _solid_itr @test Base.IteratorSize(_solid_itr) == Base.HasLength() @test Base.IteratorEltype(_solid_itr) == Base.HasEltype() - @test Base.eltype(_solid_itr) == NTuple{3, Int} + @test Base.eltype(_solid_itr) == NTuple{3,Int} @test each_solid_triangle(tri) isa DT.EachSolidTriangle _solid_tri = collect(_solid_itr) @inferred collect(_solid_itr) @@ -311,7 +321,7 @@ end _ghost_itr = each_ghost_triangle(tri) @test Base.IteratorSize(_ghost_itr) == Base.HasLength() @test Base.IteratorEltype(_ghost_itr) == Base.HasEltype() - @test Base.eltype(_ghost_itr) == NTuple{3, Int} + @test Base.eltype(_ghost_itr) == NTuple{3,Int} @test each_ghost_triangle(tri) isa DT.EachGhostTriangle _ghost_tri = collect(_ghost_itr) @inferred collect(_ghost_itr) @@ -322,9 +332,9 @@ end @test length(_ghost_tri) == DT.num_ghost_triangles(tri) == length(_ghost_itr) boundary_nodes, points = convert_boundary_points_to_indices(x, y) rng = StableRNG(9882881) - ___tri = triangulate(points; boundary_nodes, delete_ghosts = false, rng) + ___tri = triangulate(points; boundary_nodes, delete_ghosts=false, rng) A = get_area(___tri) - refine!(___tri; max_area = 1.0e-1A, rng, use_circumcenter = true) + refine!(___tri; max_area=1.0e-1A, rng, use_circumcenter=true) DT.delete_ghost_triangles!(___tri) @test collect(each_triangle(___tri)) == collect(each_solid_triangle(___tri)) @test length(collect(each_ghost_triangle(___tri))) == 0 @@ -337,10 +347,10 @@ end @testset "Edges" begin rng = StableRNG(998871) boundary_nodes, points = convert_boundary_points_to_indices(x, y) - tri = triangulate(points; rng, boundary_nodes, delete_ghosts = false) + tri = triangulate(points; rng, boundary_nodes, delete_ghosts=false) A = get_area(tri) - refine!(tri; max_area = 1.0e-1A, rng, use_circumcenter = true) - @test DT.edge_type(tri) == NTuple{2, Int} + refine!(tri; max_area=1.0e-1A, rng, use_circumcenter=true) + @test DT.edge_type(tri) == NTuple{2,Int} @inferred DT.edge_type(tri) @test DT.num_edges(tri) == length(tri.graph.edges) @inferred DT.num_edges(tri) @@ -350,7 +360,7 @@ end @test DelaunayTriangulation.each_edge(_solid_itr) == _solid_itr @test Base.IteratorSize(_solid_itr) == Base.HasLength() @test Base.IteratorEltype(_solid_itr) == Base.HasEltype() - @test Base.eltype(_solid_itr) == NTuple{2, Int} + @test Base.eltype(_solid_itr) == NTuple{2,Int} @test each_solid_edge(tri) isa DT.EachSolidEdge _solid_tri = collect(_solid_itr) @inferred collect(_solid_itr) @@ -358,7 +368,7 @@ end _ghost_itr = each_ghost_edge(tri) @test Base.IteratorSize(_ghost_itr) == Base.HasLength() @test Base.IteratorEltype(_ghost_itr) == Base.HasEltype() - @test Base.eltype(_ghost_itr) == NTuple{2, Int} + @test Base.eltype(_ghost_itr) == NTuple{2,Int} @test each_ghost_edge(tri) isa DT.EachGhostEdge _ghost_tri = collect(_ghost_itr) @inferred collect(_ghost_itr) @@ -367,9 +377,9 @@ end @test length(_ghost_tri) + length(_solid_tri) == num_edges(tri) rng = StableRNG(998871) boundary_nodes, points = convert_boundary_points_to_indices(x, y) - ___tri = triangulate(points; rng, boundary_nodes, delete_ghosts = false) + ___tri = triangulate(points; rng, boundary_nodes, delete_ghosts=false) A = get_area(___tri) - refine!(___tri; max_area = 1.0e-1A, rng, use_circumcenter = true) + refine!(___tri; max_area=1.0e-1A, rng, use_circumcenter=true) DT.delete_ghost_triangles!(___tri) @test sort(collect(filter(!DT.is_ghost_edge, each_edge(___tri)))) == sort(collect(each_solid_edge(___tri))) @test sort(collect(filter(DT.is_ghost_edge, each_edge(___tri)))) == sort(collect(each_ghost_edge(___tri))) @@ -391,7 +401,7 @@ end @test DT.num_points(tri) == length(tri.points) @inferred DT.num_points(tri) @test DT.get_point(tri, 2, 5, 6, 10) == - ntuple(i -> Tuple(tri.points[(2, 5, 6, 10)[i]]), 4) + ntuple(i -> Tuple(tri.points[(2, 5, 6, 10)[i]]), 4) @inferred DT.get_point(tri, 2, 5, 6, 10) rep = DT.get_representative_point_list(tri) rep[1] = DT.RepresentativeCoordinates(0.5, 0.3, 2) @@ -439,7 +449,7 @@ end boundary_nodes, points = convert_boundary_points_to_indices(x, y) _triq = triangulate(points; boundary_nodes, rng) A = get_area(_triq) - refine!(_triq; max_area = 1.0e-1A, rng, use_circumcenter = true, use_lens = false) + refine!(_triq; max_area=1.0e-1A, rng, use_circumcenter=true, use_lens=false) _solid_itr = each_solid_vertex(_triq) @test DelaunayTriangulation.each_vertex(_solid_itr) == _solid_itr @test Base.IteratorSize(_solid_itr) == Base.HasLength() @@ -477,7 +487,7 @@ end rng = StableRNG(887271) ___tri = triangulate(points; boundary_nodes, rng) A = get_area(___tri) - refine!(___tri; max_area = 1.0e-1A, rng, use_circumcenter = true) + refine!(___tri; max_area=1.0e-1A, rng, use_circumcenter=true) DT.delete_ghost_triangles!(___tri) @test sort(collect(filter(!DT.is_ghost_vertex, each_vertex(___tri)))) == sort(collect(each_solid_vertex(___tri))) @test sort(collect(filter(DT.is_ghost_vertex, each_vertex(___tri)))) == sort(collect(each_ghost_vertex(___tri))) @@ -505,10 +515,10 @@ end @testset "Miscellaneous" begin @test DT.integer_type(tri) == Int @test DT.number_type(tri) == Float64 - @test DT.edge_type(tri) == NTuple{2, Int} - @test DT.edges_type(tri) == Set{NTuple{2, Int}} - @test DT.triangles_type(tri) == Set{NTuple{3, Int}} - @test DT.triangle_type(tri) == NTuple{3, Int} + @test DT.edge_type(tri) == NTuple{2,Int} + @test DT.edges_type(tri) == Set{NTuple{2,Int}} + @test DT.triangles_type(tri) == Set{NTuple{3,Int}} + @test DT.triangle_type(tri) == NTuple{3,Int} @inferred DT.integer_type(tri) @inferred DT.number_type(tri) DT.clear_empty_features!(tri) @@ -535,7 +545,7 @@ end empty!(get_all_segments(_tri)) @test !DT.is_constrained(_tri) @test reverse(sort(collect(DT.all_ghost_vertices(_tri)))) == - [DT.𝒢, DT.𝒢 - 1, DT.𝒢 - 2, DT.𝒢 - 3] + [DT.𝒢, DT.𝒢 - 1, DT.𝒢 - 2, DT.𝒢 - 3] end end @@ -558,18 +568,18 @@ end bn43 = bn4[3] bn44 = bn4[4] bn5 = all_bn[5][1] - e11 = Set(((bn11[i], bn11[i + 1]) for i in 1:(length(bn11) - 1))) - e12 = Set(((bn12[i], bn12[i + 1]) for i in 1:(length(bn12) - 1))) - e13 = Set(((bn13[i], bn13[i + 1]) for i in 1:(length(bn13) - 1))) - e14 = Set(((bn14[i], bn14[i + 1]) for i in 1:(length(bn14) - 1))) - e2 = Set(((bn2[i], bn2[i + 1]) for i in 1:(length(bn2) - 1))) - e3 = Set(((bn3[i], bn3[i + 1]) for i in 1:(length(bn3) - 1))) - e41 = Set(((bn41[i], bn41[i + 1]) for i in 1:(length(bn41) - 1))) - e42 = Set(((bn42[i], bn42[i + 1]) for i in 1:(length(bn42) - 1))) - e43 = Set(((bn43[i], bn43[i + 1]) for i in 1:(length(bn43) - 1))) - e44 = Set(((bn44[i], bn44[i + 1]) for i in 1:(length(bn44) - 1))) - e5 = Set(((bn5[i], bn5[i + 1]) for i in 1:(length(bn5) - 1))) - ace = Set{NTuple{2, Int}}() + e11 = Set(((bn11[i], bn11[i+1]) for i in 1:(length(bn11)-1))) + e12 = Set(((bn12[i], bn12[i+1]) for i in 1:(length(bn12)-1))) + e13 = Set(((bn13[i], bn13[i+1]) for i in 1:(length(bn13)-1))) + e14 = Set(((bn14[i], bn14[i+1]) for i in 1:(length(bn14)-1))) + e2 = Set(((bn2[i], bn2[i+1]) for i in 1:(length(bn2)-1))) + e3 = Set(((bn3[i], bn3[i+1]) for i in 1:(length(bn3)-1))) + e41 = Set(((bn41[i], bn41[i+1]) for i in 1:(length(bn41)-1))) + e42 = Set(((bn42[i], bn42[i+1]) for i in 1:(length(bn42)-1))) + e43 = Set(((bn43[i], bn43[i+1]) for i in 1:(length(bn43)-1))) + e44 = Set(((bn44[i], bn44[i+1]) for i in 1:(length(bn44)-1))) + e5 = Set(((bn5[i], bn5[i+1]) for i in 1:(length(bn5)-1))) + ace = Set{NTuple{2,Int}}() for es in (e11, e12, e13, e14, e2, e3, e41, e42, e43, e44, e5) for e in es push!(ace, e) @@ -589,11 +599,11 @@ end bn2 = all_bn[2] bn3 = all_bn[3] bn4 = all_bn[4] - e1 = Set(((bn1[i], bn1[i + 1]) for i in 1:(length(bn1) - 1))) - e2 = Set(((bn2[i], bn2[i + 1]) for i in 1:(length(bn2) - 1))) - e3 = Set(((bn3[i], bn3[i + 1]) for i in 1:(length(bn3) - 1))) - e4 = Set(((bn4[i], bn4[i + 1]) for i in 1:(length(bn4) - 1))) - ace = Set{NTuple{2, Int}}() + e1 = Set(((bn1[i], bn1[i+1]) for i in 1:(length(bn1)-1))) + e2 = Set(((bn2[i], bn2[i+1]) for i in 1:(length(bn2)-1))) + e3 = Set(((bn3[i], bn3[i+1]) for i in 1:(length(bn3)-1))) + e4 = Set(((bn4[i], bn4[i+1]) for i in 1:(length(bn4)-1))) + ace = Set{NTuple{2,Int}}() for es in (e1, e2, e3, e4) for e in es push!(ace, e) @@ -609,8 +619,8 @@ end j = rand(1:100000, 50) all_ce = Set(((i, j) for (i, j) in zip(i, j))) bn_map = get_ghost_vertex_map(tri_3) - e = Set(((all_bn[i], all_bn[i + 1]) for i in 1:(length(all_bn) - 1))) - ace = Set{NTuple{2, Int}}() + e = Set(((all_bn[i], all_bn[i+1]) for i in 1:(length(all_bn)-1))) + ace = Set{NTuple{2,Int}}() for e in e push!(ace, e) end @@ -624,8 +634,8 @@ end j = rand(1:100000, 50) all_ce = Set(((i, j) for (i, j) in zip(i, j))) bn_map = get_ghost_vertex_map(tri_4) - e = Set(((all_bn[i], all_bn[i + 1]) for i in 1:(length(all_bn) - 1))) - ace = Set{NTuple{2, Int}}() + e = Set(((all_bn[i], all_bn[i+1]) for i in 1:(length(all_bn)-1))) + ace = Set{NTuple{2,Int}}() for e in e push!(ace, e) end @@ -636,7 +646,7 @@ end end @testset "sort_edge_by_degree" begin - tri = triangulate(rand(2, 500); delete_ghosts = false) + tri = triangulate(rand(2, 500); delete_ghosts=false) for e in each_edge(tri) new_e = DT.sort_edge_by_degree(tri, e) d1 = DT.num_neighbours(tri, e[1]) @@ -653,17 +663,17 @@ end n = 20 for _ in 1:10 n += rand(1:5) - tri1 = triangulate(12randn(2, n), delete_ghosts = false) - tri2 = triangulate(12randn(2, n), delete_ghosts = true) + tri1 = triangulate(12randn(2, n), delete_ghosts=false) + tri2 = triangulate(12randn(2, n), delete_ghosts=true) for tri in (tri1, tri2) for qi in each_solid_vertex(tri) for k in each_solid_vertex(tri) q = get_point(tri, qi) - history = DT.PointLocationHistory{NTuple{3, Int}, NTuple{2, Int}, Int}() + history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}() find_triangle( tri, q; k, - store_history = true, + store_history=true, history, ) visited_triangles = history.triangles @@ -707,7 +717,7 @@ pts = [(rand(rng), rand(rng)) for _ in 1:50] bnd_pts = [(0.3cos(θ), 0.3sin(θ)) .+ 0.5 for θ in LinRange(0, 2π - 1 / 250, 25)] bnd_id = [(51:75)..., 51] append!(pts, bnd_pts) -global tric = triangulate(pts; boundary_nodes = bnd_id, rng) +global tric = triangulate(pts; boundary_nodes=bnd_id, rng) @testset "each_segment" begin @test each_segment(tric) == each_edge(get_all_segments(tric)) @@ -737,16 +747,16 @@ end boundary_nodes, points = convert_boundary_points_to_indices(x, y) tri = triangulate(points; rng, boundary_nodes) A = get_area(tri) - refine!(tri, rng = rng, max_area = 1.0e-1A, use_circumcenter = true, use_lens = false) + refine!(tri, rng=rng, max_area=1.0e-1A, use_circumcenter=true, use_lens=false) all_bn = DT.get_all_boundary_nodes(tri) @test all_bn == Set(reduce(vcat, reduce(vcat, get_boundary_nodes(tri)))) tri2, label_map, index_map = simple_geometry() all_bn = DT.get_all_boundary_nodes(tri2) @test all_bn == Set(reduce(vcat, reduce(vcat, get_boundary_nodes(tri2)))) - tri3 = triangulate_rectangle(0, 1, 0, 1, 50, 50; delete_ghosts = false, single_boundary = false) + tri3 = triangulate_rectangle(0, 1, 0, 1, 50, 50; delete_ghosts=false, single_boundary=false) all_bn = DT.get_all_boundary_nodes(tri3) @test all_bn == Set(reduce(vcat, reduce(vcat, get_boundary_nodes(tri3)))) - tri4 = triangulate_rectangle(0, 1, 0, 1, 50, 50; delete_ghosts = false, single_boundary = true) + tri4 = triangulate_rectangle(0, 1, 0, 1, 50, 50; delete_ghosts=false, single_boundary=true) all_bn = DT.get_all_boundary_nodes(tri4) @test all_bn == Set(reduce(vcat, reduce(vcat, get_boundary_nodes(tri4)))) tri = triangulate(rand(2, 50)) @@ -759,7 +769,7 @@ end boundary_nodes, points = convert_boundary_points_to_indices(x, y) tri = triangulate(points; rng, boundary_nodes) A = get_area(tri) - refine!(tri, rng = rng, max_area = 1.0e-2A, use_circumcenter = true) + refine!(tri, rng=rng, max_area=1.0e-2A, use_circumcenter=true) for (e, (s, i)) in tri.boundary_edge_map @test DT.get_boundary_edge_map(tri, e) == (s, i) @test DT.get_boundary_edge_map(tri, e...) == (s, i) @@ -928,14 +938,14 @@ end end @testset "get_adjacent concurrency" begin # Shouldn't be an issue anymore since we removed DefaultDict, but let's keep this here anyway. The test here is simply that it doesn't error. - tri = triangulate(rand(2, 50), delete_ghosts = false) + tri = triangulate(rand(2, 50), delete_ghosts=false) Base.Threads.@threads for _ in 1:5000 get_adjacent(tri, -5, rand(1:1000)) end end @testset "has_vertex and has_ghost_vertices" begin - tri = triangulate(rand(2, 50), delete_ghosts = false) + tri = triangulate(rand(2, 50), delete_ghosts=false) @test DT.has_vertex(tri, 1) @test !DT.has_vertex(tri, 57) @test DT.has_ghost_vertices(tri) @@ -978,8 +988,8 @@ end (-2.0, 12.0), (-4.0, 12.0), (-6.0, 11.0), (-7.0, 9.0), (-6.94, 7.13), (-6.0, 5.0), (-4.0, 3.0), (-2.0, 1.0), (0.0, 0.0), ] - boundary_nodes, pts = convert_boundary_points_to_indices(boundary_points; existing_points = pts) - tri = triangulate(pts; boundary_nodes, delete_ghosts = false) + boundary_nodes, pts = convert_boundary_points_to_indices(boundary_points; existing_points=pts) + tri = triangulate(pts; boundary_nodes, delete_ghosts=false) @test DT.is_positively_oriented(tri, 1) points = [ @@ -990,7 +1000,7 @@ end segment_2 = [(14.0, 0.0), (10.0, 4.0), (4.0, 6.0), (2.0, 12.0), (0.0, 14.0)] segment_3 = [(0.0, 14.0), (0.0, 0.0)] boundary_points = [segment_1, segment_2, segment_3] - boundary_nodes, points = convert_boundary_points_to_indices(boundary_points; existing_points = points) + boundary_nodes, points = convert_boundary_points_to_indices(boundary_points; existing_points=points) tri = triangulate(points; boundary_nodes) @test DT.is_positively_oriented(tri, 1) @@ -1027,8 +1037,8 @@ end (5.6, 7.8), (5.6, 7.6), (5.6, 7.4), (6.2, 7.4), (6.0, 7.6), (7.0, 7.8), (7.0, 7.4), ] - boundary_nodes, points = convert_boundary_points_to_indices(curves; existing_points = points) - tri = triangulate(points; boundary_nodes = boundary_nodes) + boundary_nodes, points = convert_boundary_points_to_indices(curves; existing_points=points) + tri = triangulate(points; boundary_nodes=boundary_nodes) @test DT.is_positively_oriented(tri, 1) @test !DT.is_positively_oriented(tri, 2) @test !DT.is_positively_oriented(tri, 3) @@ -1061,15 +1071,15 @@ end (13.0, 19.0), (13.0, 12.0), (16.0, 8.0), (9.8, 8.0), (7.5, 6.0), (12.0, 13.0), (19.0, 15.0), ] - boundary_nodes, points = convert_boundary_points_to_indices(curves; existing_points = points) - tri = triangulate(points; boundary_nodes = boundary_nodes) + boundary_nodes, points = convert_boundary_points_to_indices(curves; existing_points=points) + tri = triangulate(points; boundary_nodes=boundary_nodes) @test DT.is_positively_oriented(tri, 1) @test !DT.is_positively_oriented(tri, 2) @test DT.is_positively_oriented(tri, 3) θ = LinRange(0, 2π, 20) |> collect θ[end] = 0 # need to make sure that 2π gives the exact same coordinates as 0 - xy = Vector{Vector{Vector{NTuple{2, Float64}}}}() + xy = Vector{Vector{Vector{NTuple{2,Float64}}}}() cx = 0.0 for i in 1:2 # Make the exterior circle @@ -1079,7 +1089,7 @@ end cx += 3.0 end boundary_nodes, points = convert_boundary_points_to_indices(xy) - tri = triangulate(points; boundary_nodes = boundary_nodes) + tri = triangulate(points; boundary_nodes=boundary_nodes) @test DT.is_positively_oriented(tri, 1) @test !DT.is_positively_oriented(tri, 2) @test DT.is_positively_oriented(tri, 3) @@ -1243,7 +1253,7 @@ end dot_4 = [[K3, L3, M3, N3, O3, P3, Q3, R3, S3, T3, U3, V3, K3]] curves = [J_curve, U_curve, L_curve, I_curve, A_curve_outline, A_curve_hole, dot_1, dot_2, dot_3, dot_4] nodes, points = convert_boundary_points_to_indices(curves) - tri = triangulate(points; boundary_nodes = nodes) + tri = triangulate(points; boundary_nodes=nodes) @test DT.is_positively_oriented(tri, 1) @test DT.is_positively_oriented(tri, 2) @test DT.is_positively_oriented(tri, 3) @@ -1272,21 +1282,21 @@ end @test tri1 == tri2 # Custom types - tri1 = triangulate(rand(2, 50); IntegerType = Int32) - tri2 = Triangulation(get_points(tri1), each_solid_triangle(tri1), get_convex_hull_vertices(tri1); IntegerType = Int32) + tri1 = triangulate(rand(2, 50); IntegerType=Int32) + tri2 = Triangulation(get_points(tri1), each_solid_triangle(tri1), get_convex_hull_vertices(tri1); IntegerType=Int32) unlock_convex_hull!(tri2) @test tri1 == tri2 @test tri1 ⊢ tri2 # Delete ghosts works properly - tri1 = triangulate(rand(2, 50), delete_ghosts = true) - tri2 = Triangulation(get_points(tri1), each_solid_triangle(tri1), get_convex_hull_vertices(tri1), delete_ghosts = true) + tri1 = triangulate(rand(2, 50), delete_ghosts=true) + tri2 = Triangulation(get_points(tri1), each_solid_triangle(tri1), get_convex_hull_vertices(tri1), delete_ghosts=true) unlock_convex_hull!(tri2) @test tri1 == tri2 # Weights work properly weights = rand(50) - tri2 = Triangulation(get_points(tri1), each_solid_triangle(tri1), get_convex_hull_vertices(tri1), weights = weights) + tri2 = Triangulation(get_points(tri1), each_solid_triangle(tri1), get_convex_hull_vertices(tri1), weights=weights) @test DT.is_weighted(tri2) && get_weights(tri2) == weights end @@ -1345,7 +1355,7 @@ end nodes, points = convert_boundary_points_to_indices(boundary_pts) push!(points, (20.0, 20.0)) rng = StableRNG(19191919) - C = Set{NTuple{2, Int}}() + C = Set{NTuple{2,Int}}() for i in 1:50 θ = 2π * rand(rng) r = 4sqrt(rand(rng)) @@ -1354,13 +1364,13 @@ end push!(points, (x, y)) push!(C, (48, 48 + i)) end - tri = triangulate(points; boundary_nodes = nodes, segments = C, rng) + tri = triangulate(points; boundary_nodes=nodes, segments=C, rng) # Vertices ns = 5_000_000 function test_vertex_sampler(tri, ns) - solid_vertices = DefaultDict{Int, Float64, Float64}(0.0) - ghost_vertices = DefaultDict{Int, Float64, Float64}(0.0) + solid_vertices = DefaultDict{Int,Float64,Float64}(0.0) + ghost_vertices = DefaultDict{Int,Float64,Float64}(0.0) for _ in 1:ns sv = rand(each_solid_vertex(tri)) gv = rand(each_ghost_vertex(tri)) @@ -1378,8 +1388,8 @@ end gm = mean(values(ghost_vertices)) @test sm ≈ 1 / DT.num_solid_vertices(tri) @test gm ≈ 1 / 2 - @test all(y -> isapprox(y, sm, rtol = 1.0e-1), values(solid_vertices)) - @test all(y -> isapprox(y, gm, rtol = 1.0e-1), values(ghost_vertices)) + @test all(y -> isapprox(y, sm, rtol=1.0e-1), values(solid_vertices)) + @test all(y -> isapprox(y, gm, rtol=1.0e-1), values(ghost_vertices)) for i in 1:250 for f in (each_vertex, each_solid_vertex, each_ghost_vertex) rng = StableRNG(i) @@ -1396,8 +1406,8 @@ end # Triangles ns = 5_000_000 function test_triangle_sampler(tri, ns) - solid_triangles = DefaultDict{NTuple{3, Int}, Float64, Float64}(0.0) - ghost_triangles = DefaultDict{NTuple{3, Int}, Float64, Float64}(0.0) + solid_triangles = DefaultDict{NTuple{3,Int},Float64,Float64}(0.0) + ghost_triangles = DefaultDict{NTuple{3,Int},Float64,Float64}(0.0) for _ in 1:ns st = rand(each_solid_triangle(tri)) gt = rand(each_ghost_triangle(tri)) @@ -1415,8 +1425,8 @@ end gm = mean(values(ghost_triangles)) @test sm ≈ 1 / DT.num_solid_triangles(tri) @test gm ≈ 1 / DT.num_ghost_triangles(tri) - @test all(y -> isapprox(y, sm, rtol = 1.0e-1), values(solid_triangles)) - @test all(y -> isapprox(y, gm, rtol = 1.0e-1), values(ghost_triangles)) + @test all(y -> isapprox(y, sm, rtol=1.0e-1), values(solid_triangles)) + @test all(y -> isapprox(y, gm, rtol=1.0e-1), values(ghost_triangles)) for i in 1:250 for f in (each_triangle, each_solid_triangle, each_ghost_triangle) rng = StableRNG(i) @@ -1433,8 +1443,8 @@ end # Edges ns = 5_000_000 function test_edge_sampler(tri, ns) - solid_edges = DefaultDict{NTuple{2, Int}, Float64, Float64}(0.0) - ghost_edges = DefaultDict{NTuple{2, Int}, Float64, Float64}(0.0) + solid_edges = DefaultDict{NTuple{2,Int},Float64,Float64}(0.0) + ghost_edges = DefaultDict{NTuple{2,Int},Float64,Float64}(0.0) for _ in 1:ns se = rand(each_solid_edge(tri)) ge = rand(each_ghost_edge(tri)) @@ -1452,8 +1462,8 @@ end gm = mean(values(ghost_edges)) @test sm ≈ 1 / DT.num_solid_edges(tri) @test gm ≈ 1 / DT.num_ghost_edges(tri) - @test all(y -> isapprox(y, sm, rtol = 1.0e-1), values(solid_edges)) - @test all(y -> isapprox(y, gm, rtol = 1.0e-1), values(ghost_edges)) + @test all(y -> isapprox(y, sm, rtol=1.0e-1), values(solid_edges)) + @test all(y -> isapprox(y, gm, rtol=1.0e-1), values(ghost_edges)) for i in 1:250 for f in (each_edge, each_solid_edge, each_ghost_edge) rng = StableRNG(i) diff --git a/test/data_structures/triangulation_cache.jl b/test/data_structures/triangulation_cache.jl index 0241d9f8a..0fb721d04 100644 --- a/test/data_structures/triangulation_cache.jl +++ b/test/data_structures/triangulation_cache.jl @@ -3,11 +3,24 @@ const DT = DelaunayTriangulation using DataStructures using StructEquality using Setfield +using AdaptivePredicates using ..DelaunayTriangulation: add_weight!, get_weight, get_weights @struct_equal DT.TriangulationCache -tri = triangulate(rand(2, 50); weights = DT.ZeroWeight()) +cache = DT.TriangulationCache() +@test isnothing(DT.get_triangulation(cache)) +@test isnothing(DT.get_triangulation_2(cache)) +@test isnothing(DT.get_marked_vertices(cache)) +@test isnothing(DT.get_interior_segments_on_hull(cache)) +@test isnothing(DT.get_fan_triangles(cache)) +@test isnothing(DT.get_incircle_cache(cache)) +@test isnothing(DT.get_orient3_cache(cache)) +@test isnothing(DT.get_insphere_cache(cache)) +@test isnothing(DT.get_surrounding_polygon(cache)) +@test typeof(cache) == DT.EmptyTriangulationCache + +tri = triangulate(rand(2, 50); weights = rand(50)) cache = DT.get_cache(tri) @test get_weights(DT.get_triangulation(DT.get_cache(tri))) === get_weights(tri) DT.unconstrained_triangulation!(DT.get_triangulation(cache)) @@ -73,3 +86,6 @@ push!(surrounding_polygon, 5) @test !isempty(surrounding_polygon) empty!(cache) @test isempty(surrounding_polygon) +@test length.(DT.get_incircle_cache(cache)) == length.(AdaptivePredicates.incircleadapt_cache(Float64)) +@test length.(DT.get_orient3_cache(cache)) == length.(AdaptivePredicates.orient3adapt_cache(Float64)) +@test length.(DT.get_insphere_cache(cache)) == length.(AdaptivePredicates.insphereexact_cache(Float64)) \ No newline at end of file diff --git a/test/helper_functions.jl b/test/helper_functions.jl index baee93cf3..abbe9994a 100644 --- a/test/helper_functions.jl +++ b/test/helper_functions.jl @@ -2299,6 +2299,22 @@ function plot_small_angle_complexes(enricher) return fig end +function Base.:(==)(cache1::DT.TriangulationCache, cache2::DT.TriangulationCache) + DT.get_triangulation(cache1) ≠ DT.get_triangulation(cache2) && return false + DT.get_triangulation_2(cache1) ≠ DT.get_triangulation_2(cache2) && return false + DT.get_marked_vertices(cache1) ≠ DT.get_marked_vertices(cache2) && return false + DT.get_interior_segments_on_hull(cache1) ≠ DT.get_interior_segments_on_hull(cache2) && return false + DT.get_surrounding_polygon(cache1) ≠ DT.get_surrounding_polygon(cache2) && return false + DT.get_fan_triangles(cache1) ≠ DT.get_fan_triangles(cache2) && return false + length(DT.get_incircle_cache(cache1)) ≠ length(DT.get_incircle_cache(cache2)) && return false # no point comparing elements + typeof(DT.get_incircle_cache(cache1)) ≠ typeof(DT.get_incircle_cache(cache2)) && return false + length(DT.get_orient3_cache(cache1)) ≠ length(DT.get_orient3_cache(cache2)) && return false # no point comparing elements + typeof(DT.get_orient3_cache(cache1)) ≠ typeof(DT.get_orient3_cache(cache2)) && return false + length(DT.get_insphere_cache(cache1)) ≠ length(DT.get_insphere_cache(cache2)) && return false # no point comparing elements + typeof(DT.get_insphere_cache(cache1)) ≠ typeof(DT.get_insphere_cache(cache2)) && return false + return true +end + using DelaunayTriangulation: test_adjacent2vertex_map_matches_adjacent_map, test_state, diff --git a/test/predicates/general.jl b/test/predicates/general.jl index ca6c035e6..957f15d71 100644 --- a/test/predicates/general.jl +++ b/test/predicates/general.jl @@ -1033,4 +1033,4 @@ end flag1 && flag2 end end -end +end \ No newline at end of file diff --git a/test/triangulation/weighted.jl b/test/triangulation/weighted.jl index 72a8375ec..e743f9765 100644 --- a/test/triangulation/weighted.jl +++ b/test/triangulation/weighted.jl @@ -54,11 +54,11 @@ end @testset "is_weighted" begin tri = Triangulation(rand(2, 10)) @test !DT.is_weighted(tri) - tri = Triangulation(rand(2, 10); weights = rand(10)) + tri = Triangulation(rand(2, 10); weights=rand(10)) @test DT.is_weighted(tri) - tri = Triangulation(rand(2, 10); weights = DT.ZeroWeight()) + tri = Triangulation(rand(2, 10); weights=DT.ZeroWeight()) @test !DT.is_weighted(tri) - tri = Triangulation(rand(2, 10); weights = zeros(Float32, 10)) + tri = Triangulation(rand(2, 10); weights=zeros(Float32, 10)) @test DT.is_weighted(tri) end @@ -135,18 +135,24 @@ end _weights[58] = weights[2] _weights[498] = weights[3] _weights[5] = weights[4] - tri = Triangulation(_points; weights = _weights) - d = DT.get_distance_to_witness_plane(tri, 5, (137, 58, 498)) + tri = Triangulation(_points; weights=_weights) + d = DT.get_distance_to_witness_plane(tri, 5, (137, 58, 498); cache=nothing) @test d ≈ -2.129523129725314 - @test DT.get_distance_to_witness_plane(tri, 5, (137, 58, 498)) ≈ - DT.get_distance_to_witness_plane(tri, 5, (58, 137, 498)) ≈ - DT.get_distance_to_witness_plane(tri, 5, (498, 58, 137)) ≈ - DT.get_distance_to_witness_plane(tri, 5, (137, 498, 58)) ≈ - DT.get_distance_to_witness_plane(tri, 5, (58, 498, 137)) ≈ - DT.get_distance_to_witness_plane(tri, 5, (498, 137, 58)) + @test DT.get_distance_to_witness_plane(tri, 5, (137, 58, 498); cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(tri, 5, (58, 137, 498); cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(tri, 5, (498, 58, 137); cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(tri, 5, (137, 498, 58); cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(tri, 5, (58, 498, 137); cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(tri, 5, (498, 137, 58); cache=DT.get_incircle_cache(tri)) + @test DT.get_distance_to_witness_plane(AdaptiveKernel(), tri, 5, (137, 58, 498), cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(FastKernel(), tri, 5, (58, 137, 498), cache=DT.get_orient3_cache(tri)) ≈ + DT.get_distance_to_witness_plane(AdaptiveKernel(), tri, 5, (498, 58, 137), cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(AdaptiveKernel(), tri, 5, (137, 498, 58), cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(ExactKernel(), tri, 5, (58, 498, 137), cache=DT.get_incircle_cache(tri)) ≈ + DT.get_distance_to_witness_plane(AdaptiveKernel(), tri, 5, (498, 137, 58), cache=DT.get_incircle_cache(tri)) _points[5] = (-3.9, 2.01) _weights[5] = -15.0 - d = DT.get_distance_to_witness_plane(tri, 5, (137, 58, 498)) + d = DT.get_distance_to_witness_plane(tri, 5, (137, 58, 498); cache=DT.get_incircle_cache(tri)) @test d ≈ 5.603226942686893 # ghost triangles: exterior points @@ -162,16 +168,16 @@ end points[2] = (6.0, 3.0) points[3] = (-7.0, 3.0) points[4] = (6.19, -1.52) - @test DT.get_distance_to_witness_plane(tri, 4, (2, 1, -1)) == -Inf + @test DT.get_distance_to_witness_plane(AdaptiveKernel(), tri, 4, (2, 1, -1)) == -Inf @test DT.get_distance_to_witness_plane(tri, 4, (1, 3, -1)) == Inf @test DT.get_distance_to_witness_plane(tri, 4, (3, 2, -1)) == Inf points[4] = (-6.0, 0.0) @test DT.get_distance_to_witness_plane(tri, 4, (2, 1, -1)) == Inf - @test DT.get_distance_to_witness_plane(tri, 4, (1, 3, -1)) == -Inf + @test DT.get_distance_to_witness_plane(tri, 4, (1, 3, -1); cache=DT.get_orient3_cache(tri)) == -Inf @test DT.get_distance_to_witness_plane(tri, 4, (3, 2, -1)) == Inf points[4] = (0.0, 6.0) @test DT.get_distance_to_witness_plane(tri, 4, (2, 1, -1)) == Inf - @test DT.get_distance_to_witness_plane(tri, 4, (1, 3, -1)) == Inf + @test DT.get_distance_to_witness_plane(ExactKernel(), tri, 4, (1, 3, -1), cache=nothing) == Inf @test DT.get_distance_to_witness_plane(tri, 4, (3, 2, -1)) == -Inf # ghost triangles: points on the solid unweighted edge @@ -181,8 +187,8 @@ end points[2] = (2.0, 3.0) points[3] = (-5.0, 5.0) points[4] = (2.0, 0.0) - @test DT.get_distance_to_witness_plane(tri, 4, (2, 1, -1)) == DT.get_distance_to_witness_plane(tri, 4, (3, 1, 2)) - @test DT.get_distance_to_witness_plane(tri, 4, (1, 3, -1)) == Inf + @test DT.get_distance_to_witness_plane(tri, 4, (2, 1, -1), cache=nothing) == DT.get_distance_to_witness_plane(tri, 4, (3, 1, 2)) + @test DT.get_distance_to_witness_plane(tri, 4, (1, 3, -1), cache=nothing) == Inf # ghost triangles: points on the solid weighted edge points[3] = (-2.0, 0.0) @@ -203,13 +209,13 @@ end for fi in 1:NUM_WEGT tri, submerged, nonsubmerged, weights = get_weighted_example(fi) for i in submerged - @test DT.is_submerged(tri, i) + @test DT.is_submerged(tri, i; cache=DT.get_incircle_cache(tri)) end for i in nonsubmerged @test !DT.is_submerged(tri, i) end for i in 1:DT.num_points(tri) - @test (DT.is_submerged(tri, i, find_triangle(tri, get_point(tri, i)))) == (i ∈ submerged) + @test (DT.is_submerged(tri, i, find_triangle(tri, get_point(tri, i)); cache=DT.get_orient3_cache(tri))) == (i ∈ submerged) end end end @@ -228,7 +234,7 @@ end for i in 3:10 for j in 3:10 tri = triangulate_rectangle(0, 10, 0, 10, i, j) - tri = triangulate(get_points(tri); weights = zeros(i * j)) + tri = triangulate(get_points(tri); weights=zeros(i * j)) @test validate_triangulation(tri) validate_statistics(tri) end @@ -249,7 +255,7 @@ end for i in 3:10 for j in 3:10 tri = triangulate_rectangle(0, 10, 0, 10, i, j) - tri = triangulate(get_points(tri); weights = 10randn() * ones(i * j)) + tri = triangulate(get_points(tri); weights=10randn() * ones(i * j)) @test validate_triangulation(tri) # Why is this failing sometimes? Is validate not branching at weighted triangulations? (i == j == 10) || validate_statistics(tri) end @@ -261,13 +267,13 @@ end tri = triangulate(rand(2, 50)) for i in 1:50 V = DT.brute_force_search_enclosing_circumcircle(tri, i) - @test !DT.is_outside(DT.point_position_relative_to_circumcircle(tri, V, i)) + @test !DT.is_outside(DT.point_position_relative_to_circumcircle(tri, V, i; cache=rand() < 1 / 2 ? nothing : DT.get_incircle_cache(tri))) end for fi in 1:NUM_WEGT tri, submerged, nonsubmerged, weights = get_weighted_example(fi) for i in eachindex(weights) - V = DT.brute_force_search_enclosing_circumcircle(tri, i) + V = DT.brute_force_search_enclosing_circumcircle(tri, i; cache=DT.get_incircle_cache(tri)) if i ∈ submerged @test V == (0, 0, 0) else @@ -348,7 +354,7 @@ end points = get_points(tri) vpoints = points[:, 1:3] vweights = weights[1:3] - rtri = triangulate(vpoints; weights = vweights) + rtri = triangulate(vpoints; weights=vweights) for j in 4:size(points, 2) add_point!(rtri, points[:, j]..., weights[j]) end diff --git a/test/utils.jl b/test/utils.jl index a7b4b9f4c..e630296ab 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,6 +1,7 @@ using ..DelaunayTriangulation const DT = DelaunayTriangulation using LinearAlgebra +using AdaptivePredicates using CairoMakie using BenchmarkTools using StableRNGs @@ -1707,3 +1708,37 @@ end @test DT.ε(Float64) == DT.ε(1.0) == sqrt(eps(Float64)) @test DT.ε(Float32) == DT.ε(1.0f0) == sqrt(eps(Float32)) end + +@testset "Fixing caches" begin + tri = triangulate(rand(2, 50)) + orient3_cache = AdaptivePredicates.orient3adapt_cache(Float64) + incircle_cache = AdaptivePredicates.incircleadapt_cache(Float64) + @test DT.validate_incircle_cache(tri, incircle_cache) + @test DT.validate_incircle_cache(tri, nothing) + @test !DT.validate_incircle_cache(tri, orient3_cache) + @test DT.validate_orient3_cache(tri, orient3_cache) + @test DT.validate_orient3_cache(tri, nothing) + @test !DT.validate_orient3_cache(tri, incircle_cache) + tri32 = triangulate(rand(Float32, 2, 50)) + orient3_cache32 = AdaptivePredicates.orient3adapt_cache(Float32) + incircle_cache32 = AdaptivePredicates.incircleadapt_cache(Float32) + @test DT.validate_incircle_cache(tri32, incircle_cache32) + @test DT.validate_incircle_cache(tri32, nothing) + @test !DT.validate_incircle_cache(tri32, orient3_cache32) + @test DT.validate_orient3_cache(tri32, orient3_cache32) + @test DT.validate_orient3_cache(tri32, nothing) + @test !DT.validate_orient3_cache(tri32, incircle_cache32) + + @test DT.fix_incircle_cache(tri, incircle_cache) === incircle_cache + @test DT.fix_incircle_cache(tri, nothing) === nothing + @test DT.fix_incircle_cache(tri, orient3_cache) === nothing + @test DT.fix_orient3_cache(tri, orient3_cache) === orient3_cache + @test DT.fix_orient3_cache(tri, nothing) === nothing + @test DT.fix_orient3_cache(tri, incircle_cache) === nothing + @test DT.fix_incircle_cache(tri32, incircle_cache32) === incircle_cache32 + @test DT.fix_incircle_cache(tri32, nothing) === nothing + @test DT.fix_incircle_cache(tri32, orient3_cache32) === nothing + @test DT.fix_orient3_cache(tri32, orient3_cache32) === orient3_cache32 + @test DT.fix_orient3_cache(tri32, nothing) === nothing + @test DT.fix_orient3_cache(tri32, incircle_cache32) === nothing +end \ No newline at end of file From 1be922e4fc111e6960db8bff32d4e5689214da8a Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 15 Sep 2024 16:18:57 +0100 Subject: [PATCH 2/3] news --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 8babe42c9..6bcc02c77 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,7 +2,7 @@ ## 1.4.0 -- Updated to AdaptivePredicates.jl v1.2, now allowing caches to be passed to the predicates involving `incircle` and `orient3`. These are only useful when using the `AdaptiveKernel()` kernel. Outside of triangulating, these caches are not passed by default, but can be provided. The functions `get_incircle_cache` and `get_orient3_cache` can be used for this purpose on a triangulation (without a triangulation, refer to AdaptivePredicate.jl's `incircleadapt_cache` and `orient3adapt_cache`). +- Updated to AdaptivePredicates.jl v1.2, now allowing caches to be passed to the predicates involving `incircle` and `orient3`. These are only useful when using the `AdaptiveKernel()` kernel. Outside of triangulating, these caches are not passed by default, but can be provided. The functions `get_incircle_cache` and `get_orient3_cache` can be used for this purpose on a triangulation (without a triangulation, refer to AdaptivePredicate.jl's `incircleadapt_cache` and `orient3adapt_cache`). See [#185](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/185). ## 1.3.1 From 52f9125c50d2c59b235b62d20a57a85fdeeded10 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 15 Sep 2024 17:33:05 +0100 Subject: [PATCH 3/3] Fix tests, add docstrings --- docs/src/extended/utils.md | 4 ++++ test/data_structures/triangulation.jl | 6 +++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/docs/src/extended/utils.md b/docs/src/extended/utils.md index 7e0072aea..a8368d1f7 100644 --- a/docs/src/extended/utils.md +++ b/docs/src/extended/utils.md @@ -33,4 +33,8 @@ DefaultAdjacentValue GhostVertex ε ∅ +fix_orient3_cache +fix_incircle_cache +validate_orient3_cache +validate_incircle_cache ``` \ No newline at end of file diff --git a/test/data_structures/triangulation.jl b/test/data_structures/triangulation.jl index 73c3a6194..ab89647e8 100644 --- a/test/data_structures/triangulation.jl +++ b/test/data_structures/triangulation.jl @@ -130,9 +130,9 @@ refine!(tri_4; max_area=1.0e-1A, use_circumcenter=true, use_lens=false) @test DT.get_ghost_vertex_ranges(tri) == tri.ghost_vertex_ranges == DT.construct_ghost_vertex_ranges(get_boundary_nodes(tri)) @test DT.get_boundary_edge_map(tri) == tri.boundary_edge_map == DT.construct_boundary_edge_map(get_boundary_nodes(tri)) @test DT.get_cache(tri) == tri.cache - @test DT.get_incircle_cache(tri) == tri.cache.incircle_cache - @test DT.get_orient3_cache(tri) == tri.cache.orient3_cache - @test DT.get_insphere_cache(tri) == tri.cache.insphere_cache + @test DT.get_incircle_cache(tri) === tri.cache.incircle_cache + @test DT.get_orient3_cache(tri) === tri.cache.orient3_cache + @test DT.get_insphere_cache(tri) === tri.cache.insphere_cache @inferred DT.get_points(tri) @inferred DT.get_triangles(tri) @inferred DT.get_boundary_nodes(tri)