Work on docs

examples/interpolation.jl

@@ -85,3 +85,56 @@ current_figure() # hide
 
 # Clearly, the spurious oscillations are strongly suppressed near the
 # boundaries.
+
+# ## Multidimensional interpolations
+#
+# Multidimensional interpolations are supported for data defined on
+# rectilinear grids.
+# This is done using
+# [tensor-product](https://en.wikipedia.org/wiki/Tensor_product) splines.
+# This means that the returned $N$-dimensional spline belongs to the space
+# spanned by the tensor product of $N$ one-dimensional spline bases.
+#
+# In other words, if we consider the two-dimensional case for simplicity, the
+# returned interpolator is of the form
+#
+# ```math
+# f(x, y) = ∑_{i = 1}^{N_x} ∑_{j = 1}^{N_y}
+# c_{ij} \, b_{i}^X(x) \, b_{j}^Y(y),
+# ```
+#
+# where the ``b_{i}^X`` and ``b_{j}^Y`` are the B-splines spanning two
+# (generally different) spline bases ``B^X`` and ``B^Y``.
+# Note that the two bases are completely independent.
+# In principle, they can have different orders, knot definitions and boundary
+# conditions.
+#
+# All of the above also applies to (arbitrary) higher dimensions ``N``.
+# See [`interpolate`](@ref) for more details.
+#
+# Let's finish with a simple example in 2D:
+
+## Define non-uniform 2D grid
+Nx, Ny = 20, 30
+rng = MersenneTwister(42)
+
+xs = range(0, 1; length = Nx) .+ 0.01 .* randn(rng, Nx)
+xs[begin] = 0
+xs[end] = 1
+
+ys = range(0, 2π; length = Ny) .+ 0.02 .* randn(rng, Ny)
+ys[begin] = 0
+ys[end] = 2π
+
+## Generate some 2D data and interpolate
+data = [exp(-x) * cos(y) + 0.01 * randn(rng) for x ∈ xs, y ∈ ys]
+
+S = interpolate((xs, ys), data, BSplineOrder(4), Natural())
+
+## Finally, let's plot the result on a finer grid
+xplot = range(0, 1; length = 4Nx)
+yplot = range(0, 2π; length = 4Ny)
+Sdata = S.(xplot, yplot')
+
+contourf(xplot, yplot, Sdata)
+current_figure() # hide

src/SplineInterpolations/SplineInterpolations.jl

@@ -215,8 +215,27 @@ For now, only the [`Natural`](@ref) boundary condition is available.
 
 See also [`interpolate!`](@ref).
 
+# Multidimensional interpolations
+
+Multidimensional (tensor-product) interpolations are supported on arbitrary
+dimensions `N`.
+For that, the following syntax should be used:
+
+ interpolate(xs, y, BSplineOrder(k), [bc = nothing])
+
+where `xs = (x1, x2, …, xN)` is a tuple of vectors containing the grid points
+along each direction, and `y` is an `N`-dimensional array containing the data to
+be interpolated.
+See further below for some examples.
+
+For now, the B-spline order and the boundary conditions are the same along all
+dimensions.
+This constraint may be generalised in the future.
+
 # Examples
 
+## One-dimensional interpolations
+
 ```jldoctest
 julia> xs = -1:0.1:1;
 
@@ -256,6 +275,34 @@ julia> (Derivative(1) * Snat)(-1)
 julia> (Derivative(2) * Snat)(-1)
 0.0
 
+```
+
+## Two-dimensional interpolations
+
+```jldoctest
+julia> x₁ = -1:0.2:1; x₂ = 0:0.1:0.8;
+
+julia> ydata = @. cospi(x₁) * sinpi(x₂');
+
+julia> summary(ydata)
+"11×9 Matrix{Float64}"
+
+julia> itp = interpolate((x₁, x₂), ydata, BSplineOrder(4), Natural())
+SplineInterpolation containing the 11×9 Spline{2, Float64}:
+ bases:
+ (1) 11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0], BCs {left => (D{2},), right => (D{2},)}
+ (2) 9-element RecombinedBSplineBasis of order 4, domain [0.0, 0.8], BCs {left => (D{2},), right => (D{2},)}
+ knots:
+ (1) [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
+ (2) [0.0, 0.0, 0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.8, 0.8, 0.8]
+ coefficients: [0.0 -0.174524 … -0.458359 -0.408184; 0.0 -0.123177 … -0.323506 -0.288093; … ; 0.0 -0.123177 … -0.323506 -0.288093; 0.0 -0.174524 … -0.458359 -0.408184]
+ interpolation points:
+ (1) -1.0:0.2:1.0
+ (2) 0.0:0.1:0.8
+
+julia> itp(0.12, 0.4242)
+0.9032397652177311
+
 ```
 """
 function interpolate(

src/Splines/Splines.jl

@@ -4,7 +4,8 @@ export
  Spline,
  Spline1D,
  coefficients,
- integral
+ integral,
+ bases
 
 export
  SplineWrapper,

src/Splines/spline.jl

@@ -197,8 +197,8 @@ Base.size(S::Spline) = size(coefficients(S))
 
 Returns type of element returned when evaluating the [`Spline`](@ref).
 """
-Base.eltype(::Type{<:Spline{N, T}}) where {N, T} = T
 Base.eltype(S::Spline) = eltype(typeof(S))
+Base.eltype(::Type{<:Spline{N, T}}) where {N, T} = T
 
 """
  ndims(::Type{<:Spline})