do -not- use Interpolations.jl

Author: AbdAlazezAhmed

Project.toml

@@ -5,7 +5,6 @@ version = "1.1.0"
 [deps]
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
@@ -31,7 +30,6 @@ BlockArrays = "0.16, 1"
 EnumX = "1"
 ForwardDiff = "0.10, 1"
 HCubature = "1.7.0"
-Interpolations = "0.16.1"
 Metis = "1.3"
 NearestNeighbors = "0.4"
 OrderedCollections = "1"

test/runtests.jl

@@ -12,7 +12,6 @@ using OrderedCollections
 using WriteVTK
 import Metis
 using HCubature: hcubature, hquadrature
-using Interpolations
 
 include("test_utils.jl")
 

test/test_l2_projection.jl

@@ -53,8 +53,7 @@ function test_projection(order, elementtype)
         elseif dim == 3
             f_res = [f((x_1, x_2, x_3)) for x_1 in qp_1D_coord, x_2 in qp_1D_coord, x_3 in qp_1D_coord]
         end
-        interp_linear = linear_interpolation(ntuple(_ -> qp_1D_coord, dim), f_res; extrapolation_bc = Interpolations.Line())
-        ae = [interp_linear(coords...) for coords in cellcoords]
+        ae = [interp_linear(qp_1D_coord, f_res, coords...) for coords in cellcoords]
     elseif order == 2
         # For a quadratic approximation the analytical solution is recovered
         ae = zeros(length(point_vars))
@@ -67,7 +66,7 @@ function test_projection(order, elementtype)
     qp_values = analytical(f_vector)
     point_vars = project(proj, qp_values, qr)
     if order == 1
-        ae = [Vec{1, Float64}((interp_linear(coords...),)) for coords in cellcoords]
+        ae = [Vec{1, Float64}((interp_linear(qp_1D_coord, f_res, coords...),)) for coords in cellcoords]
     elseif order == 2
         ae = zeros(length(point_vars))
         apply_analytical!(ae, proj.dh, :_, x -> f_vector(x)[1])
@@ -82,7 +81,7 @@ function test_projection(order, elementtype)
     point_vars = project(proj, qp_values, qr)
     point_vars_2 = project(proj, qp_values_matrix, qr)
     if order == 1
-        ae = [Tensor{2, 2, Float64}((interp_linear(coords...), 2 * interp_linear(coords...), 3 * interp_linear(coords...), 4 * interp_linear(coords...))) for coords in cellcoords]
+        ae = [Tensor{2, 2, Float64}((interp_linear(qp_1D_coord, f_res, coords...), 2 * interp_linear(qp_1D_coord, f_res, coords...), 3 * interp_linear(qp_1D_coord, f_res, coords...), 4 * interp_linear(qp_1D_coord, f_res, coords...))) for coords in cellcoords]
     elseif order == 2
         ae = zeros(4, length(point_vars))
         for i in 1:4
@@ -99,7 +98,7 @@ function test_projection(order, elementtype)
     point_vars = project(proj, qp_values, qr)
     point_vars_2 = project(proj, qp_values_matrix, qr)
     if order == 1
-        ae = [SymmetricTensor{2, 2, Float64}((interp_linear(coords...), 2 * interp_linear(coords...), 3 * interp_linear(coords...))) for coords in cellcoords]
+        ae = [SymmetricTensor{2, 2, Float64}((interp_linear(qp_1D_coord, f_res, coords...), 2 * interp_linear(qp_1D_coord, f_res, coords...), 3 * interp_linear(qp_1D_coord, f_res, coords...))) for coords in cellcoords]
     elseif order == 2
         ae = zeros(3, length(point_vars))
         for i in 1:3

test/test_utils.jl

@@ -383,3 +383,21 @@ function single_element_grid(::Type{CellType}) where {RefShape, CellType <: Ferr
     cells = [CellType(ntuple(i -> i, length(nodes)))]
     return Grid(cells, nodes)
 end
+
+interpolate_linear(t₁, t₂, y₁, y₂, t) = t₂>t₁ ? y₁*(t₂-t)/(t₂-t₁) .+ y₂*(t-t₁)/(t₂-t₁) : y₁
+# generic tensor-product interpolator using corner weights; handles any dim
+interp_linear = function (xs, ys, coords...)
+    # locate interval and compute weight for each dimension
+    is = map(c -> clamp(searchsortedlast(xs, c), 1, length(xs)-1), coords)
+    ws = [(coords[k] - xs[i])/(xs[i+1] - xs[i]) for (k,i) in enumerate(is)]
+    val = zero(eltype(ys))
+    for corner in Iterators.product((0:1 for _ in coords)...) # iterate over 2^dim corners
+        weight = one(val)
+        idxs = ntuple(k -> is[k] + corner[k], length(coords))
+        for k in 1:length(coords)
+            weight *= corner[k] == 0 ? 1 - ws[k] : ws[k]
+        end
+        val += weight * getindex(ys, idxs...)
+    end
+    return val
+end