Bases

Three concrete bases implement AbstractBSplineBasis, and which one you get is decided by the boundary condition rather than chosen by name.

using SimpleSplines

Construction

One constructor covers all three:

BSplineBasis(mesh, p)                 # clamped, the default
BSplineBasis(mesh, p, bc)             # dispatches on bc

m = UniformMesh(8, 0 .. 1)
[(bc, nameof(typeof(BSplineBasis(m, 3, bc))), nbasis(BSplineBasis(m, 3, bc)))
 for bc in (Free(), (Free(), Free()), Periodic(), :periodic,
     Dirichlet(), (Dirichlet(), Free()), Neumann(), Natural())]

8-element Vector{Tuple{Any, Symbol, Int64}}:
 (Free(), :BSplineBasis, 11)
 ((Free(), Free()), :BSplineBasis, 11)
 (Periodic(), :PeriodicBSplineBasis, 8)
 (:periodic, :PeriodicBSplineBasis, 8)
 (Dirichlet(), :RecombinedBSplineBasis, 9)
 ((Dirichlet(), Free()), :RecombinedBSplineBasis, 10)
 (Neumann(), :RecombinedBSplineBasis, 9)
 (Natural(), :RecombinedBSplineBasis, 9)

Periodic() gives a PeriodicBSplineBasis, two Free ends give the plain clamped BSplineBasis unwrapped, and anything else gives a RecombinedBSplineBasis. This is deliberate: it lets one call site select any of the three, which is exactly what a tensor product with a different condition per axis needs.

The named constructors are also available directly, and PeriodicBSplineBasis has a convenience form that builds its own uniform mesh:

BSplineBasis(mesh, p)                     BSplineBasis{T}(mesh, p)
PeriodicBSplineBasis(mesh, p)             PeriodicBSplineBasis{T}(mesh, p)
PeriodicBSplineBasis(n, p; L = 2π)        # = PeriodicBSplineBasis(UniformMesh(n, L), p)
RecombinedBSplineBasis(parent, left, right)

PeriodicBSplineBasis(16, 3) == BSplineBasis(UniformMesh(16, 2π), 3, Periodic())

true

What is required, and what is rejected

basis	requires
`BSplineBasis`	$p \ge 0$, $n \ge 1$ — a single cell carries the whole of $\mathbb{P}_p$
`PeriodicBSplineBasis`	$p \ge 0$ and $n > p$, so that a basis function does not wrap onto itself
`RecombinedBSplineBasis`	a clamped parent; each condition of order $\le p$; the two end blocks disjoint; at least one function left

for bad in (() -> BSplineBasis(UniformMesh(8, 0 .. 1), -1),
    () -> PeriodicBSplineBasis(UniformMesh(3, 0 .. 1), 3),
    () -> BSplineBasis(UniformMesh(8, 0 .. 1), 1, Natural()),
    () -> BSplineBasis(UniformMesh(1, 0 .. 1), 0, Dirichlet()),
    () -> RecombinedBSplineBasis(BSplineBasis(UniformMesh(8, 0 .. 1), 3),
        Periodic(), Dirichlet()))
    try
        bad()
    catch err
        println(err.msg, "\n")
    end
end

the degree of a B-spline basis must be non-negative, got p = -1

a periodic B-spline basis of degree 3 needs more than 3 cells, got n = 3; refine the mesh or lower the degree

the left boundary condition involves D^2, but every derivative above order 1 of a degree-1 spline vanishes identically, so the condition constrains nothing

a degree-0 basis on 1 cells has only 1 functions, too few to impose Dirichlet() on the left and Dirichlet() on the right; refine the mesh

Periodic identifies the two ends of the domain rather than constraining one, so it cannot be imposed on a clamped basis by recombination; use `PeriodicBSplineBasis(mesh, p)` for the whole axis

Accessors

