Add field gradient and hessian

Project.toml

@@ -9,8 +9,8 @@ RegionTrees = "dee08c22-ab7f-5625-9660-a9af2021b33f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-Interpolations = "0.9, 0.10, 0.11"
-RegionTrees = "0.2"
+Interpolations = "0.9, 0.10, 0.11, 0.12"
+RegionTrees = "0.2, 0.3"
 julia = "1"
 
 [extras]

src/AdaptiveDistanceFields.jl

@@ -6,15 +6,20 @@ using RegionTrees
 import RegionTrees: needs_refinement, refine_data
 using Interpolations: interpolate!,
                       extrapolate,
+                      gradient,
+                      hessian,
                       AbstractInterpolation,
                       BSpline,
                       OnGrid,
                       Linear,
                       Line
+
 using StaticArrays: SVector
 
 export AdaptiveDistanceField,
        evaluate,
+       fieldgradient,
+       fieldhessian,
        ConvexMesh
 
 include("interpolation.jl")
@@ -39,4 +44,13 @@ function (field::AdaptiveDistanceField)(point::AbstractArray)
     evaluate(field.root, point)
 end
 
+
+function fieldgradient(field::AdaptiveDistanceField,point::AbstractArray)
+    fieldgradient(field.root, point)
+end
+
+function fieldhessian(field::AdaptiveDistanceField,point::AbstractArray)
+    fieldhessian(field.root, point)
+end
+
 end

src/interpolation.jl

@@ -3,7 +3,7 @@
 end
 
 function evaluate(itp::AbstractInterpolation, point::AbstractArray)
-    itp[point...]
+    itp(point...)
 end
 
 function evaluate(cell::Cell{D}, point::AbstractArray) where {D <: AbstractInterpolation}
@@ -12,6 +12,50 @@ function evaluate(cell::Cell{D}, point::AbstractArray) where {D <: AbstractInter
 end
 
 function evaluate(interp::AbstractInterpolation, boundary::HyperRectangle, point::AbstractArray)
-    coords = (point - boundary.origin) ./ boundary.widths + 1
+    coords = (point - boundary.origin) ./ boundary.widths .+ 1
     evaluate(interp, coords)
 end
+
+
+
+
+
+@generated function fieldgradient(itp::AbstractInterpolation, point::SVector{N}) where N
+    Expr(:call, :gradient, :itp, [:(point[$i]) for i in 1:N]...)
+end
+
+function fieldgradient(cell::Cell{D}, point::AbstractArray) where {D <: AbstractInterpolation}
+    leaf = findleaf(cell, point)
+    fieldgradient(leaf.data, leaf.boundary, point)
+end
+
+function fieldgradient(interp::AbstractInterpolation, boundary::HyperRectangle, point::AbstractArray)
+    coords = (point - boundary.origin) ./ boundary.widths .+ 1
+    fieldgradient(interp, coords) ./ boundary.widths
+end
+
+function fieldgradient(itp::AbstractInterpolation, point::AbstractArray)
+    gradient(itp, point...)
+end
+
+
+
+
+
+@generated function fieldhessian(itp::AbstractInterpolation, point::SVector{N}) where N
+    Expr(:call, :hessian, :itp, [:(point[$i]) for i in 1:N]...)
+end
+
+function fieldhessian(cell::Cell{D}, point::AbstractArray) where {D <: AbstractInterpolation}
+    leaf = findleaf(cell, point)
+    fieldhessian(leaf.data, leaf.boundary, point)
+end
+
+function fieldhessian(interp::AbstractInterpolation, boundary::HyperRectangle, point::AbstractArray)
+    coords = (point - boundary.origin) ./ boundary.widths .+ 1
+    fieldhessian(interp, coords) ./ (boundary.widths*boundary.widths')
+end
+
+function fieldhessian(itp::AbstractInterpolation, point::AbstractArray)
+    hessian(itp, point...)
+end