From 5b0fef42feaf7b1334b50ae9629621adc1d8caee Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 28 Jul 2024 20:04:31 +0530 Subject: [PATCH 1/3] Add comments explaining `compute_bowyer_envelope!` --- src/interpolation/utils.jl | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/src/interpolation/utils.jl b/src/interpolation/utils.jl index 3f33f9d..fd5eeaa 100644 --- a/src/interpolation/utils.jl +++ b/src/interpolation/utils.jl @@ -1,14 +1,29 @@ function compute_bowyer_envelope!(envelope, tri::Triangulation, history::InsertionEventHistory, temp_adjacent::Adjacent, point; kwargs...) #kwargs are add_point! kwargs + #= + `add_point!` uses the Bowyer-Watson algorithm to add points. + The Bowyer-Watson algorithm removes all triangles whose circumcenters contain + the point to be inserted, creating a hole which is then retriangulated with the + new point included. + What happens here is that we record the new triangles added in `history`, and then + extract those triangles - we know that these changed after the new point was added, + so they are the "natural neighbours" of the new point. + =# empty!(history) empty!(envelope) empty!(get_adjacent(temp_adjacent)) n = num_points(tri) I = integer_type(tri) + # Note that the `peek` keyword here indicates that the original triangulation is + # never modified - but the new triangles are added to `history`. V = add_point!(tri, point; store_event_history=Val(true), event_history=history, peek=Val(true), kwargs...) all_triangles = each_added_triangle(history) if isempty(all_triangles) # This is possible for constrained triangulations + # To clarify, this means that no triangles were added. This means the point + # has no natural neighbours. return envelope, temp_adjacent, history, V end + # Add triangles to the `DelaunayTriangulation.Adjacent` dict, + # which stores the ((u, v)::Edge => w::Point` adjacency relationship. for T in all_triangles add_triangle!(temp_adjacent, T) end @@ -17,10 +32,14 @@ function compute_bowyer_envelope!(envelope, tri::Triangulation, history::Inserti v = i == I(n + 1) ? j : i # get a vertex on the envelope push!(envelope, v) for i in 2:num_triangles(all_triangles) - v = get_adjacent(temp_adjacent, I(n + 1), v) + # Get the vertex associated with the edge `(point, points[v])` such that the + # returned point (the new `v`) forms a positively oriented triangle + # with that edge. Then, add the newly found `v` to the envelope. + # To be clear, `n+1` is the index of the newly inserted vertex in the triangulation. + v = get_adjacent(temp_adjacent, I(n + 1), v) push!(envelope, v) end - push!(envelope, envelope[begin]) + push!(envelope, envelope[begin]) # Close out the envelope polygon return envelope, temp_adjacent, history, V end function compute_bowyer_envelope!(envelope, tri::Triangulation, point; kwargs...) @@ -241,4 +260,4 @@ function hessian_form(tri, u, v, w, N₀, H) # zᵢ,ⱼₖ dxᵢₖ = xₖ[1] - xᵢ[1] dyᵢₖ = xₖ[2] - xᵢ[2] return dxᵢⱼ * (dxᵢₖ * B₁₁ + dyᵢₖ * B₁₂) + dyᵢⱼ * (dxᵢₖ * B₁₂ + dyᵢₖ * B₂₂) -end \ No newline at end of file +end From b176fac8c096ec902d278cadb80fff808ba13529 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 28 Jul 2024 20:06:36 +0530 Subject: [PATCH 2/3] update styling --- src/interpolation/utils.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/interpolation/utils.jl b/src/interpolation/utils.jl index fd5eeaa..075ab96 100644 --- a/src/interpolation/utils.jl +++ b/src/interpolation/utils.jl @@ -4,6 +4,7 @@ function compute_bowyer_envelope!(envelope, tri::Triangulation, history::Inserti The Bowyer-Watson algorithm removes all triangles whose circumcenters contain the point to be inserted, creating a hole which is then retriangulated with the new point included. + What happens here is that we record the new triangles added in `history`, and then extract those triangles - we know that these changed after the new point was added, so they are the "natural neighbours" of the new point. @@ -13,8 +14,8 @@ function compute_bowyer_envelope!(envelope, tri::Triangulation, history::Inserti empty!(get_adjacent(temp_adjacent)) n = num_points(tri) I = integer_type(tri) - # Note that the `peek` keyword here indicates that the original triangulation is - # never modified - but the new triangles are added to `history`. + # Note that the `peek` keyword here indicates that the original triangulation + # is never modified - but the new triangles are added to `history`. V = add_point!(tri, point; store_event_history=Val(true), event_history=history, peek=Val(true), kwargs...) all_triangles = each_added_triangle(history) if isempty(all_triangles) # This is possible for constrained triangulations @@ -32,10 +33,12 @@ function compute_bowyer_envelope!(envelope, tri::Triangulation, history::Inserti v = i == I(n + 1) ? j : i # get a vertex on the envelope push!(envelope, v) for i in 2:num_triangles(all_triangles) - # Get the vertex associated with the edge `(point, points[v])` such that the - # returned point (the new `v`) forms a positively oriented triangle - # with that edge. Then, add the newly found `v` to the envelope. - # To be clear, `n+1` is the index of the newly inserted vertex in the triangulation. + #= + Get the vertex associated with the edge `(point, points[v])` such that the + returned point (the new `v`) forms a positively oriented triangle + with that edge. Then, add the newly found `v` to the envelope. + To be clear, `n+1` is the index of the newly inserted vertex in the triangulation. + =# v = get_adjacent(temp_adjacent, I(n + 1), v) push!(envelope, v) end From d324aa9a8bcc18bf38e3d404c4096cf0ab0950e2 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 28 Jul 2024 16:22:52 +0100 Subject: [PATCH 3/3] Thought I fixed this already --- test/doc_examples/differentiation.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/doc_examples/differentiation.jl b/test/doc_examples/differentiation.jl index 788d3df..25733d8 100644 --- a/test/doc_examples/differentiation.jl +++ b/test/doc_examples/differentiation.jl @@ -101,12 +101,12 @@ tri = triangulate(points) ∇g = generate_gradients(tri, z) fig, ε = plot_gradients(∇g, tri, f′, x, y) @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "gradient_data.png") fig -@test ε ≈ 10.2511800 +@test ε ≈ 10.2511800 rtol=1e-3 ∇gr, _ = generate_derivatives(tri, z; method=Direct()) fig, ε = plot_gradients(∇gr, tri, f′, x, y) @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "joint_gradient_data.png") fig -@test ε ≈ 7.7177915977 +@test ε ≈ 7.7177915977 rtol=1e-3 # Hessians to_mat(H::NTuple{3,Float64}) = [H[1] H[3]; H[3] H[2]] @@ -128,17 +128,17 @@ function plot_hessians(H, tri, f′′, x, y) end _, Hg = generate_derivatives(tri, z) fig, ε = plot_hessians(Hg, tri, f′′, x, y) -@test ε ≈ 42.085578794 +@test ε ≈ 42.085578794 rtol=1e-3 @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "hessian_data.png") fig _, Hg = generate_derivatives(tri, z, use_cubic_terms=false) fig, ε = plot_hessians(Hg, tri, f′′, x, y) -@test ε ≈ 35.20873081 +@test ε ≈ 35.20873081 rtol=1e-3 @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "hessian_data_no_cubic.png") fig _, Hg = generate_derivatives(tri, z, method=Iterative()) # the gradients will be generated first automatically fig, ε = plot_hessians(Hg, tri, f′′, x, y) -@test ε ≈ 39.584816 +@test ε ≈ 39.584816 rtol=1e-3 @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "hessian_data_iterative.png") fig # Away from the data sites @@ -190,7 +190,7 @@ function plot_gradients(∇g, f′, xg, yg) end fig, ε = plot_gradients(∇g, f′, xg, yg) @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "gradient_surface.png") fig -@test ε ≈ 13.1857476 +@test ε ≈ 13.1857476 rtol=1e-1 other_methods = [Sibson(), Laplace(), Nearest(), Triangle()] ∇gs = [∂(_x, _y; interpolant_method=method) for method in other_methods] @@ -236,8 +236,8 @@ figH, εH = plot_hessians(Hg, f′′, xg, yg) zlims!(figH.content[4], -25, 25) @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "hessian_surface_direct.png") figH @test_reference normpath(@__DIR__, "../..", "docs", "src", "figures", "gradient_surface_2_direct.png") fig∇ -@test εH ≈ 46.7610990050276 -@test ε∇ ≈ 9.853286514069882 +@test εH ≈ 46.7610990050276 rtol=1e-1 +@test ε∇ ≈ 9.853286514069882 rtol=1e-1 function rrmserr(z, ẑ, ∂, x, y) tri = ∂.interpolant.triangulation @@ -259,5 +259,5 @@ end εH_nohull = rrmserr(f′′.(_x, _y), to_mat.(Hg), ∂, _x, _y) ε∇_nohull = rrmserr(f′.(_x, _y), ∇g, ∂, _x, _y) -@test ε∇_nohull ≈ 7.479964687673911 -@test εH_nohull ≈ 38.884740966379056 \ No newline at end of file +@test ε∇_nohull ≈ 7.479964687673911 rtol=1e-2 +@test εH_nohull ≈ 38.884740966379056 rtol=1e-2