b = BSplineBasis(UniformMesh(8, 0 .. 1), 3)
(nbasis(b), degree(b), order(b), ncells(b), local_width(b),
    polynomial_reproduction(b), boundary(b))

(11, 3, 4, 8, 4, 3, (Free(), Free()))

call	returns
`nbasis(b)`	the dimension $N$, and the length a coefficient vector must have
`degree(b)`	the polynomial degree $p$
`order(b)`	$p + 1$
`ncells(b)`, `mesh(b)`, `domain(b)`, `domainlength(b)`, `meshwidth(b)`	forwarded to the mesh
`breakpoints(b)`	the `n+1` cell boundaries, cached in the basis
`knotvector(b)`	the knot sequence the recursion runs on
`boundary(b)`	`Periodic()`, or the pair `(left, right)`
`nodes(b)`	the Greville abscissae; `ContinuumArrays`' `grid` is the same
`nnodes(b)`	`nbasis(b)`
`local_width(b)`	how many entries `evaluate_all!` writes: $p+1$, or more near a recombined end
`polynomial_reproduction(b)`	the largest degree the span reproduces exactly, or `-1`
`basis_index(b, j)`	`j`, wrapped onto `1:N` where the basis is periodic
`eachindex(b)`, `axes(b)`, `eltype(b)`	`1:N`, `(Inclusion(domain), 1:N)`, the element type

`order` here means `p + 1`

That is the spline convention: a B-spline of order $k$ is piecewise of degree $k - 1$. It is not the meaning order carries for the bases of CompactBasisFunctions, where it is the number of basis functions. For a spline basis the two differ, and nbasis is the one you want for a dimension.

Evaluating

evaluate(b, j, x, d = 0)                # basis function j
evaluate(b, û, x, d = 0)                # the spline with coefficients û
b[x, j]     b[x, :]     b[X, j]     b[X, :]        b(x, j)
b'          # a BSplineDerivative: first order only

x may be a number or a vector; with a vector the result is a vector (or, for b[X, :], a matrix of points × functions).

xs = [0.1, 0.5, 0.9]
evaluate(b, 3, 0.3), evaluate(b, 3, 0.3, 2), size(b[xs, :]), b'[0.3, 3]

(0.03600000000000002, 38.400000000000006, (3, 11), -1.4400000000000004)

The argument order is not the same in the two spellings

evaluate takes the index first, evaluate(b, j, x, d), because the derivative order comes last. getindex and the callable form take the point first, b[x, j] and b(x, j), to match a matrix of samples whose rows are points. Both spellings are correct; mixing them up is silent whenever j and x are both plausible numbers.

Derivatives of arbitrary order

The fourth argument of evaluate is the derivative order, any non-negative integer. b' is a BSplineDerivative, which is a lazy product and first order only — the products do not compose, so a second derivative is not b''.

[evaluate(b, 3, 0.3, d) for d in 0:4]

5-element Vector{Float64}:
    0.03600000000000002
   -1.4400000000000004
   38.400000000000006
 -512.0
    0.0

Orders above p are identically zero and are returned as such rather than raised, so a loop over derivative orders needs no special case. A negative order is an error.

Outside the domain

A bounded basis evaluates to zero outside its closed domain — not an error, and not an extrapolated polynomial. That is the mathematically correct value for a compactly supported function, and it is the behaviour a particle method wants; it is also silent, so it is worth knowing.

sum(b[1.7, :]), evaluate(b, 4, -0.2)

(0.0, 0.0)

A PeriodicBSplineBasis instead accepts any real argument and reduces it onto the domain, so the result is the periodic extension:

c = BSplineBasis(UniformMesh(8, 0 .. 1), 3, Periodic())
sum(c[1.7, :]), evaluate(c, 4, 0.7) ≈ evaluate(c, 4, 0.7 + 3.0)

(1.0000000000000002, true)

At the right endpoint

Knot spans are half-open, except the topmost, which is closed on the right so that the clamped basis is interpolatory at $b$ as it should be. The value there is exactly the limit from the left.

