Fix #31

Project.toml

@@ -1,7 +1,7 @@
 name = "NaturalNeighbours"
 uuid = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.3.1"
+version = "1.3.2"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"

src/differentiation/differentiate.jl

@@ -4,6 +4,12 @@
 Differentiate a given interpolant `itp` up to degree `order` (1 or 2). The returned object is a 
 `NaturalNeighboursDifferentiator` struct, which is callable. 
 
+!!! warning "Missing vertices"
+
+    When the underlying triangulation, `tri`, has points in `get_points(tri)` that are not 
+    vertices of the triangulation itself, the associated derivatives at these points
+    will be set to zero.
+
 For calling the resulting struct, we define the following methods:
 
     (∂::NaturalNeighboursDifferentiator)(x, y, zᵢ, nc, id::Integer=1; parallel=false, method=default_diff_method(∂), kwargs...)

src/differentiation/generate.jl

@@ -137,29 +137,38 @@ end
 @inline function generate_second_order_derivatives!(∇, ℋ, method, tri, z, zrange, derivative_caches, neighbour_caches, id; alpha, use_cubic_terms, initial_gradients)
     for i in zrange
         zᵢ = z[i]
-        d_cache = derivative_caches[id]
-        n_cache = neighbour_caches[id]
-        S = get_iterated_neighbourhood(d_cache)
-        S′ = get_second_iterated_neighbourhood(d_cache)
-        if method == Direct()
-            λ, E = get_taylor_neighbourhood!(S, S′, tri, i, 2 + use_cubic_terms, n_cache)
-        elseif method == Iterative()
-            λ, E = get_taylor_neighbourhood!(S, S′, tri, i, 1, n_cache)
+        if !DelaunayTriangulation.has_vertex(tri, i)
+            ∇[i] = (zero(zᵢ), zero(zᵢ))
+            ℋ[i] = (zero(zᵢ), zero(zᵢ), zero(zᵢ))
+        else
+            d_cache = derivative_caches[id]
+            n_cache = neighbour_caches[id]
+            S = get_iterated_neighbourhood(d_cache)
+            S′ = get_second_iterated_neighbourhood(d_cache)
+            if method == Direct()
+                λ, E = get_taylor_neighbourhood!(S, S′, tri, i, 2 + use_cubic_terms, n_cache)
+            elseif method == Iterative()
+                λ, E = get_taylor_neighbourhood!(S, S′, tri, i, 1, n_cache)
+            end
+            ∇[i], ℋ[i] = generate_second_order_derivatives(method, tri, z, zᵢ, i, λ, E, derivative_caches, id; alpha, use_cubic_terms, initial_gradients)
         end
-        ∇[i], ℋ[i] = generate_second_order_derivatives(method, tri, z, zᵢ, i, λ, E, derivative_caches, id; alpha, use_cubic_terms, initial_gradients)
     end
     return nothing
 end
 
 @inline function generate_first_order_derivatives!(∇, method, tri, z, zrange, derivative_caches, neighbour_caches, id)
     for i in zrange
         zᵢ = z[i]
-        d_cache = derivative_caches[id]
-        n_cache = neighbour_caches[id]
-        S = get_iterated_neighbourhood(d_cache)
-        S′ = get_second_iterated_neighbourhood(d_cache)
-        λ, E = get_taylor_neighbourhood!(S, S′, tri, i, 1, n_cache)
-        ∇[i] = generate_first_order_derivatives(method, tri, z, zᵢ, i, λ, E, derivative_caches, id)
+        if !DelaunayTriangulation.has_vertex(tri, i)
+            ∇[i] = (zero(zᵢ), zero(zᵢ))
+        else
+            d_cache = derivative_caches[id]
+            n_cache = neighbour_caches[id]
+            S = get_iterated_neighbourhood(d_cache)
+            S′ = get_second_iterated_neighbourhood(d_cache)
+            λ, E = get_taylor_neighbourhood!(S, S′, tri, i, 1, n_cache)
+            ∇[i] = generate_first_order_derivatives(method, tri, z, zᵢ, i, λ, E, derivative_caches, id)
+        end
     end
     return nothing
 end

src/interpolation/interpolate.jl

@@ -6,6 +6,12 @@
 Construct an interpolant over the data `z` at the sites defined by the triangulation `tri` (or `points`, or `(x, y)`). See the Output 
 section for a description of how to use the interpolant `itp`.
 
+!!! warning "Missing vertices"
+
+    When the underlying triangulation, `tri`, has points in `get_points(tri)` that are not 
+    vertices of the triangulation itself, the associated derivatives (relevant only if `derivatives=true`) at these points
+    will be set to zero.
+
 # Keyword Arguments 
 - `gradient=nothing`: The gradients at the corresponding data sites of `z`. Will be generated if `isnothing(gradient)` and `derivatives==true`.
 - `hessian=nothing`: The hessians at the corresponding data sites of `z`. Will be generated if `isnothing(hessian)` and `derivatives==true`.

test/interpolation/basic_tests.jl

@@ -137,4 +137,34 @@ end
     bad_idx = identify_exterior_points(_itp_xs, _itp_ys, points, ch; tol=1e-3) # boundary effects _really_ matter...
     deleteat!(err, bad_idx)
     @test norm(err) ≈ 0 atol = 1e-2
+end
+
+
+@testset "Test that the derivatives are all zero for missing vertices" begin
+    R₁ = 0.2
+    R₂ = 1.0
+    θ = collect(LinRange(0, 2π, 100))
+    θ[end] = 0.0 # get the endpoints to match
+    x = [
+        [R₂ .* cos.(θ)], # outer first
+        [reverse(R₁ .* cos.(θ))] # then inner - reverse to get clockwise orientation
+    ]
+    y = [
+        [R₂ .* sin.(θ)], # 
+        [reverse(R₁ .* sin.(θ))]
+    ]
+    boundary_nodes, points = convert_boundary_points_to_indices(x, y)
+    tri = triangulate(points; boundary_nodes)
+    A = get_area(tri)
+    refine!(tri; max_area=1e-4A)
+    itp = interpolate(tri, ones(DelaunayTriangulation.num_points(tri)), derivatives=true, parallel=false)
+    ind = findall(DelaunayTriangulation.each_point_index(tri)) do i
+        !DelaunayTriangulation.has_vertex(tri, i)
+    end
+    for i in ind
+        ∇ = NaturalNeighbours.get_gradient(itp, i)
+        @test all(iszero, ∇)
+        H = NaturalNeighbours.get_hessian(itp, i)
+        @test all(iszero, H)
+    end
 end

test/interpolation/constrained.jl

@@ -62,6 +62,7 @@ end
     tri = triangulate(points; boundary_nodes)
     A = get_area(tri)
     D = 6.25e-4
+    refine!(tri; max_area=1e-2A)
     Tf = (x, y) -> let r = sqrt(x^2 + y^2)
         (R₂^2 - r^2) / (4D) + R₁^2 * log(r / R₂) / (2D)
     end