b[1.0, nbasis(b)], sum(b[1.0, :])

(1.0, 1.0)

The local path

A loop that needs every nonzero function at a point — a matrix assembly, a particle deposition — should not call evaluate once per index. That path runs the recursion literally, at $O(2^p)$ per value; evaluate_all! fills a caller-supplied buffer with the whole local block by de Boor's triangular scheme, in $O(p^2)$ for all p+1 together, and allocates nothing:

buf = zeros(local_width(b))
j₀ = evaluate_all!(buf, b, 0.3)
j₀, buf, sum(buf)

(3, [0.03600000000000002, 0.5386666666666667, 0.41466666666666663, 0.01066666666666666], 1.0)

all(buf[t] ≈ evaluate(b, basis_index(b, j₀ + t - 1), 0.3) for t in eachindex(buf))

true

Three rules for using it:

The buffer must have exactly local_width(b) entries. That is p+1 for a clamped or periodic basis, and possibly more for a recombined one, where a function near an end spans the union of two parent supports. A wrong length is a DimensionMismatch.
j₀ is the first index before wrapping, because the block is contiguous only in that form. On a periodic basis it can fall outside 1:N; put every index through basis_index rather than writing mod1 at the call site, since whether the wrap is needed is a property of the basis and getting it wrong on a bounded basis silently folds the two ends of the domain together.
Outside a bounded domain the buffer is zeroed and the clamped cell's j₀ is returned, so a deposition loop touches the same block it would for a point just inside and adds nothing to it.

buf2 = zeros(local_width(c))
jc = evaluate_all!(buf2, c, 0.05)             # near the periodic seam
jc, [basis_index(c, jc + t - 1) for t in eachindex(buf2)]

(-2, [6, 7, 8, 1])

The allocating form evaluate_all returns (j₀, values) with a fresh buffer, which is convenient at a call site that runs once and wasteful in a loop.

A deposition therefore reads:

coeffs = zeros(nbasis(b))
particles = [(0.13, 1.0), (0.42, 2.0), (0.97, 0.5), (1.4, 3.0)]   # the last is outside
scratch = zeros(local_width(b))
for (x, w) in particles
    j = evaluate_all!(scratch, b, x)
    for t in eachindex(scratch)
        i = basis_index(b, j + t - 1)
        1 ≤ i ≤ nbasis(b) || continue
        coeffs[i] += w * scratch[t]
    end
end
coeffs, sum(coeffs)

([0.0, 0.221184, 0.5913706666666667, 0.27481600000000006, 1.1208, 0.7762773333333332, 0.015551999999999983, 0.001152000000000003, 0.036864000000000056, 0.24249600000000016, 0.21948799999999982], 3.4999999999999996)

The sum is $1.0 + 2.0 + 0.5 = 3.5$ and not $6.5$: the partition of unity means each particle deposits its whole weight, and the one outside the domain deposited nothing.

Cell lookup

findcell(b, x) gives the cell index in 1:ncells, reading the breakpoints cached in the basis so that it allocates nothing; on a uniform mesh it is a division that does not touch them at all. local_indices(b, cell) gives the indices of the functions nonzero on a cell — cell:(cell+p) for a clamped basis, the same block wrapped for a periodic one, and a precomputed range for a recombined one.

findcell(b, 0.3), collect(local_indices(b, 3)), collect(local_indices(c, 1))

(3, [3, 4, 5, 6], [6, 7, 8, 1])

Note the periodic case: cell 1 is reached by functions 6, 7, 8, 1, wrapped.

Comparison

== compares type, degree, mesh and boundary condition; isequal additionally requires the same element type; isapprox compares the meshes approximately.

b == BSplineBasis(UniformMesh(8, 0 .. 1), 3),
b == BSplineBasis(UniformMesh(8, 0 .. 1), 3, Dirichlet()),
b == c

(true, false, false)